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.397792470904 -0.90 4.8 25.5ms
2 -8.400242797458 -2.61 -1.73 1.0 27.5ms
3 -8.400401431881 -3.80 -2.91 1.5 19.6ms
4 -8.400427755461 -4.58 -2.92 3.0 24.0ms
5 -8.400427963120 -6.68 -3.06 1.0 18.6ms
6 -8.400428150325 -6.73 -4.85 1.0 24.2ms
7 -8.400428155914 -8.25 -4.47 3.2 26.0ms
8 -8.400428156263 -9.46 -5.49 1.0 19.3ms
9 -8.400428156275 -10.91 -6.25 1.5 20.4msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397757239765 -0.89 4.8 25.6ms
2 -8.400366046872 -2.58 -1.77 0.80 2.2 18.9ms
3 -8.400419674706 -4.27 -2.95 0.80 1.0 15.6ms
4 -8.400428080287 -5.08 -3.27 0.80 2.8 20.1ms
5 -8.400428149975 -7.16 -4.30 0.80 1.0 15.8ms
6 -8.400428156065 -8.22 -5.09 0.80 2.2 24.7ms
7 -8.400428156273 -9.68 -6.22 0.80 2.5 19.9msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.869799203843 -1.06 51.5ms
2 -1.719057510066 0.41 -0.66 28.9ms
3 -4.053795415518 0.37 -0.40 44.5ms
4 -5.080269931408 0.01 -0.51 39.2ms
5 -7.098377550961 0.30 -0.55 39.3ms
6 -7.846156992087 -0.13 -1.07 43.7ms
7 -8.136620411841 -0.54 -1.31 29.0ms
8 -8.244006853318 -0.97 -1.57 29.0ms
9 -8.307607237130 -1.20 -1.70 33.9ms
10 -8.340882672324 -1.48 -2.01 29.1ms
11 -8.364561858459 -1.63 -2.17 29.0ms
12 -8.379829315697 -1.82 -2.24 29.0ms
13 -8.388949950540 -2.04 -2.12 29.0ms
14 -8.395551716420 -2.18 -2.64 33.6ms
15 -8.398318838254 -2.56 -2.40 29.0ms
16 -8.399790302132 -2.83 -2.74 29.0ms
17 -8.400149247073 -3.44 -2.81 28.9ms
18 -8.400332887340 -3.74 -3.38 33.9ms
19 -8.400378808694 -4.34 -3.21 29.0ms
20 -8.400414242477 -4.45 -3.66 29.0ms
21 -8.400419425655 -5.29 -3.62 29.0ms
22 -8.400425981403 -5.18 -4.08 33.5ms
23 -8.400426982109 -6.00 -4.12 29.1ms
24 -8.400427771711 -6.10 -4.82 29.0ms
25 -8.400427914308 -6.85 -4.47 29.0ms
26 -8.400428062717 -6.83 -5.36 34.3ms
27 -8.400428092584 -7.52 -4.68 29.7ms
28 -8.400428141765 -7.31 -5.16 29.3ms
29 -8.400428146715 -8.31 -5.23 28.9ms
30 -8.400428153637 -8.16 -5.71 29.0ms
31 -8.400428153641 -11.41 -5.74 33.8ms
32 -8.400428155477 -8.74 -6.40 29.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.397845583315 -0.90 5.0 25.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.400427979916 -1.79 536ms
2 -8.400428156277 -6.75 -4.03 370ms
3 -8.400428156277 -14.75 -7.81 123msComparison 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.6325355518889146e-6
|ρ_newton - ρ_scfv| = 5.73616892516896e-7
|ρ_newton - ρ_dm| = 3.2633679921662835e-6