Elastic constants

We compute clamped-ion elastic constants of a crystal using the algorithmic differentiation density-functional perturbation theory (AD-DFPT) approach as introduced in ^[SPH25].

We consider a crystal in its equilibrium configuration, where all atomic forces and stresses vanish. Homogeneous strains $η$ are then applied relative to this relaxed structure. The elastic constants are derived from the stress-strain relationship. In Voigt notation, the stress $\sigma$ and strain $\eta$ tensors are represented as 6-component vectors. The elastic constants $C$ are then given by the Jacobian of the stress with respect to strain, forming a $6 \times 6$ matrix

\[ C = \frac{\partial \sigma}{\partial \eta}.\]

This example computes the clamped-ion elastic tensor, keeping internal atomic positions fixed under strain. The relaxed-ion tensor includes additional corrections from internal relaxations, which can be obtained from first-order atomic displacements in DFPT (see ^[Wu2005]).

using DFTK
using PseudoPotentialData
using LinearAlgebra
using ForwardDiff
using DifferentiationInterface
using AtomsBuilder
using Unitful
using UnitfulAtomic

Computing PBE elastic constants

We start with the PBE ^[PBE1996] functional.

pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf")
a0_pbe = 10.33u"bohr"  # Equilibrium lattice constant of silicon with PBE
model0 = model_DFT(bulk(:Si; a=a0_pbe); pseudopotentials, functionals=PBE())

Ecut = recommended_cutoff(model0).Ecut
kgrid = [4, 4, 4]
tol = 1e-6

1.0e-6

The elastic_tensor postprocessing function automatically detects crystal symmetry and chooses the appropriate strain patterns to extract all independent elastic constants:

basis = PlaneWaveBasis(model0; Ecut, kgrid)
scfres = self_consistent_field(basis; tol)
(; C) = elastic_tensor(scfres)

println("C11: ", uconvert(u"GPa", C[1, 1] * u"hartree" / u"bohr"^3))
println("C12: ", uconvert(u"GPa", C[1, 2] * u"hartree" / u"bohr"^3))
println("C44: ", uconvert(u"GPa", C[4, 4] * u"hartree" / u"bohr"^3))

n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -8.453366585654                   -0.94    5.6    208ms
  2   -8.455587375260       -2.65       -1.78    1.0   99.4ms
  3   -8.455777302216       -3.72       -2.84    2.0    125ms
  4   -8.455790104820       -4.89       -3.48    3.1    156ms
  5   -8.455790183865       -7.10       -4.32    1.8    121ms
  6   -8.455790187641       -8.42       -5.22    2.4    135ms
  7   -8.455790187712      -10.15       -6.02    2.8    148ms
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -8.455790187713                   -6.70   19.2    1.16s
Solving response problem
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      17.9   7.49s  Non-interacting
   1        0       1        0.11     13.8   16.6s  
   2        0       2       -0.70     12.5   827ms  
   3        0       3       -1.71     10.9   480ms  
   4        0       4       -2.78      9.1   406ms  
   5        0       5       -4.04      6.7   330ms  
   6        0       6       -5.56      2.7   205ms  
   7        0       7       -7.73      1.0   154ms  
   8        1       1       -6.72     17.8   853ms  Restart
                                      16.8   1.83s  Final orbitals
C11: 156.51815351970552 GPa
C12: 59.57460325751817 GPa
C44: 98.61658316643093 GPa

Here are AD-DFPT results from increasing discretization parameters:

Ecut	kgrid	c11	c12	c44
18	[4, 4, 4]	156.51	59.57	98.61
18	[8, 8, 8]	153.53	56.90	100.07
24	[8, 8, 8]	153.26	56.82	99.97
24	[14, 14, 14]	153.03	56.71	100.09

For comparison, Materials Project for PBE relaxed-ion elastic constants of silicon mp-149: c11 = 153 GPa, c12 = 57 GPa, c44 = 74 GPa. Note the discrepancy in c44, which is due to us not yet including ionic relaxation in this example.

Moving to meta-GGA functionals

To make the problem a little more interesting we will now compute the elastic constants using the Laplacian-dependent r2SCAN-L functional (a deorbitalized version of r2SCAN, see ^[MT2020]), again using AD-DFPT.

model_r2scanl = model_DFT(bulk(:Si; a=a0_pbe);
                       pseudopotentials, functionals=[:mgga_x_r2scanl, :mgga_c_r2scanl])
basis_r2scanl = PlaneWaveBasis(model_r2scanl; Ecut, kgrid)
scfres_r2scanl = self_consistent_field(basis_r2scanl; tol)
(; C) = elastic_tensor(scfres_r2scanl)

(voigt_stress = [-1.9912982462690275e-5, -1.991297726981849e-5, -1.9912977269818477e-5, -4.810649983824134e-12, 0.0, 0.0], C = [0.005299351720418529 0.0020171883054278003 0.0020171883054278003 0.0 0.0 0.0; 0.0020171883054278003 0.005299351720418529 0.0020171883054278003 0.0 0.0 0.0; 0.0020171883054278003 0.0020171883054278003 0.005299351720418529 0.0 0.0 0.0; 0.0 0.0 0.0 0.0033291650782267474 0.0 0.0; 0.0 0.0 0.0 0.0 0.0033291650782267474 0.0; 0.0 0.0 0.0 0.0 0.0 0.0033291650782267474])

The r2SCAN-L elastic constants, which agree well with PBE:

println("C11: ", uconvert(u"GPa", C[1, 1] * u"hartree" / u"bohr"^3))
println("C12: ", uconvert(u"GPa", C[1, 2] * u"hartree" / u"bohr"^3))
println("C44: ", uconvert(u"GPa", C[4, 4] * u"hartree" / u"bohr"^3))

C11: 155.91231014840344 GPa
C12: 59.347728797052326 GPa
C44: 97.94741802318578 GPa

Manual AD-DFPT derivation

For reference, the following shows how elastic_tensor works under the hood for cubic crystals. The sparsity pattern of the matrix $C$ follows from crystal symmetry and is tabulated in standard references (eg. Table 9 in ^[Nye1985]). This sparsity can be used a priori to reduce the number of strain patterns that need to be probed to extract all independent components of $C$. For example, cubic crystals have only three independent elastic constants $C_{11}$, $C_{12}$ and $C_{44}$, with the pattern

\[C = \begin{pmatrix} C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{44} \\ \end{pmatrix}.\]

Thus we can just choose a suitable strain pattern $\dot{\eta} = (1, 0, 0, 1, 0, 0)^\top$, such that $C\dot{\eta} = (C_{11}, C_{12}, C_{12}, C_{44}, 0, 0)^\top$. That is, for cubic crystals like diamond silicon we obtain all independent elastic constants from a single Jacobian-vector product on the stress-strain function.

function symmetries_from_strain(model0, voigt_strain)
    lattice = DFTK.voigt_strain_to_full(voigt_strain) * model0.lattice
    model = Model(model0; lattice, symmetries=true)
    model.symmetries
end

strain_pattern = [1., 0., 0., 1., 0., 0.];  # recovers [c11, c12, c12, c44, 0, 0]

For elastic constants beyond the bulk modulus, symmetry-breaking strains are required. That is, the symmetry group of the crystal is reduced. Here we simply precompute the relevant subgroup by applying the automatic symmetry detection (spglib) to the finitely perturbed crystal.

symmetries_strain = symmetries_from_strain(model0, 0.01 * strain_pattern)


function stress_from_strain(model0, voigt_strain; symmetries, Ecut, kgrid, tol)
    lattice = DFTK.voigt_strain_to_full(voigt_strain) * model0.lattice
    model = Model(model0; lattice, symmetries)
    basis = PlaneWaveBasis(model; Ecut, kgrid)
    scfres = self_consistent_field(basis; tol)
    DFTK.full_stress_to_voigt(compute_stresses_cart(scfres))
end

stress_fn(voigt_strain) = stress_from_strain(model0, voigt_strain;
                                             symmetries=symmetries_strain,
                                             Ecut, kgrid, tol)
stress, (dstress,) = value_and_pushforward(stress_fn, AutoForwardDiff(),
                                           zeros(6), (strain_pattern,));

n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -8.453428185901                   -0.94    5.3    421ms
  2   -8.455623018462       -2.66       -1.77    1.0    168ms
  3   -8.455776463151       -3.81       -2.90    1.9    182ms
  4   -8.455790069035       -4.87       -3.28    3.0    255ms
  5   -8.455790173614       -6.98       -3.78    1.2    151ms
  6   -8.455790186670       -7.88       -4.69    1.5    172ms
  7   -8.455790187675       -9.00       -5.16    2.2    210ms
  8   -8.455790187711      -10.45       -5.93    1.9    191ms
  9   -8.455790187713      -11.61       -6.51    2.2    219ms
Solving response problem
Iter  Restart  Krydim  log10(res)  avg(CG)  Δtime   Comment
----  -------  ------  ----------  -------  ------  ---------------
                                      17.8   582ms  Non-interacting
   1        0       1        0.11     13.5   581ms  
   2        0       2       -0.70     12.4   497ms  
   3        0       3       -1.71     10.8   447ms  
   4        0       4       -2.78      8.8   389ms  
   5        0       5       -4.04      6.1   310ms  
   6        0       6       -5.56      2.3   202ms  
   7        0       7       -7.80      1.0   168ms  
   8        1       1       -6.63     17.7   787ms  Restart
                                      16.8   630ms  Final orbitals

We can inspect the stress to verify it is small (close to equilibrium):

stress

6-element StaticArraysCore.SVector{6, Float64} with indices SOneTo(6):
 -2.3120771627894893e-5
 -2.3120745836678844e-5
 -2.3120745836678817e-5
 -7.864670032483498e-10
  0.0
  0.0

The response of the stress to strain_pattern contains the elastic constants in atomic units, with the expected pattern $(c11, c12, c12, c44, 0, 0)$:

dstress

6-element StaticArraysCore.SVector{6, Float64} with indices SOneTo(6):
 0.0053199421438012
 0.0020248997477434424
 0.0020248997477434424
 0.0033519095106697974
 0.0
 0.0

Convert to GPa:

println("C11: ", uconvert(u"GPa", dstress[1] * u"hartree" / u"bohr"^3))
println("C12: ", uconvert(u"GPa", dstress[2] * u"hartree" / u"bohr"^3))
println("C44: ", uconvert(u"GPa", dstress[4] * u"hartree" / u"bohr"^3))

C11: 156.51810131794528 GPa
C12: 59.5746072624645 GPa
C44: 98.61658232710371 GPa

These results can be compared directly to finite differences of the stress_fn:

h = 1e-3
dstress_fd = (stress_fn(h * strain_pattern) - stress_fn(-h * strain_pattern)) / 2h
println("C11 (FD): ", uconvert(u"GPa", dstress_fd[1] * u"hartree" / u"bohr"^3))
println("C12 (FD): ", uconvert(u"GPa", dstress_fd[2] * u"hartree" / u"bohr"^3))
println("C44 (FD): ", uconvert(u"GPa", dstress_fd[4] * u"hartree" / u"bohr"^3))

n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -8.453447000852                   -0.94    5.3    368ms
  2   -8.455620546376       -2.66       -1.77    1.0    138ms
  3   -8.455780653223       -3.80       -2.89    1.8    190ms
  4   -8.455795166193       -4.84       -3.24    3.0    255ms
  5   -8.455795296321       -6.89       -3.69    1.2    155ms
  6   -8.455795314074       -7.75       -5.04    1.3    158ms
  7   -8.455795315193       -8.95       -5.52    3.2    276ms
  8   -8.455795315197      -11.38       -5.80    1.3    156ms
  9   -8.455795315198      -12.12       -6.86    1.1    242ms
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -8.453482927176                   -0.94    5.3    348ms
  2   -8.455624792790       -2.67       -1.78    1.0    141ms
  3   -8.455768170689       -3.84       -2.89    1.8    184ms
  4   -8.455782119567       -4.86       -3.24    3.1    254ms
  5   -8.455782235680       -6.94       -3.66    1.1    154ms
  6   -8.455782252547       -7.77       -4.79    1.1    147ms
  7   -8.455782254349       -8.74       -5.18    2.9    253ms
  8   -8.455782254369      -10.70       -5.68    1.3    163ms
  9   -8.455782254371      -11.68       -6.53    1.4    162ms
C11 (FD): 156.27960064686212 GPa
C12 (FD): 59.50474077717529 GPa
C44 (FD): 98.5442082989543 GPa

SPH25Schmitz, N. F., Ploumhans, B., & Herbst, M. F. (2026) Algorithmic differentiation for plane-wave DFT: materials design, error control and learning model parameters. npj Computational Materials 12, 6 (2026). (Supplementary material and computational scripts).
Wu2005Wu, X., Vanderbilt, D., & Hamann, D. R. (2005). Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory. Physical Review B, 72(3), 035105.
PBE1996Perdew, J. P., Burke, K., & Ernzerhof, M. (1996). Generalized Gradient Approximation Made Simple. Physical Review Letters, 77(18), 3865-3868.
MT2020Mejía-Rodríguez, D., & Trickey, S. B. (2020). Meta-GGA performance in solids at almost GGA cost. Physical Review B, 102(12), 121109(R).
Nye1985Nye, J. F. (1985). Physical Properties of Crystals. Oxford University Press. Comment: Since the elastic tensor transforms equivariantly under rotations, its numerical components depend on the chosen Cartesian coordinate frame. These tabulated patterns assume a standardized orientation of the structure with respect to conventional crystallographic axes.