Equation of State Under External Stress from Crystals to Non-Crystals

1. Introduction

The Equation of State (EOS) of a system has a very long history, with the first publication in 1662 [1]. The main purpose of it is to yield the specific fixed volume of a given material under any given external thermal and mechanical conditions. Then, the system must be in a macroscopic equilibrium state.

In fact, the rigorous EOS for any material system can be derived based on the following widely used de facto “theorem" in Statistical Physics, which reflects the first law of thermodynamics. As shown with the Equations (2.95), (3.3), and (3.129) of the statistical thermodynamics book [2] by Bellac, Mortessagne, and Batrouni, the “theorem" may be stated as:

If the infinitesimal work done by the EXTERNAL forces on a system in a macroscopic equilibrium state is written in the form of

$$dW=A_{1}d{B}_1+A_{2}d{B}_2+\dots +A_{m}d{B}_m,$$

(1)

where $B_1$, $B_2$, …, and $B_m$ are variables independent of each other, then, for any pair of the conjugate variables $A_i$ and $B_i$, one has

$$A_i=-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial lnZ}{\partial B_i}}(i=1,2,\dots ,m),$$

(2)

where Z is the partition function of the system, and $\beta =1/\left(kT\right)$ with Boltzmann constant k and temperature T.

Now let us consider such a system but under isotropic external pressure P. Then the infinitesimal work done by the external force is

$$dW=-PdV,$$

(3)

where V is the system volume. Then, based on the above “theorem", we have the EOS of the system immediately:

$$P={\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial lnZ}{\partial V}}.$$

(4)

This EOS, applying to any system under external isotropic pressure, is taught in almost any related books [2,3].

2. Crystals Under External Anisotropic Mechanical Conditions

2.1. Crystal Period Vectors and External Stress

A crystal is made of a three-dimensional periodic arrangement of the same cells concisely and tidily. Each cell, containing the same particles, is a parallelepiped. Since a crystal is essentially such an unlimited periodic structure, let us call the edge vectors of the cells as the (crystal) period vectors, shown as $\mathbf{a},\mathbf{b},\mathrm{and}\mathbf{c}$ red arrows in Figure 1.

The positions of all particles inside the cell, and the crystal period vectors are usually measured through X-ray diffraction experiments. From the pure theory point of view, the particles’ positions can also be calculated by applying the Newtonian Dynamics or Quantum Mechanics. Then what is the theory to calculate the period vectors? Furthermore, crystals may experience rather complicated external mechanical actions. For example, in order to get electric or magnetic field, an additional external pressure is applied onto piezoelectric and piezomagnetic crystals, but only in ONE direction. Here, the external pressures are not the same in all directions. Since the EOS of Eq. (4), for an isotropic pressure, accepting only one pressure value, cannot apply in this situation.

Crystals may also experience external forces containing components parallel to the crystal surfaces, as shown in Figure 2. In general, external forces can be described by the so-called stress. Like pressure (the normal force per unit area), roughly speaking, stress is the total force acting on a surface per unit area “vector". It is expressed as a second-rank tensor ($3\times 3$ matrix). The force by the stress $\mathbf{S}$ acting on the area vector $\mathbf{\sigma }$ of a surface is $\mathbf{S}·\mathbf{\sigma }$.

Still as shown in Figure 2, the crystal period vectors may change independently with each other according to the external stress $\mathbf{S}$ applied. As the period vectors determine the shape and volume of the cell, then proportional to the macroscopic shape and volume of crystals, the equation determining the period vectors is also the EOS of the crystals. It was derived based on the above “de facto `theorem ’" for general external stress $\mathbf{S}$ in our recent effort [4,5]. In the following, we will repeat it briefly.

2.2. The Derivation of the EOS

Now let us write the infinitesimal work done by the external stress $\mathbf{S}$ acting on the surfaces of a crystal cell, which is proportional to that on the surfaces of the whole crystal:

$$\begin{array}{c}\hfill dW=\left(\mathbf{S}·{\mathbf{\sigma }}_{\mathbf{a}}\right)·\mathrm{d}\mathbf{a}+\left(\mathbf{S}·{\mathbf{\sigma }}_{\mathbf{b}}\right)·\mathrm{d}\mathbf{b}+\left(\mathbf{S}·{\mathbf{\sigma }}_{\mathbf{c}}\right)·\mathrm{d}\mathbf{c},\end{array}$$

(5)

where ${\mathbf{\sigma }}_{\mathbf{a}}=\mathbf{b}\times \mathbf{c}$, ${\mathbf{\sigma }}_{\mathbf{b}}=\mathbf{c}\times \mathbf{a}$, and ${\mathbf{\sigma }}_{\mathbf{c}}=\mathbf{a}\times \mathbf{b}$ are surface area vectors.

Then, according to the above “de facto `theorem ’", we have the following equation to determine the period vectors, which is also the EOS for crystals under general external stress and temperature:

$$\mathbf{S}·{\mathbf{\sigma }}_{\mathbf{h}}=-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial lnZ}{\partial \mathbf{h}}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).$$

(6)

2.3. Tuckerman’s Internal Stress

In 2010, Tuckerman introduced the (macroscopic) internal stress ${\mathbf{P}}^{\left(\mathrm{int}\right)}$ of crystals as Equation (5.6.9) in his book [6] (Equation (5.7.9) in the 2023 version):

$${\mathbf{P}}^{\left(\mathrm{int}\right)}={\displaystyle \frac{1}{\beta V}}\sum _{\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}}{\displaystyle \frac{\partial lnZ}{\partial \mathbf{h}}}\otimes \mathbf{h}.$$

(7)

Combining it with Eq. (6) yields

$$\mathbf{S}·{\mathbf{\sigma }}_{\mathbf{h}}=-{\mathbf{P}}^{\left(\mathrm{int}\right)}·{\mathbf{\sigma }}_{\mathbf{h}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}),$$

(8)

which means

$$\mathbf{S}+{\mathbf{P}}^{\left(\mathrm{int}\right)}=0.$$

(9)

So Eq. (6) is the Macroscopic Mechanical Equilibrium Condition (MMEC): the macroscopic internal stress balances the external stress.

2.4. Isotropic External Pressure

The isotropic external pressure P is a special case of the external stress $\mathbf{S}$. That is $\mathbf{S}=-P\mathbf{I}$, with $\mathbf{I}$ being the identity tensor. Then, Eq. (6) becomes

$$P{\mathbf{\sigma }}_{\mathbf{h}}=-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial lnZ}{\partial \mathbf{h}}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}),$$

(10)

which can be further written as

$$PV={\displaystyle \frac{1}{\beta }}\left({\displaystyle \frac{\partial lnZ}{\partial \mathbf{h}}}\right)·\mathbf{h}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).$$

(11)

Since the external pressure is the same in all directions, we may assume that the crystal cell expands uniformly when the external pressure changes. This means that all the directions of the period vectors are fixed and their magnitudes are proportional to $V^{1/3}$. In other words, the period vectors depend on the volume only by (implicitly) containing the factor $V^{1/3}$. Once the factor $V^{1/3}$ is removed, the period vectors are independent of the volume V:

$${\displaystyle \frac{\partial }{\partial V}}\left({\displaystyle \frac{\mathbf{h}}{V^{\frac{1}{3}}}}\right)=0(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).$$

(12)

As a result,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \mathbf{h}}{\partial V}}& ={\displaystyle \frac{\partial }{\partial V}}\left({\displaystyle \frac{\mathbf{h}}{V^{\frac{1}{3}}}}{V}^{{\displaystyle \frac{1}{3}}}\right)=\left({\displaystyle \frac{\mathbf{h}}{V^{\frac{1}{3}}}}\right){\displaystyle \frac{\partial V^{\frac{1}{3}}}{\partial V}}\hfill \\ & ={\displaystyle \frac{1}{3V}}\mathbf{h}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(13)

Then using Eq. (13) and Eq. (11),

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial lnZ}{\partial V}}& ={\displaystyle \frac{1}{\beta }}\sum _{\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}}\left({\displaystyle \frac{\partial lnZ}{\partial \mathbf{h}}}\right)·{\displaystyle \frac{\partial \mathbf{h}}{\partial V}}\hfill \\ & ={\displaystyle \frac{1}{3V\beta }}\sum _{\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}}\left({\displaystyle \frac{\partial lnZ}{\partial \mathbf{h}}}\right)·\mathbf{h}=P.\hfill \end{array}$$

(14)

It reproduces Eq. (4). Therefore, Eq. (6) covers the special case of the isotropic external pressure.

It is better to put Eqs. (14) and (11) together as a summary for the external isotropic pressure condition:

$$\begin{array}{cc}\hfill P& ={\displaystyle \frac{1}{\beta }}\left({\displaystyle \frac{\partial lnZ}{\partial V}}\right){\displaystyle \frac{V}{V}},\hfill \\ \hfill P& ={\displaystyle \frac{1}{\beta }}\left({\displaystyle \frac{\partial lnZ}{\partial \mathbf{a}}}\right)·{\displaystyle \frac{\mathbf{a}}{V}},\hfill \\ \hfill P& ={\displaystyle \frac{1}{\beta }}\left({\displaystyle \frac{\partial lnZ}{\partial \mathbf{b}}}\right)·{\displaystyle \frac{\mathbf{b}}{V}},\hfill \\ \hfill P& ={\displaystyle \frac{1}{\beta }}\left({\displaystyle \frac{\partial lnZ}{\partial \mathbf{c}}}\right)·{\displaystyle \frac{\mathbf{c}}{V}}.\hfill \end{array}$$

(15)

In fact, all these four equations must be satisfied at the same time, and they are satisfied as long as any one of them is satisfied and the crystal expands uniformly with the external isotropic pressure P.

3. Non-Crystal System Under General External Stress

3.1. The EOS

Assuming that the basic physical properties of a non-crystal system do not depend on its macroscopic shape, the system may be regarded as only one but huge crystal “cell" without periodicity in structures, and the “period vectors" are now only the edge vectors of the “cell", but we may still call them the period vectors.

Using the same symbols, Eq. (5) is the infinitesimal work done by the external stress $\mathbf{S}$ acting on the surfaces of the system, then we also get Eq. (6) as the EOS.

Supposing that Tuckerman’s the internal stress, as in Eq.(7) also applies, the MMEC as of Eq. (9) holds as well. Furthermore, since Eqs.(10 through 15) apply, the isotropic external pressure is covered.

3.2. The Detailed MMEC in Classical Physics

In classical physics, the partition function can be factorized as:

$$Z=Z_k{Z}_u,$$

(16)

where $Z_k$ and $Z_u$ are contributed by the kinetic energy $E_k$ and the potential energy $E_p$ of the system, respectively. As in Equations (3.45)–(3.47) in the statistical reference book [2], they can be further expressed as:

$$\begin{array}{ccc}\hfill Z_k& =& {\displaystyle \frac{V^N}{N!}}\int \int \dots \int {\displaystyle \frac{1}{h^{3N}}}{e}^{-\beta E_k\left(\mathbf{p}\right)}\mathrm{d}\mathbf{p},\hfill \end{array}$$

(17)

$$\begin{array}{ccc}\hfill Z_u& =& {\displaystyle \frac{1}{V^N}}{\int }_V{\int }_V\dots {\int }_V{e}^{-\beta E_p\left(\mathbf{R}\right)}\mathrm{d}\mathbf{R},\hfill \end{array}$$

(18)

where h is the Planck constant and N is the total number of particles in the system. The integration is over all particle momenta $\mathbf{p}=\{{\mathbf{p}}_1,{\mathbf{p}}_2,\dots ,{\mathbf{p}}_N\}$ in Eq. (17), and over all particle positions $\mathbf{R}=\{{\mathbf{r}}_1,{\mathbf{r}}_2,\dots ,{\mathbf{r}}_N\}$ within the system volume for Eq. (18).

Then, Eq. (6) becomes

$$\mathbf{S}·{\mathbf{\sigma }}_{\mathbf{h}}=-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ln{Z}_k}{\partial \mathbf{h}}}-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ln{Z}_u}{\partial \mathbf{h}}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).$$

(19)

3.2.1. The Thermal Pressure

Since the integration in Eq. (17) has nothing to do with the period vectors, the first term on the right side of Eq. (19), the derivative of the partition function contributed by the kinetic energy, can be written as:

$$\begin{array}{cc}\hfill -{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ln{Z}_k}{\partial \mathbf{h}}}& =-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ln\left(V^N\right)}{\partial \mathbf{h}}}=-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{N}{V}}{\displaystyle \frac{\partial V}{\partial \mathbf{h}}}\hfill \\ & =-{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{N}{V}}{\mathbf{\sigma }}_{\mathbf{h}}=-{\displaystyle \frac{N}{V}}kT{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ & =-P_{thermal}{\mathbf{\sigma }}_{\mathbf{h}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(20)

This is the thermal pressure $P_{thermal}={\displaystyle \frac{N}{V}}kT$ on the cell surfaces, reflecting the forces in the particle collisions and/or the “force associated with the transport of momentum" introduced in our earlier work [7].

3.2.2. The Derivative of Potential Contribution with Respect to the Period Vectors

Now, let us take the derivative of the partition function contributed by the potential energy, the second term on the right side of Eq. (19), by using of Eq. (18):

$$\begin{array}{cc}& -{\displaystyle \frac{1}{\beta }}{\displaystyle \frac{\partial ln{Z}_u}{\partial \mathbf{h}}}=-{\displaystyle \frac{1}{\beta Z_u}}{\displaystyle \frac{\partial }{\partial \mathbf{h}}}{Z}_u(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(21)

Here, the derivative is with respect to the period vectors, which determine the volume V of the “cell", also the region of the integrals in Eq. (18). A careful inspection shows that Eq. (18) is the average of the function $e^{-\beta E_p\left(\mathbf{R}\right)}$, which depends on the particle positions $\mathbf{R}$, evenly and fully sampled over the “cell" space Vs. Then, as long as the points $\mathbf{R}$ are sampled exactly the same way, the calculated average is the same, even if the integrals are changed in some form. In any case, we will try to derive rigorously.

Now let us expand the position vector of each particle in the “cell", with respect to the period vectors as:

$${\mathbf{r}}_i=r_{i,\mathbf{a}}\mathbf{a}+r_{i,\mathbf{b}}\mathbf{b}+r_{i,\mathbf{c}}\mathbf{c}(i=1,\dots ,N),$$

(22)

where the scaled coordinates $r_{i,\mathbf{h}}$, in the range of $[0,1)$, can be calculated as

$$r_{i,\mathbf{h}}={\displaystyle \frac{1}{V}}{\mathbf{r}}_i·{\mathbf{\sigma }}_{\mathbf{h}}(i=1,\dots ,N;\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).$$

(23)

For simplicity, let us consider a general position vector $\mathbf{r}$ in the “cell", also expanded as $\mathbf{r}=r_{\mathbf{a}}\mathbf{a}+r_{\mathbf{b}}\mathbf{b}+r_{\mathbf{c}}\mathbf{c}$ where $r_{\mathbf{h}}={\displaystyle \frac{1}{V}}\mathbf{r}·{\mathbf{\sigma }}_{\mathbf{h}}$ with $\mathbf{h}=\mathbf{a},\mathbf{b},$ or $\mathbf{c}$. Imagine to evenly cut the “cell" into n blocks along period vector $\mathbf{a}$, m blocks along $\mathbf{b}$, and l blocks along $\mathbf{c}$. By the integral definition, for any function $f\left(\mathbf{r}\right)$, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{V}}{\int }_{V}f\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}& ={\displaystyle \frac{1}{V}}\underset{\begin{array}{c}n\to \infty \\ m\to \infty \\ l\to \infty \end{array}}{lim}\sum _{i=1}^n\sum _{j=1}^m\sum _{k=1}^{l}f\left({\mathbf{r}}_{i,j,k}\right)\Delta v_{i,j,k}\hfill \\ & =\underset{\begin{array}{c}n\to \infty \\ m\to \infty \\ l\to \infty \end{array}}{lim}\sum _{i=1}^n\sum _{j=1}^m\sum _{k=1}^{l}f\left({\mathbf{r}}_{i,j,k}\right){\displaystyle \frac{\Delta v_{i,j,k}}{V}},\hfill \end{array}$$

(24)

where $\Delta v_{i,j,k}$ is the volume of the tiny intersection of the i-th block along $\mathbf{a}$, j-th block along $\mathbf{b}$, and k-th block along $\mathbf{c}$, and ${\mathbf{r}}_{i,j,k}$ is the center position of the intersection. In fact, under the limit of $n,m,$ and l all approaching infinity, $\Delta v_{i,j,k}$ is the term $\mathrm{d}\mathbf{r}$ of the integral of the left side expression of the above equation in mathematics. Meanwhile, the expression ${\displaystyle \frac{\Delta v_{i,j,k}}{V}}$ in the above equation is the term $\mathrm{d}{r}_{\mathbf{a}}\mathrm{d}{r}_{\mathbf{b}}\mathrm{d}{r}_{\mathbf{c}}$ in the three nested integrals all from 0 to 1 in the next equation, so we have it based on the above equation and the integral definition as:

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{V}}{\int }_{V}f\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}& ={\int }_0^1{\int }_0^1{\int }_0^{1}f\left(\mathbf{r}\right)\mathrm{d}{r}_{\mathbf{a}}\mathrm{d}{r}_{\mathbf{b}}\mathrm{d}{r}_{\mathbf{c}}\hfill \\ & ={\int }_0^1{\int }_0^1{\int }_0^{1}f\left(\mathbf{r}(r_{\mathbf{a}},r_{\mathbf{b}},r_{\mathbf{c}})\right)\mathrm{d}{r}_{\mathbf{a}}\mathrm{d}{r}_{\mathbf{b}}\mathrm{d}{r}_{\mathbf{c}}.\hfill \end{array}$$

(25)

Please note that the $f\left(\mathbf{r}\right)$ remains the same here as a function of $\mathbf{r}$. However, the position $\mathbf{r}$ is not an independent variable now, but a further function of $(r_{\mathbf{a}},r_{\mathbf{b}},r_{\mathbf{c}})$, as $\mathbf{r}=r_{\mathbf{a}}\mathbf{a}+r_{\mathbf{b}}\mathbf{b}+r_{\mathbf{c}}\mathbf{c}$.

Since Eq. (25) applies to every particles, employing Eq. (22) to map every scaled coordinates to the actual particle positions $\mathbf{R}$ in the system potential energy $E_p\left(\mathbf{R}\right)$, the samplings and the average of $e^{-\beta E_p\left(\mathbf{R}\right)}$ keep unchanged, if Eq. (18) is rewritten as:

$$\begin{array}{cc}\hfill Z_u& ={\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1\dots {\int }_0^1{\int }_0^1{\int }_0^1\hfill \\ & e^{-\beta E_p\left(\mathbf{R}\right)}\hfill \\ & \mathrm{d}{r}_{1,\mathbf{a}}\mathrm{d}{r}_{1,\mathbf{b}}\mathrm{d}{r}_{1,\mathbf{c}}\mathrm{d}{r}_{2,\mathbf{a}}\mathrm{d}{r}_{2,\mathbf{b}}\mathrm{d}{r}_{2,\mathbf{c}}\dots \mathrm{d}{r}_{N,\mathbf{a}}\mathrm{d}{r}_{N,\mathbf{b}}\mathrm{d}{r}_{N,\mathbf{c}}.\hfill \end{array}$$

(26)

Since the integrals are of fixed regions now, the derivative can go inside

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial Z_u}{\partial \mathbf{h}}}& ={\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1\dots {\int }_0^1{\int }_0^1{\int }_0^1\hfill \\ & {\displaystyle \frac{\partial }{\partial \mathbf{h}}}{e}^{-\beta E_p\left(\mathbf{R}\right)}\hfill \\ & \mathrm{d}{r}_{1,\mathbf{a}}\mathrm{d}{r}_{1,\mathbf{b}}\mathrm{d}{r}_{1,\mathbf{c}}\mathrm{d}{r}_{2,\mathbf{a}}\mathrm{d}{r}_{2,\mathbf{b}}\mathrm{d}{r}_{2,\mathbf{c}}\dots \mathrm{d}{r}_{N,\mathbf{a}}\mathrm{d}{r}_{N,\mathbf{b}}\mathrm{d}{r}_{N,\mathbf{c}}\hfill \end{array}$$

(27)

$$\begin{array}{cc}& ={\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1{\int }_0^1\dots {\int }_0^1{\int }_0^1{\int }_0^1\hfill \\ & -\beta e^{-\beta E_p\left(\mathbf{R}\right)}{\displaystyle \frac{\partial }{\partial \mathbf{h}}}{E}_p\left(\mathbf{R}\right)\hfill \\ & \mathrm{d}{r}_{1,\mathbf{a}}\mathrm{d}{r}_{1,\mathbf{b}}\mathrm{d}{r}_{1,\mathbf{c}}\mathrm{d}{r}_{2,\mathbf{a}}\mathrm{d}{r}_{2,\mathbf{b}}\mathrm{d}{r}_{2,\mathbf{c}}\dots \mathrm{d}{r}_{N,\mathbf{a}}\mathrm{d}{r}_{N,\mathbf{b}}\mathrm{d}{r}_{N,\mathbf{c}}\hfill \\ & (\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

Here the particle positions $\mathbf{R}$ are not independent variables, but expressed as in Eq. (22). In a Cartesian coordinate system, the components of the period vector $\mathbf{h}$ may be denoted as $\left(h_x,h_y,h_z\right)$. Employing Eq. (22), one has

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial {\mathbf{r}}_i}{\partial h_x}}& =\left(r_{i,\mathbf{h}},0,0\right),{\displaystyle \frac{\partial {\mathbf{r}}_i}{\partial h_y}}=\left(0,r_{i,\mathbf{h}},0\right),\hfill \\ \hfill {\displaystyle \frac{\partial {\mathbf{r}}_i}{\partial h_z}}& =\left(0,0,r_{i,\mathbf{h}}\right)(i=1,\dots ,N;\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}),\hfill \end{array}$$

(28)

and, for any d-direction,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial E_p\left(\mathbf{R}\right)}{\partial h_d}}& =\sum _{i=1}^N{\displaystyle \frac{\partial E_p\left(\mathbf{R}\right)}{\partial {\mathbf{r}}_i}}·{\displaystyle \frac{\partial {\mathbf{r}}_i}{\partial h_d}}\hfill \\ & (d=x,y,z;\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(29)

Let us denote the net force acting on particle i by all other particles as :

$$\begin{array}{cc}\hfill {\mathbf{F}}_i& =-{\displaystyle \frac{\partial E_p\left(\mathbf{R}\right)}{\partial {\mathbf{r}}_i}},\hfill \end{array}$$

(30)

with components form of $\left(F_{i,x},F_{i,y},F_{i,z}\right)$, and use Eqs. (23), (28 - 30) to get:

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial E_p\left(\mathbf{R}\right)}{\partial \mathbf{h}}}& =\left({\displaystyle \frac{\partial }{\partial h_x}},{\displaystyle \frac{\partial }{\partial h_y}},{\displaystyle \frac{\partial }{\partial h_z}}\right)E_p\left(\mathbf{R}\right)\hfill \\ & =\sum _{i=1}^N\left(-F_{i,x}{r}_{i,\mathbf{h}},-F_{i,y}{r}_{i,\mathbf{h}},-F_{i,z}{r}_{i,\mathbf{h}}\right)\hfill \\ & =\sum _{i=1}^N\left(-{\mathbf{F}}_i\right)r_{i,\mathbf{h}}=\sum _{i=1}^N\left(-{\mathbf{F}}_i\right){\displaystyle \frac{{\mathbf{r}}_i·{\mathbf{\sigma }}_{\mathbf{h}}}{V}}\hfill \\ & =-{\displaystyle \frac{1}{V}}\sum _{i=1}^N\left({\mathbf{F}}_i\otimes {\mathbf{r}}_i\right)·{\mathbf{\sigma }}_{\mathbf{h}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(31)

Since the derivatives are clear now, let us go back to the “cell" spaces for the integration and use Eqs. (27 and 31), to rewrite Eq. (21) as:

$$\begin{array}{cc}& -{\displaystyle \frac{\partial ln{Z}_u}{\beta \partial \mathbf{h}}}\hfill \\ & ={\displaystyle \frac{1}{Z_u{V}^N}}{\int }_V{\int }_V\dots {\int }_V{e}^{-\beta E_p\left(\mathbf{R}\right)}\hfill \\ & \left(-{\displaystyle \frac{1}{V}}\sum _{i=1}^N{\mathbf{F}}_i\otimes {\mathbf{r}}_i·{\mathbf{\sigma }}_{\mathbf{h}}\right)\mathrm{d}\mathbf{R}\hfill \\ & (\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(32)

This is the average of Eq. (31) for all possible positions of particles, with the weight of $e^{-\beta E_p\left(\mathbf{R}\right)}$.

3.2.3. The Pair-Interaction Contributions

If only consider internal pair-interactions, Eq. (31) can be rewritten as:

$$\begin{array}{cc}& {\displaystyle \frac{\partial E_p\left(\mathbf{R}\right)}{\partial \mathbf{h}}}=-{\displaystyle \frac{1}{V}}\sum _{i=1}^N\left({\mathbf{F}}_i\otimes {\mathbf{r}}_i\right)·{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ & =-{\displaystyle \frac{1}{V}}\sum _{i=1}^N\sum _{j\ne i}^N\left({\mathbf{f}}_{j\to i}\otimes {\mathbf{r}}_i\right)·{\mathbf{\sigma }}_{\mathbf{h}}(\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}),\hfill \end{array}$$

(33)

where ${\mathbf{f}}_{j\to i}$ is the force acting on particle i by particle j. Since ${\mathbf{f}}_{j\to i}+{\mathbf{f}}_{i\to j}=0$, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial E_p\left(\mathbf{R}\right)}{\partial \mathbf{h}}}& =-{\displaystyle \frac{1}{2V}}\sum _{i=1}^N\sum _{j\ne i}^N\left({\mathbf{f}}_{j\to i}\otimes {\mathbf{r}}_i+{\mathbf{f}}_{i\to j}\otimes {\mathbf{r}}_j\right)·{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ & =-{\displaystyle \frac{1}{2V}}\sum _{i=1}^N\sum _{j\ne i}^N\left({\mathbf{f}}_{j\to i}\otimes \left({\mathbf{r}}_i-{\mathbf{r}}_j\right)\right)·{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ & =-{\displaystyle \frac{1}{V}}\sum _{i=1}^{N-1}\sum _{j>i}^N\left({\mathbf{f}}_{i\to j}\otimes \left({\mathbf{r}}_j-{\mathbf{r}}_i\right)\right)·{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ & (\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}).\hfill \end{array}$$

(34)

Let us consider all geometric planes cutting through the “cell" but parallel to ${\mathbf{\sigma }}_{\mathbf{h}}$. The total number (amount) of the planes may be defined as $\mathbf{h}·{\mathbf{\sigma }}_{\mathbf{h}}=V$, as it is proportional to the length of the period vector $\mathbf{h}$. Then, between particles j and i, there are $\left|\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)·{\mathbf{\sigma }}_{\mathbf{h}}\right|\le V$ amount of such planes. Thus, the term in the above equation

$$\begin{array}{cc}& -{\displaystyle \frac{1}{V}}\left({\mathbf{f}}_{i\to j}\otimes \left({\mathbf{r}}_j-{\mathbf{r}}_i\right)\right)·{\mathbf{\sigma }}_{\mathbf{h}}={\displaystyle \frac{\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)·{\mathbf{\sigma }}_{\mathbf{h}}}{V}}{\mathbf{f}}_{i\to j}\hfill \\ & (\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}),\hfill \end{array}$$

(35)

means the force ${\mathbf{f}}_{i\to j}$ between the particles i and j and penetrating and only penetrating the planes between them is averaged across all these planes. Please note that the force ${\mathbf{f}}_{i\to j}$ is in the direction of $\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)$, but may be not in the period vector $\mathbf{h}$ direction.

3.2.4. The Many-Body Interaction Contributions

In actual calculations, three-, four-, and/or even more-body interactions [8,9,10] may also be employed, as in molecular dynamics simulations. By definition, these interactions mean that any of the forces depend on the positions of all the participating particles, then can not essentially be replaced with some pair-interactions in general. Newton’s Third Law may be interpreted as: within a given system, any part may exert force(s) on any other part(s), but cannot accelerate the system itself as a whole, then the net of all forces and the net of all torques within any group of participating particles must be zero. For example, for an M-body interaction ($M\ge 2$), we have

$$\begin{array}{c}\hfill \sum _{k=1}^M{\mathbf{F}}_k^{\left(M\right)}=0,\end{array}$$

(36)

where ${\mathbf{F}}_k^{\left(M\right)}$ is the M-body force on particle k by the rest $M-1$ particles.

If an internal M-body interaction is employed, all the corresponding M-body potentials of various groups of M particles may be included in the system potential $E_p\left(\mathbf{R}\right)$ of Eq. (18). Then all the forces of each group additionally contribute to Eq. (32), in the form of

$$\begin{array}{cc}& \left(-{\displaystyle \frac{1}{V}}\sum _{k=1}^M{\mathbf{F}}_k^{\left(M\right)}\otimes {\mathbf{r}}_k·{\mathbf{\sigma }}_{\mathbf{h}}\right)\hfill \\ \hfill =& -{\displaystyle \frac{1}{V}}\left({\mathbf{F}}_1^{\left(M\right)}\otimes {\mathbf{r}}_1+\sum _{k=2}^M{\mathbf{F}}_k^{\left(M\right)}\otimes {\mathbf{r}}_k\right)·{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ \hfill =& -{\displaystyle \frac{1}{V}}\left(-\sum _{k=2}^M{\mathbf{F}}_k^{\left(M\right)}\otimes {\mathbf{r}}_1+\sum _{k=2}^M{\mathbf{F}}_k^{\left(M\right)}\otimes {\mathbf{r}}_k\right)·{\mathbf{\sigma }}_{\mathbf{h}}\hfill \\ \hfill =& \left(-{\displaystyle \frac{1}{V}}\sum _{k=2}^M{\mathbf{F}}_k^{\left(M\right)}\otimes \left({\mathbf{r}}_k-{\mathbf{r}}_1\right)·{\mathbf{\sigma }}_{\mathbf{h}}\right)\hfill \\ & (\mathbf{h}=\mathbf{a},\mathbf{b},\mathbf{c}),\hfill \end{array}$$

(37)

where Eq. (36) was employed.

Compared with Eq. (34), Eq. (37) may mean an M-body interaction can be instantaneously regarded as $M-1$ pair-interactions here. In each of the $M-1$ pairs, the net of the two forces is zero. However, they may not point to each other. Although they may generate torque, the total torque of the M-body interaction is zero, as required by the Newton’s Third Law. Then, the M-body interactions contribute the forces over the previously mentioned geometric planes the same way as pair-interactions do.

3.2.5. Conclusions

According to Eqs. (20, 32, 34, 35, and 37 ), the right side of Eq. (19) is the averaged internal forces over all geometric planes cutting the system but parallel to the system surface $\mathbf{h}$ and over all possible particle position distributions, while the left side is the external force on the surface $\mathbf{h}$, then Eq. (19) means the MMEC.

4. Summary

As the EOS, Eq. (6) applies to both crystals and non-crystal systems, under isotropic external pressure or anisotropic external stress and any temperature. It is also the Macroscopic Mechanical Equilibrium Condition of the system.