* Hoekstra  _et al._ (2015) H. Hoekstra, R. Herbonnet, A. Muzzin, A. Babul, A. Mahdavi, M. Viola, and M. Cacciato, MNRAS 449, 685–714 (2015).
  * Bolliet  _et al._ (2018) B. Bolliet, B. Comis, E. Komatsu, and J. F. Macías-Pérez, MNRAS 477, 4957–4967 (2018).
  * Bolliet  _et al._ (2020) B. Bolliet, T. Brinckmann, J. Chluba, and J. Lesgourgues, MNRAS 497, 1332–1347 (2020).
  * Biffi  _et al._ (2016) V. Biffi, S. Borgani, G. Murante, E. Rasia, S. Planelles, G. L. Granato, C. Ragone-Figueroa, A. M. Beck, M. Gaspari, and K. Dolag, ApJ 827, 112 (2016), arXiv:1606.02293 [astro-ph.CO] .
  * Shi and Komatsu (2014) X. Shi and E. Komatsu, MNRAS 442, 521 (2014), arXiv:1401.7657 [astro-ph.CO] .
  * Hikage  _et al._ (2019) C. Hikage, M. Oguri, T. Hamana, S. More, R. Mandelbaum, M. Takada, F. Köhlinger, H. Miyatake, A. J. Nishizawa, H. Aihara, and et al., Publications of the Astronomical Society of Japan 71 (2019), 10.1093/pasj/psz010.
  * Asgari  _et al._ (2021) M. Asgari, C.-A. Lin, B. Joachimi, B. Giblin, C. Heymans, H. Hildebrandt, A. Kannawadi, B. Stölzner, T. Tröster, J. L. van den Busch, and et al., A&A 645, A104 (2021).
  * Hilton  _et al._ (2018) M. Hilton, M. Hasselfield, C. Sifón, N. Battaglia, S. Aiola, V. Bharadwaj, J. R. Bond, S. K. Choi, D. Crichton, R. Datta, and et al., The Astrophysical Journal Supplement Series 235, 20 (2018).
  * Hilton  _et al._ (2021) M. Hilton, C. Sifón, S. Naess, M. Madhavacheril, M. Oguri, E. Rozo, E. Rykoff, T. M. C. Abbott, S. Adhikari, M. Aguena, and et al., The Astrophysical Journal Supplement Series 253, 3 (2021).
  * Hurier and Lacasa (2017) G. Hurier and F. Lacasa, A&A 604, A71 (2017).
  * Smith  _et al._ (2016) G. P. Smith, P. Mazzotta, N. Okabe, F. Ziparo, S. L. Mulroy, A. Babul, A. Finoguenov, I. G. McCarthy, M. Lieu, Y. M. Bahé, H. Bourdin, A. E. Evrard, T. Futamase, C. P. Haines, M. Jauzac, D. P. Marrone, R. Martino, P. E. May, J. E. Taylor, and K. Umetsu, MNRAS 456, L74 (2016), arXiv:1511.01919 [astro-ph.CO] .
  * Eckert  _et al._ (2019) D. Eckert, V. Ghirardini, S. Ettori, E. Rasia, V. Biffi, E. Pointecouteau, M. Rossetti, S. Molendi, F. Vazza, F. Gastaldello, M. Gaspari, S. De Grandi, S. Ghizzardi, H. Bourdin, C. Tchernin, and M. Roncarelli, A&A 621, A40 (2019), arXiv:1805.00034 [astro-ph.CO] .
  * Hill and Spergel (2014) J. C. Hill and D. N. Spergel, JCAP 2014, 030–030 (2014).
  * Ma  _et al._ (2015) Y.-Z. Ma, L. V. Waerbeke, G. Hinshaw, A. Hojjati, D. Scott, and J. Zuntz, JCAP 2015, 046–046 (2015).
  * Rotti  _et al._ (2021) A. Rotti, B. Bolliet, J. Chluba, and M. Remazeilles, MNRAS 503, 5310 (2021), arXiv:2010.07797 [astro-ph.CO] .
  * Barnes  _et al._ (2021) D. J. Barnes, M. Vogelsberger, F. A. Pearce, A.-R. Pop, R. Kannan, K. Cao, S. T. Kay, and L. Hernquist, MNRAS (2021), 10.1093/mnras/stab1276.

##  Appendix A Covariance matrix

Our full model of theory covariance, including the Gaussian and non-Gaussian
terms is shown is Eq. 34. In paper I we validated the Gaussian part of our
total covariance using simulations. We have also compared it to the jackknife
covariance estimate which partly captures the non-Gaussian contribution to the
total covariance. Our total covariance includes the contribution from poisson
fluctuations of large clusters.

In Fig. 12 we show the part of the correlation matrix using fourth source
tomographic bin. It clearly shows that due to large beam, the small scale
angular bins corresponding to $\theta<10$arcmin are more correlated in the
$\textit{Planck}\times\textsc{DES}$ part of the matrix compared to
$\textsc{ACT}\times\textsc{DES}$.

Figure 12: Correlation matrix of $\xi_{\gamma_{t}y}$ using the fourth source
bin and the two Compton-$y$ maps, binned into 20 radial bins from 2.5 to 250
arcmin.

##  Appendix B Fits with Arnaud10 model

In Fig. 13 we compare the best fits obtained from the models based on Arnaud
_et al._ (2010) with the one obtained from the Battaglia  _et al._ (2012)
model (as shown in Fig. 4). We find that all three models result in similar
goodness of fit. The PTE for A10c model is 0.02 and for A10z model is 0.0198.

Figure 13: This figure is simiar to Fig. 4 but comparing the best fit to the
data obtained from the Arnaud  _et al._ (2010) model and the Battaglia  _et
al._ (2012) model. We show the best fit for the case of fixing $\rho_{B}=0$
(dashed curve, A10c model) and free $\rho_{B}$ (dot-dashed curve, A10z model).
We refer the reader to Table 1 for details of the model parameters and priors
used in the analyses.

##  Appendix C Impact of assumed cosmological model on parameter constraints

We repeat our analysis adopting the best-fit cosmological parameters from
Aghanim  _et al._ (2020) and from the DES Year 1 analysis of (Abbott  _et al._
, 2018). The full posteriors for these two analyses are shown in Fig. 14. We
find that our results are largely insensitive to the choice of cosmological
model.

Figure 14: Constraints on the pressure profile parameters from the break model
assuming two different cosmological models.


# Calculation of elastic constants of embedded-atom-model potentials in the
NVT ensamble

Menahem Krief<EMAIL_ADDRESS>Yinon Ashkenazy  Racah Institute
of Physics, The Hebrew University, 9190401 Jerusalem, Israel

###### Abstract

A method for the calculation of elastic constants in the NVT ensamble, using
molecular dynamics (MD) simulation with a realistic many-body embedded-atom-
model (EAM) potential, is studied in detail. It is shown that in such NVT MD
simulations, the evaluation of elastic constants is robust and accurate, as it
gives the elastic tensor in a single simulation which converges using a small
number of time steps and particles. These results highlight the applicability
of this method in: (i) the calculation of local elastic constants of non-
homogeneous crystalline materials and (ii) in the calibration of interatomic
potentials, as a fast and accurate alternative to the common method of
explicit deformation, which requires a set of consistent simulations at
different conditions. The method is demonstrated for the calculation of the
elastic constants of copper in the temperature range of 0-1000K, and results
agree with the target values used for the potential calibration. The various
contributions to the values of the elastic constants, namely, the Born, stress
fluctuation and ideal gas terms, are studied as a function of temperature.

##  I Introduction

The calculation of thermo-elastic properties of materials using computer
simulations, plays a key role in understanding the response of materials to
deformation under varying conditions. In a wide range of applications,
material structure leads to local variation in the elastic response functions,
due to grain boundaries and hetero-phase interfaces. It is this local
variation of elastic properties that allow tailoring of the macroscopic
effective material properties. The development of a model for the relation
between composite material specific local property and an effective average
response function depends on understanding how local properties contribute and
interact in order to generate an observable average response. Furthermore,
local elastic properties are not accessible experimentally for a wide range of
systems, and so reliable numerical evaluation of these may play a key role in
the development of effective models for composite materials.

While it is well established that various thermal and mechanical properties
can be evaluated for various atomistic structures by atomistic simulations, it
is also well established that such evaluation requires addressing non pairwise
terms. The widely used embedded atom model (EAM)[1], is a fast, simple and
accurate method, which allows a correct description of the microscopic
interactions in crystalline materials. Molecular dynamics calculations of
adiabatic elastic constants (which are performed in the microcanonical NVE
ensamble), for metallic elements, using EAM model potentials, were performed
in the past [2, 3, 4], based on the widely used formulation of Ray et al. [5].

In this manuscript we present and analyze the feasibility, robustness and
accuracy of the calculation of elastic constants of metals under constant
temperature and volume, that is, in the canonical NVT ensamble. Using the
stress-stress formulation [6, 7, 8, 9, 4], all the components of the elastic
constant tensor are obtained in a single molecular dynamics simulation, as
opposed to the common explicit deformation method [10, 8, 11, 12, 13, 14, 15,
16], which requires several simulations under different deformations. In
addition, the method allows the evaluation of elastic properties in localized
regions within larger non-homogeneous simulation box, unlike first principle
methods which allow efficient temperature dependent calculations, but only for
a uniform system [17, 18]. Finally, it is demonstrated that NVT calculations
of elastic constants converge more rapidly in comparison to the NPT strain-
strain fluctuations formulation [19, 20, 8, 9, 21], using symplectic numerical
integrators [22, 23, 24, 25, 26, 27] and Nose-Hoover thermostat chain [28, 29,
30].

We performed calculations of the isothermal elastic constants
$C_{11},C_{12},C_{44}$ of Copper for a widely used realistic tabulated EAM
potential by Mishin et al. [31]. This potential successfully reproduces energy
and stability of several nonequilibrium configurations as well as
transformation paths between different structures. The calculations were
compared with experimental results in the temperature range 0-1000K, and a
good agreement is achieved. The manuscript is structured as follows: We start
with a short review of formalism used from thermoelasticity, then we describe
in detail how elastic constants can be calculated using a single MD simulation
in the NVT ensamble in Section 3, and conclude with Section 4 in which we
employ the described method to calculate elastic constants in Copper for a
commonly used EAM potential.

##  II Thermoelasticity

We start with a brief review of the definitions and notations of thermo-
elasticity, that will be used throughout the manuscript. Given a deformation
from a reference configuration $\boldsymbol{R}$ to a configuration
$\boldsymbol{r}=\boldsymbol{r}\left(\boldsymbol{R}\right)$, the Lagrangian
strain tensor is defined by:

$\eta_{\alpha\beta}=\frac{1}{2}\left(\frac{\partial r_{\gamma}}{\partial
R_{\alpha}}\frac{\partial r_{\gamma}}{\partial
R_{\beta}}-\delta_{\alpha\beta}\right),$ (1)

where Greek indices such as $\alpha,\beta,\gamma$ represent Cartesian
components, for which we employ Einstein summation convention. In matrix
notation, the strain tensor is related to the Jacobian (metric) tensor:

$J_{\alpha\beta}=\frac{\partial r_{\alpha}}{\partial R_{\beta}},$ (2)

by:

$\boldsymbol{\eta}=\frac{1}{2}\left(\boldsymbol{J}^{T}\boldsymbol{J}-\boldsymbol{I}\right),$
(3)

where $\boldsymbol{I}$ is the identity matrix. The difference in dimensions
due to deformation can thus be written using the strain tensor 1 via:

$\delta r^{2}-\delta R^{2}=2\eta_{\alpha\beta}\delta R_{\alpha}\delta
R_{\beta}.$ (4)

And the mechanical work due to an infinitesimal deformation with Lagrangian
strain $d\eta_{\alpha\beta}$ is given by [32, 33, 34, 35, 36, 37]:

$\delta W/V_{0}=-\tau_{\alpha\beta}d\eta_{\alpha\beta},$ (5)

where $V_{0}$ is the undeformed volume and the thermodynamic tension tensor,
also known as the second Piola-Kirchhoff stress tensor [32], is defined by:

$\boldsymbol{\tau}=\det\left(\boldsymbol{J}\right)\boldsymbol{J}^{-1}\boldsymbol{\sigma}\boldsymbol{J}^{-T},$
(6)

where $\boldsymbol{\sigma}$ is the Cauchy stress tensor, and the ’$-T$’
superscript denotes matrix transposition and inversion. We note that
$\boldsymbol{\tau}=\boldsymbol{\sigma}$ for a zero applied strain
($\boldsymbol{J=}\boldsymbol{I}$, $\boldsymbol{\eta=0}$). If the deformation
process is reversible (an assumption which excludes plasticity), then the
first law of thermodynamics for the internal energy $\mathcal{E}$ is written
as:

$d\mathcal{E}=Td\mathcal{S}+V_{0}\tau_{\alpha\beta}d\eta_{\alpha\beta},$ (7)

where $T$ denotes the temperature and $\mathcal{S}$ the entropy. And changes
in the free energy

$d\mathcal{F}=\mathcal{S}dT+V_{0}\tau_{\alpha\beta}d\eta_{\alpha\beta}.$ (8)

The thermodynamic tension tensor is the thermodynamic conjugate of the
Lagrangian strain tensor:

$\tau_{\alpha\beta}=\frac{1}{V_{0}}D_{\alpha\beta}\mathcal{F},$ (9)

where we used the symmetrical partial derivative operator:

$D_{\alpha\beta}=\frac{1}{2}\left(\frac{\partial}{\partial\eta_{\alpha\beta}}+\frac{\partial}{\partial\eta_{\beta\alpha}}\right),$
(10)

due to symmetry of the strain tensor $d\eta_{\alpha\beta}=d\eta_{\beta\alpha}$
in eq. 8. The Taylor expansion of the elastic free energy around the un-
deformed state, $\boldsymbol{\eta=0}$, is written as:

$\mathcal{F}\left(\boldsymbol{\eta}\right)/V_{0}=\mathcal{F}\left(\boldsymbol{0}\right)/V_{0}+\sigma_{\alpha\beta}\eta_{\alpha\beta}+\frac{1}{2}C_{\alpha\beta\gamma\delta}\eta_{\alpha\beta}\eta_{\gamma\delta}+...$
(11)

where:

$\sigma_{\alpha\beta}=\tau_{\alpha\beta}\left(\boldsymbol{\eta=0}\right)=\frac{1}{V_{0}}D_{\alpha\beta}\mathcal{F}\bigg{|}_{\boldsymbol{\eta}=\boldsymbol{0}},$
(12)

and the second order isothermal elastic constants are defined by:

$C_{\alpha\beta\gamma\delta}=\frac{1}{V_{0}}D_{\alpha\beta}D_{\gamma\delta}\mathcal{F}\bigg{|}_{\boldsymbol{\eta}=\boldsymbol{0}},$
(13)

where the derivatives in equations 12-13 are taken at constant $T$ and
evaluated at zero strain $\boldsymbol{\eta}=\boldsymbol{0}$.

##  III molecular dynamics and Elastic constants

In this section we briefly outline how the entire isothermal elasticity tensor
can be obtained in a single molecular dynamics simulation using the stress-
stress fluctuation method [36, 5, 38, 2, 8, 9]. This method is favorable over
the widely used explicit deformation method [8, 11, 12, 13, 14, 15, 16], which
requires multiple simulations under various deformations, in order to obtain
different components of the elasticity tensor.

Starting from the Hamiltonian of a system of $N$ particles is written as:

$\displaystyle\mathcal{H}\left(\boldsymbol{r}^{N},\boldsymbol{p}^{N}\right)$
$\displaystyle=\sum_{i}\frac{p_{i}^{2}}{2m_{i}}+\mathcal{V}\left(\boldsymbol{r}_{1},...,\boldsymbol{r}_{N}\right),$
(14)

where $m_{i}$,$\boldsymbol{r}^{N}$ and $\boldsymbol{p}^{N}$ denotes,
respectively, the masses, position and momentum vectors of the $N$ particles.
The embedded-atom-model (EAM) potential [1], is defined by a pair potential
$v=v\left(r\right)$, an embedding function $F=F\left(\rho\right)$ and a local
electron density function $\rho=\rho\left(r\right)$, so that the the potential
energy takes the form:

$\mathcal{V}\left(\boldsymbol{r}_{1},...,\boldsymbol{r}_{N}\right)=\sum_{i}F\left(\rho_{i}\right)+\sum_{i<j}v\left(r_{ij}\right),$
(15)

where $\boldsymbol{r}_{ij}=\boldsymbol{r}_{i}-\boldsymbol{r}_{j}$ is the
interatomic displacement vector and the electron density function of particle
$i$ is given by $\rho_{i}=\sum_{j\neq i}\rho\left(r_{ij}\right).$ For an EAM
potential of the form 15, the force on a particle $i$ can be written as:

$\boldsymbol{F}_{i}=\sum_{j\neq i}F_{ij}\frac{\boldsymbol{r}_{ij}}{r_{ij}},$
(16)

where:

$F_{ij}=-\left(v^{\prime}\left(r_{ij}\right)+\left[F^{\prime}\left(\rho_{i}\right)+F^{\prime}\left(\rho_{j}\right)\right]\rho^{\prime}\left(r_{ij}\right)\right).$
(17)

The instantaneous kinetic temperature is given by:

$\frac{1}{2}k_{B}T_{K}=\frac{1}{g}\sum_{i}\frac{p_{i}^{2}}{2m_{i}},$ (18)

where $k_{B}$ is Boltzmann’s constant and $g=3\left(N-1\right)$ is the number
of degrees of freedom. The instantaneous pressure tensor $P_{\alpha\beta}$ is
obtained from the virial theorem:

$\displaystyle P_{\alpha\beta}V$
$\displaystyle=\sum_{i}\left(\frac{p_{i,\alpha}p_{i,\beta}}{m_{i}}+r_{i,\alpha}F_{i,\beta}\right).$

Where $V$ is the system volume. The total pressure is given by the average of
the diagonal components:

$P=\frac{1}{3}\left(P_{11}+P_{22}+P_{33}\right).$ (19)

For an EAM potential, it follows from eq. 16 that:

$\displaystyle\sum_{i}r_{i,\alpha}F_{i,\beta}$
$\displaystyle=\sum_{i<j}F_{ij}\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}}$

so that the pressure tensor has the form:

$\displaystyle P_{\alpha\beta}V$
$\displaystyle=\sum_{i}\frac{p_{i,\alpha}p_{i,\beta}}{m_{i}}+\sum_{i<j}F_{ij}\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}}$
(20)

It can be shown [6, 7, 3, 39, 40, 8, 9], as detailed in Appendix A, that in
the NVT ensamble, the elastic constants 13 are given by:

$\displaystyle C_{\alpha\beta\gamma\delta}$ $\displaystyle=\left\langle
C_{\alpha\beta\gamma\delta}^{B}\right\rangle-\frac{V}{k_{B}T}\left[\left\langle\sigma_{\alpha\beta}^{B}\sigma_{\gamma\delta}^{B}\right\rangle-\left\langle\sigma_{\alpha\beta}^{B}\right\rangle\left\langle\sigma_{\gamma\delta}^{B}\right\rangle\right]$
$\displaystyle+\frac{Nk_{B}T}{V}\left(\delta_{\alpha\gamma}\delta_{\beta\delta}+\delta_{\alpha\delta}\delta_{\beta\gamma}\right),$
(21)

where $\left\langle\cdot\right\rangle$ represents ensemble average. The last
term in eq. 21 is a non configurational ideal gas contribution which vanishes
at zero temperature and is related to volume derivatives with respect to the
strain tensor. The first term in eq. 21, known as the Born term, is a
configurational part which is given by a canonical average of the zero-
temperature expression for the elastic constant [41, 8, 5]. The second term in
eq. 21 accounts for stress fluctuations and also vanishes at zero temperature
(this term is obtained directly from the general identity 39 detailed in
Appendix A). The Born stress tensor in eq. 21 is defined by the derivative of
the potential energy with respect to strain, evaluated at a state of zero
applied strain:

$\sigma_{\alpha\beta}^{B}=\frac{1}{V}D_{\alpha\beta}\mathcal{V}\bigg{|}_{\boldsymbol{\eta}=\boldsymbol{0}}.$
(22)

The Cauchy stress tensor is given by:

$\sigma_{\alpha\beta}=\left\langle\sigma_{\alpha\beta}^{B}\right\rangle-\frac{k_{B}TN}{V}\delta_{\alpha\beta}$
(23)

Similarly, the Born elastic constant term is given by:

$C_{\alpha\beta\gamma\delta}^{B}=\frac{1}{V}D_{\alpha\beta}D_{\gamma\delta}\mathcal{V}\bigg{|}_{\boldsymbol{\eta}=\boldsymbol{0}}.$
(24)

For an EAM potential of the form 15, it can be shown that the Born stress has
the following form:

$V\sigma_{\alpha\beta}^{B}=-\sum_{i<j}F_{ij}\frac{r_{ij,\alpha}r_{ij,\beta}}{r_{ij}},$
(25)

while the Born elastic constant takes the form [7, 2, 4, 3]:

$\displaystyle VC_{\alpha\beta\gamma\delta}^{B}$
$\displaystyle=\sum_{i<j}X_{ij}\frac{r_{ij,\alpha}r_{ij,\beta}r_{ij,\gamma}r_{ij,\delta}}{r_{ij}^{2}}$
$\displaystyle+\sum_{i}F^{\prime\prime}\left(\rho_{i}\right)g_{i,\alpha\beta}g_{i,\gamma\delta},$
(26)

where:

$\displaystyle X_{ij}$
$\displaystyle=v^{\prime\prime}\left(r_{ij}\right)-\frac{1}{r_{ij}}v^{\prime}\left(r_{ij}\right)$
$\displaystyle+\left(F^{\prime}\left(\rho_{i}\right)+F^{\prime}\left(\rho_{j}\right)\right)\left(\rho^{\prime\prime}\left(r_{ij}\right)-\frac{1}{r_{ij}}\rho^{\prime}\left(r_{ij}\right)\right),$
(27)

and:

$g_{i,\alpha\beta}=\sum_{k\neq
i}\rho^{\prime}\left(r_{ik}\right)\frac{r_{ik,\alpha}r_{ik,\beta}}{r_{ik}}.$
(28)

For completeness, a detailed derivation of equations 21-28, is given in
Appendix A.

##  IV Results

Figure 1: The thermal expansion ratio $L\left(T\right)/L\left(0\right)-1$ of
Copper, as a function of temperature, resulting from a series of NPT molecular
dynamics simulations at zero pressure in the temperature range 0-1000K. Time
dependent measures of the simulation at T=300K are presented in detail in
figures 2-3.

Figure 2: (Color online) Analysis of a molecular dynamics simulation of Copper
in the NPT ensemble at zero pressure and $T=300\text{K}$. The upper left
figure shows the instantaneous pressure (eq. 19, red solid line), the
cumulative average pressure (black dashed line) and the cumulative pressure
standard deviation. The upper middle figure shows the instantaneous kinetic
temperature (eq. 18, red solid line, left y-axis), the cumulative average
temperature (black dashed line, left y-axis), the relative error between the
current cumulative average value to the the final average value (green dotted
line) and the ratio between the cumulative standard deviation to the
cumulative average (blue solid line, right axis). Similarly, the upper right
figure shows the results for the ratio between the simulation box length and
the initial box length ($L_{0}=3.165\mathrm{\mathring{A}}$). The first
$2\times 10^{4}$ steps are plotted on a logarithmic x scale, in order to show
the initial thermalization period. The lower figures show histograms and
fitted Gaussian distributions for the values of the instantaneous pressure
(lower left), temperature (middle figure) and box size ratio (lower right).
Figure 3: (Color online) Components of the instantaneous pressure tensor (eq.
18) for the NPT simulation described in Fig. 2.

Figure 4: (Color online) Components of the instantaneous Born stress tensor
(upper figure, eq. 25) for an NVT simulation of Copper at $T=300$K for which
the volume is chosen such that the total pressure is zero (as obtained from
Fig. 1). The lower figure shows the resulting histograms and fitted Gaussian
distributions of the various Born stress tensor components.

Figure 5: (Color online) Upper figures - the Born elastic constants (eq. 26):
$C_{11}^{B}$ (left upper figure), $C_{12}^{B}$ (middle upper figure) and
$C_{44}^{B}$ (right upper figure), for an NVT simulation of Copper at $T=300$K
(as described in Fig. 4). The instantaneous value, cumulative average,
convergence error and ratio of the standard deviation to the mean are shown,
as is detailed in Fig. 2. The lower figures shows the resulting histograms and
fitted Gaussian distributions.

Figure 6: Relative error between the instantaneous conserved energy and the
initial value, for NPT (upper figure, eq. 60) and NVT (upper figure, eq. 75)
molecular dynamics simulations of Copper at $T=300$K and zero pressure (the
simulations described in figures 2-3 and 4-5, respectively).

Figure 7: (Color online) The different terms contributing to the total
isothermal elastic constants $C_{11}$ (left pane), $C_{12}$ (middle pane) and
$C_{44}$ (right pane) of Copper as a function of temperature (solid blue line)
- the Born term (first term in eq. 21, dashed orange line), Born stress
fluctuation term (minus the second term in eq. 21, dashed dotted green line)
and the ideal gas kinetic term (minus the last term in eq. 21, dotted red
line).  Figure 8: (Color online) The isothermal elastic constants $C_{11}$
(black), $C_{12}$ (red) and $C_{44}$ (blue) of Copper as a function of
temperature. Solid lines represent molecular dynamics calculations performed
in this work (see text for details), circles represent experimental values in
the temperature range 0-300K given in Ref. [42] and squares represent
experimental values in the range 300-800K given in Ref. [43] for $C_{44}$. The
blue X data point represents a calculated $C_{44}$ result given in Ref. [10].

Detailed molecular dynamics simulations are performed, in order to demonstrate
and analyze the formalism developed in the previous sections, which we employ
for the calculation of the elasticity tensor of copper.

The results presented in this work were obtained using a newly developed
molecular dynamics code. The computational model employs symplectic numerical
integrators, which preserve phase space measures of non-Hamiltonian dynamics
[22, 23, 24, 25, 26, 27], incorporates Nose-Hoover thermostat chains in the
NVT and NPT ensembles [28, 29, 30], solved usin high order Suzuki-Yoshida
decomposition [44, 45]. The equations of motions are detailed in appendices B
and C, for the NVT and NPT ensembles, respectively.

The isothermal elastic constants of Copper were calculated as a function of
temperature using a series of molecular dynamics calculations in the NVT
ensemble employing a widely used EAM potential by Mishin et al. [31]. This
potential is calibrated so that the calculated zero temperature fcc lattice
constant is $L\left(T=0\right)=3.615\mathrm{\mathring{A}}$ and the elastic
constants are given by $C_{11}=169.9\text{GPa}$, $C_{12}=122.6\text{GPa}$ and
$C_{44}=76.2\text{GPa}$ (the Voigt notation will be used throughout this
section). In the simulations presented below, a cubical system with 500 atoms
was used with their initial positions on an fcc lattice, and initial
velocities sampled from a Maxwell distribution at the appropriate temperature.
Periodic boundary conditions were used. The integration of the equations of
motion was performed with a time step of $2\text{fs}$ for $2\times 10^{5}$
steps. A Nose-Hoover thermostat chain of size $M=10$ and a relaxation time of
$\tau_{T}=50\text{fs}$ was applied.

The NVT calculations require prior evaluation of the equilibrium volume at
zero stress at the various temperatures. These volumes can be evaluated using
known thermal expansion coefficient, but we re-evaluated them using a series
of zero stress NPT calculations. These yield the lattice constant as a
function of temperature, $L\left(T\right)$. The resulting equilibrium thermal
expansion ratio, $L\left(T\right)/L\left(0\right)-1$, is shown as a function
of temperature in Fig. 1. An analysis of the thermalization and convergence of
a selected NPT simulation at $T=300$K is presented in detail in figures 2-3.
It is evident that convergence is achieved relatively fast with $N=10^{4}$ MD
steps. The instantaneous pressure tensor is shown in Fig. 3, demonstrating
that the pressure is equilibrated to zero. The NPT simulations were
initialized with an fcc lattice with the zero temperature lattice constant
$L\left(T=0\right)=3.615\mathrm{\mathring{A}}$, which is then subsequently
relaxed to the equilibrium value $L\left(T\right)$ at zero pressure (as seen
in Fig. 2). A barostat, coupled to a Nose-Hoover thermostat chain of size
$M^{\prime}=10$ with a relaxation time of $\tau_{P}=500\text{fs}$ was used.

Results of a specific NVT simulation with $T=300$K and a volume chosen such
that the total pressure is zero, are shown in figures 4-5. Fig. 4 shows the
components of the Born stress tensor (defined in eq. 25) throughout the
simulation. As expected, only the diagonal ideal gas terms are nonzero (since
the total pressure is zero, from eq. 23 we have
$\left\langle\sigma_{\alpha\beta}^{B}\right\rangle=\delta_{\alpha\beta}k_{B}TN/V$).
Fig. 5 shows the Born elastic constants $C_{11}^{B},C_{12}^{B},C_{44}^{B}$
given in eq. 26, as in Fig. 2. It is evident that the Born elastic constants
are converged relatively fast, at about $10^{4}$ MD steps.

In Fig. 6 we show the conservation of the energy like quantities for both the
NPT (given in eq. 75) and NVT calcualtions and NVT (given in eq. 60)
throughout the simulations discussed above. It is evident that a good
conservation of about $10^{-5}$ is maintained throughout the simulation, which
highlights the great advantages of measure preserving integrators of the non-
Hamltonian dynamical systems [23, 24].

In Fig. 7 the isothermal elastic constants as a function of temperature are
plotted, together with the contributions of the Born term, the Born stress
fluctuation term and the ideal gas kinetic term (see eq. 21). As expected, it
is evident that the reduction of the elastic constant with increasing
temperature is mainly due to the stress fluctuations term, which is larger at
higher temperatures. It is also evident that in the temperature range studied
here, the ideal gas contribution is negligible.

Finally, in Fig. 8, we compare the calculated isothermal elastic constants as
a function of temperature, with experimental results. In Ref. [42],
experimental results for isothermal elastic constants $C_{11},C_{12},C_{44}$
are given in the temperature range 0-300K. In Ref. [43], experimental results
for adiabatic elastic constants are given in the range 300-800K. Since the
isothermal and adiabatic elastic constants are identical for $C_{44}$ , we can
use these experimental results directly. On the other hand, the isothermal and
adiabatic elastic constants $C_{11},C_{12}$ are not identical [46, 47, 48][46,
47, 48], and were not compared here. It is evident that the agreement between
the calculated and experimental values is relatively good, especially given
the fact that the Copper potential we used [31] was calibrated to slightly
different values at $T=0$K (the experimental values of Ref. [42] are
$C_{11}=176.2\text{GPa}$, $C_{11}=124.94\text{GPa}$ and
$C_{11}=81.77\text{GPa}$). The values that were used to calibrate the
tabulated potential (at $T=0)$, are reproduced by our calculations to 5
significant digits. In addition, in Ref. [10] adiabatic elastic constants are
calculated by molecular dynamics simulations (using LAMMPS MD code [49]) at
$T=300$K, using the same EAM potential and employing the explicit deformation
method. This result for $C_{44}$ is also shown in Fig. 8, showing a very good
agreement with our calculation (less than 0.5%).

and at results calculated with the same potential (and using LAMMPS [49]), The
blue X data point represents the computational $C_{44}$ result given in Ref.
[10]

##  V Summary

The elasticity tensor of copper, modeled with a realistic tabulated EAM
potential, was calculated using a single molecular dynamics simulation in the
NVT ensamble, employing the stress-stress fluctuation formulation. It was
shown that such calculations are accurate and robust and converge within a few
thousand MD steps, with a relatively small number of particles (a few
hundreds). The calculations were performed in the temperature range of 0-1000K
and compared to experimental values, showing a good quantitative and
qualitative agreement. The various thermal contributions to the values of the
elastic constants were studied as a function of temperature.

The results suggest that this method can be applied to calculate local elastic
constants of real crystalline materials for local regions embedded within
large crystalline structure. In addition, it can improve the calibration
process of inter-atomic potentials, which typically employ the explicit
deformation method [10, 16, 8, 11, 12, 13, 14, 15], that requires the
generation of a consistent set of simulations under varying deformations, in
order to obtain the elasticity tensor. The use of a single simulation, other
than being simpler, convenient and accurate, has a lower risk for the
spontaneous appearance of defects and other phases during direct force
calculations.

## References

  * [1] Murray S Daw and Michael I Baskes.  Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals.  Physical Review B, 29(12):6443, 1984.
  * [2] Ralph J Wolf, Khalid A Mansour, Myung W Lee, and John R Ray.  Temperature dependence of elastic constants of embedded-atom models of palladium.  Physical Review B, 46(13):8027, 1992.
  * [3] T Çağın, G Dereli, M Uludoğan, and MEHMET Tomak.  Thermal and mechanical properties of some fcc transition metals.  Physical review B, 59(5):3468, 1999.
  * [4] Somchart Chantasiriwan and Frederick Milstein.  Higher-order elasticity of cubic metals in the embedded-atom method.  Physical Review B, 53(21):14080, 1996.
  * [5] John R Ray, Michael C Moody, and Aneesur Rahman.  Molecular dynamics calculation of elastic constants for a crystalline system in equilibrium.  Physical Review B, 32(2):733, 1985.
  * [6] DR Squire, AC Holt, and WG Hoover.  Isothermal elastic constants for argon. theory and monte carlo calculations.  Physica, 42(3):388–397, 1969.
  * [7] JF Lutsko.  Generalized expressions for the calculation of elastic constants by computer simulation.  Journal of applied physics, 65(8):2991–2997, 1989.
  * [8] Germain Clavier, Nicolas Desbiens, Emeric Bourasseau, Véronique Lachet, Nadège Brusselle-Dupend, and Bernard Rousseau.  Computation of elastic constants of solids using molecular simulation: comparison of constant volume and constant pressure ensemble methods.  Molecular Simulation, 43(17):1413–1422, 2017.
  * [9] Dominik Lips and Philipp Maass.  Stress-stress fluctuation formula for elastic constants in the npt ensemble.  Physical Review E, 97(5):053002, 2018.
  * [10] Seyed Moein Rassoulinejad-Mousavi, Yijin Mao, and Yuwen Zhang.  Evaluation of copper, aluminum, and nickel interatomic potentials on predicting the elastic properties.  Journal of Applied Physics, 119(24):244304, 2016.
  * [11] Michael Griebel and Jan Hamaekers.  Molecular dynamics simulations of the elastic moduli of polymer–carbon nanotube composites.  Computer methods in applied mechanics and engineering, 193(17-20):1773–1788, 2004.
  * [12] DJ Quesnel, DS Rimai, and LP DeMejo.  Elastic compliances and stiffnesses of the fcc lennard-jones solid.  Physical Review B, 48(10):6795, 1993.
  * [13] Oleg L Manevitch and Gregory C Rutledge.  Elastic properties of a single lamella of montmorillonite by molecular dynamics simulation.  The Journal of Physical Chemistry B, 108(4):1428–1435, 2004.
  * [14] Priya Vashishta, Rajiv K Kalia, Aiichiro Nakano, and José Pedro Rino.  Interaction potential for silicon carbide: A molecular dynamics study of elastic constants and vibrational density of states for crystalline and amorphous silicon carbide.  Journal of applied physics, 101(10):103515, 2007.
  * [15] Qing-Xiang Pei, Yong-Wei Zhang, and Vivek B Shenoy.  Mechanical properties of methyl functionalized graphene: a molecular dynamics study.  Nanotechnology, 21(11):115709, 2010.
  * [16] Saaketh Desai and Alejandro Strachan.  Lammps driver tool for potential calibration.  2019\.
  * [17] L Huang, Levente Vitos, SK Kwon, Börje Johansson, and Rajeev Ahuja.  Thermoelastic properties of random alloys from first-principles theory.  Physical Review B, 73(10):104203, 2006.
  * [18] Philipp Keuter, Denis Music, Volker Schnabel, Michael Stuer, and Jochen M Schneider.  From qualitative to quantitative description of the anomalous thermoelastic behavior of v, nb, ta, pd and pt.  Journal of Physics: Condensed Matter, 31(22):225402, 2019.
  * [19] M Parrinello and A Rahman.  Strain fluctuations and elastic constants.  The Journal of Chemical Physics, 76(5):2662–2666, 1982.
  * [20] Andrei A Gusev, Marcel M Zehnder, and Ulrich W Suter.  Fluctuation formula for elastic constants.  Physical Review B, 54(1):1, 1996.
  * [21] Wataru Shinoda, Motoyuki Shiga, and Masuhiro Mikami.  Rapid estimation of elastic constants by molecular dynamics simulation under constant stress.  Physical Review B, 69(13):134103, 2004.
  * [22] Glenn J Martyna, Douglas J Tobias, and Michael L Klein.  Constant pressure molecular dynamics algorithms.  The Journal of chemical physics, 101(5):4177–4189, 1994.
  * [23] Glenn J Martyna, Mark E Tuckerman, Douglas J Tobias, and Michael L Klein.  Explicit reversible integrators for extended systems dynamics.  Molecular Physics, 87(5):1117–1157, 1996.
  * [24] Mark E Tuckerman, Yi Liu, Giovanni Ciccotti, and Glenn J Martyna.  Non-hamiltonian molecular dynamics: Generalizing hamiltonian phase space principles to non-hamiltonian systems.  The Journal of Chemical Physics, 115(4):1678–1702, 2001.
  * [25] Mark E Tuckerman, José Alejandre, Roberto López-Rendón, Andrea L Jochim, and Glenn J Martyna.  A liouville-operator derived measure-preserving integrator for molecular dynamics simulations in the isothermal–isobaric ensemble.  Journal of Physics A: Mathematical and General, 39(19):5629, 2006\.
  * [26] Mark Tuckerman.  Statistical mechanics: theory and molecular simulation.  Oxford university press, 2010.
  * [27] Michael P Allen and Dominic J Tildesley.  Computer simulation of liquids.  Oxford university press, 2017.
  * [28] Glenn J Martyna, Michael L Klein, and Mark Tuckerman.  Nosé–hoover chains: The canonical ensemble via continuous dynamics.  The Journal of chemical physics, 97(4):2635–2643, 1992.
  * [29] Shuichi Nosé.  A unified formulation of the constant temperature molecular dynamics methods.  The Journal of chemical physics, 81(1):511–519, 1984.
  * [30] William G Hoover.  Canonical dynamics: Equilibrium phase-space distributions.  Physical review A, 31(3):1695, 1985.
  * [31] Yu Mishin, MJ Mehl, DA Papaconstantopoulos, AF Voter, and JD Kress.  Structural stability and lattice defects in copper: Ab initio, tight-binding, and embedded-atom calculations.  Physical Review B, 63(22):224106, 2001.
  * [32] RN Thurston.  Wave propagation in fluids and normal solids.  In Physical acoustics, pages 1–110. Elsevier, 1964.
  * [33] Duane C Wallace.  Thermoelasticity of stressed materials and comparison of various elastic constants.  Physical Review, 162(3):776, 1967.
  * [34] K Brugger.  Thermodynamic definition of higher order elastic coefficients.  Physical Review, 133(6A):A1611, 1964.
  * [35] Duane C. Wallace.  Thermodynamics of crystals.  Wiley New York, 1972.
  * [36] John R Ray and Aneesur Rahman.  Statistical ensembles and molecular dynamics studies of anisotropic solids.  The Journal of chemical physics, 80(9):4423–4428, 1984.
  * [37] KW Wojciechowski.  Constant thermodynamic tension monte carlo studies of elastic properties of a two-dimensional system of hard cyclic hexamers.  Molecular Physics, 61(5):1247–1258, 1987.
  * [38] John R Ray.  Elastic constants and statistical ensembles in molecular dynamics.  Computer physics reports, 8(3):109–151, 1988.
  * [39] Kevin Van Workum, Kenji Yoshimoto, Juan J de Pablo, and Jack F Douglas.  Isothermal stress and elasticity tensors for ions and point dipoles using ewald summations.  Physical Review E, 71(6):061102, 2005.
  * [40] J-L Barrat.  Microscopic elasticity of complex systems.  In Computer Simulations in Condensed Matter Systems: From Materials to Chemical Biology Volume 2, pages 287–307. Springer, 2006.
  * [41] Max Born and Kun Huang.  Dynamical theory of crystal lattices.  Clarendon press, 1954.
  * [42] WC Overton Jr and John Gaffney.  Temperature variation of the elastic constants of cubic elements. i. copper.  Physical Review, 98(4):969, 1955.
  * [43] YA Chang and L Himmel.  Temperature dependence of the elastic constants of cu, ag, and au above room temperature.  Journal of Applied Physics, 37(9):3567–3572, 1966.
  * [44] Haruo Yoshida.  Construction of higher order symplectic integrators.  Physics letters A, 150(5-7):262–268, 1990.
  * [45] Masuo Suzuki.  General theory of fractal path integrals with applications to many-body theories and statistical physics.  Journal of Mathematical Physics, 32(2):400–407, 1991.
  * [46] Y Wang, JJ Wang, H Zhang, VR Manga, SL Shang, LQ Chen, and ZK Liu.  A first-principles approach to finite temperature elastic constants.  Journal of Physics: Condensed Matter, 22(22):225404, 2010.
  * [47] O Gülseren and RE Cohen.  High-pressure thermoelasticity of body-centered-cubic tantalum.  Physical Review B, 65(6):064103, 2002.
  * [48] Hieu H Pham, Michael E Williams, Patrick Mahaffey, Miladin Radovic, Raymundo Arroyave, and Tahir Cagin.  Finite-temperature elasticity of fcc al: Atomistic simulations and ultrasonic measurements.  Physical Review B, 84(6):064101, 2011.
  * [49] Steve Plimpton.  Fast parallel algorithms for short-range molecular dynamics.  Journal of computational physics, 117(1):1–19, 1995.

##  Appendix A Elastic constants in the canonical ensemble

In the canonical ensemble, the free energy is given by:

$\mathcal{F}=-k_{B}T\ln Z,$ (29)

where the canonical partition function is given by:

$\displaystyle Z$ $\displaystyle=\frac{1}{h^{3N}}\int
d^{3}\boldsymbol{p}^{N}\int_{V}d^{3}\boldsymbol{r}^{N}e^{-\beta\mathcal{H}\left(\boldsymbol{r}^{N},\boldsymbol{p}^{N}\right)}$
$\displaystyle=\frac{1}{\prod_{i}^{N}\Lambda_{i}^{3}}\int_{V}d^{3}\boldsymbol{r}^{N}e^{-\beta\mathcal{V}\left(\boldsymbol{r}^{N}\right)},$
(30)

where $\beta=1/k_{B}T$ and
$\Lambda_{i}=\left(\frac{2\pi\hbar^{2}\beta}{m_{i}}\right)^{\frac{1}{2}}$ is
the thermal de-Broglie wavelength. We note that in order calculate strain
derivatives (i.e. as in equations 12-13), the
$\frac{1}{\prod_{i}^{N}\Lambda_{i}^{3}}$ factor can be dropped, so that for
the derivations presented in this appendix, we re-denote:

$Z=\int_{V}d^{3}\boldsymbol{r}^{N}e^{-\beta\mathcal{V}\left(\boldsymbol{r}^{N}\right)}.$
(31)

Consider the canonical average of a configurational operator
$A=A\left(\boldsymbol{r}^{N}\right)$:

$\left\langle
A\right\rangle=\frac{1}{Z}\int_{V}d^{3}\boldsymbol{r}^{N}A\left(\boldsymbol{r}^{N}\right)e^{-\beta\mathcal{V}\left(\boldsymbol{r}^{N}\right)}.$
(32)

In order to evaluate the strain derivative $\frac{\partial\left\langle
A\right\rangle}{\partial\eta_{\alpha\beta}}$, we write the integral in terms
of the reference configuration coordinates $\boldsymbol{R}$, since the strain
derivative cannot be brought inside the integral, because the deformed volume
$V$ depends on $\eta$. Therefore, we write:

$\int_{V}d^{3}\boldsymbol{r}^{N}A\left(\boldsymbol{r}^{N}\right)e^{-\beta\mathcal{V}\left(\boldsymbol{r}^{N}\right)}=\int_{V_{0}}d^{3}\boldsymbol{R}^{N}\left(\det\boldsymbol{J}\right)^{N}Ae^{-\beta\mathcal{V}},$
(33)

where the integrand in the RHS is understood to be evaluated at
$\boldsymbol{r}\left(\boldsymbol{R}\right)$. We can now write the strain
derivative:

$\displaystyle\frac{\partial}{\partial\eta_{\alpha\beta}}\int_{V}d^{3}\boldsymbol{r}^{N}A\left(\boldsymbol{r}^{N}\right)e^{-\beta\mathcal{V}\left(\boldsymbol{r}^{N}\right)}=$
$\displaystyle\int_{V_{0}}d^{3}\boldsymbol{R}^{N}\frac{\partial}{\partial\eta_{\alpha\beta}}\left[\left(\det\boldsymbol{J}\right)^{N}Ae^{-\beta\mathcal{V}}\right],$
(34)

the integrand reads:

$\displaystyle\frac{\partial}{\partial\eta_{\alpha\beta}}\left[\left(\det\boldsymbol{J}\right)^{N}Ae^{-\beta\mathcal{V}}\right]=$
$\displaystyle
e^{-\beta\mathcal{V}}\left(\det\boldsymbol{J}\right)^{N}\left[A\frac{N}{\det\left(\boldsymbol{J}\right)}\frac{\partial\det\left(\boldsymbol{J}\right)}{\partial\eta_{\alpha\beta}}+\frac{\partial
A}{\partial\eta_{\alpha\beta}}-\beta
A\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right].$

Using the well known rule for the derivative of a determinant:

$\frac{\partial\det\left(\boldsymbol{M}\right)}{\partial
M_{\alpha\beta}}=\det\left(\boldsymbol{M}\right)M_{\alpha\beta}^{-T},$ (35)

and eq. 3, it is readily shown that:

$\frac{1}{\det\left(\boldsymbol{J}\right)}\frac{\partial\det\left(\boldsymbol{J}\right)}{\partial\eta_{\alpha\beta}}=\frac{1}{\det\left(\boldsymbol{J}\right)}D_{\alpha\beta}\det\left(\boldsymbol{J}\right)=\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\alpha\beta}^{-T}.$
(36)

As a result, eq. 34 reads:

$\displaystyle\frac{\partial}{\partial\eta_{\alpha\beta}}\left(\int_{V}d^{3}\boldsymbol{r}^{N}A\left(\boldsymbol{r}^{N}\right)e^{-\beta\mathcal{V}\left(\boldsymbol{r}^{N}\right)}\right)=$
$\displaystyle\int_{V_{0}}d^{3}\boldsymbol{R}^{N}\left(\det\boldsymbol{J}\right)^{N}\left[A\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\alpha\beta}^{-T}+\frac{\partial
A}{\partial\eta_{\alpha\beta}}-\beta
A\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right]e^{-\beta\mathcal{V}}$
$\displaystyle=\int_{V}d^{3}\boldsymbol{r}^{N}\left[A\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\alpha\beta}^{-T}+\frac{\partial
A}{\partial\eta_{\alpha\beta}}-\beta
A\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right]e^{-\beta\mathcal{V}}.$
(37)

Using this for the particular case $A\equiv 1$, one finds:

$\displaystyle\frac{1}{Z}\frac{\partial Z}{\partial\eta_{\alpha\beta}}$
$\displaystyle=\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\alpha\beta}^{-T}-\beta\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right\rangle.$
(38)

Differentiation of eq. 32 using the derivative product rule and equations
37-38, results in the general strain derivative rule:

$\frac{\partial\left\langle
A\right\rangle}{\partial\eta_{\alpha\beta}}=\left\langle\frac{\partial
A}{\partial\eta_{\alpha\beta}}\right\rangle-\beta\left[\left\langle
A\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right\rangle-\left\langle
A\right\rangle\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right\rangle\right].$
(39)

For an EAM potential of the form 15, we can write:

$\displaystyle\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}$
$\displaystyle=\sum_{i}\frac{\partial\mathcal{V}}{\partial\boldsymbol{r}_{i}}\cdot\frac{\partial\boldsymbol{r}_{i}}{\partial\eta_{\alpha\beta}}=-\sum_{i}F_{i,\gamma}\frac{\partial
r_{i,\gamma}}{\partial\eta_{\alpha\beta}}$ $\displaystyle=-\sum_{i}\sum_{j\neq
i}F_{ij}\frac{r_{ij,\gamma}}{r_{ij}}\frac{\partial
r_{i,\gamma}}{\partial\eta_{\alpha\beta}}$
$\displaystyle=-\frac{1}{2}\sum_{i}\sum_{j\neq
i}\left(F_{ij}\frac{r_{ij,\gamma}}{r_{ij}}\frac{\partial
r_{i,\gamma}}{\partial\eta_{\alpha\beta}}+F_{ji}\frac{r_{ji,\gamma}}{r_{ji}}\frac{\partial
r_{j,\gamma}}{\partial\eta_{\alpha\beta}}\right)$
$\displaystyle=-\frac{1}{2}\sum_{i}\sum_{j\neq
i}F_{ij}\frac{r_{ij,\gamma}}{r_{ij}}\frac{\partial
r_{ij,\gamma}}{\partial\eta_{\alpha\beta}}=-\frac{1}{2}\sum_{i<j}\frac{F_{ij}}{r_{ij}}\frac{\partial
r_{ij}^{2}}{\partial\eta_{\alpha\beta}},$

where in the last term, the summation is over ordered pairs. From eq. 4 it is
evident that:

$\frac{\partial r^{2}}{\partial\eta_{\alpha\beta}}=2R_{\alpha}R_{\beta},$ (40)

and, more generally, for any radial function $f=f\left(r\right)$, we have:

$D_{\alpha\beta}f=\frac{\partial f}{\partial\eta_{\alpha\beta}}=\frac{\partial
f}{\partial\eta_{\beta\alpha}}=\frac{df}{dr}\frac{R_{\alpha}R_{\beta}}{r}.$
(41)

Hence, we can finally write:

$\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}=\frac{\partial\mathcal{V}}{\partial\eta_{\beta\alpha}}=D_{\alpha\beta}\mathcal{V}=-\sum_{i<j}F_{ij}\frac{R_{ij,\alpha}R_{ij,\beta}}{r_{ij}}.$
(42)

From eq. 38 it follows that:

$\displaystyle\frac{\partial\mathcal{F}}{\partial\eta_{\alpha\beta}}=$
$\displaystyle\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right\rangle-
k_{B}TN\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\alpha\beta}^{-T},$
(43)

which is a symmetric tensor. Hence, eq. 43, when evaluated at zero strain,
proves equations 23 and 25.

Next, in order to evaluate the isothermal elastic constant defined by eq. 13,
consider the second derivative tensor:

$D_{\alpha\beta}D_{\gamma\delta}\mathcal{F}=D_{\alpha\beta}\left[\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\gamma\delta}}\right\rangle-
k_{B}TN\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\gamma\delta}^{-T}\right].$
(44)

In order to evaluate the second term, we use the following identity for the
derivative of a matrix inverse:

$\frac{\partial M_{\delta\gamma}^{-1}}{\partial
M_{\alpha\beta}}=-M_{\delta\alpha}^{-1}M_{\beta\gamma}^{-1},$ (45)

which results in the relation:

$\displaystyle
D_{\alpha\beta}\left[\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\gamma\delta}^{-T}\right]$
$\displaystyle=-\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\delta\alpha}^{-1}\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\beta\gamma}^{-1}$
$\displaystyle-\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\delta\beta}^{-1}\left(2\boldsymbol{\eta}+\boldsymbol{I}\right)_{\alpha\gamma}^{-1}.$
(46)

In order to evaluate the first term in 44, we use the identity 39:

$\displaystyle\frac{\partial}{\partial\eta_{\alpha\beta}}\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\gamma\delta}}\right\rangle=$
$\displaystyle\left\langle\frac{\partial^{2}\mathcal{V}}{\partial\eta_{\alpha\beta}\partial\eta_{\gamma\delta}}\right\rangle-\beta\left[\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\gamma\delta}}\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right\rangle-\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\gamma\delta}}\right\rangle\left\langle\frac{\partial\mathcal{V}}{\partial\eta_{\alpha\beta}}\right\rangle\right]$
$\displaystyle=\left\langle\frac{\partial^{2}\mathcal{V}}{\partial\eta_{\alpha\beta}\partial\eta_{\gamma\delta}}\right\rangle-\frac{V^{2}}{k_{B}T}\left[\left\langle\sigma_{\alpha\beta}^{B}\sigma_{\gamma\delta}^{B}\right\rangle-\left\langle\sigma_{\alpha\beta}^{B}\right\rangle\left\langle\sigma_{\gamma\delta}^{B}\right\rangle\right].$
(47)

We note that the second term is a symmetric tensor with respect to
$\alpha\leftrightarrow\beta$ and $\gamma\leftrightarrow\delta$. In order to
evaluate the first term, we differentiate eq. 42:

$\displaystyle\frac{\partial^{2}\mathcal{V}}{\partial\eta_{\alpha\beta}\partial\eta_{\gamma\delta}}=\frac{\partial}{\partial\eta_{\alpha\beta}}\left(-\sum_{i<j}F_{ij}\frac{R_{ij,\gamma}R_{ij,\delta}}{r_{ij}}\right)$
$\displaystyle=-\sum_{i<j}\left[\frac{R_{ij,\gamma}R_{ij,\delta}}{r_{ij}}\frac{\partial
F_{ij}}{\partial\eta_{\alpha\beta}}+R_{ij,\gamma}R_{ij,\delta}F_{ij}\frac{\partial}{\partial\eta_{\alpha\beta}}\left(\frac{1}{r_{ij}}\right)\right].$
(48)

Using the identity 41, we get:

$\displaystyle\frac{\partial}{\partial\eta_{\alpha\beta}}\left(\frac{1}{r}\right)$
$\displaystyle=-\frac{R_{\alpha}R_{\beta}}{r^{3}},$ (49)

and using eq. 17 and the identity 41 gives:

$\displaystyle\frac{\partial F_{ij}}{\partial\eta_{\alpha\beta}}$
$\displaystyle=\left(v^{\prime\prime}\left(r_{ij}\right)+\left[F^{\prime}\left(\rho_{i}\right)+F^{\prime}\left(\rho_{j}\right)\right]\rho^{\prime\prime}\left(r_{ij}\right)\right)\frac{R_{ij,\alpha}R_{ij,\beta}}{r_{ij}}$
$\displaystyle+\rho^{\prime}\left(r_{ij}\right)\frac{\partial}{\partial\eta_{\alpha\beta}}\left[F^{\prime}\left(\rho_{i}\right)+F^{\prime}\left(\rho_{j}\right)\right].$
(50)

In order to evaluate the last term in eq. 50, we write:

$\displaystyle\frac{\partial}{\partial\eta_{\alpha\beta}}F^{\prime}\left(\rho_{i}\right)$
$\displaystyle=F^{\prime\prime}\left(\rho_{i}\right)\frac{\partial\rho_{i}}{\partial\eta_{\alpha\beta}}=F^{\prime\prime}\left(\rho_{i}\right)\sum_{k\neq
i}\frac{\partial\rho}{\partial\eta_{\alpha\beta}}\bigg{|}_{r_{ik}}$
$\displaystyle=F^{\prime\prime}\left(\rho_{i}\right)\sum_{k\neq
i}\rho^{\prime}\left(r_{ik}\right)\frac{R_{ik,\alpha}R_{ik,\beta}}{r_{ik}^{2}},$
(51)

where in the last step the identity 41 was used again. Plugging equations
49-50 back in eq. 48 gives:

$\displaystyle\frac{\partial^{2}\mathcal{V}}{\partial\eta_{\alpha\beta}\partial\eta_{\gamma\delta}}=D_{\alpha\beta}D_{\gamma\delta}\mathcal{V}=\sum_{i<j}X_{ij}\frac{R_{ij,\alpha}R_{ij,\beta}R_{ij,\gamma}R_{ij,\delta}}{r_{ij}^{2}}$
$\displaystyle+\sum_{i}\left(\sum_{k\neq
i}\frac{R_{ik,\alpha}R_{ik,\beta}}{r_{ik}^{2}}\rho^{\prime}\left(r_{ik}\right)\right)\left(\sum_{k\neq
i}\frac{R_{ik,\alpha}R_{ik,\beta}}{r_{ik}^{2}}\rho^{\prime}\left(r_{ik}\right)\right),$
(52)

where $X_{ij}$ is defined by eq. 27, and we have used the fact that the
resulting expression defines a symmetric tensor with respect to
$\alpha\leftrightarrow\beta$ and $\gamma\leftrightarrow\delta$.

Finally, when evaluated at zero strain ($\boldsymbol{\eta}=\boldsymbol{0}$,
$\boldsymbol{R}=\boldsymbol{r}$), the combination of equations 44, 46, 47 and
52 proves the relations 21 and 26-28.

##  Appendix B Equations of motion in the NVT ensemble

In this appendix we list the equations of motion of the widely used Nose-
Hoover-Chains method [28, 26, 27], which was used in the calculations
performed in this work. We use a chain of $M$ thermostats with coordinates
$\eta_{j}$ and momenta $p_{\eta_{j}}$. The equations of motions for the
$i=1...N$ particles are:

$\frac{d\boldsymbol{r}_{i}}{dt}=\frac{\boldsymbol{p}_{i}}{m_{i}},$ (53)
$\frac{d\boldsymbol{p}_{i}}{dt}=\boldsymbol{F}_{i}-\left(\frac{p_{\eta_{1}}}{Q_{1}}\right)\boldsymbol{p}_{i}.$
(54)

The equations of motion for the $j=1...M$ thermostat variables are given by:

$\frac{d\eta_{j}}{dt}=\frac{p_{\eta_{j}}}{Q_{j}},$ (55)

$\frac{dp_{\eta_{j}}}{dt}=\begin{cases}G_{j}-\left(\frac{p_{\eta_{j+1}}}{Q_{j+1}}\right)p_{\eta_{j}},&j=1...M-1\\\
G_{M},&j=M\end{cases}$ (56)

where the forces are:

$G_{j}=\begin{cases}\sum_{i}\frac{p_{i}^{2}}{m_{i}}-gk_{B}T,&j=1\\\
\frac{p_{\eta_{j-1}}^{2}}{Q_{j-1}}-k_{B}T,&j=2...M\end{cases}$ (57)

with the number of degrees of freedom:

$g=3\left(N-1\right).$ (58)

The thermostat “masses” $Q_{j}$ can be written in terms of a thermostat
relaxation timescale $\tau_{T}$ [28]:

$Q_{j}=\begin{cases}gk_{B}T\tau_{T}^{2},&j=1\\\
k_{B}T\tau_{T}^{2},&j=2...M\end{cases}$ (59)

The non-Hamiltonian system 53-56 has the conserved quantity [28, 26]:

$E_{\text{cons}}=\mathcal{H}+\sum_{j=1}^{M}\frac{p_{\eta_{j}}^{2}}{2Q_{j}}+gk_{B}T\eta_{1}+\sum_{j=2}^{M}k_{B}T\eta_{j},$
(60)

where $\mathcal{H}$ is the true Hamiltonian given by eq. 14.

##  Appendix C Equations of motion in the NPT ensemble

In this appendix we list the equations of motion in the widely used Nose-
Hoover-Chains method in the isothermal-isobaric (NPT) ensemble [22, 25, 26,
27], which was used in the calculations performed in this work. As in appendix
B, we use a chain of $M$ thermostats whose coordinates and momenta are
$\eta_{j}$ and $p_{\eta_{j}}$, which are coupled the the equations of motion
of the $N$ particles. In the NPT ensemble, the volume $V=V\left(t\right)$ is
treated as a dynamical variable via the generalized coordinate:

$\epsilon\left(t\right)=\frac{1}{3}\ln\left(\frac{V\left(t\right)}{V\left(0\right)}\right),$
(61)

whose momentum $p_{\epsilon}$ is coupled to an additional chain of
$M^{\prime}$ thermostats with coordinates $\eta_{j}^{\prime}$ and momenta
$p_{\eta_{j}}^{\prime}$. The equations of motions for the $i=1...N$ particles
are:

$\frac{d\boldsymbol{r}_{i}}{dt}=\frac{\boldsymbol{p}_{i}}{m_{i}}+\frac{p_{\epsilon}}{W}\boldsymbol{r}_{i},$
(62)
$\frac{d\boldsymbol{p}_{i}}{dt}=\boldsymbol{F}_{i}-\left(1+\frac{3}{g}\right)\left(\frac{p_{\epsilon}}{W}\right)\boldsymbol{p}_{i}-\left(\frac{p_{\eta_{1}}}{Q_{1}}\right)\boldsymbol{p}_{i},$
(63)

The equations of motion for the barostat variables are:

$\frac{d\epsilon}{dt}=\frac{p_{\epsilon}}{W},$ (64)

$\frac{dp_{\epsilon}}{dt}=3\left(P_{inst}V-PV\right)+\frac{3}{g}\sum_{i=1}^{N}\frac{p_{i}^{2}}{m_{i}}-\left(\frac{p_{\eta_{1}^{\prime}}}{Q_{1}^{\prime}}\right)p_{\epsilon},$
(65)

where $P$ is the applied external pressure and the instantaneous pressure is
given by the virial theorem:

$P_{inst}V=\frac{1}{3}\sum_{i}\frac{p_{i}^{2}}{m_{i}}+\frac{1}{3}\sum_{i<j}F_{ij}r_{ij}.$
(66)

The equations of motion for the $j=1...M$ particle-coupled thermostat
variables are given by:

$\frac{d\eta_{j}}{dt}=\frac{p_{\eta_{j}}}{Q_{j}},$ (67)

$\frac{dp_{\eta_{j}}}{dt}=\begin{cases}G_{j}-\left(\frac{p_{\eta_{j+1}}}{Q_{j+1}}\right)p_{\eta_{j}},&j=1...M-1\\\
G_{M},&j=M\end{cases}$ (68)

where the forces are:

$G_{j}=\begin{cases}\sum_{i}\frac{p_{i}^{2}}{m_{i}}-gk_{B}T,&j=1\\\
\frac{p_{\eta_{j-1}}^{2}}{Q_{j-1}}-k_{B}T,&j=2...M\end{cases}$ (69)

The equations of motion for the $j=1...M^{\prime}$ barostat-coupled thermostat
variables are given by:

$\frac{d\eta_{j}^{\prime}}{dt}=\frac{p_{\eta_{j}^{\prime}}}{Q_{j}^{\prime}},$
(70)

$\frac{dp_{\eta_{j}^{\prime}}}{dt}=\begin{cases}G_{j}^{\prime}-\left(\frac{p_{\eta_{j+1}^{\prime}}}{Q_{j+1}^{\prime}}\right)p_{\eta_{j}^{\prime}},&j=1...M^{\prime}-1\\\
G_{M^{\prime}}^{\prime},&j=M^{\prime}\end{cases}$ (71)

where the forces are:

$G_{j}^{\prime}=\begin{cases}\frac{p_{\epsilon}^{2}}{W}-k_{B}T,&j=1\\\
\frac{p_{\eta_{j-1}^{\prime}}^{2}}{Q_{j-1}^{\prime}}-k_{B}T,&j=2...M^{\prime}\end{cases}$
(72)

The particle-coupled thermostat masses $Q_{j}$ are given in terms of a
thermostat relaxation timescale $\tau_{T}$ as in 59, while the barostat mass
$W$ and the barostat-coupled thermostat masses $Q_{j}^{\prime}$ are written in
terms of a barostat relaxation timescale $\tau_{p}$ as [28, 23]:

$W=gk_{B}T\tau_{p}^{2},$ (73) $Q_{j}^{\prime}=k_{B}T\tau_{p}^{2}.$ (74)

The non-Hamiltonian system 62-71 has the conserved quantity [22, 25, 26]:

$\displaystyle E_{\text{cons}}$
$\displaystyle=\mathcal{H}+PV+\frac{p_{\epsilon}^{2}}{2W}+\sum_{j=1}^{M}\frac{p_{\eta_{j}}^{2}}{2Q_{j}}+\sum_{j=1}^{M^{\prime}}\frac{p_{\eta_{j}^{\prime}}^{2}}{2Q_{j}^{\prime}}$
$\displaystyle+gk_{B}T\eta_{1}+\sum_{j=2}^{M}k_{B}T\eta_{j}+\sum_{j=1}^{M^{\prime}}k_{B}T\eta_{j}^{\prime}.$
(75)