# 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.868354622367 NaN 1.88e-01 2.9 2 -7.872614925719 -4.26e-03 2.85e-02 1.0 3 -7.872845946845 -2.31e-04 2.57e-03 1.2 4 -7.872877511842 -3.16e-05 1.16e-03 2.4 5 -7.872878044605 -5.33e-07 5.82e-04 1.0 6 -7.872878206595 -1.62e-07 1.41e-05 1.0 7 -7.872878213103 -6.51e-09 6.08e-06 2.8 5.131100 seconds (1.21 M allocations: 612.371 MiB, 3.87% 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.868309908741 NaN 1.88e-01 3.6 2 -7.872647101844 -4.34e-03 2.89e-02 1.0 3 -7.872846898754 -2.00e-04 2.46e-03 1.1 4 -7.872877605032 -3.07e-05 1.20e-03 2.4 5 -7.872878009782 -4.05e-07 7.06e-04 1.0 6 -7.872878207075 -1.97e-07 6.42e-05 1.0 7 -7.872878213008 -5.93e-09 8.54e-06 1.8 0.924420 seconds (197.08 k allocations: 112.294 MiB, 2.84% 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.ρ.real - scfres_nosym.ρ.real),
norm(values(scfres_sym.energies) .- values(scfres_nosym.energies)))
```

(4.410269542873316e-6, 2.7795709514438916e-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.124368 -0.124369 0.0951535 0.0951526 0.176059 0.176059 0.231688 0.231687 0.366277 0.366276 0.397684 0.397684 0.405286 0.405286