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.397805908113 -0.90 5.0 27.4ms
2 -8.400251294577 -2.61 -1.73 1.0 50.9ms
3 -8.400404203757 -3.82 -2.94 1.5 21.4ms
4 -8.400427778812 -4.63 -2.91 3.0 25.9ms
5 -8.400427993688 -6.67 -3.09 1.0 20.3ms
6 -8.400428150572 -6.80 -4.94 1.0 20.1ms
7 -8.400428155954 -8.27 -4.48 3.5 36.0ms
8 -8.400428156261 -9.51 -5.23 1.2 20.7ms
9 -8.400428156275 -10.85 -6.47 1.0 20.6msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397805578481 -0.90 5.2 35.0ms
2 -8.400385316731 -2.59 -1.78 0.80 2.0 19.1ms
3 -8.400422720763 -4.43 -3.01 0.80 1.0 16.0ms
4 -8.400428100384 -5.27 -3.43 0.80 2.5 20.8ms
5 -8.400428154099 -7.27 -4.60 0.80 1.5 17.6ms
6 -8.400428156243 -8.67 -5.99 0.80 2.2 20.5ms
7 -8.400428156277 -10.48 -5.96 0.80 3.5 31.7ms
8 -8.400428156277 -13.05 -7.40 0.80 1.0 16.7msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.775605526982 -1.04 53.6ms
2 -1.730247044336 0.40 -0.67 29.7ms
3 -3.783187989242 0.31 -0.41 47.9ms
4 -4.818446786771 0.02 -0.43 41.2ms
5 -6.681735688021 0.27 -0.58 41.0ms
6 -7.564880777112 -0.05 -0.79 46.5ms
7 -7.785080579520 -0.66 -1.22 29.8ms
8 -8.223099938303 -0.36 -1.51 29.8ms
9 -8.302718611457 -1.10 -1.77 29.9ms
10 -8.334198139995 -1.50 -1.94 36.9ms
11 -8.355535677473 -1.67 -2.51 30.0ms
12 -8.381194643363 -1.59 -2.42 30.1ms
13 -8.388700031717 -2.12 -2.73 29.8ms
14 -8.393094200855 -2.36 -2.91 36.3ms
15 -8.396950942706 -2.41 -2.84 30.0ms
16 -8.397938153637 -3.01 -3.19 40.8ms
17 -8.398358054314 -3.38 -3.17 46.6ms
18 -8.399063246555 -3.15 -3.09 30.2ms
19 -8.399821838791 -3.12 -3.31 30.6ms
20 -8.400191096201 -3.43 -3.63 30.5ms
21 -8.400329134992 -3.86 -3.62 29.8ms
22 -8.400382806771 -4.27 -3.83 35.8ms
23 -8.400406255921 -4.63 -4.08 29.7ms
24 -8.400421494718 -4.82 -4.51 29.4ms
25 -8.400423586431 -5.68 -4.17 29.5ms
26 -8.400426723156 -5.50 -4.40 36.6ms
27 -8.400427153156 -6.37 -4.65 30.1ms
28 -8.400427838909 -6.16 -4.68 30.0ms
29 -8.400427977121 -6.86 -5.08 29.9ms
30 -8.400428096030 -6.92 -5.36 36.2ms
31 -8.400428122624 -7.58 -5.70 30.2ms
32 -8.400428146072 -7.63 -5.79 29.8ms
33 -8.400428148635 -8.59 -6.08 29.6msNewton 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.397855807909 -0.90 5.0 32.9msRemove 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.400427984616 -1.78 566ms
2 -8.400428156277 -6.77 -4.03 380ms
3 -8.400428156277 -14.45 -7.83 102msComparison 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| = 9.865262464816923e-7
|ρ_newton - ρ_scfv| = 2.6138953650540496e-8
|ρ_newton - ρ_dm| = 3.040292291772617e-6