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.847813849318 -0.70 4.5 29.5ms
2 -7.852412540530 -2.34 -1.53 1.0 17.0ms
3 -7.852610829942 -3.70 -2.55 1.0 16.7ms
4 -7.852645464017 -4.46 -2.78 2.2 21.4ms
5 -7.852646095994 -6.20 -2.88 1.0 16.9ms
6 -7.852646666224 -6.24 -3.83 1.0 17.2ms
7 -7.852646684709 -7.73 -4.60 1.2 17.9ms
8 -7.852646686686 -8.70 -5.03 2.0 21.0ms
9 -7.852646686726 -10.41 -5.78 1.2 18.1ms
10 -7.852646686728 -11.57 -5.64 1.8 20.6ms
11 -7.852646686730 -11.86 -6.14 1.0 29.8msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.847760240034 -0.70 5.0 30.5ms
2 -7.852559190284 -2.32 -1.62 0.80 2.0 19.7ms
3 -7.852641752770 -4.08 -2.68 0.80 1.0 15.7ms
4 -7.852646509485 -5.32 -3.43 0.80 1.8 18.6ms
5 -7.852646681417 -6.76 -4.39 0.80 2.0 19.0ms
6 -7.852646686614 -8.28 -4.86 0.80 2.5 21.7ms
7 -7.852646686725 -9.95 -5.52 0.80 1.0 15.9ms
8 -7.852646686729 -11.38 -6.64 0.80 1.8 18.3msDirect 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.289714437899 -1.04 44.4ms
2 -1.354926437320 0.42 -0.70 24.8ms
3 -3.656180602273 0.36 -0.47 32.8ms
4 -4.712803391798 0.02 -0.58 33.0ms
5 -6.607905445474 0.28 -0.68 32.8ms
6 -7.333718723848 -0.14 -1.12 32.8ms
7 -7.644586616196 -0.51 -1.41 24.8ms
8 -7.753806486603 -0.96 -1.86 24.7ms
9 -7.797048807441 -1.36 -2.15 24.8ms
10 -7.819300213693 -1.65 -2.14 24.7ms
11 -7.835211614059 -1.80 -2.18 24.9ms
12 -7.842470279422 -2.14 -2.41 24.7ms
13 -7.846763713019 -2.37 -2.60 24.7ms
14 -7.848553883071 -2.75 -3.10 24.7ms
15 -7.850565083726 -2.70 -3.21 24.8ms
16 -7.851866136754 -2.89 -3.36 24.6ms
17 -7.852403449189 -3.27 -3.25 24.7ms
18 -7.852570070248 -3.78 -3.31 24.7ms
19 -7.852621400661 -4.29 -4.12 24.9ms
20 -7.852637708265 -4.79 -4.13 24.7ms
21 -7.852643623282 -5.23 -4.31 24.7ms
22 -7.852645665172 -5.69 -4.89 24.7ms
23 -7.852646287425 -6.21 -5.05 24.9ms
24 -7.852646527888 -6.62 -5.12 24.8ms
25 -7.852646623050 -7.02 -5.59 24.7ms
26 -7.852646659174 -7.44 -5.58 24.7ms
27 -7.852646677709 -7.73 -5.78 25.0ms
28 -7.852646684627 -8.16 -6.09 24.7ms
29 -7.852646686014 -8.86 -6.44 24.8ms
30 -7.852646686462 -9.35 -6.64 24.8ms
31 -7.852646686633 -9.77 -6.80 24.8ms
32 -7.852646686689 -10.25 -7.08 32.2ms
33 -7.852646686717 -10.55 -7.20 25.5ms
34 -7.852646686726 -11.03 -7.69 24.9ms
35 -7.852646686729 -11.59 -7.87 24.7ms
36 -7.852646686729 -12.08 -7.84 24.7ms
37 -7.852646686730 -12.57 -8.05 24.7ms
38 -7.852646686730 -13.01 -8.34 31.4ms
39 -7.852646686730 -13.60 -8.87 24.8ms
40 -7.852646686730 -14.15 -8.98 24.9ms
41 -7.852646686730 -14.27 -8.98 24.8ms
42 -7.852646686730 -14.45 -9.54 24.8ms
┌ 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.847774083254 -0.70 4.8 39.0msRemove 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.852645844027 -1.63 383ms
2 -7.852646686730 -6.07 -3.69 289ms
3 -7.852646686730 -13.21 -7.20 119msComparison 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| = 7.246569562664095e-7
|ρ_newton - ρ_scfv| = 6.003049770832736e-7
|ρ_newton - ρ_dm| = 2.052852744478766e-9