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.847784922373 -0.70 4.8 35.9ms
2 -7.852412836877 -2.33 -1.53 1.0 18.5ms
3 -7.852611196162 -3.70 -2.54 1.0 18.2ms
4 -7.852645346415 -4.47 -2.78 2.0 23.4ms
5 -7.852646081798 -6.13 -2.87 1.0 18.4ms
6 -7.852646668236 -6.23 -3.95 1.0 18.3ms
7 -7.852646686048 -7.75 -4.43 1.8 22.3ms
8 -7.852646686694 -9.19 -5.43 1.2 19.8ms
9 -7.852646686727 -10.48 -5.56 2.5 25.6ms
10 -7.852646686729 -11.58 -5.95 1.0 19.3ms
11 -7.852646686730 -12.34 -6.59 1.0 18.8msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.847781266583 -0.70 4.8 36.1ms
2 -7.852561457636 -2.32 -1.62 0.80 2.0 22.3ms
3 -7.852641075024 -4.10 -2.71 0.80 1.0 17.1ms
4 -7.852646501951 -5.27 -3.38 0.80 1.8 21.0ms
5 -7.852646681880 -6.74 -4.42 0.80 1.8 21.0ms
6 -7.852646686619 -8.32 -4.87 0.80 2.2 24.0ms
7 -7.852646686726 -9.97 -5.51 0.80 1.2 18.5ms
8 -7.852646686729 -11.47 -6.62 0.80 1.5 19.6msDirect 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.286355930195 -1.01 57.0ms
2 -1.299846290199 0.41 -0.64 25.8ms
3 -3.573127666365 0.36 -0.43 33.5ms
4 -4.536687730223 -0.02 -0.57 33.5ms
5 -6.378414191374 0.27 -0.67 33.4ms
6 -7.112116295715 -0.13 -1.07 33.5ms
7 -7.435383518221 -0.49 -1.49 25.3ms
8 -7.634760662205 -0.70 -1.92 25.5ms
9 -7.709576554902 -1.13 -2.17 25.5ms
10 -7.738598312409 -1.54 -2.18 25.3ms
11 -7.759474704592 -1.68 -2.10 26.2ms
12 -7.781226858798 -1.66 -2.08 25.3ms
13 -7.809504486942 -1.55 -2.09 25.3ms
14 -7.827450132110 -1.75 -2.28 25.2ms
15 -7.838075909678 -1.97 -2.50 30.5ms
16 -7.846000975481 -2.10 -2.86 34.3ms
17 -7.850571298907 -2.34 -2.92 27.1ms
18 -7.852072079404 -2.82 -2.97 25.5ms
19 -7.852467760838 -3.40 -3.15 25.4ms
20 -7.852599963932 -3.88 -3.50 25.4ms
21 -7.852636332571 -4.44 -4.06 25.3ms
22 -7.852643497006 -5.14 -4.13 25.2ms
23 -7.852645747267 -5.65 -4.41 25.3ms
24 -7.852646331541 -6.23 -5.17 25.3ms
25 -7.852646549487 -6.66 -5.07 25.3ms
26 -7.852646639926 -7.04 -5.68 25.4ms
27 -7.852646667546 -7.56 -5.86 25.5ms
28 -7.852646680096 -7.90 -5.85 25.5ms
29 -7.852646684011 -8.41 -6.49 25.9ms
30 -7.852646685608 -8.80 -6.12 25.8ms
31 -7.852646686297 -9.16 -6.61 25.6ms
32 -7.852646686553 -9.59 -6.68 25.6ms
33 -7.852646686670 -9.93 -6.66 25.7ms
34 -7.852646686712 -10.38 -7.25 25.5ms
35 -7.852646686724 -10.92 -7.43 25.6ms
36 -7.852646686728 -11.42 -7.30 25.4ms
37 -7.852646686729 -11.83 -7.56 25.6ms
38 -7.852646686730 -12.27 -8.25 25.7ms
39 -7.852646686730 -12.76 -7.93 25.6ms
40 -7.852646686730 -13.30 -8.22 25.5ms
41 -7.852646686730 -14.05 -9.08 25.6ms
42 -7.852646686730 -14.35 -9.13 25.7ms
43 -7.852646686730 -14.75 -9.21 33.7ms
44 -7.852646686730 -14.75 -9.17 25.4ms
45 -7.852646686730 + -Inf -9.96 33.4ms
46 -7.852646686730 + -Inf -9.78 33.5ms
┌ 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.847775126146 -0.70 4.8 35.8msRemove 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.852645843584 -1.63 390ms
2 -7.852646686730 -6.07 -3.69 289ms
3 -7.852646686730 -13.23 -7.20 112msComparison 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| = 2.4563028441990545e-7
|ρ_newton - ρ_scfv| = 5.821154799600917e-7
|ρ_newton - ρ_dm| = 5.998030954344661e-10