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.844854079310 -0.70 4.5 32.9ms
2 -7.848889916730 -2.39 -1.53 1.0 22.7ms
3 -7.849103116539 -3.67 -2.49 1.5 24.0ms
4 -7.849135141387 -4.49 -2.80 2.2 27.2ms
5 -7.849135786450 -6.19 -3.02 1.0 22.5ms
6 -7.849136042943 -6.59 -4.04 1.0 22.5ms
7 -7.849136053996 -7.96 -4.71 1.8 25.5ms
8 -7.849136054365 -9.43 -5.23 1.5 25.0ms
9 -7.849136054391 -10.58 -5.87 1.5 24.6ms
10 -7.849136054392 -11.79 -6.98 1.8 25.2msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.844860932941 -0.71 4.5 32.4ms
2 -7.849061120600 -2.38 -1.65 0.80 2.2 23.9ms
3 -7.849130024556 -4.16 -2.75 0.80 1.0 19.2ms
4 -7.849135905358 -5.23 -3.32 0.80 1.8 22.7ms
5 -7.849136045735 -6.85 -4.37 0.80 1.8 21.7ms
6 -7.849136054307 -8.07 -5.03 0.80 2.2 24.3ms
7 -7.849136054388 -10.09 -5.57 0.80 1.2 20.5ms
8 -7.849136054392 -11.44 -6.58 0.80 1.5 21.2msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.505037091176 -0.97 68.4ms
2 -1.734460128162 0.51 -0.75 43.5ms
3 -3.728855616682 0.30 -0.46 86.5ms
4 -4.908276172778 0.07 -0.59 49.1ms
5 -6.586999893947 0.22 -0.69 49.1ms
6 -7.304325093490 -0.14 -1.09 49.0ms
7 -7.554989648191 -0.60 -1.38 36.3ms
8 -7.630549333544 -1.12 -1.73 35.9ms
9 -7.659649357154 -1.54 -2.10 36.1ms
10 -7.690865928891 -1.51 -1.91 36.0ms
11 -7.702959887071 -1.92 -1.81 48.3ms
12 -7.718953030361 -1.80 -2.15 48.5ms
13 -7.750408112618 -1.50 -2.45 48.4ms
14 -7.760431831604 -2.00 -2.29 35.8ms
15 -7.790606600313 -1.52 -2.66 36.0ms
16 -7.798112673770 -2.12 -2.69 36.0ms
17 -7.816870130746 -1.73 -2.94 35.9ms
18 -7.827746403272 -1.96 -2.70 35.6ms
19 -7.841074443501 -1.88 -2.41 36.7ms
20 -7.845903212507 -2.32 -2.70 40.9ms
21 -7.848111992392 -2.66 -2.89 42.6ms
22 -7.848968482542 -3.07 -3.18 42.5ms
23 -7.849084232257 -3.94 -3.90 42.5ms
24 -7.849117077970 -4.48 -4.21 42.3ms
25 -7.849130778122 -4.86 -4.33 41.9ms
26 -7.849134051770 -5.48 -4.74 36.1ms
27 -7.849135388577 -5.87 -4.90 36.1ms
28 -7.849135879871 -6.31 -5.44 36.3ms
29 -7.849135989216 -6.96 -5.34 36.6ms
30 -7.849136031512 -7.37 -5.57 36.0ms
31 -7.849136044956 -7.87 -5.98 35.9ms
32 -7.849136050579 -8.25 -5.86 35.8ms
33 -7.849136053240 -8.58 -6.07 36.0msNewton 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.844878866373 -0.71 4.8 33.1msRemove 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.849135378826 -1.65 510ms
2 -7.849136054392 -6.17 -3.72 383ms
3 -7.849136054393 -13.24 -7.22 152msComparison 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| = 2.084936328419631e-7
|ρ_newton - ρ_scfv| = 4.2278999501898224e-7
|ρ_newton - ρ_dm| = 2.156356993353712e-6