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
system = attach_psp(bulk(:Si); Si="hgh/lda/Si-q4")
model = model_DFT(system; functionals=PBEsol())
basis = PlaneWaveBasis(model; Ecut=5, kgrid=[3, 3, 3])
# Convergence we desire in the density
tol = 1e-61.0e-6Density-based self-consistent field
scfres_scf = self_consistent_field(basis; tol);n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -7.844826470045 -0.71 4.5 35.1ms
2 -7.848885785438 -2.39 -1.53 1.0 24.8ms
3 -7.849099650509 -3.67 -2.49 1.2 25.5ms
4 -7.849134992103 -4.45 -2.78 2.5 30.3ms
5 -7.849135710717 -6.14 -2.96 1.0 32.9ms
6 -7.849136041012 -6.48 -4.02 1.0 25.1ms
7 -7.849136054049 -7.88 -4.83 2.0 28.8ms
8 -7.849136054369 -9.49 -5.15 2.0 28.7ms
9 -7.849136054390 -10.68 -5.67 1.2 32.0ms
10 -7.849136054392 -11.68 -6.25 1.0 26.1msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.844893190179 -0.71 4.8 35.7ms
2 -7.849060851015 -2.38 -1.66 0.80 2.0 32.0ms
3 -7.849129354398 -4.16 -2.75 0.80 1.0 22.1ms
4 -7.849135909440 -5.18 -3.31 0.80 2.0 25.6ms
5 -7.849136044987 -6.87 -4.36 0.80 1.2 22.8ms
6 -7.849136054314 -8.03 -5.12 0.80 2.0 25.9ms
7 -7.849136054390 -10.12 -5.61 0.80 1.5 23.8ms
8 -7.849136054392 -11.61 -6.89 0.80 1.2 29.0msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.410593455830 -1.01 75.8ms
2 -1.347811410028 0.44 -0.68 43.6ms
3 -3.658104485535 0.36 -0.45 65.6ms
4 -4.940100256223 0.11 -0.59 66.0ms
5 -6.592509357714 0.22 -0.75 63.2ms
6 -7.250146423404 -0.18 -1.23 58.5ms
7 -7.560667326600 -0.51 -1.44 43.5ms
8 -7.690448332350 -0.89 -1.86 43.0ms
9 -7.740532613015 -1.30 -1.98 48.8ms
10 -7.773722994413 -1.48 -2.25 43.4ms
11 -7.802668970872 -1.54 -2.02 43.4ms
12 -7.822044618640 -1.71 -2.25 43.3ms
13 -7.836134829832 -1.85 -2.30 48.0ms
14 -7.840591976793 -2.35 -2.57 43.3ms
15 -7.841273153663 -3.17 -3.17 43.6ms
16 -7.841671136660 -3.40 -2.88 43.3ms
17 -7.845233860945 -2.45 -2.91 48.6ms
18 -7.847749961520 -2.60 -2.99 43.4ms
19 -7.848546148032 -3.10 -3.18 43.2ms
20 -7.848938412796 -3.41 -3.45 43.5ms
21 -7.849077510415 -3.86 -3.82 48.3ms
22 -7.849116793517 -4.41 -3.88 43.7ms
23 -7.849130691152 -4.86 -4.55 43.5ms
24 -7.849134422032 -5.43 -4.71 43.3ms
25 -7.849135464470 -5.98 -4.92 49.1ms
26 -7.849135855273 -6.41 -5.02 46.0ms
27 -7.849135987604 -6.88 -5.65 48.7ms
28 -7.849136034937 -7.32 -5.44 44.1ms
29 -7.849136047497 -7.90 -6.39 48.4msNewton 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 -7.844881523071 -0.71 4.8 35.4msRemove 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 -7.849135371178 -1.65 572ms
2 -7.849136054392 -6.17 -3.72 430ms
3 -7.849136054393 -13.21 -7.21 168msComparison 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| = 4.068993874732651e-7
|ρ_newton - ρ_scfv| = 1.716931260669828e-7
|ρ_newton - ρ_dm| = 1.21327835625516e-6