review comments

Author: niklasschmitz

examples/elastic_constants.jl

@@ -1,13 +1,14 @@
 # # 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 compute *clamped-ion* elastic constants of a crystal using
+# the algorithmic differentiation density-functional perturbation theory (AD-DFPT) approach
+# as introduced in [^SPH25].
 # 
 # [^SPH25]: 
 #     Schmitz, N. F., Ploumhans, B., & Herbst, M. F. (2025)
 #     *Algorithmic differentiation for plane-wave DFT: materials design, error control and learning model parameters.*
 #     [arXiv:2509.07785](https://arxiv.org/abs/2509.07785)
-
+#
 # 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.
@@ -24,8 +25,8 @@
 # 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
+# For example, cubic crystals have only three independent elastic constants
+# $C_{11}$, $C_{12}$ and $C_{44}$, with the pattern
 # ```math
 # C = \begin{pmatrix}
 #   C_{11} & C_{12} & C_{12} & 0      & 0      & 0 \\
@@ -84,7 +85,7 @@ function symmetries_from_strain(model0, voigt_strain)
     model.symmetries
 end
 
-strain_pattern = [1., 0., 0., 1., 0., 0.];  # should yield [c11, c12, c12, c44, 0, 0] for cubic crystal
+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.
@@ -101,22 +102,31 @@ function stress_from_strain(model0, voigt_strain; symmetries, Ecut, kgrid, 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,));
-@show stress dstress;
+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,));
+
+# We can inspect the stress to verify it is small (close to equilibrium):
+stress
 
-c11 = uconvert(u"GPa", dstress[1] * u"hartree" / u"bohr"^3)
-c12 = uconvert(u"GPa", dstress[2] * u"hartree" / u"bohr"^3)
-c44 = uconvert(u"GPa", dstress[4] * u"hartree" / u"bohr"^3)
-@show c11 c12 c44;
+# 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
 
-# These results can be compared directly to finite differences of the stress-strain relation:
+# 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))
+
+# 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
-c11_fd = uconvert(u"GPa", dstress_fd[1] * u"hartree" / u"bohr"^3)
-c12_fd = uconvert(u"GPa", dstress_fd[2] * u"hartree" / u"bohr"^3)
-c44_fd = uconvert(u"GPa", dstress_fd[4] * u"hartree" / u"bohr"^3)
-@show c11_fd c12_fd c44_fd;
+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))
+
 
 # Here are AD-DFPT results from increasing discretization parameters:
 #
@@ -127,7 +137,9 @@ c44_fd = uconvert(u"GPa", dstress_fd[4] * u"hartree" / u"bohr"^3)
 # | 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](https://next-gen.materialsproject.org/materials/mp-149):
+# For comparison, Materials Project for PBE *relaxed-ion* elastic constants of
+# silicon [mp-149](https://next-gen.materialsproject.org/materials/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.
+# Note the discrepancy in c44, which is due to us not yet including
+# ionic relaxation in this example.