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.oncvpsp3.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.397823823660 -0.90 5.8 36.3ms
2 -8.400250696577 -2.61 -1.74 1.0 26.9ms
3 -8.400398671478 -3.83 -2.93 1.2 33.1ms
4 -8.400428650524 -4.52 -2.91 3.2 26.3ms
5 -8.400428737663 -7.06 -2.96 1.0 20.0ms
6 -8.400429016894 -6.55 -4.91 1.0 19.7ms
7 -8.400429023880 -8.16 -4.36 3.5 33.5ms
8 -8.400429024380 -9.30 -4.89 2.0 23.4ms
9 -8.400429024419 -10.40 -6.26 1.0 20.1msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397767475238 -0.90 5.2 33.5ms
2 -8.400381673860 -2.58 -1.79 0.80 2.0 19.2ms
3 -8.400422852281 -4.39 -3.02 0.80 1.0 16.2ms
4 -8.400428967103 -5.21 -3.43 0.80 2.5 20.9ms
5 -8.400429019821 -7.28 -4.59 0.80 1.0 16.4ms
6 -8.400429024393 -8.34 -5.71 0.80 2.8 21.2ms
7 -8.400429024420 -10.56 -6.09 0.80 2.8 29.2msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.645081890766 -1.08 53.8ms
2 -1.582445044509 0.35 -0.67 29.7ms
3 -3.875017961513 0.36 -0.39 45.6ms
4 -4.774994147164 -0.05 -0.40 40.3ms
5 -6.794104993463 0.31 -0.53 40.1ms
6 -7.661148785690 -0.06 -0.80 46.5ms
7 -7.933006307087 -0.57 -1.32 30.4ms
8 -8.159287645716 -0.65 -1.76 29.7ms
9 -8.218646751820 -1.23 -1.85 29.8ms
10 -8.271935761035 -1.27 -1.95 35.0ms
11 -8.309388835252 -1.43 -1.95 30.0ms
12 -8.341101171990 -1.50 -2.13 29.7ms
13 -8.369868627760 -1.54 -2.41 30.8ms
14 -8.389404214142 -1.71 -2.26 34.6ms
15 -8.396948094341 -2.12 -2.70 30.0ms
16 -8.398939137135 -2.70 -2.76 29.6ms
17 -8.399894246026 -3.02 -3.04 30.0ms
18 -8.400147761804 -3.60 -3.13 34.7ms
19 -8.400329791378 -3.74 -3.19 32.9ms
20 -8.400395016582 -4.19 -4.12 33.6ms
21 -8.400413509306 -4.73 -3.78 30.1ms
22 -8.400421360320 -5.11 -4.38 34.8ms
23 -8.400423733725 -5.62 -4.16 30.0ms
24 -8.400427369677 -5.44 -4.56 29.8ms
25 -8.400428213141 -6.07 -4.52 29.4ms
26 -8.400428728278 -6.29 -4.74 34.7ms
27 -8.400428825966 -7.01 -4.66 29.9ms
28 -8.400428917051 -7.04 -5.41 29.9ms
29 -8.400428967768 -7.29 -4.93 29.5ms
30 -8.400429002840 -7.46 -5.51 34.6ms
31 -8.400429013191 -7.99 -5.40 29.8ms
32 -8.400429020025 -8.17 -5.96 29.8ms
33 -8.400429022585 -8.59 -5.62 29.5ms
34 -8.400429023496 -9.04 -6.15 36.7msNewton 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.397819474771 -0.90 5.2 27.2msRemove 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.400428845121 -1.78 559ms
2 -8.400429024420 -6.75 -4.02 362ms
3 -8.400429024420 -14.45 -7.80 113msComparison 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.77658025917077e-7
|ρ_newton - ρ_scfv| = 3.364858177986376e-7
|ρ_newton - ρ_dm| = 2.632320403273696e-6