Analytic elastic coefficients in molecular calculations: Finite strain, nonaffine displacements, and many-body interatomic potentials

Elastic moduli are among the most fundamental and important properties of solid materials, which is why they are routinely characterized in both experiments and simulations. While conceptually simple, the treatment of elastic coefficients is complicated by two factors not yet having been concurrently discussed: finite-strain and nonaffine, internal displacements. Here, we revisit the theory behind zero-temperature, finite-strain elastic moduli and extend it to explicitly consider nonaffine displacements. We further present analytical expressions for second-order derivatives of the potential energy for two-body and generic many-body interatomic potentials, such as cluster and empirical bond-order potentials. Specifically, we revisit the elastic constants of silicon, silicon carbide, and silicon dioxide under hydrostatic compression and dilatation. Based on existing and recent results, we outline the effect of multiaxial stress states as opposed to volumetric deformation on the limits of stability of their crystalline lattices.

Figure 1

Thermodynamics of elastic deformation. The external stress $\underline{\mathrm{\Pi }}$ is a property of a hypothetical ideal embedding medium. Equilibrium is the case where embedding medium and system pressure are equal, $\underline{\mathrm{\Pi }}={\underline{\sigma }}^{\eta }$. (a) Our reference configuration is not necessarily at equilibrium. (b) Equilibrium is obtained at some finite strain ${\eta }^{\text{eq}}$ where bath and system pressure are equal. (c) We make small deformations around this equilibrium condition to probe the elastic properties of the system at constant stress $\underline{\mathrm{\Pi }}$.

Figure 2

Mechanical work required to deform an RVE along a simple, exemplary closed deformation path while maintaining either the Cauchy stress $\underline{\sigma }$ or the Piola–Kirchhoff stress $\underline{\mathrm{\Pi }}$ constant. At constant $\underline{\sigma }$ (blue), the work around this cycle is not zero while it is zero at constant $\underline{\mathrm{\Pi }}$ (orange). The reason for this is that when extending the area $\mathring{A}$ to $A_1$, the effective force on this area increases at constant $\underline{\sigma }$ but is constant at constant $\underline{\mathrm{\Pi }}$. This path dependence implies that the state of the bath has changed after this cycle at constant $\underline{\sigma }$.

Figure 3

Numerical and analytical calculations of the elastic coefficients of silicon carbide as a function of pressure. Results are shown for the TIII interatomic potential with a fixed bonding topology.

Figure 4

Pressure-dependent elastic coefficients at zero temperature for silicon (a)–(d), silicon carbide (e)–(h), and silicon dioxide (i)–(l). DFT results for silicon are taken from Ref. [87], for silicon carbide from Ref. [88], and for silicon dioxide from Ref. [89]. $C_{ij}$: Elastic coefficients with nonaffine contributions. $c_{ij}$: Elastic coefficients without nonaffine contributions.

Figure 5

Effect of a fixed bonding topology versus a variable topology on the pressure dependence of $C_{44}$ for silicon (Kum) (a) and silicon carbide (EA) (b). The results deviate in regions where neighbors sit within the cutoff distance.

Figure 6

Eigenvalues of the elastic tensor $\underline{\underline{C}}$, denoted by $M_1, M_2, M_3$ for cubic solids and $N_1, N_2, N_3, N_4$ for silicon dioxide, and smallest eigenvalue ${\lambda }_{\text{min}}$ of the Hessian $\mathbb{H}$ as a function of pressure $P$. Shown are results for silicon (a)–(d), silicon carbide (e)–(h), and silicon dioxide (i)–(l). The eigenvalues are normalized to their value at zero hydrostatic pressure. The systems lose mechanical stability when an eigenvalue becomes negative.

Figure 7

Smallest eigenvalue ${\gamma }^{\text{min}}$ of the elastic tensor ${\underline{\underline{C}}}^{\text{sym}}$ and smallest eigenvalue ${\lambda }^{\text{min}}$ of the Hessian $\mathbb{H}$ of diamond cubic silicon. The multiaxial stress state is triaxial compression, ${\sigma }_{xx}={\sigma }_{yy}\ne {\sigma }_{zz}$, and ${\sigma }_{xy}={\sigma }_{xz}={\sigma }_{yz}=0$.