# 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
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.1521640
AtomicLocal         -2.1284083
AtomicNonlocal      1.5890982
Ewald               -8.4004648
PspCorrection       -0.2948928
Hartree             0.5469838
Xc                  -2.3964736

total               -7.931993463235

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, [:gga_x_pbe, :gga_c_pbe]) model_DFT(lattice, atoms, :lda_xc_teter93) where the latter is equivalent to model_LDA. Specifying no functional is the reduced Hartree-Fock model: model_DFT(lattice, atoms, [])
• 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:
Ecut                 : 15.0 Ha
fft_size             : (27, 27, 27)
kgrid type           : Monkhorst-Pack
kgrid                : [4, 4, 4]
kshift               : [0.5, 0.5, 0.5]
num. irred. kpoints  : 10

Discretized Model(lda_xc_teter93, 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)

num. electrons       : 8
spin polarization    : none
temperature          : 0 Ha

terms                : Kinetic()
AtomicLocal()
AtomicNonlocal()
Ewald()
PspCorrection()
Hartree()
Xc(:lda_xc_teter93)

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 G_to_r and r_to_G 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))
0.9999999999999998

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 = G_to_r(basis, basis.kpoints[1], ψtest)
sum(abs2.(ψreal)) * basis.dvol
0.9999999999999998

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

basis.kpoints
10-element Vector{Kpoint{Float64}}:
Kpoint{Float64}(1, [0.125, 0.125, 0.125], [0.07654952859624253, 0.07654952859624253, 0.07654952859624253], [1, 2, 3, 4, 5, 6, 23, 24, 25, 26  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 428, 2108 => 231, 17437 => 488, 1546 => 209, 16655 => 427, 699 => 78, 18089 => 539, 673 => 69, 3705 => 388, 17417 => 487…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [0.375, 0.125, 0.125], [-0.07654952859624253, 0.2296485857887276, 0.2296485857887276], [1, 2, 3, 4, 5, 23, 24, 25, 26, 27  …  19658, 19659, 19660, 19661, 19678, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 427, 2108 => 234, 17437 => 490, 16655 => 426, 699 => 77, 18089 => 539, 673 => 69, 3705 => 388, 17417 => 489, 18116 => 545…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [-5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-6, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [-0.375, 0.125, 0.125], [0.38274764298121267, -0.2296485857887276, -0.2296485857887276], [1, 2, 3, 4, 5, 6, 24, 25, 26, 27  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 429, 2108 => 231, 1546 => 208, 16655 => 428, 699 => 75, 18089 => 544, 673 => 66, 16797 => 468, 3705 => 388, 17417 => 492…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [-0.125, 0.125, 0.125], [0.2296485857887276, -0.07654952859624253, -0.07654952859624253], [1, 2, 3, 4, 5, 6, 24, 25, 26, 27  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 432, 2108 => 231, 17437 => 495, 1546 => 209, 16655 => 431, 699 => 76, 18089 => 546, 673 => 67, 16797 => 470, 3705 => 389…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [0.375, 0.375, 0.125], [0.07654952859624253, 0.07654952859624253, 0.38274764298121267], [1, 2, 3, 4, 5, 23, 24, 25, 26, 27  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 429, 2108 => 231, 17437 => 489, 16655 => 428, 699 => 76, 18089 => 537, 19546 => 704, 673 => 67, 3705 => 385, 17417 => 488…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [-5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [-0.375, 0.375, 0.125], [0.5358467001736977, -0.38274764298121267, -0.07654952859624253], [1, 2, 3, 4, 5, 6, 24, 25, 26, 27  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 429, 2108 => 232, 1546 => 211, 16655 => 428, 699 => 76, 18089 => 538, 19546 => 702, 673 => 67, 16797 => 467, 3705 => 387…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [-0.125, 0.375, 0.125], [0.38274764298121267, -0.2296485857887276, 0.07654952859624253], [1, 2, 3, 4, 5, 6, 24, 25, 26, 27  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 425, 2108 => 232, 17437 => 487, 1546 => 210, 16655 => 424, 699 => 75, 18089 => 535, 19546 => 702, 673 => 66, 3705 => 387…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [-0.375, -0.375, 0.125], [0.07654952859624253, 0.07654952859624253, -0.5358467001736977], [1, 2, 3, 4, 5, 6, 24, 25, 26, 27  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 430, 2108 => 236, 1546 => 212, 699 => 72, 673 => 64, 16797 => 467, 3705 => 385, 17417 => 487, 18116 => 542, 112 => 42…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [0.375, 0.375, 0.375], [0.2296485857887276, 0.2296485857887276, 0.2296485857887276], [1, 2, 3, 4, 5, 6, 23, 24, 25, 26  …  19657, 19658, 19659, 19660, 19661, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 423, 2108 => 233, 17437 => 493, 16655 => 422, 699 => 80, 18089 => 545, 19546 => 717, 673 => 71, 18116 => 552, 17417 => 492…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0]  …  [0, -1, -1], [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])
Kpoint{Float64}(1, [-0.375, 0.375, 0.375], [0.6889457573661828, -0.2296485857887276, -0.2296485857887276], [1, 2, 3, 4, 5, 6, 24, 25, 26, 27  …  19658, 19659, 19660, 19661, 19662, 19679, 19680, 19681, 19682, 19683], Dict(16656 => 413, 2108 => 228, 1546 => 208, 16655 => 412, 699 => 75, 18089 => 528, 673 => 67, 16797 => 452, 3705 => 372, 17417 => 477…), StaticArrays.SVector{3, Int64}[[0, 0, 0], [1, 0, 0], [2, 0, 0], [3, 0, 0], [4, 0, 0], [5, 0, 0], [-4, 0, 0], [-3, 0, 0], [-2, 0, 0], [-1, 0, 0]  …  [1, -1, -1], [2, -1, -1], [3, -1, -1], [4, -1, -1], [5, -1, -1], [-5, -1, -1], [-4, -1, -1], [-3, -1, -1], [-2, -1, -1], [-1, -1, -1]])

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]
10-element Vector{Int64}:
743
748
753
757
751
748
747
747
763
744
ik = 1
G_vectors(basis, basis.kpoints[ik])[1:4]
4-element Vector{StaticArrays.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)
(19683, 19683)
collect(G_vectors(basis))[1:4]
4-element Vector{StaticArrays.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))
19683
collect(r_vectors(basis))[1:4]
4-element Vector{StaticArrays.SVector{3, Float64}}:
[0.0, 0.0, 0.0]
[0.037037037037037035, 0.0, 0.0]
[0.07407407407407407, 0.0, 0.0]
[0.1111111111111111, 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.(r_to_G(basis, scfres.ρ)[:]);
yscale=:log10, ylims=(1e-12, 1), label="", xlabel="Energy", ylabel="|ρ|^2")