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.397768532432 -0.90 5.2 27.3ms
2 -8.400233699668 -2.61 -1.73 1.0 28.6ms
3 -8.400402689831 -3.77 -2.93 1.8 21.0ms
4 -8.400427787689 -4.60 -2.94 3.2 25.7ms
5 -8.400427965696 -6.75 -3.07 1.0 19.5ms
6 -8.400428149657 -6.74 -4.61 1.0 19.4ms
7 -8.400428155838 -8.21 -4.44 2.8 31.1ms
8 -8.400428156226 -9.41 -4.91 1.0 19.5ms
9 -8.400428156274 -10.31 -6.56 1.0 19.6msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397870696068 -0.90 5.2 34.1ms
2 -8.400383743001 -2.60 -1.79 0.80 2.0 19.0ms
3 -8.400422322189 -4.41 -3.04 0.80 1.0 15.7ms
4 -8.400428107597 -5.24 -3.41 0.80 2.2 19.8ms
5 -8.400428152830 -7.34 -4.60 0.80 1.0 15.8ms
6 -8.400428156243 -8.47 -5.66 0.80 2.5 20.3ms
7 -8.400428156277 -10.48 -6.14 0.80 3.0 22.1msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.365032036701 -1.05 59.5ms
2 -1.983377260075 0.52 -0.64 29.8ms
3 -4.678214008827 0.43 -0.41 40.2ms
4 -6.349259178802 0.22 -0.51 45.4ms
5 -7.600674646341 0.10 -0.77 40.0ms
6 -7.995413014429 -0.40 -1.42 29.2ms
7 -8.174763646439 -0.75 -1.57 35.0ms
8 -8.244307905716 -1.16 -1.77 29.3ms
9 -8.295511626795 -1.29 -2.00 29.8ms
10 -8.315677164875 -1.70 -2.13 29.5ms
11 -8.350855230388 -1.45 -1.95 29.9ms
12 -8.371231678507 -1.69 -2.17 35.7ms
13 -8.389037704450 -1.75 -2.63 29.5ms
14 -8.395719713110 -2.18 -2.63 29.7ms
15 -8.398381357096 -2.57 -2.81 30.1ms
16 -8.399355410478 -3.01 -3.30 35.9ms
17 -8.400087986417 -3.14 -3.58 29.7ms
18 -8.400247433165 -3.80 -3.29 29.5ms
19 -8.400375527903 -3.89 -3.47 29.6ms
20 -8.400403860803 -4.55 -4.00 35.4ms
21 -8.400420066216 -4.79 -4.16 29.9ms
22 -8.400424864249 -5.32 -4.45 29.9ms
23 -8.400426991483 -5.67 -4.76 29.9ms
24 -8.400427659799 -6.18 -4.52 29.7ms
25 -8.400427967024 -6.51 -5.34 36.2ms
26 -8.400428054337 -7.06 -5.09 29.8ms
27 -8.400428126033 -7.14 -5.67 30.0ms
28 -8.400428140764 -7.83 -5.40 29.7ms
29 -8.400428152030 -7.95 -5.75 35.9ms
30 -8.400428154235 -8.66 -5.69 30.2ms
31 -8.400428155617 -8.86 -6.44 30.1msNewton 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.397802581192 -0.90 5.2 27.6msRemove 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.400427978560 -1.79 582ms
2 -8.400428156277 -6.75 -4.03 368ms
3 -8.400428156277 -14.27 -7.81 131msComparison 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| = 3.8256921367762754e-7
|ρ_newton - ρ_scfv| = 3.005435864491837e-7
|ρ_newton - ρ_dm| = 2.2031824056874472e-6