Calculation of elastic constants of embedded-atom-model potentials in the $NVT$ ensemble

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 because 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 nonhomogeneous crystalline materials and (ii) 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–1000 K, 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.

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–1000 K. Time-dependent measures of the simulation at $T=300$ K are presented in detail in Figs. 2 and 3.

Figure 2

Analysis of a molecular-dynamics simulation of Copper in the $NPT$ ensemble at zero pressure and $T=300$ 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 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\AA$). 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

Components of the instantaneous pressure tensor [Eq. (18)] for the $NPT$ simulation described in Fig. 2.

Figure 4

Components of the instantaneous Born stress tensor [upper figure, Eq. (25)] for a $NVT$ simulation of Copper at $T=300\mathrm{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

Upper figures show the Born elastic constants [Eq. (26)]: $C_{11}^B$ (upper-left figure), $C_{12}^B$ (middle-upper figure), and $C_{44}^B$ (upper-right figure) for a $NVT$ simulation of Copper at $T=300\mathrm{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 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. (B8)] and $NVT$ [upper figure, Eq. (C15)] molecular-dynamics simulations of Copper at $T=300\mathrm{K}$ and zero pressure (the simulations described in Figs. 2, 3 and 4, 5, respectively).

Figure 7

The different terms contributing to the total isothermal elastic constants $C_{11}$ (left panel), $C_{12}$ (middle panel), and $C_{44}$ (right panel) of Copper as a function of temperature (solid blue line): the Born term [first term in Eq. (21), dashed orange line], the 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

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–300 K given in Ref. [42] and squares represent experimental values in the range 300–800 K given in Ref. [43] for $C_{44}$. The blue $\times$ data point represents a calculated $C_{44}$ result given in Ref. [10].