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.846940822172 -0.70 4.8
2 -7.852336526202 -2.27 -1.53 1.0 38.7ms
3 -7.852616947394 -3.55 -2.56 1.2 18.0ms
4 -7.852646013551 -4.54 -2.89 2.5 23.1ms
5 -7.852646509208 -6.30 -3.19 1.0 17.2ms
6 -7.852646680349 -6.77 -4.07 1.2 18.0ms
7 -7.852646686409 -8.22 -5.22 1.8 20.5ms
8 -7.852646686724 -9.50 -5.49 2.2 22.2ms
9 -7.852646686729 -11.29 -5.91 1.2 18.2ms
10 -7.852646686730 -12.20 -6.34 1.2 18.7msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.846876925013 -0.70 4.5
2 -7.852528531221 -2.25 -1.63 0.80 2.0 19.7ms
3 -7.852637238439 -3.96 -2.71 0.80 1.0 15.8ms
4 -7.852646514414 -5.03 -3.29 0.80 2.2 21.1ms
5 -7.852646674111 -6.80 -4.15 0.80 1.0 16.0ms
6 -7.852646686551 -7.91 -4.80 0.80 1.8 19.2ms
7 -7.852646686724 -9.76 -5.53 0.80 1.5 18.2ms
8 -7.852646686730 -11.29 -6.59 0.80 1.8 18.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);Iter Function value Gradient norm
0 1.394672e+01 3.544632e+00
* time: 0.009553909301757812
1 1.385849e+00 1.794337e+00
* time: 0.02740788459777832
2 -1.185393e+00 2.190993e+00
* time: 0.04525589942932129
3 -3.617802e+00 1.854679e+00
* time: 0.07098102569580078
4 -4.977092e+00 1.641886e+00
* time: 0.09677791595458984
5 -6.758541e+00 8.713809e-01
* time: 0.12250304222106934
6 -6.848992e+00 1.647052e+00
* time: 0.1402430534362793
7 -7.472111e+00 1.300939e+00
* time: 0.15796804428100586
8 -7.486335e+00 2.475845e+00
* time: 0.17573809623718262
9 -7.621276e+00 2.282987e+00
* time: 0.19350790977478027
10 -7.639091e+00 1.951443e+00
* time: 0.2192380428314209
11 -7.646306e+00 1.845222e+00
* time: 0.2449660301208496
12 -7.667894e+00 1.432112e+00
* time: 0.270658016204834
13 -7.683349e+00 1.207230e+00
* time: 0.29653096199035645
14 -7.715872e+00 6.314722e-01
* time: 0.32306909561157227
15 -7.729993e+00 1.486678e+00
* time: 0.34087705612182617
16 -7.760390e+00 7.062637e-01
* time: 0.35875606536865234
17 -7.772689e+00 6.727363e-01
* time: 0.3845078945159912
18 -7.793503e+00 6.439758e-01
* time: 0.40236496925354004
19 -7.817270e+00 4.881956e-01
* time: 0.42012906074523926
20 -7.830550e+00 3.453619e-01
* time: 0.4378840923309326
21 -7.840620e+00 1.092378e-01
* time: 0.4557220935821533
22 -7.845939e+00 5.563220e-02
* time: 0.4920768737792969
23 -7.851151e+00 3.876657e-02
* time: 0.5178570747375488
24 -7.852210e+00 1.980195e-02
* time: 0.5382630825042725
25 -7.852483e+00 1.456035e-02
* time: 0.5561549663543701
26 -7.852611e+00 6.933888e-03
* time: 0.5739340782165527
27 -7.852638e+00 3.404811e-03
* time: 0.591731071472168
28 -7.852643e+00 2.032716e-03
* time: 0.6094698905944824
29 -7.852645e+00 1.024032e-03
* time: 0.6271910667419434
30 -7.852646e+00 4.012938e-04
* time: 0.6449639797210693
31 -7.852646e+00 3.457974e-04
* time: 0.6627469062805176
32 -7.852647e+00 2.393512e-04
* time: 0.6805789470672607
33 -7.852647e+00 1.441701e-04
* time: 0.6984009742736816
34 -7.852647e+00 5.113364e-05
* time: 0.7161760330200195
35 -7.852647e+00 4.267315e-05
* time: 0.7339038848876953
36 -7.852647e+00 3.439727e-05
* time: 0.7516989707946777
37 -7.852647e+00 2.289903e-05
* time: 0.7695050239562988
38 -7.852647e+00 7.868150e-06
* time: 0.7872960567474365
39 -7.852647e+00 5.090627e-06
* time: 0.8050329685211182
40 -7.852647e+00 3.425136e-06
* time: 0.8227748870849609
41 -7.852647e+00 1.263145e-06
* time: 0.8405370712280273
42 -7.852647e+00 9.913267e-07
* time: 0.8583879470825195
43 -7.852647e+00 6.137473e-07
* time: 0.876147985458374
44 -7.852647e+00 3.478600e-07
* time: 0.8939189910888672
45 -7.852647e+00 2.894924e-07
* time: 0.9116950035095215
46 -7.852647e+00 9.339685e-08
* time: 0.9294009208679199
47 -7.852647e+00 7.359717e-08
* time: 0.9471628665924072
48 -7.852647e+00 4.499015e-08
* time: 0.9649350643157959
49 -7.852647e+00 1.524489e-08
* time: 0.9826650619506836
50 -7.852647e+00 1.132541e-08
* time: 1.0005860328674316
51 -7.852647e+00 9.684323e-09
* time: 1.0342068672180176
52 -7.852647e+00 7.387184e-09
* time: 1.059967041015625Newton 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.846883795509 -0.70 4.8Remove 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.852645919875 -1.64
2 -7.852646686730 -6.12 -3.71 290ms
3 -7.852646686730 -13.30 -7.25 121msComparison 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.3302220192353144e-7
|ρ_newton - ρ_scfv| = 2.567017146265263e-7
|ρ_newton - ρ_dm| = 6.109249844477483e-10