1 Introduction

In recent years, predicting crystal structures has become increasingly important and feasible [1,2,3,4,5,6,7,8,9,10]. For this purpose, even before crystals are synthesized, equations are required to determine the discrete particle positions and the period vectors (cell edge vectors). Since particles (atoms, ions, electrons) inside the crystal obey either Newton’s second law or the Schrödinger equation, the only unknown is an explicit equation governing the period vectors, especially when the crystal is under general external stress and temperature. In 1980, by extending Andersen’s idea [11], Parrinello and Rahman proposed their theory regarding this fundamental physics problem in molecular dynamics simulations with the periodic boundary condition being applied on the system [12, 13]. Since then, much more effort has been devoted to it based on the Lagrangian, or the Hamiltonian, or the Newtonian dynamics [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]. However, since these works were carried out in the framework of classical physics, further modeling is needed to apply them to quantum systems. Alternatively, cell volumes and period vectors were determined by minimizing the Gibbs free energy for crystals under external pressure and temperature [32,33,34,35]. Here, we derive a new equation for the period vectors, starting with a rigorous derivation of the work done by arbitrary external stress on the “central” crystal cell. Since it is based on the principles of statistical physics, the new equation applies to both classical and quantum systems immediately. Furthermore, it is consistent and can be combined with previously derived corresponding one in classical Newtonian dynamics. The existing theory for crystals under external pressure and temperature is covered as a special case.

In 2010, Tuckerman presented an expression for the internal stress of crystals in his book on statistical molecular simulations [36], and applied it successfully in a solid silicon system with the mechanical equilibrium condition [37]. With this internal stress, our newly derived equation for the period vectors also turns out to be the mechanical equilibrium condition. Hence, the internal stress expression and our equation are mutually consistent.

Since it determines the period vectors for crystals under arbitrarily given external stress and temperature, the new equation is also the equation of state of the crystals. It may therefore be used to predict crystal structures and to study structural phase transitions and various crystal expansions induced by the external stress and/or temperature. It may play a key role in the study of piezoelectric and piezomagnetic materials, as the period vectors change due to external stress. For linear elastic crystals, it takes the microscopic form of the generalized Hooke’s law, but with a temperature dependence. In this article, we address all these aspects individually, followed by a one-dimensional simplified example and a summary.

2 Equation derivation for period vectors in a crystal under external stress and temperature

Consider a crystal under arbitrarily given external stress \({\mathbf {S}}\), a second-rank tensor (\(3\times 3\) matrix), and temperature T. From a microscopic point of view, the crystal consists of an unlimited periodic arrangement of identical cells in the three-dimensional space. We may study the system by focusing on a “center” cell interacting with all other cells. Equivalently, external forces that act on the surfaces of the macroscopic crystal bulk can be described as the action on the surfaces of the “center” cell. To derive the equation for the period (cell edge) vectors: \({\mathbf {h}}={\mathbf {a}}\), \({\mathbf {b}}\), \({\mathbf {c}}\), which form a right-handed system, we start with the formulation of the work done by the external stress on the “center” cell. In order to do it rigorously, we need to bear in mind that in the definition of work, the subject of the displacement and the object being acted on by the corresponding force must be exactly the same physical object.

Fig. 1
figure 1

If the positions of the particles in a crystal are kept fixed, all the cells can be collectively translated anywhere. The green frame is one example of the new cell placement, and the lime frame of dashed lines is another (though only one layer of cells is shown here). As shown in the lower part, the red particle in the bottom-left corner of the amplified original cell on the left, is in the top-right corner of the amplified new green-framed cell on the right, viewed through the front surfaces

Full size image

As shown in Fig. 1, while the positions of the particles inside a crystal are fixed, all the cells can be translated anywhere together. The green and the lime frames in the figure are two new placements. In physics, all the translations are equivalent to each other, as is their average. Alternatively, if the cells are kept fixed, the particles may be collectively translated anywhere. As a result, while all relative positions between the particles are kept constant, individual particles may be placed anywhere in the cell. An average of all these cell movements corresponds to a situation where the fixed cells are made of a continuous medium. Therefore, any cell is now an actual object with matter in it everywhere, instead of, as normally described, a vacuum region of three-dimensional space containing some discrete particles. In fact, any instance of all fixed cells with all fixed particle positions is one of the indistinguishable microscopic states, discussed in Section 4 of our previous work [31]. Furthermore, the cell surfaces are physical entities with an infinitesimal fraction of the fixed cell mass, rather than pure geometric planes. Let us denote the area vectors of these surfaces as \(\sigma _{{\mathbf {a}}}={\mathbf {b}}\times {\mathbf {c}}\), \(\sigma _{{\mathbf {b}}}={\mathbf {c}}\times {\mathbf {a}}\), and \(\sigma _{{\mathbf {c}}}={\mathbf {a}}\times {\mathbf {b}}\) with respect to the period vectors \({\mathbf {h}}={\mathbf {a}}\), \({\mathbf {b}}\), \({\mathbf {c}}\) respectively. Then, the force, described by the external stress \({\mathbf {S}}\), acting on the physical surface \(\sigma _{{\mathbf {h}}}\) of the cell, is a real force and expressed as \({\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}\). The displacement of the physical cell surface \(\sigma _{{\mathbf {h}}}\) is \(d{\mathbf {h}}\) (= \(d{\mathbf {a}}\), \(d{\mathbf {b}}\), \(d{\mathbf {c}}\)) as it represents a physical quantity. Then according to the definition, the work done by the external stress on the crystal “center” cell is

$$\begin{aligned} dW =({\mathbf {S}}\cdot \sigma _{{\mathbf {a}}}) \cdot d{\mathbf {a}}+ ({\mathbf {S}}\cdot \sigma _{{\mathbf {b}}}) \cdot d{\mathbf {b}}+ ({\mathbf {S}}\cdot \sigma _{{\mathbf {c}}}) \cdot d{\mathbf {c}}. \end{aligned}$$
(1)

Since, as shown in Eq. (1), the displacement \(d{\mathbf {h}}\) is the conjugate variable of the force \({\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}\) acting on the physical cell surface \(\sigma _{{\mathbf {h}}}\), based on the principle of statistical physicsFootnote 1, we arrive at

$$\begin{aligned} {\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}=-\frac{1}{\beta } \frac{\partial \ln Z}{\partial {\mathbf {h}}} \ \ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}), \end{aligned}$$
(2)

where \(\beta =1/(kT)\), and k and Z are the Boltzmann constant and the system partition function, respectively. The partition function is a function of the system energy and the temperature. Equation (2) determines the period vectors for a crystal under external stress \({\mathbf {S}}\) and temperature T if all other variables are known.

3 Existing statistical theory for external pressure as a special case

In statistical physics, the theory for crystals under external pressure has been established for a long time [38]. Let us re-derive the basic equations for this special case from the above formulae of general external stress to ensure their consistency.

For an external pressure P, the stress reduces to: \({\mathbf {S}}= -P {\mathbf {I}}\), where \({\mathbf {I}}\) is the identity tensor. Then, Eq. (1) becomes

$$\begin{aligned} dW =-(P\sigma _{{\mathbf {a}}}) \cdot d{\mathbf {a}}- (P\sigma _{{\mathbf {b}}}) \cdot d{\mathbf {b}}- (P\sigma _{{\mathbf {c}}}) \cdot d{\mathbf {c}}. \end{aligned}$$
(3)

Since \(d{\mathbf {h}}\) is the conjugate variable of the force \(-P\sigma _{{\mathbf {h}}}\) acting on the cell surface \(\sigma _{{\mathbf {h}}}\), in view of the above principle of statistical physics, we have

$$\begin{aligned} P \sigma _{{\mathbf {h}}}=\frac{1}{\beta } \frac{\partial \ln Z}{\partial {\mathbf {h}}} \ \ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(4)

The right-hand side of Eq. (4) has the same form as in Eq. (2), except for the sign. However, this equation is for crystals under external pressure, while Eq. (2) is for external stress.

The cell volume \(V=({\mathbf {a}}\times {\mathbf {b}}) \cdot {\mathbf {c}}\), so it follows that

$$\begin{aligned} dV =\sigma _{{\mathbf {a}}}\cdot d{\mathbf {a}}+ \sigma _{{\mathbf {b}}}\cdot d{\mathbf {b}}+ \sigma _{{\mathbf {c}}}\cdot d{\mathbf {c}}. \end{aligned}$$
(5)

Therefore, Eq. (3) can be rewritten as in the familiar form

$$\begin{aligned} dW =-PdV. \end{aligned}$$
(6)

As the cell volume V is the conjugate variable of the pressure P, based on the principle of statistical physics, we have

$$\begin{aligned} P =\frac{1}{\beta } \frac{\partial \ln Z}{\partial V}. \end{aligned}$$
(7)

This is a standard equation, for instance, as in Equation (2.96) from the reference book [38]. In essence, it is the equation of state of crystals under external pressure at equilibrium. For a given external pressure and temperature, the crystal cell volume can be calculated from this equation. This was extensively studied and presented in the book by Anderson [39]. However, Eq. (7) is not explicitly about the period vectors. Combining Eqs. (7) and (4) yields:

$$\begin{aligned} P V= & {} \frac{V}{\beta } \frac{\partial \ln Z}{\partial V} =\frac{{\mathbf {a}}}{\beta } \cdot \frac{\partial \ln Z}{\partial {\mathbf {a}}} =\frac{{\mathbf {b}}}{\beta } \cdot \frac{\partial \ln Z}{\partial {\mathbf {b}}} =\frac{{\mathbf {c}}}{\beta } \cdot \frac{\partial \ln Z}{\partial {\mathbf {c}}}. \end{aligned}$$
(8)

This means that the crystal cell shape must keep certain symmetries under external pressure. Under such circumstances, Eq. (4) for any specific period vector, e.g. \({\mathbf {h}}={\mathbf {b}}\), must be equivalent to Eq. (7). However, since neither Eq. (4) nor Eq. (7) applies, there is no similar cell shape restriction as shown in Eq. (8), on crystals under general external stress.

Then, as expected, Eq. (1) reduces to Eq. (3) and further to Eq. (6), whereas Eq. (2) reduces to Eq. (4) and further to Eq. (7) for the situation of external pressure.

4 Internal stress and mechanical equilibrium condition

In Equation (5.6.9) of his 2010 book [36], Tuckerman expresses the crystal internal stress, also a second-rank tensor. This expression can also be written in the following form

$$\begin{aligned} {\mathbf {P}}^{\text {(int)}}= \frac{1}{\beta V} \sum _{{\mathbf {h}}={\mathbf {a}},{\mathbf {b}},{\mathbf {c}}} \frac{\partial \ln Z}{\partial {\mathbf {h}}} \otimes {\mathbf {h}}. \end{aligned}$$
(9)

Let us perform three dot products on the right side of Eq. (9) with the cell surface area vectors (\(\sigma _{{\mathbf {a}}}\), \(\sigma _{{\mathbf {b}}}\), \(\sigma _{{\mathbf {c}}}\)) separately and apply the relationship \( {\mathbf {h}}\cdot \sigma _{{\mathbf {x}}}= V \delta _{{\mathbf {h}},{\mathbf {x}}}\), to obtain

$$\begin{aligned} {\mathbf {P}}^{\text {(int)}}\cdot \sigma _{{\mathbf {h}}}= \frac{1}{\beta } \frac{\partial \ln Z}{\partial {\mathbf {h}}} \ \ \ ({\mathbf {h}}={\mathbf {a}},{\mathbf {b}},{\mathbf {c}}). \end{aligned}$$
(10)

Substituting Eq. (10) into Eq. (2), we have

$$\begin{aligned} {\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}= - {\mathbf {P}}^{\text {(int)}}\cdot \sigma _{{\mathbf {h}}}\ \ \ ({\mathbf {h}}={\mathbf {a}},{\mathbf {b}},{\mathbf {c}}), \end{aligned}$$
(11)

which is equivalent to

$$\begin{aligned} {\mathbf {S}}+ {\mathbf {P}}^{\text {(int)}}= 0 , \end{aligned}$$
(12)

as the three cell surface area vectors are not in the same plane.

Since Eq. (12) means that the internal and the external stresses balance each other, Eq. (2) is the equilibrium condition, where the internal stress is expressed by Eq. (9), for the crystal under external stress and temperature, from mechanical point of view. If Eq. (2) cannot be satisfied, the system disintegrates. For example, if the temperature is increased to a point where Eq. (2) cannot hold, the crystal melts. It also means that Eq. (2) and the internal stress expression are mutually consistent.

5 Equation of state and system expansion

In fact, Eq. (2) involves the period vectors, the position vectors of all the particles in the “center” cell, the external stress, and the temperature. To solve the whole crystal system, Eq. (2) has to be combined with Newton’s second law or the Schrödinger equation for the motion of all particles in the cell. All particles’ motion can be solved by Newton’s second law and/or the Schrödinger equation, for a set of given period vectors as parameters. Equation (2) can then be used to solve the period vectors for given external stress, temperature, and the solved particle motions. This process should be repeated until all variables converge. Since Newton’s second law and the Schrödinger equation apply generally, we may regard this process as a means to solve the period vectors by Eq. (2) for external conditions. Therefore, Eq. (2) is the equation of state for crystals under external stress and temperature in terms of the period vectors, instead of the cell volume.

As expressed by the period vectors, various crystal expansions caused by the change of the external stress and/or temperature may be calculated by solving Eq. (2) under all external conditions. Alternatively, for instance, if the temperature is fixed but the external stress is changed, we may take partial derivatives of Eq. (2) with respect to the changed stress components and solve the isothermal expansion in terms of the external stress.

As another example, let us consider the “isobaric” thermal expansion, in which the external stress is fixed but the temperature is changed. A particle’s motion inside a crystal is usually separated into the motion of equilibrium position and harmonic oscillation around it. The harmonic oscillation may be represented by its frequency \(\omega \) (or frequencies), which is a function of the period vectors and all the particle equilibrium position vectors in the cell. For simplicity, suppose that there is only one particle in the cell, then its equilibrium position vector may be set to zero without loss of generality. The partition function Z can be considered as a function of the frequency \(\omega \) and the temperature. When the temperature is changed, the right-hand side of Eq. (2), \(-kT {\partial \ln Z(\omega , T)}/{\partial {\mathbf {h}}}\), changes. This causes a change in \({\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}\). Since the external stress \({\mathbf {S}}\) is fixed, \(\sigma _{{\mathbf {h}}}\) changes, which adjusts the period vectors. This means that the harmonic oscillation causes crystal isobaric thermal expansion.

As mentioned above, for the special case of external pressure, the cell shape should keep certain symmetries, then the inside particle equilibrium positions should also show the symmetries with respect to the cell. If no structural phase transition happens, the crystal can only expand or contract uniformly. This means that the particle equilibrium positions and the period vectors must change proportionally. Therefore, the particle equilibrium positions in terms of the period vectors must stay the same. Specifically, the directions of the period vectors may be set fixed, and their lengths be set proportional to \(V^{1/3}\). Furthermore, all particle equilibrium positions should be functions of the period vectors, and therefore of V. As a result, Eq. (2) only involves the crystal cell volume, the external pressure, and the temperature variables, as in the familiar form of the state equation.

6 Work and energy in quasi-equilibrium process

Let us consider an initial instant when the system is in an equilibrium state with the period vectors \({\mathbf {a}}_i\), \({\mathbf {b}}_i\), \({\mathbf {c}}_i\), under given external stress and temperature. If the external conditions vary only very slowly, then the system changes accordingly and reaches a new equilibrium final state with period vectors \({\mathbf {a}}_f\), \({\mathbf {b}}_f\), \({\mathbf {c}}_f\). If for each moment of this process the system is in an equilibrium state satisfying Eq. (2), the whole process is a quasi-equilibrium one. Although almost all variables (external stress, temperature, period vectors, cell surface vectors, and all particle position vectors) change, we get the work done by the external stress for the whole process by integration of Eq. (1) from the initial to the final state:

$$\begin{aligned} W =\int _{{\mathbf {a}}_i}^{{\mathbf {a}}_f} ({\mathbf {S}}\cdot \sigma _{{\mathbf {a}}}) \cdot d{\mathbf {a}}+ \int _{{\mathbf {b}}_i}^{{\mathbf {b}}_f} ({\mathbf {S}}\cdot \sigma _{{\mathbf {b}}}) \cdot d{\mathbf {b}}+ \int _{{\mathbf {c}}_i}^{{\mathbf {c}}_f} ({\mathbf {S}}\cdot \sigma _{{\mathbf {c}}}) \cdot d{\mathbf {c}}. \end{aligned}$$
(13)

The internal energy of the system is the total energy of the center cell \({E_t}\). Further denoting the internal energies of the initial and the final states as \(E_{t,i}\) and \(E_{t,f}\), respectively, the total heat absorbed by the system in the process can be expressed as

$$\begin{aligned} Q= E_{t,f} - E_{t,i} - W. \end{aligned}$$
(14)

7 Discussion

7.1 Piezoelectricity and piezomagnetism

The fascinating piezoelectric and piezomagnetic phenomena are induced when the external pressure on a crystal in one direction is changed while those in the other directions are kept constant. Since the external stress is changed, the period vectors also change according to the equation of state (2). Therefore, employing Eq. (2) to obtain the changed crystal structure accurately will prove helpful for the study of these phenomena.

7.2 Comparison with elasticity theory and generalized Hooke’s law

Since the theory of elasticity [40] also deals with the action of external stress on crystals, the work here should be compared with it.

In the first place, the elasticity theory mainly studies macroscopic properties of crystals by regarding them as a continuous medium, and therefore does not involve crystal period vectors explicitly. This work treats cells as a continuous medium only when formulating the work done by the external stress on the crystal cell, and derives the equation for determining the period vectors as the sole objective. Let us suppose elasticity theory also studies microscopic structures of crystals in the following.

Secondly, the elasticity theory normally employs a reference state that is not being acted on by any external forces, for which all the particle positions and the period vectors should be assumed to be known. This work determines the state under external stress without relying on other states.

Thirdly, the elasticity theory introduces and uses the concept strain to describe crystal deformation caused by external stress with respect to the reference state. Structural phase transitions may not be easily or straightforwardly described. This work does not use the strain concept, but regards all the period vectors and particle positions as independent variables, and as a result, straightforwardly describes any new crystal structures created by the change of the external stress and/or temperature.

Finally, the elasticity theory employs the generalized Hooke’s law as a principle, in which the stress and strain are assumed linearly related by elastic constants, which are usually determined by experiment. This work does not assume any analytical relationship between period vectors and external stress. Actually, Eq. (2) is the temperature-dependent relationship between them, but with no additional constants introduced, supposing everything in Eq. (2) can be at least calculated numerically. Then, Eq. (2) is the microscopic form of the generalized Hooke’s law in a situation where the period vectors and the external stress approximately change proportionally in reality. The elastic constants in the generalized Hooke’s law can be calculated by solving Eq. (2) for a series of values of the external stress, for a given temperature.

7.3 Consistency with previous work in classical physics

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

$$\begin{aligned} Z=Z_k Z_u, \end{aligned}$$
(15)

where \(Z_k\) and \(Z_u\) are contributed from the particles’ kinetic energy \(E_k\) and the cell potential energy \(E_p\), respectively. As in Equations (3.45)–(3.47) in the statistical reference book [38], they may be expressed as

$$\begin{aligned} Z_k= & {} \frac{V^N}{N!} \int \int \cdots \int \frac{1}{h^{3N}} e^{-\beta E_k({{\mathbf {p}}})} d {{\mathbf {p}}}, \end{aligned}$$
(16)
$$\begin{aligned} Z_u= & {} \frac{1}{V^N} \int _V \int _V \cdots \int _V e^{-\beta E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})} d {{\mathbf {R}}}, \end{aligned}$$
(17)

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

Since the integration in Eq. (16) has nothing to do with the period vectors, the derivative is straightforward:

$$\begin{aligned} \frac{\partial \ln Z_k}{\partial {\mathbf {h}}} = \frac{N}{V} \frac{\partial V}{\partial {\mathbf {h}}} = \frac{N}{V} \sigma _{{\mathbf {h}}}\ \ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(18)

However, Eq. (17) is a little more complicated, as the cell potential energy \(E_p\) depends on all period vectors and particle positions \({\mathbf {R}}\). For each particle in the “center” cell, the position vector may be expanded with respect to the period vectors as:

$$\begin{aligned} {\mathbf {r}}_i=r_{i,{\mathbf {a}}} {\mathbf {a}}+ r_{i,{\mathbf {b}}} {\mathbf {b}}+ r_{i,{\mathbf {c}}} {\mathbf {c}}\ \ (i=1, \cdots , N), \end{aligned}$$
(19)

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

$$\begin{aligned} r_{i,{\mathbf {h}}} = \frac{1}{V} {\mathbf {r}}_i \cdot \sigma _{{\mathbf {h}}}\ \ \ \ (i=1, \cdots , N; \ {\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(20)

The following derivation is similar to the one in page 135 of the reference book [38]. Let us first carry out the integration over \({\mathbf {r}}_i\) in the cell, per unit volume:

$$\begin{aligned} \frac{1}{V} \int _V \cdots d {{\mathbf {r}}_i}= & {} \int _V \cdots \frac{d {{\mathbf {r}}_i}}{V} = \int _0^1 \int _0^1 \int _0^1 \cdots d r_{i,{\mathbf {a}}} d r_{i,{\mathbf {b}}} d r_{i,{\mathbf {c}}} \quad (i=1, \cdots , N). \end{aligned}$$
(21)

Then using Eq. (17), the derivatives of \(Z_u\) with respect to the period vectors are just the ones of the cell potential \(E_p\) inside the integration:

$$\begin{aligned} \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 \cdots \int _0^1 \int _0^1 \int _0^1 \frac{\partial }{\partial {\mathbf {h}}} e^{-\beta E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})} \nonumber \\&d r_{1,{\mathbf {a}}} d r_{1,{\mathbf {b}}} d r_{1,{\mathbf {c}}} d r_{2,{\mathbf {a}}} d r_{2,{\mathbf {b}}} d r_{2,{\mathbf {c}}} \cdots d r_{N,{\mathbf {a}}} d r_{N,{\mathbf {b}}} d r_{N,{\mathbf {c}}} \quad ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(22)

We must now consider the dependency of \({\mathbf {R}}\) on the period vectors in Eq. (19). For this, the derivatives should be separated:

$$\begin{aligned} \frac{\partial E_p}{\partial {\mathbf {h}}}= & {} \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial {\mathbf {h}}} \Bigr |_{\mathbf {R}}+ \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial {\mathbf {h}}} \Bigr |_{{\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}} \quad ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(23)

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. (19), one has

$$\begin{aligned} \frac{\partial {\mathbf {r}}_i}{\partial h_x}= & {} \left( r_{i,{\mathbf {h}}}, 0,0 \right) , \ \ \frac{\partial {\mathbf {r}}_i}{\partial h_y} = \left( 0, r_{i,{\mathbf {h}}}, 0 \right) , \ \ \frac{\partial {\mathbf {r}}_i}{\partial h_z} = \left( 0,0, r_{i,{\mathbf {h}}} \right) \nonumber \\&\qquad \qquad (i=1, \cdots , N; \ {\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}), \end{aligned}$$
(24)

and, for any d-direction,

$$\begin{aligned} \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial h_d} \Bigr |_{{\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}} = \sum _{i=1}^N \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial {\mathbf {r}}_i} \cdot \frac{\partial {\mathbf {r}}_i}{\partial h_d} \ \ \ \ (d=x,y,z; \ {\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(25)

We denote the net force acting on particle i from all other particles in any cell as \({\mathbf {F}}_i = - \partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}}) / {\partial {\mathbf {r}}_i}\) and use Eqs. (20), (24), and (25) to obtain the last term of Eq. (23):

$$\begin{aligned} \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial {\mathbf {h}}} \Bigr |_{{\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}}= & {} \left( \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial h_x}, \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial h_y}, \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial h_z} \right) \Bigr |_{{\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}} \nonumber \\= & {} \sum _{i=1}^N - {\mathbf {F}}_i r_{i,{\mathbf {h}}} = \sum _{i=1}^N - {\mathbf {F}}_i \frac{{\mathbf {r}}_i \cdot \sigma _{{\mathbf {h}}}}{V} = - \frac{1}{V} \sum _{i=1}^N \left( {\mathbf {F}}_i \otimes {\mathbf {r}}_i \right) \cdot \sigma _{{\mathbf {h}}}\ \ \ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \nonumber \\ \end{aligned}$$
(26)

Combining Eqs. (22), (23), and (26) yields

$$\begin{aligned} - \frac{\partial \ln Z_u}{\beta \partial {\mathbf {h}}}= & {} - \frac{1}{\beta Z_u} \frac{\partial Z_u}{\partial {\mathbf {h}}} = - \frac{1}{\beta Z_u V^N} \int _V \int _V \cdots \int _V \frac{\partial }{\partial {\mathbf {h}}} e^{-\beta E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})} d {{\mathbf {R}}} \nonumber \\= & {} \frac{1}{Z_u V^N} \int _V \int _V \cdots \int _V e^{-\beta E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})} \nonumber \\&\left( \frac{\partial E_p({\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}, {{\mathbf {R}}})}{\partial {\mathbf {h}}} \Bigr |_{\mathbf {R}}- \frac{1}{V} \sum _{i=1}^N {\mathbf {F}}_i \otimes {\mathbf {r}}_i \cdot \sigma _{{\mathbf {h}}}\right) d {{\mathbf {R}}} \quad ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(27)

Substituting Eqs. (18) and (27) into Eq. (2), we arrive at

$$\begin{aligned} {\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}=- \frac{1}{V}NkT \sigma _{{\mathbf {h}}}- \frac{1}{\beta } \frac{\partial \ln Z_u}{\partial {\mathbf {h}}} \ \ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}). \end{aligned}$$
(28)

For “crystals” only containing ideal gases under external pressure, Eq. (28) becomes the ideal gas law: \(PV=NkT\), which is also the equilibrium condition from mechanical point of view.

Now let us consider equilibrium states where

$$\begin{aligned} {\ddot{{\mathbf {h}}}}= & {} 0 \ \ \ \ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}), \end{aligned}$$
(29)
$$\begin{aligned} \ddot{\mathbf {r}}_i= & {} 0 \ \ \ \ \ (i=1, \cdots , N). \end{aligned}$$
(30)

Then, the previously derived equation for the period vectors based on the Newtonian dynamics: Eq. (27) in paper [31] becomes

$$\begin{aligned} {\mathbf {S}}\cdot \sigma _{{\mathbf {h}}}= & {} - \frac{2}{3V}\sum _{i=1}^N \frac{1}{2} m_i\left| {\dot{{\mathbf {r}}}}_i\right| ^2 \sigma _{{\mathbf {h}}}+ \frac{\partial E_{p}}{\partial {\mathbf {h}}} \Bigr |_{\mathbf {R}}- \frac{1}{V} \sum _{i=1}^N \left( {\mathbf {F}}_i \otimes {\mathbf {r}}_i \right) \cdot \sigma _{{\mathbf {h}}}\quad ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}), \end{aligned}$$
(31)

where Eqs. (9), (17), (19), (25), and (26) in the paper [31] were employed, with \(m_i\) being the mass of particle i in the cell. The last term \(-kT \partial \ln Z_u / \partial {\mathbf {h}}\) in Eq. (28) as shown in Eq. (27) is the averaged form of the last two terms \( \left( \partial E_{p}/ \partial {\mathbf {h}}\right) |_{\mathbf {R}}\) and \( - \sum _{i=1}^N {\mathbf {F}}_i \otimes {\mathbf {r}}_i \cdot \sigma _{{\mathbf {h}}}/V\) of Eq. (31). They essentially share the same physics.

Since the kinetic energy terms in the two equations are normally also regarded as the same, Eq. (28) derived here through statistical physics, and Eq. (31) derived based on the Newtonian dynamics previously [31] are mutually consistent in classical physics.

7.4 Combination with previous work

The previously derived Eq. (27) in paper [31] for the period vectors based on the Newtonian dynamics is in the form of dynamical equation, then can be used as an “algorithm” for solving the equilibrium states. Let us replace its internal stress with Tuckerman’s expression, Eq. (9), to get

$$\begin{aligned} \alpha _{{\mathbf {h}},{\mathbf {h}}} {\ddot{{\mathbf {h}}}}=\left( {\mathbf {S}}+ {\mathbf {P}}^{\text {(int)}}\right) \cdot \sigma _{{\mathbf {h}}}\ \ ({\mathbf {h}}={\mathbf {a}}, {\mathbf {b}}, {\mathbf {c}}), \end{aligned}$$
(32)

where \(\alpha _{{\mathbf {h}},{\mathbf {h}}}\) is an equivalent mass, which can be chosen to be any positive real number. Eq. (32) can be used as an “algorithm” to iteratively solve Eq. (2) in both classical physics and quantum physics. It tells us how the period vectors should be changed if the internal and the external stresses do not balance each other.

In the Born-Oppenheimer approximation, ions and electrons in crystals are usually solved separately by Newtonian dynamics and quantum mechanics. The partition function of the whole system may then be factorized accordingly, meaning that Eq. (2) is a combination of classical and quantum physics.

8 A one-dimensional simplified model system

In 2010, Ma and Tuckerman applied Eq. (2) (Equations (1) and (2) in their paper [37]) in a solid silicon system with a success. However, since it is based on the Density Functional Theory, such kind of calculation cannot be done quickly from scratch, even if more detailed algorithms for the internal stress are available [41]. Instead, we present a much simpler one-dimensional example from classical physics here.

8.1 The model

We assume our “crystal” to consist of particles in a straight line. The period vector is just a single scalar variable, denoted as \(a>0\), which is also the cell volume, as the corresponding cell surface area vector is assumed to be a unit vector along the line. Let \(f_{ext}\) be the external scalar force, described by the external stress and acting on this “surface” along the line. If \(f_{ext}>0\), it acts outwards, i.e. it is a pulling force; if \(f_{ext}<0\), it acts inwards, i.e. a pressing force.

Let us further assume that there is only one atom in each cell. Considering the periodicity of the system, the net force on any atom from all the others is always zero, so that the atom in the center cell may be placed at the coordinate origin.

For simplicity, we model the kinetic energy of the center cell (i.e. the atom) as

$$\begin{aligned} E_k=\frac{3}{2}kT . \end{aligned}$$
(33)

As the interaction potential between any pair of atoms, we use the Lennard–Jones (L–J) 12-6 potential:

$$\begin{aligned} \varphi ^{(L-J)} \left( r \right) = 4 \epsilon \left[ { \left( \frac{\lambda }{r} \right) }^{12} - { \left( \frac{\lambda }{r} \right) }^{6} \right] , \end{aligned}$$
(34)

where r is the distance between the atoms. The parameters we chose here are \(\epsilon =3.500 \times 10^{-20}\hbox { J}\), and \(\lambda = 2.800 \, \AA \). The force exerted on an atom at r by the atom at the origin is

$$\begin{aligned} f^{(L-J)} \left( r \right)= & {} -\frac{d}{d r} \varphi ^{(L-J)} \left( r \right) = \frac{4 \epsilon }{r} \left[ { 12 \left( \frac{\lambda }{r} \right) }^{12} - { 6 \left( \frac{\lambda }{r} \right) }^{6} \right] . \end{aligned}$$
(35)

Based on Eq. (34), since only one atom is in each cell, the total cell potential is half of the total potential between this cell and the others, as each belongs to both of the two cells:

$$\begin{aligned} E_p^{(L-J)} \left( a \right)= & {} \frac{1}{2} \sum _{j=-\infty } ^{\infty (j \ne 0)} \ \varphi ^{(L-J)} \left( j a \right) = \sum _{j=1} ^\infty \ {4 \epsilon } \left[ { \left( \frac{\lambda }{j a} \right) }^{12} - { \left( \frac{\lambda }{j a} \right) }^{6} \right] . \end{aligned}$$
(36)

Therefore, the total energy of the center cell is

$$\begin{aligned} E_t= E_k + E_p^ {(L-J)} \left( a \right) . \end{aligned}$$
(37)

Substituting Eqs. (33) and (36) into Eq. (31) yields the equation for the period:

$$\begin{aligned} f_{ext} + \frac{kT}{a} + F_{L\rightarrow R} \left( a \right) = 0 , \end{aligned}$$
(38)

where

$$\begin{aligned} F_{L\rightarrow R} \left( a \right)= & {} - \frac{d}{da} E_p^{(L-J)} \left( a \right) = \frac{1}{2} \sum _{ j=-\infty } ^{\infty (j \ne 0)} \ j f^{(L-J)} \left( j a \right) \nonumber \\= & {} \sum _{ j=1} ^\infty \ \frac{4 \epsilon }{a} \left[ { 12 \left( \frac{\lambda }{ j a} \right) }^{12} - { 6 \left( \frac{\lambda }{ j a} \right) }^{6} \right] , \end{aligned}$$
(39)

which is the net force applied on the right half of the system by the left half, wherever the system is bisected.

Fig. 2
figure 2

The red line is the force acting on an atom by the atom in the center cell, as in Eq. (35); the blue line is the net force acting on the right half of the system by the left half when the system is bisected, as in Eq. (39). The two lines coincide almost completely as the data are very close. The inside figure shows a magnification near the minimum

Full size image

Figure 2 shows the force \(f^{(L-J)} \left( r \right) \) as a function of inter-atomic distance r as in Eq. (35) (red line), and the net force \(F_{L\rightarrow R} \left( a \right) \) as a function of system period a, as in Eq. (39) (blue line). As two atoms get close, they repel each other. This repulsion becomes stronger as the distance gets smaller, with a singularity at zero. When moving further apart, they attract each other. This attraction weakens and goes to zero asymptotically. The force between two neighbouring atoms is the first term of the total force \(F_{L\rightarrow R} \left( a \right) \), as shown in Eq. (39). The similarity between the two curves indicates that the rest of the terms in Eq. (39) decrease very rapidly.

Figure 2 also shows that the force \(F_{L\rightarrow R} \left( a \right) \) exhibits a minimum. Taking the derivative of the force with respect to the period, setting it to zero, and solving the equation

$$\begin{aligned}&\frac{d}{d a} F_{L\rightarrow R}\left( a \right) = \sum _{ j=1} ^\infty \ - \frac{4 \epsilon }{a^2} \left[ { 12 \times 13 \left( \frac{\lambda }{ j a} \right) }^{12} - { 6 \times 7 \left( \frac{\lambda }{ j a} \right) }^{6} \right] =0 , \end{aligned}$$
(40)

we obtain

$$\begin{aligned} a_{min}= \lambda \left( \frac{ \sum _{ j=1}^\infty 12 \times 13 j^{-12} }{ \sum _{ j=1}^\infty 6 \times 7 j^{-6} }\right) ^{\frac{1}{6}} =3.102 \AA . \end{aligned}$$
(41)

When \(a=a_{min}\), \(F_{L\rightarrow R}\left( a \right) \) has the unique minimum value of \(F_{L\rightarrow R}\left( a_{min} \right) = -5.769 \times 10^{-10}\hbox { N}\).

When the period a runs from zero to \(a_{min}\), the force \(F_{L\rightarrow R}\left( a \right) \) decreases monotonically, and when it runs from \(a_{min}\) to positive infinity, the force \(F_{L\rightarrow R}\left( a \right) \) increases monotonically approaching zero. As a matter of fact, the system is not stable in the region \([a_{min}, \infty ) \). For simplicity, let us consider zero temperature for the region first. Since \(F_{L\rightarrow R}<0\), Eq. (38) becomes

$$\begin{aligned} f_{ext}=-F_{L\rightarrow R}>0 , \end{aligned}$$
(42)

which means that the external force has a pulling effect. Suppose that the system is in an equilibrium state in the region. If any fluctuation factor results in \(f_{ext} > \left| F_{L\rightarrow R} \right| \) (even just a little bit), the atoms move apart, and the period increases (this may also be observed based on Eq. (32)). Thus, according to Fig. 2, \(\left| F_{L\rightarrow R} \right| \) decreases, causing an imbalance between \(f_{ext} \) and \( \left| F_{L\rightarrow R} \right| \). This further causes the atoms to move further apart, meaning the period increases. The value \(\left| F_{L\rightarrow R} \right| \) decreases more, causing greater imbalance between \(f_{ext} \) and \( \left| F_{L\rightarrow R} \right| \), etc., until the system breaks apart. If the fluctuation factor results in \(f_{ext} < \left| F_{L\rightarrow R} \right| \), a similar effect occurs; however, rather than breaking the system, the period becomes smaller and leaves the region \([a_{min}, \infty ) \). Then, an appropriate position in the region \((0, a_{min}] \) sets up a new equilibrium state. For nonzero temperature, as the term kT/a is always positive as in Eq. (38), it cannot make the system more stable in the region \([a_{min}, \infty )\). Therefore, let us restrict our calculation to the region \((0, a_{min}] \) in further computations.

Since constant volume is not of much interest, let us focus on constant temperature and constant external force.

8.2 Constant external temperature

Fig. 3
figure 3

System period a as a function of the external force \(f_{ext}\) at constant temperatures \(0\hbox {K}\) (red), \(1000\hbox {K}\) (blue), and \(5000\hbox {K}\) (green)

Full size image

Figure 3 shows changing system period as a function of the external force, for three constant external temperatures. While the pressing force may take any negative value, since the force \(F_{L\rightarrow R}(a)\) has the unique minimum at \(a_{min}\), the pulling force is limited to a maximum based on Eq. (38):

$$\begin{aligned} f_{ext,max,pull}= & {} - F_{L\rightarrow R} \left( a_{min} \right) - \frac{kT}{a_{min}}\nonumber \\&= 5.769 \times 10^{-10}\hbox { N} - \left( 4.452 \times 10^{-14}\hbox { N K}^{-1}\right) T. \end{aligned}$$
(43)

For external forces greater than this, the system disintegrates, i.e. the curves exhibit a singularity.

For external forces ranging from zero to \(f_{ext,e} = 1.0 \times 10^{-10}\hbox { N}\), let us suppose that the period and the external force change approximately linearly. The elastic constants of Hooke’s law, \(\tau \), may be calculated as in Table 1, where \(a_0\) and \(a_1\) are the periods corresponding to zero external force and \(f_{ext,e}\) respectively, and the change of the external force is \(\varDelta f_{ext} =f_{ext,e} - 0 =\) \(1.0 \times 10^{-10}\hbox { N}\).

Table 1 Elastic constants: \(\tau \)
Full size table
Fig. 4
figure 4

Work done on the system by the external force in the quasi-equilibrium process from the state with zero external force to the state where the period is \(a_{min}\), for three constant temperatures

Full size image

Then, the Hooke’s law can be expressed as

$$\begin{aligned} f_{ext} \approx \tau \frac{a- a_0}{a_0} . \end{aligned}$$
(44)

In fact, the relationship between the period and the external force may also be calculated by taking the derivative of Eq. (38) with respect to the period.

Fig. 5
figure 5

Heat absorbed by the system in the quasi-equilibrium process from the state with zero external force to the state where the period is \(a_{min}\), for three constant temperatures

Full size image

Employing Eq. (13) and choosing the state with the external force being zero as the initial one, and the state with the period being \(a_{min}\) as the final one, we may calculate the work done by the external force for the quasi-equilibrium process under the three constant temperatures. The results are shown in Fig. 4. As computed from Eqs. (14) and (37), Fig. 5 shows the heat absorbed by the system during this process. For zero temperature, no heat is needed by the system from outside (red curve). Consider an object moving in a potential field \(E_p^{(L-J)} \left( r \right) \), where r is the distance from the coordinate origin, as shown in Eq. (36). This object experiences a force of \(- dE_p^{(L-J)} \left( r \right) /dr\) from the field. According to the first equality in Eq. (39), that force is exactly the force \(F_{L\rightarrow R} \left( a \right) \) for \(r=a\). Furthermore, considering the force balance in Eq. (42), the quasi-equilibrium process in which the external force increases the system period is physically equivalent to the process of raising an object in a gravitational field with a force balancing the object’s weight from one stationary state to another. In the latter process, the energy needed for the potential change is exactly compensated by the work done by the force raising it, without any additional energy exchange. This is the situation for the zero-temperature curve in Fig. 5.

For nonzero constant temperatures, the system kinetic energy (Eq. (33)) is fixed, so that there is no change between the initial and the final states, and no heat is needed from the kinetic-energy perspective. Furthermore, Eq. (38) means that the difference between the external force and the force \(F_{L\rightarrow R} \left( a \right) \) is kT/a. In Fig. 5, the period a is in the range of (2.8, 3.2) so that kT/a is almost fixed as well. When an object is raised by an external force in a gravitational field of fixed difference from its weight, the additional energy needed to compensate the change in potential over the work done by the external force is proportional to the change of the object’s height/altitude. For the same reason, we see the lines in Fig. 5 are almost straight for positive temperatures. The heat may be exchanged between the system and an external heat reservoir through inelastic collisions between the atoms of the system and particles of the external reservoir around the system surfaces, assuming there is a temperature difference between them.

Fig. 6
figure 6

Under three constant external forces respectively, the period is changed by the temperature

Full size image

8.3 Constant external force

Figure 6 shows the period as a function of the temperature for fixed values for the external force. Based on Eq. (38), there is the upper limit

$$\begin{aligned} T_{max}= & {} - \frac{a_{min}}{k} \left( f_{ext} + F_{L\rightarrow R} \left( a_{min} \right) \right) \nonumber \\= & {} -2.247 \times 10^{13}\hbox { K N}^{-1} \left( f_{ext} -5.769 \times 10^{-10} \hbox { N }\right) , \end{aligned}$$
(45)

for the temperature. If the temperature is higher than \(T_{max}\), the system cannot exist, as Eq. (38) or Eq. (2) fails. This means \(T_{max}\) is the melting point. In Fig. 6, all lines for various given external forces end at \(T_{max}\). The figure also shows that the weaker the pulling (positive) force or the stronger the (negative) pressing force, the higher the melting point.

Fig. 7
figure 7

Relative coefficient of thermal expansion with temperature under the three given external forces

Full size image

Taking the derivative of Eq. (38) with respect to the period a under constant external force, we obtain the relative thermal expansion coefficient:

$$\begin{aligned} \zeta = \frac{1}{a} \frac{da}{dT} = \frac{k}{a^2} \left( \frac{kT}{a^2} - \frac{d}{da} F_{L\rightarrow R} \left( a \right) \right) ^{-1}, \end{aligned}$$
(46)

which is shown in Fig. 7 for the three given external forces. In Fig. 2, the force \(F_{L\rightarrow R}(a)\) curve is the flattest around the minimum in the region \((0, a_{min})\). Essentially, the temperature plays a very similar role as the force \(F_{L\rightarrow R}(a)\) plays in Eq. (38). Then, Fig. 7 presents the maximum values of the expansion coefficient, meaning the system is the most sensitive to temperature near the melting point.

The work done by the external force and the heat needed by the system for the quasi-equilibrium process under given external forces may also be calculated the same way as for the constant-temperature process. They are not presented here, as not much new physics is involved.

The FORTRAN 90 source code used for this work can be downloaded from

https://github.com/LiuGangKingston/one-dimensional

or obtained by email request from the author.

9 Summary

By formulating the work done by the external forces on the system and applying the principles of statistical physics, we derived Eq. (2) for the determination of the period vectors of a crystal under any external stress and temperature. It applies to both classical physics and quantum physics. The existing statistical theory for external pressure is covered as a special case. Equation (2) is also consistent and can be combined with the previously derived equation for the period vectors in the framework of the Newtonian dynamics.

Equation (2) is not only the mechanical equilibrium condition, but also the equation of state for a crystal under external stress and temperature, expressed in terms of the period vectors. It can be used to predict crystal structures and to study structural phase transitions and any kind of crystal expansions. Furthermore, Eq. (2) is also the microscopic temperature-dependent form of the generalized Hooke’s law for linear elastic crystals, and can be used to calculate the corresponding elastic constants for given temperatures. This stands to play an important role in studies of piezoelectric and piezomagnetic phenomena that are caused by an increase of the external pressure on crystals in specific direction.

Finally, we illustrated the maximum external pulling force before the system disintegrates, the melting temperature, the changes in the crystal period induced by the external force and temperature, Hooke’s law, and some thermodynamic properties and processes using a simple one-dimensional model.