# Data structures

In this section we assume a calculation of silicon LDA model in the setup described in Tutorial.

`Model`

datastructure

The physical model to be solved is defined by the `Model`

datastructure. It contains the unit cell, number of electrons, atoms, type of spin polarization and temperature. Each atom has an atomic type (`Element`

) specifying the number of valence electrons and the potential (or pseudopotential) it creates with respect to the electrons. The `Model`

structure also contains the list of energy terms defining the energy functional to be minimised during the SCF. For the silicon example above, we used an LDA model, which consists of the following terms^{[2]}:

`typeof.(model.term_types)`

```
7-element Vector{DataType}:
Kinetic{BlowupIdentity}
AtomicLocal
AtomicNonlocal
Ewald
PspCorrection
Hartree
Xc
```

DFTK computes energies for all terms of the model individually, which are available in `scfres.energies`

:

`scfres.energies`

```
Energy breakdown (in Ha):
Kinetic 3.1739355
AtomicLocal -2.1467809
AtomicNonlocal 1.5858684
Ewald -8.4004648
PspCorrection -0.2948928
Hartree 0.5586701
Xc -2.4032005
total -7.926865082194
```

For now the following energy terms are available in DFTK:

- Kinetic energy
- Local potential energy, either given by analytic potentials or specified by the type of atoms.
- Nonlocal potential energy, for norm-conserving pseudopotentials
- Nuclei energies (Ewald or pseudopotential correction)
- Hartree energy
- Exchange-correlation energy
- Power nonlinearities (useful for Gross-Pitaevskii type models)
- Magnetic field energy
- Entropy term

Custom types can be added if needed. For examples see the definition of the above terms in the `src/terms`

directory.

By mixing and matching these terms, the user can create custom models not limited to DFT. Convenience constructors are provided for common cases:

`model_LDA`

: LDA model using the Teter parametrisation`model_DFT`

: Assemble a DFT model using any of the LDA or GGA functionals of the libxc library, for example:`model_DFT(lattice, atoms, positions, [:gga_x_pbe, :gga_c_pbe]) model_DFT(lattice, atoms, positions, :lda_xc_teter93)`

where the latter is equivalent to`model_LDA`

. Specifying no functional is the reduced Hartree-Fock model:`model_DFT(lattice, atoms, positions, [])`

`model_atomic`

: A linear model, which contains no electron-electron interaction (neither Hartree nor XC term).

`PlaneWaveBasis`

and plane-wave discretisations

The `PlaneWaveBasis`

datastructure handles the discretization of a given `Model`

in a plane-wave basis. In plane-wave methods the discretization is twofold: Once the $k$-point grid, which determines the sampling *inside* the Brillouin zone and on top of that a finite plane-wave grid to discretise the lattice-periodic functions. The former aspect is controlled by the `kgrid`

argument of `PlaneWaveBasis`

, the latter is controlled by the cutoff energy parameter `Ecut`

:

`PlaneWaveBasis(model; Ecut, kgrid)`

```
PlaneWaveBasis discretization:
architecture : DFTK.CPU()
num. mpi processes : 1
num. julia threads : 1
num. DFTK threads : 1
num. blas threads : 2
num. fft threads : 1
Ecut : 15.0 Ha
fft_size : (30, 30, 30), 27000 total points
kgrid : MonkhorstPack([4, 4, 4])
num. red. kpoints : 64
num. irred. kpoints : 8
Discretized Model(lda_x+lda_c_pw, 3D):
lattice (in Bohr) : [0 , 5.13 , 5.13 ]
[5.13 , 0 , 5.13 ]
[5.13 , 5.13 , 0 ]
unit cell volume : 270.01 Bohr³
atoms : Si₂
atom potentials : ElementPsp(Si; psp="hgh/lda/si-q4")
ElementPsp(Si; psp="hgh/lda/si-q4")
num. electrons : 8
spin polarization : none
temperature : 0 Ha
terms : Kinetic()
AtomicLocal()
AtomicNonlocal()
Ewald(nothing)
PspCorrection()
Hartree()
Xc(lda_x, lda_c_pw)
```

The `PlaneWaveBasis`

by default uses symmetry to reduce the number of `k`

-points explicitly treated. For details see Crystal symmetries.

As mentioned, the periodic parts of Bloch waves are expanded in a set of normalized plane waves $e_G$:

\[\begin{aligned} \psi_{k}(x) &= e^{i k \cdot x} u_{k}(x)\\ &= \sum_{G \in \mathcal R^{*}} c_{G} e^{i k \cdot x} e_{G}(x) \end{aligned}\]

where $\mathcal R^*$ is the set of reciprocal lattice vectors. The $c_{{G}}$ are $\ell^{2}$-normalized. The summation is truncated to a "spherical", $k$-dependent basis set

\[ S_{k} = \left\{G \in \mathcal R^{*} \,\middle|\, \frac 1 2 |k+ G|^{2} \le E_\text{cut}\right\}\]

where $E_\text{cut}$ is the cutoff energy.

Densities involve terms like $|\psi_{k}|^{2} = |u_{k}|^{2}$ and therefore products $e_{-{G}} e_{{G}'}$ for ${G}, {G}'$ in $S_{k}$. To represent these we use a "cubic", $k$-independent basis set large enough to contain the set $\{{G}-G' \,|\, G, G' \in S_{k}\}$. We can obtain the coefficients of densities on the $e_{G}$ basis by a convolution, which can be performed efficiently with FFTs (see `ifft`

and `fft`

functions). Potentials are discretized on this same set.

The normalization conventions used in the code is that quantities stored in reciprocal space are coefficients in the $e_{G}$ basis, and quantities stored in real space use real physical values. This means for instance that wavefunctions in the real space grid are normalized as $\frac{|\Omega|}{N} \sum_{r} |\psi(r)|^{2} = 1$ where $N$ is the number of grid points.

For example let us check the normalization of the first eigenfunction at the first $k$-point in reciprocal space:

```
ψtest = scfres.ψ[1][:, 1]
sum(abs2.(ψtest))
```

`1.0000000000000002`

We now perform an IFFT to get ψ in real space. The $k$-point has to be passed because ψ is expressed on the $k$-dependent basis. Again the function is normalised:

```
ψreal = ifft(basis, basis.kpoints[1], ψtest)
sum(abs2.(ψreal)) * basis.dvol
```

`1.0000000000000002`

The list of $k$ points of the basis can be obtained with `basis.kpoints`

.

`basis.kpoints`

```
8-element Vector{Kpoint{Float64, Vector{StaticArraysCore.SVector{3, Int64}}}}:
KPoint([ 0, 0, 0], spin = 1, num. G vectors = 725)
KPoint([ 0.25, 0, 0], spin = 1, num. G vectors = 754)
KPoint([ -0.5, 0, 0], spin = 1, num. G vectors = 754)
KPoint([ 0.25, 0.25, 0], spin = 1, num. G vectors = 729)
KPoint([ -0.5, 0.25, 0], spin = 1, num. G vectors = 748)
KPoint([ -0.25, 0.25, 0], spin = 1, num. G vectors = 754)
KPoint([ -0.5, -0.5, 0], spin = 1, num. G vectors = 740)
KPoint([ -0.25, -0.5, 0.25], spin = 1, num. G vectors = 744)
```

The $G$ vectors of the "spherical", $k$-dependent grid can be obtained with `G_vectors`

:

`[length(G_vectors(basis, kpoint)) for kpoint in basis.kpoints]`

```
8-element Vector{Int64}:
725
754
754
729
748
754
740
744
```

```
ik = 1
G_vectors(basis, basis.kpoints[ik])[1:4]
```

```
4-element Vector{StaticArraysCore.SVector{3, Int64}}:
[0, 0, 0]
[1, 0, 0]
[2, 0, 0]
[3, 0, 0]
```

The list of $G$ vectors (Fourier modes) of the "cubic", $k$-independent basis set can be obtained similarly with `G_vectors(basis)`

.

`length(G_vectors(basis)), prod(basis.fft_size)`

`(27000, 27000)`

`collect(G_vectors(basis))[1:4]`

```
4-element Vector{StaticArraysCore.SVector{3, Int64}}:
[0, 0, 0]
[1, 0, 0]
[2, 0, 0]
[3, 0, 0]
```

Analogously the list of $r$ vectors (real-space grid) can be obtained with `r_vectors(basis)`

:

`length(r_vectors(basis))`

`27000`

`collect(r_vectors(basis))[1:4]`

```
4-element Vector{StaticArraysCore.SVector{3, Float64}}:
[0.0, 0.0, 0.0]
[0.03333333333333333, 0.0, 0.0]
[0.06666666666666667, 0.0, 0.0]
[0.1, 0.0, 0.0]
```

## Accessing Bloch waves and densities

Wavefunctions are stored in an array `scfres.ψ`

as `ψ[ik][iG, iband]`

where `ik`

is the index of the $k$-point (in `basis.kpoints`

), `iG`

is the index of the plane wave (in `G_vectors(basis, basis.kpoints[ik])`

) and `iband`

is the index of the band. Densities are stored in real space, as a 4-dimensional array (the third being the spin component).

```
rvecs = collect(r_vectors(basis))[:, 1, 1] # slice along the x axis
x = [r[1] for r in rvecs] # only keep the x coordinate
plot(x, scfres.ρ[:, 1, 1, 1], label="", xlabel="x", ylabel="ρ", marker=2)
```

```
G_energies = [sum(abs2.(model.recip_lattice * G)) ./ 2 for G in G_vectors(basis)][:]
scatter(G_energies, abs.(fft(basis, scfres.ρ)[:]);
yscale=:log10, ylims=(1e-12, 1), label="", xlabel="Energy", ylabel="|ρ|")
```

Note that the density has no components on wavevectors above a certain energy, because the wavefunctions are limited to $\frac 1 2|k+G|^2 ≤ E_{\rm cut}$.

- 2If you are not familiar with Julia syntax,
`typeof.(model.term_types)`

is equivalent to`[typeof(t) for t in model.term_types]`

.