First-Principles Calculation of Third-Order Elastic Constants via Numerical Differentiation of the Second Piola-Kirchhoff Stress Tensor

A general method is presented to calculate from first principles the full set of third-order elastic constants of a material of arbitrary symmetry. The method here illustrated relies on a plane-wave density functional theory scheme to calculate the Cauchy stress and the numerical differentiation of the second Piola-Kirchhoff stress tensor to evaluate the elastic constants. It is shown that finite difference formulas lead to a cancellation of the finite basis set errors, whereas simple solutions are proposed to eliminate numerical errors arising from the use of Fourier interpolation techniques. Applications to diamond, silicon, aluminum, magnesium, graphene, and a graphane conformer give results in excellent agreement with both experiments and previous calculations based on fitting energy density curves, demonstrating both the accuracy and generality of our new methodology to investigate nonlinear elastic behaviors of materials.

Figure 1

Energy density $U$ [Eq. (3)] of diamond versus normal uniaxial (Lagrangian) strain ${ϵ}_1$. Inset: Component $P_1$ of the 2nd-PK stress tensor versus ${ϵ}_1$. Colored disks show the results obtained from DFT calculations [19] using an energy cutoff of 300 Ry. Solid lines are guides to the eye. To calculate $C_{11}^{(2)}$ and $C_{111}^{(3)}$, we use the value of $P_1$ at the unstressed state (large black disk) and those obtained at either ${ϵ}_1=\pm 0.005$ or ${ϵ}_1=\pm 0.010$ (circles).

Figure 2

$C_{144}^{(3)}$ of diamond obtained from Eq. (7) by using values of $P_1$ calculated by DFT at increasing energy cutoffs. $C_{144}^{(3)}$ is calculated by using the value of $P_1$ at an unstressed state and values of $P_1$ of crystals accommodating a shear strain of ${ϵ}_4=\pm 0.005$. The blue (black) solid line shows the results obtained by (not) including in the Cauchy stress the Pulay corrective terms of Eq. (9). The red dashed line shows the results obtained by replacing $P_1$ at the unstressed state in Eq. (7) with the value obtained by extrapolation from the values of $P_1$ at ${ϵ}_4=\pm 0.005$ and $\pm 0.010$.

Figure 3

Component $P_1$ of the 2nd-PK stress tensor of a diamond crystal accommodating a shear strain ${ϵ}_4$. Disks show results obtained from DFT calculations carried out using energy cutoffs of 100 (blue and black) and 200 Ry (red), with atoms in the unstressed primitive unit cell having fractional coordinates (0,0,0) and (0.25,0.25,0.25) (blue) and $(x,y,z)$ and $(0.25+x,0.25+y,0.25+z)$ (black and red), where $x$, $y$, $z$ are random numbers in the interval (0,1). Insets: Schematic representations of unstressed (blue and black) and shear strained (red) unit cells, with lattice coordinates that are (blue) or are not (black and red) aligned with the real-space grid of points used to represent wave functions in plane-wave-based DFT calculations.