Comparison of DFT solvers
We compare four different approaches for solving the DFT minimisation problem, namely a density-based SCF, a potential-based SCF, direct minimisation and Newton.
First we setup our problem
using AtomsBuilder
using DFTK
using LinearAlgebra
using PseudoPotentialData
pseudopotentials = PseudoFamily("dojo.nc.sr.pbesol.v0_4_1.standard.upf")
model = model_DFT(bulk(:Si); functionals=PBEsol(), pseudopotentials)
basis = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3])
# Convergence we desire in the density
tol = 1e-6
1.0e-6
Density-based self-consistent field
scfres_scf = self_consistent_field(basis; tol);
n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.397812342863 -0.90 5.2 27.9ms
2 -8.400245356082 -2.61 -1.74 1.0 19.3ms
3 -8.400405365200 -3.80 -2.94 1.5 20.1ms
4 -8.400427856883 -4.65 -2.98 3.0 34.3ms
5 -8.400427955866 -7.00 -3.05 1.0 19.0ms
6 -8.400428149838 -6.71 -4.66 1.0 18.9ms
7 -8.400428155886 -8.22 -4.45 2.8 24.2ms
8 -8.400428156252 -9.44 -5.14 1.0 26.2ms
9 -8.400428156275 -10.64 -6.45 1.0 19.7ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397813621133 -0.90 5.2 27.5ms
2 -8.400386200899 -2.59 -1.79 0.80 2.0 19.4ms
3 -8.400422947101 -4.43 -3.02 0.80 1.0 46.2ms
4 -8.400428101604 -5.29 -3.47 0.80 2.5 21.2ms
5 -8.400428153761 -7.28 -4.69 0.80 1.5 17.3ms
6 -8.400428156257 -8.60 -5.84 0.80 2.5 20.3ms
7 -8.400428156277 -10.71 -6.29 0.80 3.0 22.1ms
Direct minimization
scfres_dm = direct_minimization(basis; tol);
n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.919451469661 -1.03 58.3ms
2 -1.677858818021 0.41 -0.65 28.8ms
3 -3.747092432909 0.32 -0.37 38.8ms
4 -4.334474487935 -0.23 -0.44 44.0ms
5 -6.359246507447 0.31 -0.47 39.2ms
6 -7.130156016926 -0.11 -0.70 39.4ms
7 -7.854761104643 -0.14 -0.93 44.7ms
8 -8.003162964487 -0.83 -1.48 28.7ms
9 -8.143751525957 -0.85 -2.06 28.9ms
10 -8.200082744663 -1.25 -1.92 33.7ms
11 -8.248960387334 -1.31 -2.14 29.0ms
12 -8.272892444389 -1.62 -2.26 28.8ms
13 -8.292469693735 -1.71 -2.21 28.7ms
14 -8.315498130415 -1.64 -2.17 28.7ms
15 -8.348688120333 -1.48 -2.08 34.9ms
16 -8.383177369085 -1.46 -2.00 28.7ms
17 -8.395098669857 -1.92 -2.29 28.8ms
18 -8.397400039851 -2.64 -2.93 28.7ms
19 -8.399426976794 -2.69 -3.15 33.7ms
20 -8.399942836387 -3.29 -3.33 28.8ms
21 -8.400277156899 -3.48 -3.28 28.6ms
22 -8.400328383263 -4.29 -3.64 28.7ms
23 -8.400406322208 -4.11 -3.60 33.8ms
24 -8.400417394703 -4.96 -4.17 28.9ms
25 -8.400425423165 -5.10 -4.41 29.0ms
26 -8.400426101124 -6.17 -4.80 28.9ms
27 -8.400427297183 -5.92 -4.60 33.7ms
28 -8.400427786686 -6.31 -5.04 29.4ms
29 -8.400428004881 -6.66 -4.76 29.4ms
30 -8.400428077999 -7.14 -5.31 29.3ms
31 -8.400428121595 -7.36 -5.16 34.9ms
32 -8.400428143526 -7.66 -5.58 29.1ms
33 -8.400428150935 -8.13 -5.60 28.9ms
34 -8.400428154231 -8.48 -6.27 28.7ms
Newton algorithm
Start not too far from the solution to ensure convergence: We run first a very crude SCF to get close and then switch to Newton.
scfres_start = self_consistent_field(basis; tol=0.5);
n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.397842246617 -0.90 5.2 26.9ms
Remove the virtual orbitals (which Newton cannot treat yet)
ψ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ
scfres_newton = newton(basis, ψ; tol);
n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 -8.400427983409 -1.78 546ms
2 -8.400428156277 -6.76 -4.02 375ms
3 -8.400428156277 -14.75 -7.83 95.4ms
Comparison of results
println("|ρ_newton - ρ_scf| = ", norm(scfres_newton.ρ - scfres_scf.ρ))
println("|ρ_newton - ρ_scfv| = ", norm(scfres_newton.ρ - scfres_scfv.ρ))
println("|ρ_newton - ρ_dm| = ", norm(scfres_newton.ρ - scfres_dm.ρ))
|ρ_newton - ρ_scf| = 7.127859140011735e-7
|ρ_newton - ρ_scfv| = 4.7364703607356785e-8
|ρ_newton - ρ_dm| = 3.211491780471981e-6