# 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.864822577779 -0.72 3.6 351ms
2 -7.868987894933 -2.38 -1.54 1.0 152ms
3 -7.869172945149 -3.73 -2.58 1.0 197ms
4 -7.869208850165 -4.44 -2.86 2.4 254ms
5 -7.869209467288 -6.21 -3.03 1.0 153ms
6 -7.869209804122 -6.47 -4.61 1.0 153ms
7 -7.869209817287 -7.88 -4.62 2.5 271ms
8 -7.869209817518 -9.64 -5.63 1.0 156ms
9 -7.869209817527 -11.05 -5.54 2.1 259ms
10 -7.869209817530 -11.49 -6.21 1.0 175ms
2.128302 seconds (2.48 M allocations: 601.878 MiB, 3.22% 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.865080711593 -0.72 4.4 58.9ms
2 -7.869035285294 -2.40 -1.54 1.0 27.4ms
3 -7.869184338641 -3.83 -2.59 1.0 27.3ms
4 -7.869208928577 -4.61 -2.90 2.1 37.7ms
5 -7.869209560490 -6.20 -3.11 1.1 28.4ms
6 -7.869209803417 -6.61 -4.02 1.0 27.9ms
7 -7.869209817174 -7.86 -4.97 1.8 34.8ms
8 -7.869209817506 -9.48 -5.01 1.9 37.5ms
9 -7.869209817509 -11.51 -5.05 1.0 28.4ms
10 -7.869209817528 -10.72 -5.53 1.0 27.8ms
11 -7.869209817531 -11.59 -7.05 1.0 28.1ms
0.368754 seconds (372.57 k allocations: 122.515 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))`

`(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.8882434054443205e-7, 8.517978391272548e-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.0
0.0 0.25
0.0 0.0
```

The eigenvalues match also:

`[scfres_sym.eigenvalues[ikpt_irred] scfres_nosym.eigenvalues[ikpt_red]]`

```
7×2 Matrix{Float64}:
-0.14075 -0.14075
0.120309 0.120308
0.233237 0.233237
0.233237 0.233237
0.34283 0.34283
0.391626 0.391626
0.391626 0.391626
```