# Computation of the thermal elastic constants for arbitrary manybody potentials in LAMMPS using the stress-fluctuation formalism

Germain Clavier, Aidan P Thompson

► To cite this version:

Germain Clavier, Aidan P Thompson. Computation of the thermal elastic constants for arbitrary manybody potentials in LAMMPS using the stress-fluctuation formalism. Computer Physics Communications, 2023, 286, pp.108674. 10.1016/j.cpc.2023.108674. hal-04411734

HAL Id: hal-04411734

https://hal.science/hal-04411734v1

Submitted on 26 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Computation of the thermal elastic constants for arbitrary manybody potentials in LAMMPS using the stress-fluctuation formalism ☆

<figure><img src="image_2.png" /></figure>

Germain Clavier<sup>a, *</sup>, Aidan P. Thompson<sup>b</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Article history:**

Received 19 September 2022

Received in revised form 10 January 2023

Accepted 17 January 2023

Available online 20 January 2023

**Keywords:**

Molecular simulation

Elastic constants

Stress-stress fluctuations

LAMMPS

## ABSTRACT

This paper describes the implementation of the stress-fluctuation technique into the LAMMPS code to compute the anisotropic thermal elastic constants tensor of materials. The implementation provides both methods for computing the analytical fluctuation expressions and also a generic numerical derivative method. The former makes the extension to new potentials straightforward, as it requires writing code only for the second derivatives of each energy term w.r.t. distance, angle, etc. The latter provides a generic interface to compute an accurate approximation of the elastic constants for any potential already implemented in LAMMPS. We show how both methods compare with the direct deformation computation in several test cases and discuss the implementation advantages and limitations.

© 2023 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license

(http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Computation of elastic constants of material through molecular simulation has been a research topic since early on in the field. Early works from Squire et al. [1], Parrinello and Rahman [2] or Ray [3] tried to derive several ways to compute the elastic constants of materials from atomistic simulations of crystalline models. With the important increase of computational power in the last forty years and the development of more sophisticated molecular models, molecular modeling has taken an important place in understanding materials' properties. This can be illustrated by the wide variety of materials that are now studied using molecular modeling, ranging from alloys to polymer materials, including complex crystals structures, glassy unordered materials, or materials in extreme pressure and temperature conditions. Accurate computation of materials properties through molecular modeling thus became increasingly important. However, in the case of elastic constants, the spread of refined generic methods suffered from several difficulties. On the one hand, methods simply relying on volume fluctuations in constant stress ensembles (NσT) are dependent on the computational method [4]: the dynamics of the pressure bath can affect the covariance of the strain tensor components, which are the values used to compute elastic constants

at constant pressure. On the other hand, in the NVT ensemble, the necessity of complex analytical derivations made the use of the general formulas complicated past pair potential interactions. A general framework had to wait the work of Lutsko [5,6], and following works still had to find workarounds for cases including 3 bodies potential or complex derivation for complex potentials such as AIREBO [7,8]. The values computed this way seemed however, less dependent on the simulation techniques and parameters [4].

In this paper, we present the integration of generic NVT elastic constant computation framework in the LAMMPS simulation software. To our knowledge, this is the first publicly available integration of such method directly in a molecular simulation package, despite several uses in the literature. Compared to other codes, this integration can use both analytical derivations and numerical derivative to evaluate elastic constants. The analytical derivations actually need forcefield specific terms which need to be implemented for all potential interactions formulas, but provide exact contribution to the elastic constants and are faster to compute than using numerical derivative. In the mean time, numerical derivative technique allows to double check new analytical implementations and replace them for potential forms that are too complicated for analytical derivation. The paper will start with a general reminder of the underlying formulas and will use it as a base to explain the actual implementations. We will then present examples showing the agreement between both methods and compare them with elastic constants obtained through direct deformation simulations. Finally we will present future works and extension of the method in the LAMMPS code.

<footnote><sup>a</sup> Eindhoven University of Technology, Dep. of Mech. Eng., P.O. Box 513, 5600MB, Eindhoven, the Netherlands</footnote>

<footnote><sup>b</sup> Sandia National Laboratories, Albuquerque, NM 87185, USA</footnote>

<footnote>☆ The review of this paper was arranged by Prof. W. Jong.</footnote>

<footnote>* Corresponding author.
E-mail address: g.m.g.c.clavier@tue.nl (G. Clavier).</footnote>
## 2. Computation of elastic constants in the NVT ensemble

The stress tensor $\sigma$ and the elastic constant tensor $C$ are the first and second derivatives of the free energy with respect to the strain tensor $\epsilon$. In the NVT ensemble, this reduces to the following equations for the stress [5,8]:

$$ \sigma_{\alpha\beta} = \langle \sigma_{\alpha\beta}^B \rangle - \rho k_B T \delta_{\alpha\beta} \quad (1) $$

and for the elastic constant tensor

$$ C_{\alpha\beta\mu\nu} = \langle C_{\alpha\beta\mu\nu}^B \rangle - \frac{V}{k_B T} \left[ \langle \sigma_{\alpha\beta}^B \sigma_{\mu\nu}^B \rangle - \langle \sigma_{\alpha\beta}^B \rangle \langle \sigma_{\mu\nu}^B \rangle + \rho k_B T (\delta_{\alpha\mu} \delta_{\beta\nu} + \delta_{\alpha\nu} \delta_{\beta\mu}) \right] \quad (2) $$

in this equation $\sigma_{\alpha\beta}^B$ and $C_{\alpha\beta\mu\nu}^B$ are called the Born stress tensor and the Born matrix. These are the first and second derivatives of the potential energy $U$ w.r.t. strain, given by:

$$ \sigma_{\alpha\beta}^B = \frac{1}{V} \frac{\partial U}{\partial \epsilon_{\alpha\beta}} \quad (3) $$

and

$$ C_{\alpha\beta\mu\nu}^B = \frac{1}{V} \frac{\partial^2 U}{\partial \epsilon_{\alpha\beta} \partial \epsilon_{\mu\nu}}. \quad (4) $$

A detailed computation of these terms can be found in the work of Van Workum et al. [7,8] or Vogiatzis et al. [9]. Here $T$ is the temperature, $V$ the volume, $\rho = N/V$ is the number density of particles, $N$ is the number of particles, $k_B$ the Boltzmann constant and $\delta_{\alpha\beta}$ is the Kronecker identity tensor. The Born stress tensor appearing in Eq. (1) is equivalent to the virial contribution to the microscopic stress tensor generally used in molecular simulation [10]. Similarly the kinetic part of the microscopic stress tensor averages to the second term of the right-hand side of Eq. (1). Because the Born stress term is already computed in most molecular simulation codes, the second term of the right-hand side of Eq. (2) can be easily computed in a post-process step. However, the first term of Eq. (4), the Born matrix, is more problematic and is not expressible in terms of standard quantities. This motivated the development of general LAMMPS implementation of the Born matrix calculation. The two following sections go into the details of the two complementary methods we developed to compute the Born matrix.

### 2.1. Analytical derivative method

The internal energy of molecular systems can be split into two contributions: the potential energy and the kinetic energy. The kinetic energy depends only on the particle momenta and its contribution to stress and elastic constants average out in the last terms of Eqs. (1) and (2). The potential energy part depends on the interaction potential and the particle coordinates. We can write the interaction potential as $U(q)$, where $q$ represents the 3N coordinate of the atoms corresponding to a particular point in configuration space. Equivalently, we can consider every internal coordinate involved in every energy term, such as inter-atomic distances for two-body terms, bond angles for three-body terms, etc. Representing the generalized coordinates for a single energy term as $\chi$, we can say that $U(q) = \sum_{\chi} u(\chi)$ where $\chi$ is the set of generalized coordinates related to the energy term $u(\chi)$. Using the chain rule, Eq. (3) and Eq. (4) can be rewritten as [8]:

$$ \sigma_{\alpha\beta}^B = \frac{1}{V} \sum_{\chi} \frac{\partial u(\chi)}{\partial \chi} \frac{\partial \chi}{\partial \epsilon_{\alpha\beta}}, \quad (5) $$

and

$$ C_{\alpha\beta\mu\nu}^B = \frac{1}{V} \sum_{\chi} \left[ \frac{\partial^2 u(\chi)}{\partial \chi^2} \frac{\partial \chi}{\partial \epsilon_{\alpha\beta}} \frac{\partial \chi}{\partial \epsilon_{\mu\nu}} + \frac{\partial u(\chi)}{\partial \chi} \frac{\partial^2 \chi}{\partial \epsilon_{\alpha\beta} \partial \epsilon_{\mu\nu}} \right]. \quad (6) $$

The expression for the Born stress, Eq. (5), is equivalent to the group form of the virial already implemented in LAMMPS [10].

The derivation for pair potentials is straightforward and, for the stress, gives back the general formula using relative particle coordinates and pairwise force (see eq 15 of [4]). Explicit expressions for the Born matrix, Eq. (6) for the specific cases of angular and dihedral potentials have been published previously [8,9]. However, the expressions have not been implemented in any commonly available molecular simulation packages before now.

The key idea of the general implementation presented here uses the fact that the derivatives appearing in the Born matrix are products of two independent factors. One is the derivative of the energy with respect to a given generalized coordinate (angle, dihedral, ...). The other is a configurational factor which depends only on the relative positions of the particles. On a computational level, this allows the two factors to be handled separately. A general function computes the configurational factor for every interaction. The energy derivative factor is calculated by a potential-specific function. In the current LAMMPS implementation, each potential term is coded as a separate C++ class derived from general Pair, Bond or Angle parent classes. Each child class only needs to provide the energy derivative factor. The new Compute derived class for the Born matrix computes all the configurational factors, multiplies both terms and sums over all interactions. This greatly reduces the effort required to implement the Born matrix for any new potential.

### 2.2. Numerical derivative method

Yoshimoto et al. [11] derived a general formula for three-body interaction potentials. This employed small displacements of atoms to compute local contributions to a specific form of the Born matrix. Inspired by this result, Zhen et al. [12] rewrote the Born matrix as follows:

$$ C_{\alpha\beta\mu\nu}^B = \delta_{\alpha\nu} \sigma_{\beta\mu}^B + \delta_{\alpha\mu} \sigma_{\beta\nu}^B + D_{\alpha\beta\mu\nu} \quad (7) $$

with

$$ D_{\alpha\beta\mu\nu} = - \frac{d\sigma_{\mu\nu}^B}{\epsilon_{\alpha\beta}} = - \lim_{e_{\alpha\beta} \to 0} \frac{\sigma_{\mu\nu}^{B,1} - \sigma_{\mu\nu}^{B,0}}{e_{\alpha\beta}}. \quad (8) $$

In Eq. (8), $e_{\alpha\beta}$ is a small strain value, $\sigma_{\mu\nu}^{B,0}$ is the stress in a reference state and $\sigma_{\mu\nu}^{B,1}$ is the stress obtained after applying the strain $\epsilon_{\alpha\beta}$ of amplitude $e_{\alpha\beta}$ to the reference state. As shown in reference [12], $\sigma_{\mu\nu}^{B,1}$ can be efficiently calculated in LAMMPS by applying an affine deformation to the atomic coordinates and recalculating the forces in the normal way. An explicit box deformation is not required, which eliminates the need for parallel communication, as well as eliminating any dependence on periodic boundary conditions.

This method allows the computation of elastic constants for arbitrary complex potential interactions without the need for additional code. This comes at the price of increased computational cost. The implementation of the numerical derivative in Eq. (8) uses a 3-point centered difference scheme. This requires the computation of two additional sets of forces for each of the 6 strain directions. This can make the computation longer than using the analytical derivative method. Concerning the implementation, the
compute born/matrix command handles the entire calculation, but it must be provided with a specific instance of the virial stress computation from LAMMPS. Details can be found in the command documentation.

In the next section, both methods will be compared for a simple Lennard-Jones system crystal. As will be shown, they both agree with one another and with standard direct deformation simulations for the computation of elastic constants. For bond angle and dihedral angle potential terms, we will show with a simple polyethylene crystal model that both methods also give the same results, which provides a stringent validation of the implementations. Finally we will show how the numerical derivative method compares with the explicit deformation method for a cubic crystal with more complex manybody potentials, including the MEAM and ADP potentials for copper and the Stilling-Weber potential for silicon, and for a wurtzite structure of zinc oxide using a Tersoff-style potential.

# 3. Results and discussion

## 3.1. General procedure

In the subsequent sections, the simulation procedure is as follows, except where stated otherwise. The atomic system is generated in a cubic crystalline configuration. Before any production simulation, the system is equilibrated in the anisotropic NPT ensemble at 0.1 MPa using a Nose-Hoover chain thermostat and anisotropic barostat with 10 chain particles each and with respective coupling constants $\tau_T = 100$ fs and $\tau_P = 1000$ fs. We set the timestep to $dt = 1$ fs. The thermostat settings for all subsequent simulations in the NVT ensemble are the same as for the NPT ensemble. The temperature of the simulations depends on the simulated material.

A first NPT equilibration phase of 100 ps is performed followed by an NPT acquisition phase of 1 ns. During the acquisition phase, the running average of the volume $\langle V \rangle$ is computed using volume on every timestep. At the end, the shape of the box is set as a cube of length $l^0 = \sqrt[3]{\langle V \rangle}$. This system is the reference system for both the numerical derivative method and the stress-stress fluctuation computation.

From then on, the stress tensor and Born matrix are sampled in a 1 ns long NVT simulation, using either the analytical (if possible) or numerical derivative methods. The running average of the Born matrix is sampled every timestep on the fly. For numerical stability, the fluctuation term (second term of Eq. (2)) is not computed on the fly but at the end of the simulation. The instantaneous value on every timestep is used to compute the covariance terms.

For the direct deformation simulations, the procedure is the same as in previous work [4] but a series of NVT simulations are performed with two kinds of deformation applied to the simulation cell

* For uniaxial strains, the box is deformed by changing the x-dimension to $l_x = \xi l_x^0$.

* For shear strains, the box is deformed by displacing the period lattice vector initially aligned with in the y-direction by $xy = \xi l_x^0$.

In the case of uniaxial strains, $\xi$ values range from -0.005 to 0.005 with a step of 0.001 are used. For the shear strains, only the positive values of $\xi$ are used. In both cases, $\xi$ is the strain value. Each NVT simulation of the deformed boxes begins with an equilibration phase of 100 ps followed by an acquisition phase of 1 ns.

<figure><img src="image_4.png" /><figcaption>Fig. 1. Stress strain curve of the Argon crystal at 60 K with regard to their respective strain. $\sigma_{11}$ and $\sigma_{22}$ are plotted with regard to $\epsilon_{11}$ while $\sigma_{12}$ is plotted with regard to $\epsilon_{12}$.</figcaption></figure>

<caption>Table 1</caption>

Anisotropic elastic constants [GPa] for argon crystal computed from the 3 different methods at 60 K. In the case of direct deformation, the error bar is taken as the square root of the variance of the estimated parameters. In the case of stress fluctuation, it is the standard error from the 3 symmetrical values of the C matrix.

<table><thead><tr><th></th><th>Deformation</th><th>Analytical</th><th>Numerical</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>2.19 ± 0.03</td><td>2.21 ± 0.01</td><td>2.23 ± 0.01</td></tr><tr><td>C<sub>12</sub></td><td>1.22 ± 0.02</td><td>1.23 ± 0.01</td><td>1.22 ± 0.01</td></tr><tr><td>C<sub>44</sub></td><td>1.37 ± 0.02</td><td>1.42 ± 0.01</td><td>1.41 ± 0.01</td></tr></tbody></table>

## 3.2. Lennard-Jones argon crystal

The first model is an argon FCC crystal. The objective of this section is to show that the direct deformation, the analytical, and the numerical stress fluctuation methods all give statistically equivalent for a simple reference model. The simulation cell is a 5 × 5 × 5 unit cell structure (500 atoms) with an initial lattice constant of 0.5392 nm. The atomic mass is set to 39.95 g mol⁻¹, and the atoms interact through a Lennard-Jones potential with parameters $\epsilon = 119.8 k_B$ and $\sigma = 0.3405$ nm. The cutoff distance is set to 1.2 nm. The simulations are performed at a temperature $T = 60$ K. Due to the cubic symmetry of the FCC crystal structure, there are only three independent values in the C tensor: $C_{11}$, $C_{12}$ and $C_{44}$.

After NPT equilibration, the average box length is 2.7047 nm which is equivalent to a crystal lattice of 0.5409 nm. All three dimensions of the simulation box are then resized to this average length for NVT simulations. The average pressure measured in the reference configuration in the NVT ensemble is $P = -0.0453$ MPa ± 1.76, where the error estimate is the standard deviation about the mean. As described, the system is then on the one hand subjected to direct deformation simulation, and on the other hand simulated in the NVT ensemble to compute the elastic constants directly. A plot of the stress values $\sigma_{11}$, $\sigma_{22}$ and $\sigma_{12}$ with regard to strain is shown in Fig. 1 and the results of all methods are given in Table 1.

The agreement of the three methods is within the error-bars in all cases.
<caption>**Table 2**</caption>

Anisotropic elastic constants [GPa] for polyethylene crystal at 300 K computed with three different methods. Due to the anisotropy of the C tensor, no error-bar is provided for the stress-fluctuation computations, but the agreement between the method is comparable to the case of Argon. Also note the high error from direct deformation involving stress along z axis.

<table><thead><tr><th></th><th>Deformation</th><th>Analytical</th><th>Numerical</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>6.32 ± 0.17</td><td>8.54</td><td>8.40</td></tr><tr><td>C<sub>22</sub></td><td>6.08 ± 0.15</td><td>8.05</td><td>7.91</td></tr><tr><td>C<sub>33</sub></td><td>172 ± 1.04</td><td>163.92</td><td>163.64</td></tr><tr><td>C<sub>12</sub></td><td>5.85 ± 0.15</td><td>6.77</td><td>6.72</td></tr><tr><td>C<sub>13</sub></td><td>4.12 ± 0.53</td><td>5.66</td><td>5.70</td></tr><tr><td>C<sub>23</sub></td><td>4.51 ± 0.58</td><td>5.69</td><td>5.73</td></tr></tbody></table>

We also performed the computation of the elastic constants in the NVE ensemble using analytical stress-fluctuations. The results obtained are $C_{11} = 2.59 \pm 0.01$, $C_{12} = 1.61 \pm 0.01$ and $C_{44} = 1.43 \pm 0.01$. The difference of 0.40 GPa with the $C_{11}$ and $C_{12}$ values from the NVT ensemble is close to the estimated value of 0.45 GPa estimated by Hoover et al. [13] and the value computed from the experimental results of Dobbs and Go [4,14].

### 3.3. Molecular terms

Molecular forcefields can include terms involving multiple bonded consecutive atoms (or three atoms all bonded to a fourth one in the case of improper angles). In these cases the potential energy and the forces are generally some functions depending on the bonds angle or dihedral angle. Following Lutsko. Van Workum et al. [7,8] gave analytical derivations of the elastic constant contribution based only on the atoms coordinates at a given time. As for pair interactions (bonded and non-bonded), the analytical derivation must be done for every potential form.

In the following, we use a molecular system of crystalline polyethylene with the OPLS [15] united atom (OPLS-UA) model, as defined by the Moltemplate software [16]. Our interest is in assessing the validity of the approach for molecular systems, rather than accurately describing the mechanical properties of the material. The use of OPLS-UA and long alkane molecules combines simple analytical forms with simple and documented molecular behavior. The system is composed of periodic molecules along the z-axis, where $CH_2$ groups are represented through united atoms. The initial configuration uses cell coordinates from Bunn [17,18] with $5 \times 7 \times 25$ cells respectively along the x, y and z axes, for a total of 3500 particles (the final system is made of 70 molecules each composed of 50 atoms).

The molecular structure of crystalline polyethylene has an orthorhombic symmetry and the cell axes are aligned with the Cartesian frame. As such, 9 $C_{ij}$ values are independent: $C_{11}$, $C_{22}$, $C_{33}$, $C_{12}$, $C_{13}$, $C_{23}$, $C_{44}$, $C_{55}$, and $C_{66}$. It is known that, without proper optimization, united atoms models have lower shear responses than all-atoms models and experimental values [19,20]. They also exhibit unexpected rotator phases at relatively low temperature in alkane crystals [21]. A more problematic behavior is an unphysical displacement of the molecules along the c-axis during the simulations. While the crystalline structure is preserved with the periodic boundary conditions, this unrealistic behavior can affect the simulation results, especially the stress tensor, and this behavior can be enhanced or lowered by strains. For simplicity, we limit our comparison to the $C_{11}$, $C_{22}$, $C_{33}$, $C_{12}$, $C_{13}$ and $C_{23}$ terms, while omitting $C_{44}$, $C_{55}$, and $C_{66}$. The results are presented in Table 2.

As mentioned, united atoms potentials are softer compared to all-atom potentials. This results in computed values systematically lower than experimental estimations [22]. But the order of magnitude of the results is in agreement with similar work using comparable potentials as well as the discrepancy between direct deformation computation and stress-stress fluctuations [18]. The

<caption>**Table 3**</caption>

Anisotropic elastic constants from copper crystals [GPa]computed at 300 K and silicon crystals at 1477 K. Note that the discrepancy between silicon experimental and computed elastic constants is the same as found by Kluge and Ray [26].

<table><thead><tr><th></th><th></th><th>Deformation</th><th>Numerical</th><th>Experimental</th></tr></thead><tbody><tr><td rowspan="3">Cu (ADP)</td><td>C<sub>11</sub></td><td>165.76 ± 0.31</td><td>164.72 ± 0.03</td><td>168.40 [28]</td></tr><tr><td>C<sub>12</sub></td><td>121.51 ± 0.32</td><td>120.63 ± 0.05</td><td>121.40</td></tr><tr><td>C<sub>44</sub></td><td>71.82 ± 0.26</td><td>71.81 ± 0.01</td><td>75.40</td></tr><tr><td rowspan="4">Cu (MEAM)</td><td>C<sub>11</sub></td><td>164.05 ± 0.42</td><td>162.82 ± 0.06</td><td>—</td></tr><tr><td>C<sub>12</sub></td><td>114.49 ± 0.34</td><td>114.09 ± 0.05</td><td>—</td></tr><tr><td>C<sub>44</sub></td><td>69.82 ± 0.59</td><td>68.96 ± 0.04</td><td>—</td></tr><tr><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Si (SW)</td><td>C<sub>11</sub></td><td>134.56 ± 0.24</td><td>132.98 ± 0.04</td><td>148.0 [26]</td></tr><tr><td>C<sub>12</sub></td><td>74.76 ± 0.25</td><td>73.81 ± 0.03</td><td>57.5</td></tr><tr><td>C<sub>44</sub></td><td>46.03 ± 0.28</td><td>43.24 ± 1.21</td><td>70.0</td></tr></tbody></table>

expected symmetry is also retrieved. Yet, a non zero displacement is also seen along the molecules axis during simulations, which can affect results of both methods in different magnitude. For such systems, the use of a more accurate, but more complex potential would compare significantly better with experimental estimations [23].

### 3.4. Many body potentials

The numerical derivative technique becomes very useful to use with many body potentials. Indeed, these potentials take into account local structures surrounding atoms be it through additional angular forces, embedded energy or magnetic terms. While possible if needed, the analytical derivation of Born matrix for such potential can become complicated to derive and implement in a generic way in a simulation code. As shown by Zhen and Chu [12], the numerical derivative method provides an elegant workaround which, as shown for Lennard-Jones crystal, gives similar results with regard to the analytical derivation.

We tested two different copper potentials using different analytical forms, MEAM [24] and ADP [25], as well as the standard Stillinger-Weber potential for Si crystals [26]. The simulation process is the same as described in previous sections. The results are given in Table 3. The potentials were provided through the OpenKim interface [27].

This example shows the reliability of the method at higher temperature using complex potentials. While analytical derivations are not available for the more elaborated models, the numerical derivative methods implemented offers a good replacement.

As a final example, we present the case of zinc oxide (ZnO) in the wurtzite structure. The wurtzite structure has a hexagonal primitive cell of 4 atoms whose parameters are $a=b=0.324$ nm, $c=0.522$ nm, $\alpha=\beta=90^{\circ}$ and $\gamma=120^{\circ}$. The resulting material has six independent elastic constant, which are $C_{11}=C_{22}=C_{33}$, $C_{12}=C_{13}=C_{23}$, $C_{44}=C_{55}$ and $C_{66}$. The simulation box x axis was aligned along the cell's [100] axis and the box z axis along the cell's [001] axis. This makes the box y axis aligned with the [120] direction. An initial structure was made from the primitive hexagonal cell reproduced 10 times along the x axis and along the y axis and 5 times along the z axis. The resulting final system is a simulation box containing 2000 atoms and is transformed in an orthorhombic box of dimensions $l_x = 3.2285$ nm, $l_y = 2.7959$ nm and $l_z = 2.6357$ nm. The model used is a modified Tersoff-style bond-order potential by Albe et al. [29] and adapted to zinc oxide by Erhart et al. [30]. Again, the potential was obtained through the OpenKim interface [27]. Due to the higher number of independent elastic constants, we made a full series of deformations along the three axes of the simulation box and along the three shear directions to get the full elastic constant tensor and we also computed the values using our stress fluctuation implementation. The simulation cell was equilibrated and deformed as described
<table><caption>**Table 4**
Anisotropic elastic constants from wurtzite [GPa] computed at 300 K. The error bars are standard error to the mean from 3 different simulations sets.</caption><thead><tr><th></th><th>Exp [32]</th><th>Erhart et al. [30]</th><th>Deformation</th><th>Numerical</th></tr></thead><tbody><tr><td>C<sub>11</sub></td><td>207</td><td>212</td><td>198.90 ± 0.06</td><td>197.5 ± 1.0</td></tr><tr><td>C<sub>33</sub></td><td>210</td><td>219</td><td>208.33 ± 0.05</td><td>206.4 ± 1.8</td></tr><tr><td>C<sub>12</sub></td><td>121</td><td>116</td><td>140.46 ± 0.03</td><td>139.3 ± 0.6</td></tr><tr><td>C<sub>13</sub></td><td>106</td><td>109</td><td>132.46 ± 0.06</td><td>132.1 ± 1.2</td></tr><tr><td>C<sub>44</sub></td><td>43</td><td>43</td><td>27.08 ± 0.05</td><td>26.0 ± 0.4</td></tr><tr><td>C<sub>66</sub></td><td>45</td><td>48</td><td>29.71 ± 0.03</td><td>28.0 ± 1.1</td></tr></tbody></table>

in Section 3.2, but at 300 K and 0.1 MPa with anisotropic barostat along the three x, y and z axes for NPT equilibration. For the stress-fluctuation procedure, a 1 ns acquisition run resulted in some discrepancies of 6 to 10 GPa in the results compared with deformations. A longer 2 ns acquisition was performed to allow for better convergence in the stress covariance term without computing the Born values. This is because as stated by Zhen et al. [12] the Born term converges much faster than the stress covariance. This means that after Born term convergence, which is computationally expensive, the simulations can be run longer independently, computing only the virial stress as needed. The results are presented in Table 4. It is worth noting that the elastic constants from Ref [30] were computed using molecular statics from equilibrium configurations [31] at zero temperature for both wurtzite and zincblende structures. But, since the potential was fitted using experimental values measured at room temperature [32], we can expect significant differences when using it to compute finite temperature elastic constant through molecular dynamics. This is actually the case but we found again a remarkable agreement between our deformation simulations and the stress-fluctuation computation. The values order is also the same as expected from Erhart's et al. comments regarding their fit with regard to the elastic response of the material ($C_{33} > C_{11}$ but $C_{13} < C_{12}$). This shows the applicability of the method to non orthogonal crystal structures but, unfortunately, that this model of ZnO might not be suited for quantitative elastic constant computation through finite temperature molecular dynamics.

As a last note, to give a general idea of the gain in term of computation time, in the case of ZnO, the 2 ns long simulations of were 1 h 38 min long, while each of the three uniaxial strain deformation were 1 h 50 min long on average and the three shear strains 50 min long on average. All the simulations were carried on similar nodes using 32 processors.

## 4. Conclusion

The previous examples show the agreement between LAMMPS integration of the stress-fluctuation and the typical direct deformation methods for the computation of the elastic constant of both atomic and molecular solids. It is also demonstrated that the analytical derivation and numerical derivative methods yield similar results in all possible cases. Yet some limitations are still to be considered. First, the current implementation only takes into account the analytical derivation presented in the literature, that is pair interaction, bond angles and dihedral angles. Long-range coulombic interaction and improper dihedral angles contributions are not yet implemented. Also only a limited amount of forcefields have their derivations implemented, but apart from the case of complex non-bonded potentials (such as AIREBO, ADP, MEAM or SW), the implementation of the derivations are straightforward and is an ongoing process. It is to be noted that, as the analytical derivations exist in the literature for example for AIREBO [7] and MEAM [33], the problem is essentially to insure accurate and fast implementation of such potentials, which might not use standard computation algorithms. Plus, while the numerical derivative method can com-

pensate the lack of analytical derivation, this comes with a severe cost in term of computation time. Indeed, the numerical derivative requires the computation of forces 12 times per evaluation (2 times for each of the 6 elementary strain directions). This resulted in a relative increase of 38% and 147% of the computation time in the argon and the polyethylene cases respectively.

We can also consider several follow-up for this work, enabled by the current implementation. For example, the methods could be extended to compute separate atomic contributions, which would be of interest to study heterogeneous or glassy systems local mechanical properties. One could also think about extensions of the formalism to the other ensembles such as the NPT ensemble as suggested by Lips and Maas [34].

## CRediT authorship contribution statement

Germain Clavier: Software, Formal analysis, Writing
Aidan P. Thompson: Software, Validation, Writing

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data availability

No data was used for the research described in the article.

## Acknowledgements

This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525.

## References

[1] D. Squire, A. Holt, W. Hoover, Physica 42 (3) (1969) 388-397.

[2] M. Parrinello, A. Rahman, J. Chem. Phys. 76 (5) (1982) 2662-2666.

[3] J.R. Ray, Comput. Phys. Rep. 8 (3) (1988) 109-151.

[4] G. Clavier, N. Desbiens, E. Bourasseau, V. Lachet, N. Brusselle-Dupend, B. Rousseau, Mol. Simul. 43 (17) (2017) 1413-1422.

[5] J. Lutsko, J. Appl. Phys. 64 (3) (1988) 1152-1154.

[6] J. Lutsko, J. Appl. Phys. 65 (8) (1989) 2991-2997.

[7] K. Van Workum, K. Yoshimoto, J.J. de Pablo, J.F. Douglas, Phys. Rev. E 71 (6) (2005) 061102.

[8] K. Van Workum, G. Gao, J.D. Schall, J.A. Harrison, J. Chem. Phys. 125 (14) (2006) 144506.

[9] G.G. Vogiatzis, L.C. van Breemen, D.N. Theodorou, M. Hütter, Comput. Phys. Commun. 249 (2020) 107008.

[10] A.P. Thompson, S.J. Plimpton, W. Mattson, J. Chem. Phys. 131 (15) (2009) 154107.

[11] K. Yoshimoto, G.J. Papakonstantopoulos, J.F. Lutsko, J.J. de Pablo, Phys. Rev. B 71 (18) (2005) 184108.

[12] Y. Zhen, C. Chu, Comput. Phys. Commun. 183 (2) (2012) 261-265.

[13] W. Hoover, A. Holt, D. Squire, Physica 44 (3) (1969) 437-443.

[14] E. Dobbs, G.O. Jones, Rep. Prog. Phys. 20 (1) (1957) 516.

[15] W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (45) (1996) 11225-11236, https://doi.org/10.1021/ja9621760.

[16] A.I. Jewett, D. Stelter, J. Lambert, S.M. Saladi, O.M. Roscioni, M. Ricci, L. Autin, M. Maritan, S.M. Bashusqeh, T. Keyes, R.T. Dame, J.-E. Shea, G.J. Jensen, D.S. Goodsell, J. Mol. Biol. 433 (11) (2021) 166841.

[17] C.W. Bunn, Trans. Faraday Soc. 35 (1939) 482-491, https://doi.org/10.1039/TF9393500482.

[18] G. Clavier, Étude à l'échelle moléculaire des propriétés mécaniques des polymères semi-cristallins, These, Université Paris Saclay (COmUE), https://tel.archives-ouvertes.fr/tel-01682585, Nov. 2017.
[19] C. Nieto-Draghi, P. Ungerer, B. Rousseau, J. Chem. Phys. 125 (4) (2006) 044517.

[20] C. Nieto-Draghi, P. Bonnaud, P. Ungerer, J. Phys. Chem. C 111 (43) (2007) 15686-15699.

[21] P. Yi, G.C. Rutledge, J. Chem. Phys. 131 (13) (2009) 134902.

[22] M. Matsuo, C. Sawatari, Macromolecules 19 (7) (1986) 2036-2040, https://doi.org/10.1021/ma00161a042.

[23] E. Roguet, K. Akhan, N. Brusselle-Dupend, V. Le Corre, M. Sidhom, L. Cangemi, M. Moreaud, G. Clavier, V. Lachet, B. Rousseau, Comput. Mater. Sci. 167 (2019) 77-84, https://doi.org/10.1016/j.commatsci.2019.05.006.

[24] B. Jelinek, S. Groh, M.F. Horstemeyer, J. Houze, S.G. Kim, G.J. Wagner, A. Moitra, M.I. Baskes, Phys. Rev. B 85 (2012) 245102, https://doi.org/10.1103/PhysRevB.85.245102, https://link.aps.org/doi/10.1103/PhysRevB.85.245102.

[25] F. Apostol, Y. Mishin, Phys. Rev. B 83 (2011) 054116, https://doi.org/10.1103/PhysRevB.83.054116, https://link.aps.org/doi/10.1103/PhysRevB.83.054116.

[26] M.D. Kluge, J.R. Ray, A. Rahman, J. Chem. Phys. 85 (7) (1986) 4028-4031.

[27] E.B. Tadmor, R.S. Elliott, J.P. Sethna, R.E. Miller, C.A. Becker, JOM 63 (7) (2011) 17.

[28] J.L.H. Richard W. Hertzberg, Richard P. Vinci, Deformation and Fracture Mechanics of Engineering Materials, 5th Edition, Wiley, 2012.

[29] K. Albe, K. Nordlund, R.S. Averback, Phys. Rev. B 65 (2002) 195124, https://doi.org/10.1103/PhysRevB.65.195124, https://link.aps.org/doi/10.1103/PhysRevB.65.195124.

[30] P. Erhart, N. Juslin, O. Goy, K. Nordlund, R. Müller, K. Albe, J. Phys. Condens. Matter 18 (29) (2006) 6585.

[31] J. Nord, K. Albe, P. Erhart, K. Nordlund, J. Phys. Condens. Matter 15 (32) (2003) 5649.

[32] T. Bateman, J. Appl. Phys. 33 (11) (1962) 3309-3312.

[33] M. Krief, Y. Ashkenazy, Phys. Rev. E 103 (2021) 063307, https://doi.org/10.1103/PhysRevE.103.063307, https://link.aps.org/doi/10.1103/PhysRevE.103.063307.

[34] D. Lips, P. Maass, Phys. Rev. E 97 (2018) 053002, https://doi.org/10.1103/PhysRevE.97.053002, https://link.aps.org/doi/10.1103/PhysRevE.97.053002.