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-61.0e-6Density-based self-consistent field
scfres_scf = self_consistent_field(basis; tol);n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.397823282473 -0.90 5.0 27.1ms
2 -8.400246991447 -2.62 -1.73 1.0 30.6ms
3 -8.400398717033 -3.82 -2.90 1.5 20.4ms
4 -8.400427783043 -4.54 -2.90 3.0 25.3ms
5 -8.400427938664 -6.81 -3.01 1.0 19.6ms
6 -8.400428149377 -6.68 -4.76 1.0 26.1ms
7 -8.400428155842 -8.19 -4.41 3.0 25.9ms
8 -8.400428156259 -9.38 -5.27 1.0 19.7ms
9 -8.400428156275 -10.80 -6.28 1.0 19.6msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397780067398 -0.90 5.0 27.4ms
2 -8.400384598445 -2.58 -1.79 0.80 2.2 25.9ms
3 -8.400423032705 -4.42 -3.02 0.80 1.0 22.4ms
4 -8.400428112098 -5.29 -3.41 0.80 2.8 29.3ms
5 -8.400428154513 -7.37 -4.43 0.80 1.5 18.5ms
6 -8.400428156194 -8.77 -5.57 0.80 2.0 27.4ms
7 -8.400428156276 -10.08 -6.10 0.80 2.5 21.7msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.170147278802 -0.99 52.5ms
2 -1.408260481398 0.41 -0.65 29.7ms
3 -3.859164187223 0.39 -0.34 46.4ms
4 -4.726920549268 -0.06 -0.37 40.4ms
5 -6.965936136752 0.35 -0.44 45.6ms
6 -7.785622452166 -0.09 -0.85 39.8ms
7 -8.131002673265 -0.46 -1.23 29.2ms
8 -8.244219803286 -0.95 -1.51 29.5ms
9 -8.330581993500 -1.06 -1.70 35.3ms
10 -8.364604586949 -1.47 -2.14 29.6ms
11 -8.382858660496 -1.74 -2.15 29.5ms
12 -8.388365750008 -2.26 -2.43 29.3ms
13 -8.392875411500 -2.35 -3.07 35.0ms
14 -8.396460177440 -2.45 -2.91 29.6ms
15 -8.398684564435 -2.65 -2.91 29.8ms
16 -8.399341940360 -3.18 -3.32 30.4ms
17 -8.399942065139 -3.22 -3.33 29.8ms
18 -8.400164818743 -3.65 -3.34 36.1ms
19 -8.400323687737 -3.80 -3.66 29.8ms
20 -8.400385292745 -4.21 -3.73 29.7ms
21 -8.400411933313 -4.57 -3.98 29.8ms
22 -8.400420084299 -5.09 -3.88 35.1ms
23 -8.400424108507 -5.40 -4.17 30.0ms
24 -8.400425568807 -5.84 -4.23 32.5ms
25 -8.400427188636 -5.79 -4.32 30.2ms
26 -8.400427867257 -6.17 -4.74 36.0ms
27 -8.400427998822 -6.88 -5.06 29.9ms
28 -8.400428073397 -7.13 -5.00 29.7ms
29 -8.400428110685 -7.43 -6.19 30.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 -8.397791071625 -0.90 5.2 27.8msRemove 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.400427985670 -1.79 558ms
2 -8.400428156277 -6.77 -4.03 403ms
3 -8.400428156277 -14.45 -7.83 115msComparison 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| = 1.6082831575002699e-6
|ρ_newton - ρ_scfv| = 1.4561100317111284e-7
|ρ_newton - ρ_dm| = 1.2709707668745447e-5