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.847680505586 -0.70 4.2 27.9ms
2 -7.852388703731 -2.33 -1.53 1.0 16.8ms
3 -7.852608855028 -3.66 -2.54 1.2 17.5ms
4 -7.852645473561 -4.44 -2.78 2.5 21.8ms
5 -7.852646130446 -6.18 -2.89 1.2 17.3ms
6 -7.852646668229 -6.27 -3.87 1.0 16.8ms
7 -7.852646685960 -7.75 -4.60 2.0 20.0ms
8 -7.852646686689 -9.14 -5.15 1.5 18.2ms
9 -7.852646686726 -10.43 -6.28 1.2 18.1msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.847668809469 -0.70 4.8 29.0ms
2 -7.852559081664 -2.31 -1.62 0.80 2.2 19.9ms
3 -7.852641245248 -4.09 -2.69 0.80 1.0 15.4ms
4 -7.852646500486 -5.28 -3.39 0.80 1.8 18.2ms
5 -7.852646681385 -6.74 -4.38 0.80 2.0 18.7ms
6 -7.852646686645 -8.28 -4.89 0.80 2.5 21.2ms
7 -7.852646686726 -10.09 -5.58 0.80 1.0 15.6ms
8 -7.852646686730 -11.43 -6.64 0.80 1.8 17.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.314028346966 -0.99 54.8ms
2 -1.726867674529 0.48 -0.67 24.9ms
3 -3.891611559688 0.34 -0.44 32.9ms
4 -5.223477079509 0.12 -0.59 32.6ms
5 -6.880634702019 0.22 -0.77 32.7ms
6 -7.459678585182 -0.24 -1.25 32.6ms
7 -7.680838395671 -0.66 -1.40 24.6ms
8 -7.764821658843 -1.08 -1.81 24.5ms
9 -7.801693397895 -1.43 -2.16 24.6ms
10 -7.822336088811 -1.69 -2.21 23.9ms
11 -7.834292578520 -1.92 -2.40 24.5ms
12 -7.842117181193 -2.11 -2.50 24.7ms
13 -7.845702512644 -2.45 -2.44 24.6ms
14 -7.848138041284 -2.61 -2.73 24.6ms
15 -7.850450452343 -2.64 -2.98 24.5ms
16 -7.851734560129 -2.89 -3.02 24.4ms
17 -7.852270840983 -3.27 -3.53 24.3ms
18 -7.852515722900 -3.61 -3.61 24.3ms
19 -7.852609278459 -4.03 -3.66 24.6ms
20 -7.852635461079 -4.58 -4.06 24.6ms
21 -7.852642940022 -5.13 -3.97 24.7ms
22 -7.852645447891 -5.60 -4.18 24.6ms
23 -7.852646173663 -6.14 -5.02 24.5ms
24 -7.852646456613 -6.55 -4.71 24.5ms
25 -7.852646583400 -6.90 -5.08 24.6ms
26 -7.852646646298 -7.20 -5.83 24.7ms
27 -7.852646675741 -7.53 -5.51 24.6ms
28 -7.852646683689 -8.10 -5.63 24.6ms
29 -7.852646685817 -8.67 -6.38 24.7ms
30 -7.852646686451 -9.20 -6.13 24.7ms
31 -7.852646686635 -9.73 -6.80 24.8ms
32 -7.852646686700 -10.19 -7.21 24.7ms
33 -7.852646686721 -10.67 -7.41 24.7ms
34 -7.852646686727 -11.24 -7.51 24.7ms
35 -7.852646686728 -11.74 -8.02 24.6ms
36 -7.852646686729 -12.12 -7.89 24.6ms
37 -7.852646686730 -12.39 -8.29 24.6ms
38 -7.852646686730 -12.73 -8.23 24.6ms
39 -7.852646686730 -13.30 -8.70 24.6ms
40 -7.852646686730 -13.80 -8.95 24.6ms
41 -7.852646686730 -14.27 -9.03 24.6ms
42 -7.852646686730 -14.27 -9.09 24.5ms
43 -7.852646686730 -14.27 -9.19 30.8ms
44 -7.852646686730 + -Inf -9.61 101ms
45 -7.852646686730 + -Inf -8.49 32.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.847779003699 -0.70 4.5 29.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.852645852050 -1.63 404ms
2 -7.852646686730 -6.08 -3.69 287ms
3 -7.852646686730 -13.23 -7.20 111msComparison 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.677551900730093e-7
|ρ_newton - ρ_scfv| = 1.1440261958962013e-7
|ρ_newton - ρ_dm| = 3.957898427853214e-9