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.397803491316 -0.90 5.0 25.4ms
2 -8.400247842464 -2.61 -1.74 1.0 25.0ms
3 -8.400400158381 -3.82 -2.93 1.2 19.5ms
4 -8.400427820361 -4.56 -2.93 3.5 24.3ms
5 -8.400427843918 -7.63 -2.93 1.0 18.9ms
6 -8.400428143517 -6.52 -4.85 1.0 23.6ms
7 -8.400428151659 -8.09 -4.35 3.2 25.0ms
8 -8.400428152188 -9.28 -5.32 1.0 19.1ms
9 -8.400428152206 -10.73 -6.24 1.2 19.6ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397783895886 -0.90 4.8 25.6ms
2 -8.400375079640 -2.59 -1.77 0.80 2.0 18.7ms
3 -8.400422579348 -4.32 -2.98 0.80 1.0 16.2ms
4 -8.400428097197 -5.26 -3.47 0.80 2.2 19.6ms
5 -8.400428146339 -7.31 -5.07 0.80 1.2 21.6ms
6 -8.400428152201 -8.23 -5.11 0.80 3.5 23.0ms
7 -8.400428152208 -11.10 -6.09 0.80 1.0 16.6ms
Direct minimization
scfres_dm = direct_minimization(basis; tol);
┌ Warning: x_tol is deprecated. Use x_abstol or x_reltol instead. The provided value (-1) will be used as x_abstol.
└ @ Optim ~/.julia/packages/Optim/7krni/src/types.jl:110
┌ Warning: f_tol is deprecated. Use f_abstol or f_reltol instead. The provided value (-1) will be used as f_reltol.
└ @ Optim ~/.julia/packages/Optim/7krni/src/types.jl:120
n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.770751616314 -1.08 57.8ms
2 -1.974847677590 0.44 -0.68 37.3ms
3 -4.549990948652 0.41 -0.41 45.9ms
4 -6.296444665171 0.24 -0.52 43.8ms
5 -7.678488592807 0.14 -0.82 47.6ms
6 -8.120425408081 -0.35 -1.23 32.4ms
7 -8.301331615857 -0.74 -1.64 32.3ms
8 -8.348200075642 -1.33 -1.99 36.7ms
9 -8.379034791541 -1.51 -2.28 32.8ms
10 -8.387102024781 -2.09 -2.30 32.8ms
11 -8.394604125935 -2.12 -2.60 32.6ms
12 -8.396683147532 -2.68 -3.02 36.0ms
13 -8.398741335371 -2.69 -3.25 32.4ms
14 -8.399477915917 -3.13 -3.11 32.3ms
15 -8.400015511962 -3.27 -3.33 32.3ms
16 -8.400231743397 -3.67 -3.72 35.6ms
17 -8.400352002041 -3.92 -3.86 36.4ms
18 -8.400398112254 -4.34 -4.24 33.1ms
19 -8.400413354026 -4.82 -4.46 32.4ms
20 -8.400421459434 -5.09 -4.35 35.8ms
21 -8.400425237909 -5.42 -4.87 32.4ms
22 -8.400426784577 -5.81 -4.78 32.4ms
23 -8.400427509456 -6.14 -5.03 32.7ms
24 -8.400427915859 -6.39 -5.37 36.0ms
25 -8.400428042437 -6.90 -5.06 32.5ms
26 -8.400428113421 -7.15 -5.25 32.3ms
27 -8.400428125904 -7.90 -5.91 32.4ms
28 -8.400428144267 -7.74 -5.66 35.7ms
29 -8.400428147806 -8.45 -6.07 32.2ms
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.397814918866 -0.90 5.2 26.0ms
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.400427982521 -1.78 601ms
2 -8.400428152209 -6.77 -4.03 380ms
3 -8.400428152209 -14.45 -7.84 102ms
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| = 2.241678877515628e-6
|ρ_newton - ρ_scfv| = 3.821740981839904e-7
|ρ_newton - ρ_dm| = 2.1309095011256974e-6