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 DFTK
using LinearAlgebra
a = 10.26 # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
[1 0 1.];
[1 1 0.]]
Si = ElementPsp(:Si; psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]
model = model_LDA(lattice, atoms, positions)
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 -7.847756544273 -0.70 4.5 35.3ms
2 -7.852403528783 -2.33 -1.53 1.0 19.0ms
3 -7.852610395495 -3.68 -2.55 1.2 20.0ms
4 -7.852645296118 -4.46 -2.77 2.0 24.5ms
5 -7.852646045339 -6.13 -2.86 1.2 20.0ms
6 -7.852646666423 -6.21 -3.94 1.0 18.6ms
7 -7.852646686208 -7.70 -4.51 2.0 23.2ms
8 -7.852646686705 -9.30 -5.46 1.2 19.6ms
9 -7.852646686724 -10.72 -5.41 2.2 24.5ms
10 -7.852646686729 -11.27 -6.01 1.0 19.0msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.847742880307 -0.70 4.8 36.1ms
2 -7.852560037552 -2.32 -1.62 0.80 2.0 22.5ms
3 -7.852640557742 -4.09 -2.71 0.80 1.0 17.1ms
4 -7.852646523033 -5.22 -3.36 0.80 2.0 22.0ms
5 -7.852646681379 -6.80 -4.35 0.80 1.5 19.8ms
6 -7.852646686618 -8.28 -4.83 0.80 2.2 42.5ms
7 -7.852646686726 -9.97 -5.51 0.80 1.0 17.7ms
8 -7.852646686729 -11.41 -6.80 0.80 1.8 20.9msDirect minimization
Note: Unlike the other algorithms, tolerance for this one is in the energy, thus we square the density tolerance value to be roughly equivalent.
scfres_dm = direct_minimization(basis; tol=tol^2);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.339642711599 -1.02 45.7ms
2 -1.695465957873 0.48 -0.63 25.4ms
3 -3.455251673248 0.25 -0.43 33.5ms
4 -4.660286048631 0.08 -0.55 33.2ms
5 -6.290902913251 0.21 -0.64 33.5ms
6 -7.155006548665 -0.06 -0.97 33.4ms
7 -7.429532076731 -0.56 -1.52 25.2ms
8 -7.669741929177 -0.62 -1.80 25.3ms
9 -7.728812920525 -1.23 -1.95 25.2ms
10 -7.777966819068 -1.31 -1.96 25.3ms
11 -7.811904316156 -1.47 -2.07 25.4ms
12 -7.831034187211 -1.72 -2.38 25.4ms
13 -7.845039616951 -1.85 -2.23 25.2ms
14 -7.847391979402 -2.63 -2.98 25.3ms
15 -7.848966658762 -2.80 -2.73 25.3ms
16 -7.850857721567 -2.72 -3.40 25.3ms
17 -7.852024719079 -2.93 -3.21 25.6ms
18 -7.852427997472 -3.39 -3.71 25.3ms
19 -7.852573322360 -3.84 -3.91 25.2ms
20 -7.852621243407 -4.32 -4.10 25.1ms
21 -7.852637839899 -4.78 -4.35 25.4ms
22 -7.852643701397 -5.23 -4.64 25.3ms
23 -7.852645528771 -5.74 -4.81 25.2ms
24 -7.852646258490 -6.14 -5.08 25.4ms
25 -7.852646525632 -6.57 -5.32 25.5ms
26 -7.852646632049 -6.97 -5.25 25.6ms
27 -7.852646671435 -7.40 -5.37 26.0ms
28 -7.852646682255 -7.97 -6.05 25.5ms
29 -7.852646685380 -8.51 -5.92 25.5ms
30 -7.852646686260 -9.06 -6.01 25.4ms
31 -7.852646686566 -9.51 -6.31 25.5ms
32 -7.852646686682 -9.94 -6.72 25.4ms
33 -7.852646686715 -10.48 -6.90 25.5ms
34 -7.852646686725 -11.01 -7.06 25.4ms
35 -7.852646686728 -11.52 -7.69 25.3ms
36 -7.852646686729 -11.86 -7.85 25.4ms
37 -7.852646686729 -12.24 -8.03 25.3ms
38 -7.852646686730 -12.63 -8.33 25.3ms
39 -7.852646686730 -13.01 -8.76 25.3ms
40 -7.852646686730 -13.36 -8.82 25.4ms
41 -7.852646686730 -13.88 -9.20 25.6ms
42 -7.852646686730 -14.45 -9.16 25.4ms
43 -7.852646686730 + -Inf -9.59 25.3ms
44 -7.852646686730 -14.45 -9.53 25.3ms
45 -7.852646686730 + -Inf -9.76 25.4ms
46 -7.852646686730 + -Inf -10.15 33.6ms
┌ Warning: DM not converged.
└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/scf_callbacks.jl:60Newton 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 -7.847789314938 -0.70 4.8 36.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 -7.852645851205 -1.63 403ms
2 -7.852646686730 -6.08 -3.69 285ms
3 -7.852646686730 -13.21 -7.20 143msComparison 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.84685640649203e-7
|ρ_newton - ρ_scfv| = 6.792426351774256e-7
|ρ_newton - ρ_dm| = 9.527468030590699e-10