Crystal symmetries
Theory
In this discussion we will only describe the situation for a monoatomic crystal $\mathcal C \subset \mathbb R^3$, the extension being easy. A symmetry of the crystal is an orthogonal matrix $W$ and a real-space vector $w$ such that
\[W \mathcal{C} + w = \mathcal{C}.\]
The symmetries where $W = 1$ and $w$ is a lattice vector are always assumed and ignored in the following.
We can define a corresponding unitary operator $U$ on $L^2(\mathbb R^3)$ with action
\[ (Uu)(x) = u\left( W x + w \right).\]
We assume that the atomic potentials are radial and that any self-consistent potential also respects this symmetry, so that $U$ commutes with the Hamiltonian.
This operator acts on a plane-wave as
\[\begin{aligned} (U e^{iq\cdot x}) (x) &= e^{iq \cdot w} e^{i (W^T q) x}\\ &= e^{- i(S q) \cdot \tau } e^{i (S q) \cdot x} \end{aligned}\]
where we set
\[\begin{aligned} S &= W^{T}\\ \tau &= -W^{-1}w. \end{aligned}\]
It follows that the Fourier transform satisfies
\[\widehat{Uu}(q) = e^{- iq \cdot \tau} \widehat u(S^{-1} q)\]
(all of these equations being also valid in reduced coordinates). In particular, if $e^{ik\cdot x} u_{k}(x)$ is an eigenfunction, then by decomposing $u_k$ over plane-waves $e^{i G \cdot x}$ one can see that $e^{i(S^T k) \cdot x} (U u_k)(x)$ is also an eigenfunction: we can choose
\[u_{Sk} = U u_k.\]
This is used to reduce the computations needed. For a uniform sampling of the Brillouin zone (the reducible $k$-points), 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. Let $\mathcal{K}_\text{reducible}$ denote the reducible $k$-points sampling the Brillouin zone, $\mathcal{S}$ be the group of all crystal symmetries that leave this BZ mesh invariant ($\mathcal{S}\mathcal{K}_\text{reducible} = \mathcal{K}_\text{reducible}$) and $\mathcal{K}$ be the irreducible $k$-points obtained from $\mathcal{K}_\text{reducible}$ using the symmetries $\mathcal{S}$. Clearly
\[\mathcal{K}_\text{red} = \{Sk \, | \, S \in \mathcal{S}, k \in \mathcal{K}\}.\]
Let $Q$ be a $k$-dependent quantity to sum (for instance, energies, densities, forces, etc). $Q$ transforms in a particular way under symmetries: $Q(Sk) = S(Q(k))$ where the (linear) action of $S$ on $Q$ depends on the particular $Q$.
\[\begin{aligned} \sum_{k \in \mathcal{K}_\text{red}} Q(k) &= \sum_{k \in \mathcal{K}} \ \sum_{S \text{ with } Sk \in \mathcal{K}_\text{red}} S(Q(k)) \\ &= \sum_{k \in \mathcal{K}} \frac{1}{N_{\mathcal{S},k}} \sum_{S \in \mathcal{S}} S(Q(k))\\ &= \frac{1}{N_{\mathcal{S}}} \sum_{S \in \mathcal{S}} \left(\sum_{k \in \mathcal{K}} \frac{N_\mathcal{S}}{N_{\mathcal{S},k}} Q(k) \right) \end{aligned}\]
Here, $N_\mathcal{S} = |\mathcal{S}|$ is the total number of symmetry operations and $N_{\mathcal{S},k}$ denotes the number of operations such that leave $k$ invariant. The latter operations form a subgroup of the group of all symmetry operations, sometimes called the small/little group of $k$. The factor $\frac{N_\mathcal{S}}{N_{S,k}}$, also equal to the ratio of number of reducible points encoded by this particular irreducible $k$ to the total number of reducible points, determines the weight of each irreducible $k$ point.
Example
Let us demonstrate this in practice. We consider silicon, setup appropriately in the lattice, atoms and positions objects as in Tutorial and to reach a fast execution, we take a small Ecut of 5 and a [4, 4, 4] Monkhorst-Pack grid. First we perform the DFT calculation disabling symmetry handling
model_nosym = model_DFT(lattice, atoms, positions; functionals=LDA(), symmetries=false)
basis_nosym = PlaneWaveBasis(model_nosym; Ecut=5, kgrid=[4, 4, 4])
scfres_nosym = @time self_consistent_field(basis_nosym, tol=1e-6)n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -7.864804325128 -0.72 3.6 326ms
2 -7.868982466382 -2.38 -1.54 1.0 144ms
3 -7.869173534157 -3.72 -2.58 1.0 142ms
4 -7.869208882023 -4.45 -2.86 2.4 243ms
5 -7.869209489694 -6.22 -3.05 1.0 140ms
6 -7.869209804823 -6.50 -4.67 1.0 149ms
7 -7.869209817370 -7.90 -4.72 2.5 267ms
8 -7.869209817521 -9.82 -5.66 1.0 139ms
9 -7.869209817529 -11.10 -5.61 2.0 203ms
10 -7.869209817531 -11.84 -6.65 1.0 141ms
1.901616 seconds (1.75 M allocations: 487.834 MiB, 4.36% gc time)and then redo it using symmetry (the default):
model_sym = model_DFT(lattice, atoms, positions; functionals=LDA())
basis_sym = PlaneWaveBasis(model_sym; Ecut=5, kgrid=[4, 4, 4])
scfres_sym = @time self_consistent_field(basis_sym, tol=1e-6)n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -7.865038838251 -0.72 4.5 67.9ms
2 -7.869030569584 -2.40 -1.54 1.0 22.8ms
3 -7.869183043963 -3.82 -2.59 1.1 23.5ms
4 -7.869209110335 -4.58 -2.89 2.6 33.8ms
5 -7.869209528597 -6.38 -3.06 1.0 23.0ms
6 -7.869209806191 -6.56 -4.05 1.0 22.9ms
7 -7.869209817241 -7.96 -4.98 1.8 28.1ms
8 -7.869209817520 -9.55 -5.24 2.1 32.2ms
9 -7.869209817528 -11.10 -5.73 1.0 41.6ms
10 -7.869209817529 -12.33 -5.62 1.1 23.9ms
11 -7.869209817530 -11.79 -6.27 1.0 29.8ms
0.355571 seconds (287.95 k allocations: 107.510 MiB, 10.26% gc time)Clearly both yield the same energy but the version employing symmetry is faster, since less $k$-points are explicitly treated:
(length(basis_sym.kpoints), length(basis_nosym.kpoints))(8, 64)Both SCFs would even agree in the convergence history if exact diagonalization was used for the eigensolver in each step of both SCFs. But since DFTK adjusts this diagtol value adaptively during the SCF to increase performance, a slightly different history is obtained. Try adding the keyword argument determine_diagtol=(args...; kwargs...) -> 1e-8 in each SCF call to fix the diagonalization tolerance to be 1e-8 for all SCF steps, which will result in an almost identical convergence history.
We can also explicitly verify both methods to yield the same density:
(norm(scfres_sym.ρ - scfres_nosym.ρ),
norm(values(scfres_sym.energies) .- values(scfres_nosym.energies)))(3.582497924570466e-7, 8.03015132598893e-7)The symmetries can be used to map reducible to irreducible points:
ikpt_red = rand(1:length(basis_nosym.kpoints))
# find a (non-unique) corresponding irreducible point in basis_nosym,
# and the symmetry that relates them
ikpt_irred, symop = DFTK.unfold_mapping(basis_sym, basis_nosym.kpoints[ikpt_red])
[basis_sym.kpoints[ikpt_irred].coordinate symop.S * basis_nosym.kpoints[ikpt_red].coordinate]3×2 StaticArraysCore.SMatrix{3, 2, Float64, 6} with indices SOneTo(3)×SOneTo(2):
-0.25 0.75
0.25 0.25
0.0 0.0The eigenvalues match also:
[scfres_sym.eigenvalues[ikpt_irred] scfres_nosym.eigenvalues[ikpt_red]]7×2 Matrix{Float64}:
-0.0935298 -0.0935298
0.0674043 0.0674044
0.122319 0.122319
0.212029 0.212029
0.355694 0.355694
0.438697 0.438726
0.454756 0.454751