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.1734770 
    AtomicLocal         -2.1476191
    AtomicNonlocal      1.5857327 
    Ewald               -8.3994719
    PspCorrection       -0.2947882
    Hartree             0.5588153 
    Xc                  -2.4030074

    total               -7.926861673456

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_DFT: Assemble a DFT model using any of the LDA or GGA functionals of the libxc library, for example: model_DFT(lattice, atoms, positions; functionals=LDA()) model_DFT(lattice, atoms, positions; functionals=[:lda_x, :lda_c_pw]) model_DFT(lattice, atoms, positions; functionals=[:gga_x_pbe, :gga_c_pbe]) For common functional combinations DFTK additionally offers shorthands. E.g. in the above example specifying LDA expands to [:lda_x, :lda_c_pw], such that the first two examples are identical. Note, that specifying no functional is the reduced Hartree-Fock model: model_DFT(lattice, atoms, positions; functionals=[])
• 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=15, kgrid=[4, 4, 4])

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             : (32, 32, 32), 32768 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.13061   , 5.13061   ]
                               [5.13061   , 0         , 5.13061   ]
                               [5.13061   , 5.13061   , 0         ]
        unit cell volume     : 270.11 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.0000000000000007

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.0000000000000009

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 =   749)
 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}:
 749
 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)

(32768, 32768)

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))

32768

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

4-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [0.0, 0.0, 0.0]
 [0.03125, 0.0, 0.0]
 [0.0625, 0.0, 0.0]
 [0.09375, 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}$.