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 reasonable base line value for
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
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:
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
Ecuts are likely needed.
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 $k$-points employed in each dimension of the Brillouin zone …
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
… and plot the obtained convergence:
using Plots plot(result.nkpts, result.errors, dpi=300, lw=3, m=:o, yaxis=:log, xlabel="k-grid", ylabel="energy absolute error")
We continue to do the convergence in Ecut using the suggested $k$-point grid.
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
… and plot it:
plot(result.Ecuts, result.errors, dpi=300, lw=3, m=:o, yaxis=:log, xlabel="Ecut", ylabel="energy absolute error")
Repeating the above exercise for more realistic settings, namely …
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