Performing a convergence study

This example shows how to perform a convergence study to find an appropriate discretisation parameters for the Brillouin zone (kgrid) and kinetic energy cutoff (Ecut), such that the simulation results are converged to a desired accuracy tolerance.

Such a convergence study is generally performed by starting with a reasonalbe base line value for kgrid and Ecut and then increasing these parameters (i.e. using finer discretisations) until a desired property (such as the energy) changes less than the tolerance.

This procedure must be performed for each discretisation parameter. Beyond the Ecut and the kgrid also convergence in the smearing temperature or other numerical parameters should be checked. For simplicity we will neglect this aspect in this example and concentrate on Ecut and kgrid. Moreover we will restrict ourselves to using the same number of $k$-points in each dimension of the Brillouin zone.

As the objective of this study we consider bulk platinum. For running the SCF conveniently we define a function:

In [1]:
using DFTK
using LinearAlgebra
using Statistics

function run_scf(; a=5.0, Ecut, nkpt, tol)
    atoms    = [ElementPsp(:Pt, psp = load_psp("hgh/lda/Pt-q10"))]
    position = [zeros(3)]
    lattice  = a * Matrix(I, 3, 3)

    model  = model_LDA(lattice, atoms, position; temperature=1e-2)
    basis  = PlaneWaveBasis(model; Ecut, kgrid=(nkpt, nkpt, nkpt))
    println("nkpt = $nkpt Ecut = $Ecut")
    self_consistent_field(basis; tol)
end;

Moreover we define some parameters. To make the calculations run fast for the automatic generation of this documentation we target only a convergence to 1e-2. In practice smaller tolerances (and thus larger upper bounds for nkpts and Ecuts are likely needed.

In [2]:
tol   = 1e-2      # Tolerance to which we target to converge
nkpts = 1:7       # K-point range checked for convergence
Ecuts = 10:2:24;  # Energy cutoff range checked for convergence

As the first step we converge in the number of kpoints employed in each dimension of the Brillouin zone.

In [3]:
function converge_kgrid(nkpts; Ecut, tol)
    energies = [run_scf(; nkpt, tol=tol/10, Ecut).energies.total for nkpt in nkpts]
    errors = abs.(energies[1:end-1] .- energies[end])
    iconv = findfirst(errors .< tol)
    (; nkpts=nkpts[1:end-1], errors, nkpt_conv=nkpts[iconv])
end
result = converge_kgrid(nkpts; Ecut=mean(Ecuts), tol)
nkpt_conv = result.nkpt_conv
nkpt = 1 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -26.49741906817                   -0.22    8.0
  2   -26.58836237620       -1.04       -0.62    1.0
  3   -26.61280068980       -1.61       -1.40    2.0
  4   -26.61326646593       -3.33       -2.04    2.0
nkpt = 2 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.79178984950                   -0.09    5.2
  2   -26.23084159405       -0.36       -0.69    2.0
  3   -26.23811363561       -2.14       -1.32    2.0
  4   -26.23848448641       -3.43       -2.31    1.8
nkpt = 3 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.78310889219                   -0.09    4.8
  2   -26.23614628831       -0.34       -0.77    2.0
  3   -26.25073967638       -1.84       -1.64    2.5
  4   -26.25103354166       -3.53       -2.13    2.0
nkpt = 4 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.90974905416                   -0.11    4.7
  2   -26.29041841880       -0.42       -0.75    2.1
  3   -26.30835878292       -1.75       -1.74    2.3
  4   -26.30842667771       -4.17       -2.69    2.1
nkpt = 5 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.90349837347                   -0.11    3.7
  2   -26.26204072277       -0.45       -0.69    2.0
  3   -26.28543461720       -1.63       -1.63    2.1
  4   -26.28570510050       -3.57       -2.21    2.2
nkpt = 6 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.87456844697                   -0.10    4.8
  2   -26.26994940171       -0.40       -0.74    2.1
  3   -26.28809805537       -1.74       -1.71    2.2
  4   -26.28818805094       -4.05       -2.57    2.0
nkpt = 7 Ecut = 17.0
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.89717243209                   -0.11    3.4
  2   -26.27487981066       -0.42       -0.73    2.0
  3   -26.29415697635       -1.71       -1.75    2.2
  4   -26.29420805436       -4.29       -2.65    2.1
Out[3]:
5

... and plot the obtained convergence:

In [4]:
using Plots
plot(result.nkpts, result.errors, dpi=300, lw=3, m=:o, yaxis=:log,
     xlabel="k-grid", ylabel="energy absolute error")
Out[4]:

We continue to do the convergence in Ecut using the suggested k-point grid.

In [5]:
function converge_Ecut(Ecuts; nkpt, tol)
    energies = [run_scf(; nkpt, tol=tol/100, Ecut).energies.total for Ecut in Ecuts]
    errors = abs.(energies[1:end-1] .- energies[end])
    iconv = findfirst(errors .< tol)
    (; Ecuts=Ecuts[1:end-1], errors, Ecut_conv=Ecuts[iconv])
end
result = converge_Ecut(Ecuts; nkpt=nkpt_conv, tol)
Ecut_conv = result.Ecut_conv
nkpt = 5 Ecut = 10
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.57958080618                   -0.16    3.4
  2   -25.77480454892       -0.71       -0.74    1.9
  3   -25.78626211028       -1.94       -1.87    2.0
  4   -25.78631880670       -4.25       -2.81    2.0
nkpt = 5 Ecut = 12
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.78766408098                   -0.12    3.5
  2   -26.07365187343       -0.54       -0.70    2.0
  3   -26.09343910867       -1.70       -1.70    2.0
  4   -26.09374098797       -3.52       -2.29    2.2
  5   -26.09375538806       -4.84       -2.82    2.1
nkpt = 5 Ecut = 14
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.86882912494                   -0.11    3.6
  2   -26.20473737205       -0.47       -0.69    2.0
  3   -26.22672180482       -1.66       -1.66    2.1
  4   -26.22700483042       -3.55       -2.24    2.2
  5   -26.22702786299       -4.64       -2.82    2.4
nkpt = 5 Ecut = 16
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.89833708665                   -0.11    3.7
  2   -26.25250685839       -0.45       -0.69    2.0
  3   -26.27560211915       -1.64       -1.64    2.1
  4   -26.27586992265       -3.57       -2.22    2.2
  5   -26.27589483703       -4.60       -2.87    2.5
nkpt = 5 Ecut = 18
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.90621377534                   -0.11    3.8
  2   -26.26805157266       -0.44       -0.70    2.0
  3   -26.29103419590       -1.64       -1.64    2.2
  4   -26.29128156854       -3.61       -2.22    2.2
  5   -26.29130575794       -4.62       -2.87    2.3
nkpt = 5 Ecut = 20
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.90867209840                   -0.11    3.9
  2   -26.27264237845       -0.44       -0.70    2.0
  3   -26.29535676157       -1.64       -1.65    2.2
  4   -26.29558543326       -3.64       -2.24    2.2
  5   -26.29560740647       -4.66       -2.87    2.5
nkpt = 5 Ecut = 22
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.90876810644                   -0.11    3.8
  2   -26.27344653728       -0.44       -0.70    2.0
  3   -26.29618671572       -1.64       -1.64    2.1
  4   -26.29641445640       -3.64       -2.23    2.3
  5   -26.29643720946       -4.64       -2.89    2.5
nkpt = 5 Ecut = 24
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -25.90854155019                   -0.11    3.9
  2   -26.27374209523       -0.44       -0.70    2.0
  3   -26.29627546620       -1.65       -1.65    2.2
  4   -26.29649181343       -3.66       -2.25    2.3
  5   -26.29651326593       -4.67       -2.91    2.5
Out[5]:
18

... and plot it:

In [6]:
plot(result.Ecuts, result.errors, dpi=300, lw=3, m=:o, yaxis=:log,
     xlabel="Ecut", ylabel="energy absolute error")
Out[6]:

A more realistic example.

Repeating the above exercise for more realistic settings, namely ...

In [7]:
tol   = 1e-4  # Tolerance to which we target to converge
nkpts = 1:20  # K-point range checked for convergence
Ecuts = 20:1:50;

...one obtains the following two plots for the convergence in kpoints and Ecut.