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=kgrid, use_symmetry=false)
scfres_nosym = @time self_consistent_field(basis_nosym, tol=1e-8)

n     Energy            Eₙ-Eₙ₋₁     ρout-ρin   Diag
---   ---------------   ---------   --------   ----
  1   -7.868391602124         NaN   1.88e-01    2.9
  2   -7.872765950294   -4.37e-03   2.85e-02    1.9
  3   -7.872865299416   -9.93e-05   2.52e-03    1.0
  4   -7.872877863981   -1.26e-05   7.23e-04    2.1
  5   -7.872878169332   -3.05e-07   2.19e-04    1.2
  6   -7.872878211209   -4.19e-08   1.45e-05    1.2
  7   -7.872878213111   -1.90e-09   3.43e-06    2.2
  3.686766 seconds (1.10 M allocations: 576.247 MiB, 2.96% gc time)

and then redo it using symmetry (the default):

basis_sym = PlaneWaveBasis(model, Ecut; kgrid=kgrid)
scfres_sym = @time self_consistent_field(basis_sym, tol=1e-8)

n     Energy            Eₙ-Eₙ₋₁     ρout-ρin   Diag
---   ---------------   ---------   --------   ----
  1   -7.868203285517         NaN   1.88e-01    3.5
  2   -7.872767794516   -4.56e-03   2.87e-02    1.9
  3   -7.872866028605   -9.82e-05   2.58e-03    1.0
  4   -7.872877871864   -1.18e-05   7.03e-04    2.0
  5   -7.872878175695   -3.04e-07   2.07e-04    1.3
  6   -7.872878211689   -3.60e-08   7.44e-06    1.1
  7   -7.872878213117   -1.43e-09   2.31e-06    2.6
  0.728844 seconds (193.66 k allocations: 111.112 MiB)

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.ρ.real - scfres_nosym.ρ.real),
 norm(values(scfres_sym.energies) .- values(scfres_nosym.energies)))

(4.259199839071749e-6, 3.6821645913498963e-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 Array{Tuple{StaticArrays.SArray{Tuple{3,3},Int64,2,9},StaticArrays.SArray{Tuple{3},Float64,1,3}},1}:
 ([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 Array{Float64,2}:
 -0.12437    -0.12437
  0.0951518   0.095152
  0.176058    0.176058
  0.231686    0.231686
  0.366275    0.366276
  0.397683    0.397684
  0.405285    0.405285