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.397818572677 -0.90 5.5 27.2ms
2 -8.400242521428 -2.62 -1.74 1.0 28.4ms
3 -8.400405848455 -3.79 -2.97 1.5 20.3ms
4 -8.400428719961 -4.64 -2.98 3.2 25.1ms
5 -8.400428788690 -7.16 -3.02 1.0 25.0ms
6 -8.400429017610 -6.64 -4.52 1.0 19.4ms
7 -8.400429023929 -8.20 -4.42 2.5 23.6ms
8 -8.400429024381 -9.35 -5.02 1.0 23.4ms
9 -8.400429024418 -10.43 -6.70 1.2 26.1msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397807802645 -0.90 5.2 27.1ms
2 -8.400386459407 -2.59 -1.78 0.80 2.0 18.7ms
3 -8.400423751817 -4.43 -3.00 0.80 1.0 21.6ms
4 -8.400428971403 -5.28 -3.43 0.80 2.2 20.0ms
5 -8.400429020369 -7.31 -4.64 0.80 1.0 15.9ms
6 -8.400429024391 -8.40 -5.62 0.80 2.8 20.9ms
7 -8.400429024420 -10.55 -6.06 0.80 2.5 20.5msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.873710274638 -1.12 57.9ms
2 -1.939697762866 0.45 -0.68 29.1ms
3 -4.029117176812 0.32 -0.42 43.7ms
4 -4.949784836629 -0.04 -0.48 58.8ms
5 -6.893081744060 0.29 -0.52 39.7ms
6 -7.651819735263 -0.12 -0.91 45.0ms
7 -7.926321084876 -0.56 -1.29 29.2ms
8 -8.199925622069 -0.56 -1.71 34.4ms
9 -8.248954679408 -1.31 -1.76 41.3ms
10 -8.282728502588 -1.47 -1.61 35.2ms
11 -8.326400136859 -1.36 -1.58 39.7ms
12 -8.365887863726 -1.40 -2.01 44.3ms
13 -8.378230169664 -1.91 -2.39 29.5ms
14 -8.389584415811 -1.94 -2.39 29.3ms
15 -8.393596950358 -2.40 -2.83 29.2ms
16 -8.396448595836 -2.54 -2.89 34.3ms
17 -8.398578795306 -2.67 -2.74 32.1ms
18 -8.399086282481 -3.29 -2.58 39.6ms
19 -8.399453431244 -3.44 -3.10 33.3ms
20 -8.399821733346 -3.43 -3.38 31.0ms
21 -8.400235611506 -3.38 -3.16 29.3ms
22 -8.400345178199 -3.96 -4.37 29.1ms
23 -8.400393550675 -4.32 -3.53 34.6ms
24 -8.400415130026 -4.67 -3.96 29.5ms
25 -8.400425891700 -4.97 -5.05 29.2ms
26 -8.400427515738 -5.79 -4.72 29.2ms
27 -8.400428331609 -6.09 -5.09 33.7ms
28 -8.400428653700 -6.49 -5.11 29.2ms
29 -8.400428898016 -6.61 -5.30 29.1ms
30 -8.400428970130 -7.14 -5.29 33.9ms
31 -8.400429007019 -7.43 -5.61 29.4ms
32 -8.400429018247 -7.95 -5.73 29.2ms
33 -8.400429021569 -8.48 -5.89 29.4ms
34 -8.400429023365 -8.75 -6.00 57.1ms
35 -8.400429023780 -9.38 -6.26 29.5msNewton 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.397738140049 -0.90 5.5 30.5msRemove 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.400428856737 -1.78 545ms
2 -8.400429024420 -6.78 -4.04 354ms
3 -8.400429024420 + -Inf -7.84 129msComparison 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| = 4.067498356477956e-7
|ρ_newton - ρ_scfv| = 2.7616369834993597e-7
|ρ_newton - ρ_dm| = 3.919434565591946e-7