Isothermal and adiabatic elastic constants from virial fluctuations

A conventional way to treat periodic boundary conditions is to assume that the system of interest which is placed in the central parallelepiped (or triclinic) box is surrounded by an infinite number of translated identical image boxes that fill up the space [15, 16, 17, 18]. It is convenient to specify these boxes with vector ${\bf n}$ given by

where ${\bf a}$, ${\bf b}$, and ${\bf c}$ are the vectors specifying the three edges of the central box and $\xi_{a}$, $\xi_{b}$, and $\xi_{c}$ can take any integer values. The central box corresponds to ${\bf n}=0$. Then, the classical Hamiltonian of a system of $N$ atoms subject to periodic boundary conditions can be written as

Here

is the kinetic energy, in which $p_{i\alpha}$ is the $\alpha$th component of the momentum of atom $i$ with mass $m_{i}$ and $U\left(\{{\bf x}_{{\bf n}i}\}\right)$ is the system potential energy, which depends on the set of atomic coordinates $\{{\bf x}_{{\bf n}i}\}$ both for the central box and the image boxes. These coordinates are given by

where ${\bf x}_{i}$ is the coordinate vector of atom $i$ in the central box. For example, in the case of pair interactions, the potential energy can be written as [15, 18]

Here the sum over ${\bf n}$ extends over all image boxes. However, for potentials with a cutoff and under the minimum image convention only the central box and its 26 nearest-neighbor image boxes contribute to the sum over ${\bf n}$.

An alternative but equivalent form of the potential energy for a periodic system that is useful for numerical applications can be written in terms of the atomic-group potential energies as was done in Ref. [18]. The system potential energy is usually given by a sum of two-atom, three-atom, four-atom or, in general, few-atom interaction terms. Each of these few-atom terms is commonly referred to as a group [18] and can be specified with an index $k$. The number of groups, the number of atoms in each group, and the number of groups in which an atom participates depends on the total potential energy but otherwise are completely arbitrary. The potential energy of the group $k$ is denoted by $u_{k}(\{{\bf x}_{i}^{k}\})$, where $\{{\bf x}_{i}^{k}\}$ is the set of coordinates for group $k$ with the index $i$ running from 1 to $N_{k}$, the total number of atom in the group $k$. (Note that the same ${\bf x}_{i}$ for atom $i$ appears as ${\bf x}_{i}^{k}$ with a different $k$ for every group in which atom $i$ participates.) The Hamiltonian of a periodic system in the atomic-group form can be written as

where the sum over $k$ runs over all groups associated with the central box. Thompson et al. [18] discuss in detail how such groups are defined and selected. Both Eq. (2) and Eq. (6) forms of the Hamiltonian are used below to derive expressions for the elastic constants.

Elastic constants describe system response to the imposed strain and, in general, depend on a particular form of strain used. Consider a parallelepiped-shaped sample of solid material specified by three basis vectors starting from the common vertex. Any homogeneous transformation of this sample can be represented by a real $3\times 3$ matrix M that transforms vectors ${\bf x}$ in the unstrained system to vectors ${\bf x}^{\prime}$ in the strained system

The matrix ${\bf M}$ has nine independent components of which six represent pure strain and three represent rotation or reflection. A convenient way to separate strain from rotation or reflection is to use the polar decomposition of the matrix ${\bf M}$ wherein ${\bf M}$ is uniquely written as [19]

where ${\bf O}$ is an orthogonal matrix and ${\bf S}$ is a symmetric matrix. Thus, the action of matrix ${\bf M}$ can be viewed as inducing a pure strain with the matrix ${\bf S}$ followed by a pure rotation or reflection, given by the matrix ${\bf O}$, in which the strain state remains unchanged. Clearly, if ${\bf M}$ is symmetric then ${\bf M}={\bf S}$ and the matrix ${\bf M}$ itself induces only pure strain with no rotation or reflection. The difference between the matrix ${\bf S}$ and the identity matrix ${\bf I}$

arguably represents the simplest definition of strain in that it is given by a linear deviation of the matrix ${\bf S}$ from the identity matrix. From Eq. (9) the matrix S is written in terms of $\bm{\mathcal{E}}$ as

Using Eq. (8) the stress $\bm{\mathcal{E}}$ can be re-expressed through ${\bf M}$ as

where superscript T denotes the matrix transpose. The strain $\bm{\mathcal{E}}$ is sometimes referred to as the Biot strain [20].

A measure of strain most often used in the literature on elastic constants is the Lagrangian strain $\bm{\mathcal{H}}$ given by

for which

In this work we choose to define elastic constants using $\bm{\mathcal{E}}$ rather than $\bm{\mathcal{H}}$ for the following reasons. First, it is more natural to study system response to linear changes in geometry as given by Eq. (10) rather than the nonlinear ones as given by Eq. (13). Second, because of the linearity of Eq. (10) the use of $\bm{\mathcal{E}}$ leads to simpler general expressions for the elastic constants. Third, there exists a one-to-one correspondence between the elastic constants defined using $\bm{\mathcal{E}}$ and those defined using $\bm{\mathcal{H}}$ (see Appendix A). Thus one set of constants can always be converted to the other one if necessary. (Note that for pair potentials the use of $\bm{\mathcal{H}}$ rather than $\bm{\mathcal{E}}$ leads to slightly simpler microscopic expressions for elastic constants. This is generally not true for more complex potentials.)

The second-order isothermal elastic constants can be defined using Helmholtz free energy $A$ as [21]

Here $V$ is the system volume and $\varepsilon_{\alpha\beta}$ denotes the component of the symmetric strain tensor $\bm{\mathcal{E}}$; subscripts $T$ and $\bm{\mathcal{E}}=0$ mean that the derivatives are evaluated at constant temperature and for $\bm{\mathcal{E}}=0$. Since the tensor $\bm{\mathcal{E}}$ is symmetric, the derivative with respect to $\varepsilon_{\alpha\beta}$ in Eq. (14) means $(\partial/\partial\varepsilon_{\alpha\beta}+\partial/\partial\varepsilon_{\beta\alpha})/2$.

The Helmholtz free energy in Eq. (14) is given by

where

Here the integration extends over the whole phase space. The dependence of $A$ on strain comes from the dependence of the integral in Eq. (16) on integration limits in coordinate space. A standard technique to avoid the mathematically inconvenient strain-dependence of the integration limits is to transfer this dependence to the Hamiltonian explicitly while holding the integration limits fixed. This can be achieved by using a suitable variable transformation. This approach was first used by Bogoliubov to derive the microscopic expression for pressure as a derivative of $A$ with respect to volume [22]. It was later generalized to express the dependence of $A$ on the full strain tensor [2, 4, 6].

Here we use the following strain-dependent canonical transformation of coordinates and momenta of the Hamiltonian (2)

Note that the transformation of coordinates affects both the atoms in the central box and the image boxes. Inserting $\widetilde{\bf p}_{i}$, $\widetilde{\bf x}_{{\bf n}i}$ in the Hamiltonian (2) leads to the strain-dependent transformed Hamiltonian $\widetilde{H}=\widetilde{H}_{0}+\widetilde{U}$, and the strain-dependent Helmholtz free energy $A$ for which the derivatives with respect to the strain components can now be evaluated. Using Eqs. (14, 15, 16) we obtains

where

is the so-called fluctuation term,

is the kinetic term, and

is the Born term. The derivatives in Eqs. (19, 20, 21) are evaluated at constant $T$ and for $\bm{\mathcal{E}}=0$ as in Eq. (14) with the corresponding subscripts dropped for brevity; the brackets denote averaging over the canonical ensemble.

The derivatives of $\widetilde{H}_{0}$ and $\widetilde{U}$ with respect to strain components in Eqs. (19, 20, 21) can be rewritten in terms of derivatives with respect to momenta and coordinates using Eq. (17). Calculation of these derivatives is rather straightforward for $\widetilde{U}$ because of the linear dependence of $\widetilde{\bf x}_{{\bf n}i}$ on $\bm{\mathcal{E}}$. The derivatives of $\widetilde{H}_{0}$ can also be calculated easily if one recalls that $\sum_{\alpha}\widetilde{p}_{i\alpha}^{2}$ is a dot product of vector $\widetilde{\bf p}_{i}$ with itself and

Calculating the derivatives explicitly in Eq. (19) we obtain for the fluctuation term

where

is the microscopic stress tensor. The second line of the last equation is the symmetrized form of the virial tensor (divided by volume) for a system with periodic boundaries which was derived by Thompson et al. [18] without using strain derivatives explicitly. This form of the virial was referred to as the atom form in Ref. [18].

Using Eqs. (20,22) the kinetic term is calculated to be

Note that the kinetic term (25) is $\frac{3}{4}$ of the kinetic term defined using the Lagrangian strain [12]. We will use this fact in subsection III.2 when comparing the elastic constants defined using the Biot and Lagrangian strains.

The Born term is

Here the subscript sym is introduced to shorten the corresponding expression. It means that the expression in brackets is the average of four terms, namely the term shown and the three terms obtained from it by the following subscript changes $(\alpha\leftrightarrow\beta)$, $(\mu\leftrightarrow\nu)$, and $(\alpha\leftrightarrow\beta,\mu\leftrightarrow\nu)$. The last equation represents the atom form of the Born term for a system with periodic boundaries.

For practical applications it is useful to derive the equivalent atomic-group forms of Eqs. (24,26). These can be obtained from the atomic-group form of the Hamiltonian (6) by using the canonical transformation

and repeating the steps that lead to Eqs. (24,26). This gives

where

is the symmetrized form of the atomic-group virial tensor for a system with periodic boundaries which was also derived in Ref. [18] without using strain derivatives explicitly. Thus, the atomic-group virial $W_{\alpha\beta}^{k}$ given by Eq. (29) is expressed in terms of coordinate components of atoms in a given group and the partial force components due to the group potential energy $u_{k}$.

For the atomic-group form of the Born term we obtain

Equation (30) can be viewed as an extension of Eq. (29) for the group-form of the virial for periodic systems: the virial is obtained from the first derivatives of $\widetilde{U}$ with respect to strain components whereas Eq. (30) is derived using the second derivatives. Equation (30) can be used directly to calculate the Born term contribution to the elastic constants: once the potential energy in the group form is expressed in Cartesian coordinates the second derivatives in Eq. (30) can be calculated. These expressions can be complicated for potentials that involve few-particle interactions.

Here we show that the right-hand side of Eq. (30) can be reduced to a form that involves only first derivatives by using integration by parts with respect to one of the derivatives. The basic steps of the derivation are similar to those used in the derivation of the equipartition theorem [23, 24]. A similar approach was used in Ref. [25] to express an ensemble-averaged matrix of second derivatives in terms of the force-force covariance matrix. The details of the derivation are given in Appendix B. The final expression, which is the central result of this work, is

where the virial $\widehat{W}_{\alpha\beta}^{k}$ is given by

with

Comparison of $\widehat{W}_{\alpha\beta}^{k}$ given by Eq. (32) and $W_{\alpha\beta}^{k}$ given by Eq. (29) shows that $\widehat{W}_{\alpha\beta}^{k}$ depends on the same atomic coordinate components but the partial forces are replaced with the total force component for a given atom. The second line of Eq. (32) involves products of the arithmetic average of coordinate components of all atoms in the group and the sum of force components for atoms in the group; these terms ensure that $\widehat{W}_{\alpha\beta}^{k}$ is translationally invariant. Translational invariance of $\widehat{W}_{\alpha\beta}^{k}$ in Eq. (32) can be verified by displacing the system along any of the Cartesian directions and using Eqs. (33) together with the translational invariance of the total force components. Thus, the Born term given by Eq. (31) depends on the first derivatives of either atomic-group or total potential energies but does not contain the second derivatives. Because of this fact, Eq. (31) can be of interest for numerical calculation of elastic constants, particularly for complex potentials that involve few-particles interaction terms.

Adiabatic elastic constants are defined using system internal energy $E$ as [21]

where the derivatives are now evaluated at constant entropy $S$ and for $\bm{\mathcal{E}}=0$. It was shown in Refs. [26, 7, 6] that the adiabatic elastic constants $C_{\alpha\beta\mu\nu}^{S}$ can be calculated using the same expressions as for the isothermal ones but with all canonical-ensemble averages replaced with averages over the microcanonical ensemble. Thus, Eqs. (18,19,20,21,23,25,26,30) remain valid for adiabatic constants with the superscript $S$ replacing the superscript $T$ and the brackets now denoting microcanonical averaging. We show in Appendix C that Eq. (31) also remains valid for the adiabatic case (with averages performed over the microcanonical ensemble). Thus, adiabatic elastic constants can also be calculated using atomic-group virials.

In the case of argon isothermal elastic constants were calculated for several values of pressure and temperature. The results for adiabatic constants for argon are qualitatively similar and will not be reported here. Solid argon forms a face-centered cubic (fcc) lattice. There are three independent elastic constants for cubic crystals [27], which are commonly chosen to be $C_{1111}^{T},C_{1122}^{T}$, and $C_{1212}^{T}$ whose values we report below.

The atomic interactions in argon were modeled with the Lennard-Jones potential

where $r$ is the distance between the two atoms, the parameters $\sigma$ and $\epsilon$ were set to be $3.4$ Å and $1.67\times 10^{-21}$ J, respectively [28, 27]. The cubic simulation cell consisting of $12\times 12\times 12$ conventional four-atom fcc unit cells was used. The cutoff value for the potential was set at the relatively high value of 28 Å to minimize any possible cutoff effects. A time step of 0.5 fs was used. Atomic coordinates were recorded every 100 fs. These coordinates were used in Eqs. (30) and (31) to calculate two forms of the Born term as functions of time. The stress tensor components for the fluctuation term were calculated using LAMMPS. Elastic constants were calculated at zero pressure and at $P=1$ GPa. For each of the two pressures three temperature values were considered. These temperature values corresponded, approximately, to 2%, 40%, and 90% of the argon melting temperature at a given pressure, which we will refer to as the low, moderate and high temperature, respectively. The melting temperature of argon is approximately 84 K at zero pressure [29] and 250 K at 1 GPa [30]. Thus, the temperatures of 2 K, 33 K, and 75 K for zero pressure and 5 K, 100 K, and 225 K for the 1 GPa pressure were considered. The fluctuation terms, the kinetic term, and the two forms of the Born term were calculated from 1 ns equilibrium NVT trajectory for each set of pressure and temperature values. As an example, the behavior of these terms as functions of time for $C_{1111}^{T}$ is shown in Fig. 1 for the case of zero pressure and $T=33$ K.

Tables 1 and 2 list the elastic constants along with the kinetic, fluctuation, and Born terms for each set of pressure and temperature values. In the case of the Born terms two values are listed: the first line shows the values obtained from Eq. (30) and the second line gives the values obtained from Eq. (31). The two values for the elastic constants in Tables 1 and 2 are calculated using the two forms of the Born term with the same kinetic and fluctuation terms. The uncertainties in Tables 1 and 2 correspond to one-sigma standard error. They were calculated using statistical analysis for the time-dependent correlated data [31].

One can see from Tables 1 and 2 that elastic constants calculated using two forms of the Born term are indeed the same within uncertainties for all pressure and temperature values considered. Note that in all cases the Born term calculated using Eq. (31) converges slower than the one obtained from the original expression (30), as reflected by the corresponding uncertainties in Tables 1 and 2. (See also panel (b) of Fig. 1). However, for moderate and high temperatures the Born term in the form of Eq. (31) converges faster than the fluctuation term (cf. Fig. 1), whereas for low temperatures it converges slower than the fluctuation term. Thus, for moderate and high temperatures the overall convergence of the elastic constants calculated using the Born term given by Eq. (31) is approximately the same as that for the elastic constants obtained using the Born term given by Eq. (30). For low temperatures the elastic constants calculated using Eq. (31) converge slower.

In the case of silicon, adiabatic rather than isothermal elastic constants were calculated. This was done for the following two reasons: to show that the new expression for the Born term can be also be applied to obtain adiabatic constants and to compare the calculated values to earlier MD-based adiabatic results for elastic constants of silicon [32, 11]. Similarly to fcc argon, silicon is a cubic crystal. Thus, only $C^{S}_{1111}$, $C^{S}_{1122}$, and $C^{S}_{1212}$ needed to be calculated. The system was studied at zero pressure and the following four temperatures: 300 K, 888 K, 1164 K, and 1477 K. The latter three temperatures corresponded to the ones considered in Refs. [32, 11].

The atomic interactions in silicon were modeled with the Stillinger-Weber potential [33] which consists of two- and three-atom interaction terms. A cubic simulation cell of 216 silicon atoms was used. All force-field and other simulation parameters reported in Refs. [32, 11] were kept except for the total number of simulation steps. References [32] and [11] used 150000 and 1 million steps, respectively. We used 2 million steps for the three higher temperatures and 6 million steps for 300 K to obtain good convergence of the Born terms obtained using Eq. (31).

Kinetic, fluctuation, and Born terms were calculated using microcanonical averaging, i. e. NVE trajectories, for each of the four temperatures. Two forms of the Born terms were calculated: the original one using Eq. (30) and the new form based on Eq. (31). For each of these two forms both two-atom and three-atom contributions to the total Born terms were calculated. These results are summarized in Table 3. The two values for the elastic constants reported in Table 3 were calculated using the two forms of the total Born term with the same kinetic and fluctuation terms.

One can see that elastic constants calculated using two forms of the Born term agree with each other within uncertainties for all temperatures considered. This remains true for the two-atom, three-atom, and total Born terms calculated using the two approaches. For all cases the Born terms calculated using Eq. (31) converges slower than the ones obtained from the original expression (30), as reflected by the corresponding uncertainties. Comparison of uncertainties for the fluctuation terms and the Born terms calculated from Eq. (31) shows that $B_{1111}^{S}$ and $B_{1122}^{S}$ converge slower and $B_{1212}^{S}$ much faster than the corresponding fluctuation terms for all four temperatures.

The results of Ref. [11] for the elastic constants of silicon along with the constitutive Born, kinetic, and fluctuation terms for 888 K, 1164 K, and 1477 K are also listed in Table 3. Also shown are the elastic constants reported in Ref. [32] for the same three temperatures. When comparing our results to the previously published data we need to bear in mind that the elastic constants reported in Refs. [32, 11] were calculated using the Lagrangian strain whereas our results were based of the Biot strain. Fortunately, the transformation between the Lagrangian-strain and Biot-strain elastic constants (and their constitutive terms) is straightforward. It follows from Eq. (40) that at zero pressure the adiabatic elastic constants obtained using the Lagrangian and Biot strains are identical, thus they can be compared directly. However, this is not true for the kinetic and Born terms for the two strains considered separately. It was noted following Eq. (25) that the Biot kinetic term is $\frac{3}{4}$ of the Lagrangian one. The fluctuation terms for the two strain measures are identical because the stress tensors obtained using the two strains are identical. This implies that the Biot Born term is equal the Lagrangian Born term plus one fourth of the Lagrangian kinetic term. Thus, the original Lagrangian results of Ref. [11] for the kinetic and total Born terms were converted to the Biot form using the recipe just above when reporting them in Table 3.

One can see that the total Born terms of Ref. [11] and both forms of the total Born terms calculated by us are generally in good agreement: for most cases they are identical within uncertainties. Elastic constants calculated by us also compare well to the elastic constants from Refs. [32] and [11]. The larger differences for $C_{1212}^{S}$ at 1164 K and 1477 K are still within uncertainties and are due to the very large uncertainties of the corresponding fluctuation terms.