|
| 1 | +# # Elastic constants |
| 2 | + |
| 3 | +# We compute *clamped-ion* elastic constants of a crystal using |
| 4 | +# the algorithmic differentiation density-functional perturbation theory (AD-DFPT) approach |
| 5 | +# as introduced in [^SPH25]. |
| 6 | +# |
| 7 | +# [^SPH25]: |
| 8 | +# Schmitz, N. F., Ploumhans, B., & Herbst, M. F. (2025) |
| 9 | +# *Algorithmic differentiation for plane-wave DFT: materials design, error control and learning model parameters.* |
| 10 | +# [arXiv:2509.07785](https://arxiv.org/abs/2509.07785) |
| 11 | +# |
| 12 | +# We consider a crystal in its equilibrium configuration, where all atomic forces |
| 13 | +# and stresses vanish. Homogeneous strains $η$ are then applied |
| 14 | +# relative to this relaxed structure. |
| 15 | +# The elastic constants are derived from the stress-strain relationship. |
| 16 | +# In [Voigt notation](https://en.wikipedia.org/wiki/Voigt_notation), |
| 17 | +# the stress $\sigma$ and strain $\eta$ tensors are represented as 6-component vectors. |
| 18 | +# The elastic constants $C$ are then given by |
| 19 | +# the Jacobian of the stress with respect to strain, forming a $6 \times 6$ matrix |
| 20 | +# ```math |
| 21 | +# C = \frac{\partial \sigma}{\partial \eta}. |
| 22 | +# ``` |
| 23 | +# |
| 24 | +# The sparsity pattern of the matrix $C$ follows from crystal symmetry |
| 25 | +# and is tabulated in standard references (eg. Table 9 in [^Nye1985]). |
| 26 | +# This sparsity can be used a priori to reduce the number of strain patterns |
| 27 | +# that need to be probed to extract all independent components of $C$. |
| 28 | +# For example, cubic crystals have only three independent elastic constants |
| 29 | +# $C_{11}$, $C_{12}$ and $C_{44}$, with the pattern |
| 30 | +# ```math |
| 31 | +# C = \begin{pmatrix} |
| 32 | +# C_{11} & C_{12} & C_{12} & 0 & 0 & 0 \\ |
| 33 | +# C_{12} & C_{11} & C_{12} & 0 & 0 & 0 \\ |
| 34 | +# C_{12} & C_{12} & C_{11} & 0 & 0 & 0 \\ |
| 35 | +# 0 & 0 & 0 & C_{44} & 0 & 0 \\ |
| 36 | +# 0 & 0 & 0 & 0 & C_{44} & 0 \\ |
| 37 | +# 0 & 0 & 0 & 0 & 0 & C_{44} \\ |
| 38 | +# \end{pmatrix}. |
| 39 | +# ``` |
| 40 | +# Thus we can just choose a suitable strain pattern $\dot{\eta} = (1, 0, 0, 1, 0, 0)^\top$, |
| 41 | +# such that $C\dot{\eta} = (C_{11}, C_{12}, C_{12}, C_{44}, 0, 0)^\top$. That is, |
| 42 | +# for cubic crystals like diamond silicon we obtain all independent elastic |
| 43 | +# constants from a single Jacobian-vector product on the stress-strain function. |
| 44 | +# |
| 45 | +# [^Nye1985]: |
| 46 | +# Nye, J. F. (1985). |
| 47 | +# *Physical Properties of Crystals*. Oxford University Press. |
| 48 | +# Comment: Since the elastic tensor transforms equivariantly under rotations, |
| 49 | +# its numerical components depend on the chosen Cartesian coordinate frame. |
| 50 | +# These tabulated patterns assume a standardized orientation of the structure |
| 51 | +# with respect to conventional crystallographic axes. |
| 52 | +# |
| 53 | +# This example computes the *clamped-ion* elastic tensor, keeping internal |
| 54 | +# atomic positions fixed under strain. The *relaxed-ion* tensor includes |
| 55 | +# additional corrections from internal relaxations, which can be obtained |
| 56 | +# from first-order atomic displacements in DFPT (see [^Wu2005]). |
| 57 | +# |
| 58 | +# [^Wu2005]: |
| 59 | +# Wu, X., Vanderbilt, D., & Hamann, D. R. (2005). |
| 60 | +# *Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory.* |
| 61 | +# [Physical Review B, 72(3), 035105](https://doi.org/10.1103/PhysRevB.72.035105). |
| 62 | + |
| 63 | + |
| 64 | +using DFTK |
| 65 | +using PseudoPotentialData |
| 66 | +using LinearAlgebra |
| 67 | +using ForwardDiff |
| 68 | +using DifferentiationInterface |
| 69 | +using AtomsBuilder |
| 70 | +using Unitful |
| 71 | +using UnitfulAtomic |
| 72 | + |
| 73 | + |
| 74 | +pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_4_1.standard.upf") |
| 75 | +a0_pbe = 10.33u"bohr" # Equilibrium lattice constant of silicon with PBE |
| 76 | +model0 = model_DFT(bulk(:Si; a=a0_pbe); pseudopotentials, functionals=PBE()) |
| 77 | + |
| 78 | +Ecut = recommended_cutoff(model0).Ecut |
| 79 | +kgrid = [4, 4, 4] |
| 80 | +tol = 1e-6 |
| 81 | + |
| 82 | +function symmetries_from_strain(model0, voigt_strain) |
| 83 | + lattice = DFTK.voigt_strain_to_full(voigt_strain) * model0.lattice |
| 84 | + model = Model(model0; lattice, symmetries=true) |
| 85 | + model.symmetries |
| 86 | +end |
| 87 | + |
| 88 | +strain_pattern = [1., 0., 0., 1., 0., 0.]; # recovers [c11, c12, c12, c44, 0, 0] |
| 89 | + |
| 90 | +# For elastic constants beyond the bulk modulus, symmetry-breaking strains |
| 91 | +# are required. That is, the symmetry group of the crystal is reduced. |
| 92 | +# Here we simply precompute the relevant subgroup by applying the automatic |
| 93 | +# symmetry detection (spglib) to the finitely perturbed crystal. |
| 94 | +symmetries_strain = symmetries_from_strain(model0, 0.01 * strain_pattern) |
| 95 | + |
| 96 | + |
| 97 | +function stress_from_strain(model0, voigt_strain; symmetries, Ecut, kgrid, tol) |
| 98 | + lattice = DFTK.voigt_strain_to_full(voigt_strain) * model0.lattice |
| 99 | + model = Model(model0; lattice, symmetries) |
| 100 | + basis = PlaneWaveBasis(model; Ecut, kgrid) |
| 101 | + scfres = self_consistent_field(basis; tol) |
| 102 | + DFTK.full_stress_to_voigt(compute_stresses_cart(scfres)) |
| 103 | +end |
| 104 | + |
| 105 | +stress_fn(voigt_strain) = stress_from_strain(model0, voigt_strain; |
| 106 | + symmetries=symmetries_strain, |
| 107 | + Ecut, kgrid, tol) |
| 108 | +stress, (dstress,) = value_and_pushforward(stress_fn, AutoForwardDiff(), |
| 109 | + zeros(6), (strain_pattern,)); |
| 110 | + |
| 111 | +# We can inspect the stress to verify it is small (close to equilibrium): |
| 112 | +stress |
| 113 | + |
| 114 | +# The response of the stress to `strain_pattern` contains the elastic constants |
| 115 | +# in atomic units, with the expected pattern $(c11, c12, c12, c44, 0, 0)$: |
| 116 | +dstress |
| 117 | + |
| 118 | +# Convert to GPa: |
| 119 | +println("C11: ", uconvert(u"GPa", dstress[1] * u"hartree" / u"bohr"^3)) |
| 120 | +println("C12: ", uconvert(u"GPa", dstress[2] * u"hartree" / u"bohr"^3)) |
| 121 | +println("C44: ", uconvert(u"GPa", dstress[4] * u"hartree" / u"bohr"^3)) |
| 122 | + |
| 123 | +# These results can be compared directly to finite differences of the stress_fn: |
| 124 | +h = 1e-3 |
| 125 | +dstress_fd = (stress_fn(h * strain_pattern) - stress_fn(-h * strain_pattern)) / 2h |
| 126 | +println("C11 (FD): ", uconvert(u"GPa", dstress_fd[1] * u"hartree" / u"bohr"^3)) |
| 127 | +println("C12 (FD): ", uconvert(u"GPa", dstress_fd[2] * u"hartree" / u"bohr"^3)) |
| 128 | +println("C44 (FD): ", uconvert(u"GPa", dstress_fd[4] * u"hartree" / u"bohr"^3)) |
| 129 | + |
| 130 | + |
| 131 | +# Here are AD-DFPT results from increasing discretization parameters: |
| 132 | +# |
| 133 | +# | Ecut | kgrid | c11 | c12 | c44 | |
| 134 | +# |------|---------------|-------:|------:|-------:| |
| 135 | +# | 18 | [4, 4, 4] | 156.51 | 59.57 | 98.61 | |
| 136 | +# | 18 | [8, 8, 8] | 153.53 | 56.90 | 100.07 | |
| 137 | +# | 24 | [8, 8, 8] | 153.26 | 56.82 | 99.97 | |
| 138 | +# | 24 | [14, 14, 14] | 153.03 | 56.71 | 100.09 | |
| 139 | +# |
| 140 | +# For comparison, Materials Project for PBE *relaxed-ion* elastic constants of |
| 141 | +# silicon [mp-149](https://next-gen.materialsproject.org/materials/mp-149): |
| 142 | +# c11 = 153 GPa, c12 = 57 GPa, c44 = 74 GPa. |
| 143 | +# Note the discrepancy in c44, which is due to us not yet including |
| 144 | +# ionic relaxation in this example. |
| 145 | + |
0 commit comments