# Tutorial

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 density-functional theory chapter for some introductory material on the topic.

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.

```
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.]]
```

3×3 Matrix{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(Å,), 𝐋, nothing}}}: 0.0 Å 2.7155 Å 2.7155 Å 2.7155 Å 0.0 Å 2.7155 Å 2.7155 Å 2.7155 Å 0.0 Å

By default, all numbers passed as arguments are assumed to be in atomic units. Quantities such as temperature, energy cutoffs, lattice vectors, and the k-point grid spacing can optionally be annotated with Unitful units, which are automatically converted to the atomic units used internally. For more details, see the Unitful package documentation and the UnitfulAtomic.jl package.

```
# Load HGH pseudopotential for Silicon
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
# Specify type and positions of atoms
atoms = [Si => [ones(3)/8, -ones(3)/8]]
# 2. Select model and basis
model = model_LDA(lattice, atoms)
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=kgrid)
# 3. Run the SCF procedure to obtain the ground state
scfres = self_consistent_field(basis, tol=1e-8);
```

n Energy Eₙ-Eₙ₋₁ ρout-ρin Diag --- --------------- --------- -------- ---- 1 -7.905180225069 NaN 1.95e-01 4.0 2 -7.909819131969 -4.64e-03 2.99e-02 1.0 3 -7.910051971089 -2.33e-04 3.02e-03 3.1 4 -7.910052831210 -8.60e-07 4.68e-04 2.4 5 -7.910052852099 -2.09e-08 4.21e-05 2.0 6 -7.910052854032 -1.93e-09 4.96e-06 3.1

That's it! Now you can get various quantities from the result of the SCF. For instance, the different components of the energy:

`scfres.energies`

Energy breakdown: Kinetic 3.0795156 AtomicLocal -2.1806378 AtomicNonlocal 1.7339394 Ewald -8.3979253 PspCorrection -0.2946254 Hartree 0.5417571 Xc -2.3920764 total -7.910052854032

Eigenvalues:

`hcat(scfres.eigenvalues...)`

7×10 Matrix{Float64}: -0.170182 -0.131801 -0.0883283 … -0.0562607 -0.11494 -0.0700168 0.201343 0.0909036 0.0122924 0.0111144 0.0420612 0.0176494 0.249295 0.174773 0.176137 0.132964 0.220114 0.11233 0.249295 0.231428 0.20237 0.161043 0.220114 0.190456 0.350986 0.360027 0.340134 0.291807 0.320727 0.327376 0.36997 0.395897 0.389483 … 0.331814 0.38819 0.460234 0.36997 0.401676 0.412474 0.565538 0.38819 0.462717

`eigenvalues`

is an array (indexed by kpoints) 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 kpoints). There are 7 eigenvalues per kpoint 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 `n_ep_extra`

argument to `self_consistent_field`

).

We can check the occupations ...

`hcat(scfres.occupation...)`

7×10 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 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 0.0 0.0 0.0 0.0 0.0 0.0

... and density, where we use that the density objects in DFTK are indexed as ρ[iσ, ix, iy, iz], i.e. first in the spin component and then in the 3-dimensional real-space grid.

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

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

`plot_bandstructure(scfres, kline_density=5)`

or get the cartesian forces (in Hartree / Bohr)

`compute_forces_cart(scfres)[1] # Select silicon forces`

2-element Vector{StaticArrays.SVector{3, Float64}}: [0.0004485948737865541, 0.00044860053216400395, 0.000448598626599144] [-0.00044860053216473395, -0.00044859487378597014, -0.00044859862659913735]

The `[1]`

extracts the forces for the first kind of atoms, i.e. `Si`

(silicon) in the setup of the `atoms`

list of step 1 above. As expected, they are almost zero in this highly symmetric configuration.

## Where to go from here

Take a look at the example index to continue exploring DFTK.