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.

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.

using DFTK
using Plots

# 1. Define lattice and atomic positions
a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]

# 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 in Hartree
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.905285006889         NaN   1.95e-01    4.2
  2   -7.910014743956   -4.73e-03   2.97e-02    2.0
  3   -7.910099133021   -8.44e-05   3.23e-03    1.0
  4   -7.910108969387   -9.84e-06   6.25e-04    1.9
  5   -7.910109400868   -4.31e-07   8.79e-05    1.6
  6   -7.910109413823   -1.30e-08   1.46e-05    1.7
  7   -7.910109414107   -2.83e-10   1.32e-06    2.0

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.0807728 
    AtomicLocal         -2.1784490
    AtomicNonlocal      1.7341054 
    Ewald               -8.4004648
    PspCorrection       -0.2948928
    Hartree             0.5413908 
    Xc                  -2.3925718

    total               -7.910109414107

Eigenvalues:

hcat(scfres.eigenvalues...)
7×10 Array{Float64,2}:
 -0.170033  -0.13162    -0.0881023  …  -0.0560208  -0.114744  -0.0697773
  0.201708   0.0912095   0.0125263      0.0114072   0.042318   0.0179176
  0.249664   0.175086    0.176465       0.133226    0.220466   0.112584 
  0.249664   0.231788    0.202704       0.16134     0.220466   0.190778 
  0.351459   0.360507    0.340562       0.292133    0.321216   0.32774  
  0.370379   0.396308    0.389908   …   0.332169    0.388572   0.4608   
  0.370379   0.402093    0.41287        0.566123    0.388572   0.463173 

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 Array{Float64,2}:
 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:

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.ρ.real[:, 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, unit=:eV)

or forces

forces(scfres)[1]  # Select silicon forces
2-element Array{StaticArrays.SArray{Tuple{3},Float64,1,3},1}:
 [0.004623068410215786, 0.004623023330046665, 0.004623003797512429]   
 [-0.004623023330047861, -0.004623068410216881, -0.004623003797512874]

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.