Tutorial

This document provides an overview of the structure of the code and how to access basic information about calculations. If you have not installed DFTK yet, please refer to the Installation instructions.

Basic familiarity with the concepts of plane-wave density functional theory is assumed in this tutorial. Take a look at the Introductory resources for some broader introductory material. There is also a Mathematical Tutorial.

Setting up a first calculation

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
using PseudoPotentialData

# 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.]];

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.

We use a pseudodojo pseudopotential (see PseudoPotentialData for more details on PseudoFamily):

pd_lda_family = PseudoFamily("dojo.nc.sr.lda.v0_4_1.standard.upf")
Si = ElementPsp(:Si, pd_lda_family)

# Specify type and positions of atoms
atoms     = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]

2-element Vector{Vector{Float64}}:
 [0.125, 0.125, 0.125]
 [-0.125, -0.125, -0.125]

Note that DFTK supports a few other ways to supply atomistic structures, see for example the sections on AtomsBase integration and Input and output formats for details.

# 2. Select model and basis
model = model_DFT(lattice, atoms, positions; functionals=LDA())

kgrid = KgridSpacing(0.3 / u"bohr")  # Regular k-point grid (Monkhorst-Pack grid)
#                                      with spacing 0.3/bohr between k-points
# kgrid = [4, 4, 4]                    Alternative: Number of k-points per dimension
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);

n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime
---   ---------------   ---------   ---------   ----   ------
  1   -8.505317530262                   -0.93    5.1   45.8ms
  2   -8.508287871130       -2.53       -1.79    1.0   23.2ms
  3   -8.508467219248       -3.75       -2.93    1.8   26.8ms
  4   -8.508483073209       -4.80       -3.25    3.0   34.8ms
  5   -8.508483196435       -6.91       -3.63    1.2   60.0ms
  6   -8.508483225650       -7.53       -4.88    1.2   25.2ms
  7   -8.508483227616       -8.71       -5.23    3.1   35.8ms

Inspecting the results

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 (in Ha):
    Kinetic             3.0841945 
    AtomicLocal         -2.3554808
    AtomicNonlocal      1.3116670 
    Ewald               -8.3979253
    PspCorrection       0.3948680 
    Hartree             0.5559191 
    Xc                  -3.1017258

    total               -8.508483227616

Eigenvalues:

stack(scfres.eigenvalues)

7×8 Matrix{Float64}:
 -0.264798  -0.235352   -0.179074   …  -0.189309   -0.112072   -0.106306
  0.174324   0.0295318  -0.0825371     -0.0265553  -0.112072   -0.106306
  0.174324   0.146148    0.129858       0.0347393   0.0686555   0.0309537
  0.174324   0.146148    0.129858       0.12521     0.0686555   0.0309537
  0.267151   0.247827    0.229557       0.263041    0.198144    0.329876
  0.267151   0.302405    0.297036   …   0.349551    0.198144    0.330058
  0.267152   0.302405    0.297036       0.362159    0.541895    0.356808

eigenvalues is an array (indexed by k-points) of arrays (indexed by eigenvalue number).

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

stack(scfres.occupation)

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

... and density, where we use that the density objects in DFTK are indexed as ρ[ix, iy, iz, iσ], i.e. first in the 3-dimensional real-space grid and then in 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)

We can also perform various postprocessing steps: We can get the Cartesian forces (in Hartree / Bohr):

compute_forces_cart(scfres)

2-element Vector{StaticArraysCore.SVector{3, Float64}}:
 [-2.1182491126014674e-14, -2.0941200786420585e-14, -2.0847997753675167e-14]
 [2.1421405818800503e-14, 2.0964802108625105e-14, 2.1109219370338126e-14]

As expected, they are numerically zero in this highly symmetric configuration. We could also compute a band structure: we compute the bands along a k-point line (determined automatically) and plot the result:

bands1 = compute_bands(scfres; kline_density=10)
plot_bandstructure(bands1)

Next we want to plot a density of states. We can use the scfres directly, but the resulting DOS is quite sharp:

plot_dos(scfres; temperature=1e-3, smearing=Smearing.FermiDirac())

To get a better result we first increase the kgrid to get a better discretisation of the Brillouin zone, then re-do the plot:

bands2 = compute_bands(scfres, MonkhorstPack(6, 6, 6))
plot_dos(bands2; temperature=1e-3, smearing=Smearing.FermiDirac())

Note, that some other codes would refer to the functionality we provide with compute_bands as "performing a NSCF calculation".