DFTK is a Julia package for playing with plane-wave density-functional theory algorithms. In its basic formulation it solves periodic Kohn-Sham equations.

This document provides an overview of the structure of the code and how to access basic information about calculations. Basic familiarity with the concepts of plane-wave density functional theory is assumed throughout. Feel free to take a look at the Periodic problems or the Introductory resources chapters for some introductory material on the topic.

!!! note "Convergence parameters in the documentation" We use rough parameters in order to be able to automatically generate this documentation very quickly. Therefore results are far from converged. Tighter thresholds and larger grids should be used for more realistic results.

For our discussion we will use the classic example of computing the LDA ground state of the silicon crystal. Performing such a calculation roughly proceeds in three steps.

In [1]:

```
using DFTK
using Plots
using Unitful
using UnitfulAtomic
# 1. Define lattice and atomic positions
a = 5.431u"angstrom" # Silicon lattice constant
lattice = a / 2 * [[0 1 1.]; # Silicon lattice vectors
[1 0 1.]; # specified column by column
[1 1 0.]];
```

In [2]:

```
# Load HGH pseudopotential for Silicon
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
# Specify type and positions of atoms
atoms = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]
# 2. Select model and basis
model = model_LDA(lattice, atoms, positions)
kgrid = [4, 4, 4] # k-point grid (Regular Monkhorst-Pack grid)
Ecut = 7 # kinetic energy cutoff
# Ecut = 190.5u"eV" # Could also use eV or other energy-compatible units
basis = PlaneWaveBasis(model; Ecut, kgrid)
# Note the implicit passing of keyword arguments here:
# this is equivalent to PlaneWaveBasis(model; Ecut=Ecut, kgrid=kgrid)
# 3. Run the SCF procedure to obtain the ground state
scfres = self_consistent_field(basis, tol=1e-5);
```

In [3]:

```
scfres.energies
```

Out[3]:

Energy breakdown (in Ha): Kinetic 3.1020977 AtomicLocal -2.1987839 AtomicNonlocal 1.7296072 Ewald -8.3979253 PspCorrection -0.2946254 Hartree 0.5530398 Xc -2.3986215 total -7.905211531382

Eigenvalues:

In [4]:

```
hcat(scfres.eigenvalues...)
```

Out[4]:

7×8 Matrix{Float64}: -0.176942 -0.14744 -0.0911692 … -0.101219 -0.023977 -0.018408 0.261073 0.116915 0.00482507 0.0611643 -0.023977 -0.018408 0.261073 0.23299 0.216733 0.121636 0.155532 0.117747 0.261073 0.23299 0.216733 0.212134 0.155532 0.117747 0.354532 0.335109 0.317102 0.350436 0.285692 0.417258 0.354532 0.389829 0.384601 … 0.436925 0.285692 0.417516 0.354533 0.389829 0.384601 0.44929 0.627662 0.443806

`eigenvalues`

is an array (indexed by k-points) of arrays (indexed by
eigenvalue number). The "splatting" operation `...`

calls `hcat`

with all the inner arrays as arguments, which collects them into a
matrix.

The resulting matrix is 7 (number of computed eigenvalues) by 8
(number of irreducible k-points). There are 7 eigenvalues per
k-point because there are 4 occupied states in the system (4 valence
electrons per silicon atom, two atoms per unit cell, and paired
spins), and the eigensolver gives itself some breathing room by
computing some extra states (see the `bands`

argument to
`self_consistent_field`

as well as the `AdaptiveBands`

documentation).
There are only 8 k-points (instead of 4x4x4) because symmetry has been used to reduce the
amount of computations to just the irreducible k-points (see
Crystal symmetries
for details).

We can check the occupations ...

In [5]:

```
hcat(scfres.occupation...)
```

Out[5]:

7×8 Matrix{Float64}: 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

In [6]:

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

Out[6]:

We can also perform various postprocessing steps: for instance compute a band structure

In [7]:

```
plot_bandstructure(scfres; kline_density=10)
```

Computing bands along kpath: Γ -> X -> U and K -> Γ -> L -> W -> X Diagonalising Hamiltonian kblocks: 100%|████████████████| Time: 0:00:01

Out[7]:

or get the cartesian forces (in Hartree / Bohr)

In [8]:

```
compute_forces_cart(scfres)
```

Out[8]:

2-element Vector{StaticArraysCore.SVector{3, Float64}}: [-3.5336251758888616e-16, -1.7807673505274671e-16, -4.694043518787138e-16] [-4.890859093161763e-16, 2.5304492174484446e-17, -1.1342162044248493e-16]

As expected, they are numerically zero in this highly symmetric configuration.