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 $\widetilde{S}$ and a real-space vector $\widetilde{\tau}$ such that
\[\widetilde{S} \mathcal{C} + \widetilde{\tau} = \mathcal{C}.\]
The symmetries where $\widetilde{S} = 1$ and $\widetilde{\tau}$ 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( \widetilde{S} x + \widetilde{\tau} \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 \widetilde{\tau}} e^{i (\widetilde{S}^T q) x}\\ &= e^{- i(S q) \cdot \tau } e^{i (S q) \cdot x} \end{aligned}\]
where we set
\[\begin{aligned} S &= \widetilde{S}^{T}\\ \tau &= -\widetilde{S}^{-1}\widetilde{\tau}. \end{aligned}\]
(these equations being also valid in reduced coordinates).
It follows that the Fourier transform satisfies
\[\widehat{Uu}(q) = e^{- iq \cdot \tau} \widehat u(S^{-1} q)\]
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.
Example
Let us demonstrate this in practice. We consider silicon, setup appropriately in the lattice and atoms 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 = model_LDA(lattice, atoms)
basis_nosym = PlaneWaveBasis(model; Ecut, kgrid, use_symmetry=false)
scfres_nosym = @time self_consistent_field(basis_nosym, tol=1e-8)n Energy Eₙ-Eₙ₋₁ ρout-ρin α Diag --- --------------- --------- -------- ---- ---- 1 -7.868422459138 NaN 1.88e-01 0.80 2.9 2 -7.872625462555 -4.20e-03 2.85e-02 0.80 1.0 3 -7.872877492267 -2.52e-04 2.59e-03 0.80 2.8 4 -7.872878197932 -7.06e-07 4.24e-04 0.80 2.4 5 -7.872878211874 -1.39e-08 2.97e-05 0.80 1.1 6 -7.872878213106 -1.23e-09 5.06e-06 0.80 3.0 5.519218 seconds (1.58 M allocations: 556.791 MiB, 2.89% gc time, 2.93% compilation time)
and then redo it using symmetry (the default):
basis_sym = PlaneWaveBasis(model; Ecut, kgrid)
scfres_sym = @time self_consistent_field(basis_sym, tol=1e-8)n Energy Eₙ-Eₙ₋₁ ρout-ρin α Diag --- --------------- --------- -------- ---- ---- 1 -7.868305338837 NaN 1.88e-01 0.80 3.7 2 -7.872644471338 -4.34e-03 2.89e-02 0.80 1.0 3 -7.872877290873 -2.33e-04 2.54e-03 0.80 2.8 4 -7.872878194676 -9.04e-07 4.46e-04 0.80 2.6 5 -7.872878211574 -1.69e-08 4.08e-05 0.80 1.5 6 -7.872878213101 -1.53e-09 4.99e-06 0.80 2.7 0.975307 seconds (220.84 k allocations: 102.534 MiB, 4.16% 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))(10, 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)))(9.820816713944777e-6, 8.051028090281091e-6)
To demonstrate the mapping between k-points due to symmetry, we pick an arbitrary k-point in the irreducible Brillouin zone:
ikpt_irred = 2
kpt_irred_coord = basis_sym.kpoints[ikpt_irred].coordinate
basis_sym.ksymops[ikpt_irred]6-element Vector{Tuple{StaticArrays.SMatrix{3, 3, Int64, 9}, StaticArrays.SVector{3, Float64}}}:
([1 0 0; 0 1 0; 0 0 1], [0.0, 0.0, 0.0])
([0 0 1; 1 0 0; 0 1 0], [0.0, 0.0, 0.0])
([0 1 0; 0 0 1; 1 0 0], [0.0, 0.0, 0.0])
([0 0 -1; 0 -1 0; -1 0 0], [0.0, 0.0, 0.0])
([0 -1 0; -1 0 0; 0 0 -1], [0.0, 0.0, 0.0])
([-1 0 0; 0 0 -1; 0 -1 0], [0.0, 0.0, 0.0])This is a list of all symmetries operations $(S, \tau)$ that can be used to map this irreducible $k$-point to reducible $k$-points. Let's pick the third symmetry operation of this $k$-point and check.
S, τ = basis_sym.ksymops[ikpt_irred][3]
kpt_red_coord = S * basis_sym.kpoints[ikpt_irred].coordinate
ikpt_red = findfirst(kcoord -> kcoord ≈ kpt_red_coord,
[k.coordinate for k in basis_nosym.kpoints])
[scfres_sym.eigenvalues[ikpt_irred] scfres_nosym.eigenvalues[ikpt_red]]7×2 Matrix{Float64}:
-0.12437 -0.12437
0.0951514 0.0951519
0.176057 0.176058
0.231686 0.231686
0.366275 0.366275
0.397683 0.397684
0.405285 0.405285