Implement stresses without unfold_bz

docs/src/advanced/symmetries.md

@@ -51,6 +51,25 @@ one can find a reduced set of ``k``-points
 (the *irreducible ``k``-points*) such that the eigenvectors at the
 reducible ``k``-points can be deduced from those at the irreducible ``k``-points.
 
+### Symmetrization
+Quantities that are calculated by summing over the reducible ``k`` points can be
+calculated by first summing over the irreducible ``k`` points and then symmetrizing.
+In this subsection, we denote by ``S`` the combined rotation and fractional transformation.
+If ``S`` is the symmetry of the system, ``f(Sk) = S(f(k))`` holds, where ``S(f)`` is the
+symmetry transformation of ``f`` which is a linear function and depends on the quantity
+``f``. Then, we find
+```math
+\sum_{k\ \mathrm{reducible}} f(k)
+= \sum_{k\ \mathrm{irreducible}} \sum_{S\text{ that maps $k$ to reducible $k$}} S(f(k))
+= \sum_{k\ \mathrm{irreducible}} \frac{N_S}{N_{S,k}} \sum_{S} S(f(k))
+= S \left( \sum_{k\ \mathrm{irreducible}} \frac{N_S}{N_{S,k}} f(k) \right).
+```
+Here, ``N_S`` and ``N_{S,k}`` are the total number of symmetry operations and the
+number of operations such that ``k=Sk``, respectively. The latter operations form
+a subgroup of the group of all symmetry operations. This subgroup is often called
+the "small/little group of ``k``".
+The factor ``\frac{N_S}{N_{S,k}}`` determines the weight of each ``k`` point.
+
 ## Example
 ```@setup symmetries
 using DFTK

src/postprocess/stresses.jl

@@ -4,7 +4,6 @@ Compute the stresses (= 1/Vol dE/d(M*lattice), taken at M=I) of an obtained SCF
 """
 @timing function compute_stresses(scfres)
     # TODO optimize by only computing derivatives wrt 6 independent parameters
-    scfres = unfold_bz(scfres)
     # compute the Hellmann-Feynman energy (with fixed ψ/occ/ρ)
     function HF_energy(lattice)
         T = eltype(lattice)
@@ -31,5 +30,6 @@ Compute the stresses (= 1/Vol dE/d(M*lattice), taken at M=I) of an obtained SCF
     end
     L = scfres.basis.model.lattice
     Ω = scfres.basis.model.unit_cell_volume
-    ForwardDiff.gradient(M -> HF_energy((I+M) * L), zero(L)) / Ω
+    stresses_not_symmetrized = ForwardDiff.gradient(M -> HF_energy((I+M) * L), zero(L)) / Ω
+    symmetrize_stresses(L, scfres.basis.symmetries, stresses_not_symmetrized)
 end

src/symmetry.jl

@@ -40,6 +40,8 @@
 # - The set of symmetry operations that we use to reduce the
 #   reducible Brillouin zone (RBZ) to the irreducible (IBZ) (basis.ksymops)
 
+# See https://juliamolsim.github.io/DFTK.jl/stable/advanced/symmetries for details.
+
 @doc raw"""
 Return the ``k``-point symmetry operations associated to a lattice and atoms.
 """
@@ -314,3 +316,14 @@ function unfold_bz(scfres)
     new_scfres = (; basis=basis_unfolded, ψ, ham, eigenvalues, occupation)
     merge(scfres, new_scfres)
 end
+
+# symmetrize the stress tensor which is a rank-2 contravariant tensor in reduced coordinates
+function symmetrize_stresses(lattice, symmetries, stresses::Mat3)
+    stresses_symmetrized = zero(stresses)
+    for (S, τ) in symmetries
+        S_reduced = inv(lattice) * S * lattice
+        stresses_symmetrized += S_reduced' * stresses * S_reduced
+    end
+    stresses_symmetrized /= length(symmetries)
+    stresses_symmetrized
+end

test/runtests.jl

@@ -91,6 +91,7 @@ Random.seed!(0)
         include("energies_guess_density.jl")
         include("compute_density.jl")
         include("forces.jl")
+        include("stresses.jl")
     end
 
     if "all" in TAGS
@@ -107,9 +108,5 @@ Random.seed!(0)
         include("aqua.jl")
     end
 
-    if "all" in TAGS
-        mpi_nprocs() == 1 && include("stresses.jl")
-    end
-
     ("example" in TAGS) && include("runexamples.jl")
 end

test/stresses.jl

@@ -2,6 +2,7 @@ using Test
 using DFTK
 using ForwardDiff
 import FiniteDiff
+using MPI
 include("testcases.jl")
 
 # Hellmann-Feynman stress
@@ -12,7 +13,7 @@ include("testcases.jl")
         Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
         atoms = [Si => silicon.positions]
         model = model_DFT(lattice, atoms, [:lda_x, :lda_c_vwn]; symmetries)
-        kgrid = [2, 2, 1]
+        kgrid = [3, 3, 1]
         Ecut = 7
         PlaneWaveBasis(model; Ecut, kgrid)
     end
@@ -42,7 +43,7 @@ include("testcases.jl")
     stresses = compute_stresses(scfres)
     @test isapprox(stresses, compute_stresses(scfres_nosym), atol=1e-10)
 
-    dir = randn(3, 3)
+    dir = MPI.bcast(randn(3, 3), 0, MPI.COMM_WORLD)
 
     dE_stresses = dot(dir, stresses) * scfres.basis.model.unit_cell_volume
     ref_recompute = FiniteDiff.finite_difference_derivative(0.0) do ε