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
tol = 1e-12
is_converged = DFTK.ScfConvergenceDensity(tol);Density-based self-consistent field
scfres_scf = self_consistent_field(basis; is_converged);n Energy log10(ΔE) log10(Δρ) Diag
--- --------------- --------- --------- ----
1 -7.846879465649 -0.70 4.5
2 -7.852323288389 -2.26 -1.53 1.0
3 -7.852646227872 -3.49 -2.52 3.0
4 -7.852646678225 -6.35 -3.37 2.5
5 -7.852646686016 -8.11 -4.70 1.2
6 -7.852646686726 -9.15 -5.31 3.5
7 -7.852646686730 -11.44 -6.22 1.2
8 -7.852646686730 -12.74 -7.63 2.5
9 -7.852646686730 + -15.05 -7.75 3.0
10 -7.852646686730 -15.05 -8.70 2.0
11 -7.852646686730 + -Inf -8.88 1.0
12 -7.852646686730 -15.05 -10.61 2.0
13 -7.852646686730 + -15.05 -10.50 2.0
14 -7.852646686730 -15.05 -11.65 1.5
15 -7.852646686730 + -15.05 -12.72 2.8Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; is_converged);n Energy log10(ΔE) log10(Δρ) α Diag
--- --------------- --------- --------- ---- ----
1 -7.846855597255 -0.70 4.8
2 -7.852525576107 -2.25 -1.64 0.80 2.0
3 -7.852634641447 -3.96 -2.74 0.80 1.0
4 -7.852646498626 -4.93 -3.22 0.80 2.2
5 -7.852646672112 -6.76 -4.02 0.80 1.0
6 -7.852646686240 -7.85 -4.69 0.80 1.5
7 -7.852646686718 -9.32 -5.60 0.80 1.8
8 -7.852646686729 -10.95 -6.76 0.80 1.8
9 -7.852646686730 -12.43 -7.49 0.80 2.8
10 -7.852646686730 -14.27 -7.96 0.80 2.0
11 -7.852646686730 -15.05 -9.13 0.80 1.0
12 -7.852646686730 + -15.05 -9.50 0.80 2.8
13 -7.852646686730 + -Inf -10.63 0.80 1.0
14 -7.852646686730 -14.75 -10.62 0.80 2.8
15 -7.852646686730 + -14.45 -12.05 0.80 1.0Direct minimization
scfres_dm = direct_minimization(basis; tol);Iter Function value Gradient norm
0 1.362834e+01 3.671831e+00
* time: 0.013882160186767578
1 1.041927e+00 2.087498e+00
* time: 0.04104304313659668
2 -1.838201e+00 2.461604e+00
* time: 0.06888198852539062
3 -3.889363e+00 2.130812e+00
* time: 0.11160898208618164
4 -5.223937e+00 1.673251e+00
* time: 0.1499011516571045
5 -6.797038e+00 1.250613e+00
* time: 0.18706011772155762
6 -7.347661e+00 4.849194e-01
* time: 0.22365903854370117
7 -7.532242e+00 1.615324e-01
* time: 0.249406099319458
8 -7.625804e+00 2.244724e-01
* time: 0.27483415603637695
9 -7.690670e+00 1.126843e-01
* time: 0.30014514923095703
10 -7.737975e+00 1.402996e-01
* time: 0.32523202896118164
11 -7.787923e+00 1.372769e-01
* time: 0.3504910469055176
12 -7.820711e+00 8.621622e-02
* time: 0.3774580955505371
13 -7.844840e+00 3.536448e-02
* time: 0.4040100574493408
14 -7.849804e+00 3.987997e-02
* time: 0.4292171001434326
15 -7.851889e+00 2.367369e-02
* time: 0.4546520709991455
16 -7.852439e+00 7.460525e-03
* time: 0.48195505142211914
17 -7.852581e+00 8.377179e-03
* time: 0.5072400569915771
18 -7.852610e+00 4.217875e-03
* time: 0.5326650142669678
19 -7.852631e+00 1.868054e-03
* time: 0.5580241680145264
20 -7.852642e+00 1.823901e-03
* time: 0.5841591358184814
21 -7.852645e+00 1.021985e-03
* time: 0.6095001697540283
22 -7.852646e+00 3.768700e-04
* time: 0.6392440795898438
23 -7.852647e+00 3.033489e-04
* time: 0.6668670177459717
24 -7.852647e+00 1.611514e-04
* time: 0.6950380802154541
25 -7.852647e+00 1.158842e-04
* time: 0.7224111557006836
26 -7.852647e+00 3.677708e-05
* time: 0.7546300888061523
27 -7.852647e+00 2.004913e-05
* time: 0.7818961143493652
28 -7.852647e+00 1.334430e-05
* time: 0.8101541996002197
29 -7.852647e+00 6.533916e-06
* time: 0.8373889923095703
30 -7.852647e+00 7.269142e-06
* time: 0.8644170761108398
31 -7.852647e+00 2.925066e-06
* time: 0.9419059753417969
32 -7.852647e+00 1.671245e-06
* time: 0.969520092010498
33 -7.852647e+00 1.091535e-06
* time: 0.9981181621551514
34 -7.852647e+00 6.211008e-07
* time: 1.0241641998291016
35 -7.852647e+00 5.391051e-07
* time: 1.0497171878814697
36 -7.852647e+00 3.728997e-07
* time: 1.0750751495361328
37 -7.852647e+00 9.735543e-08
* time: 1.102931022644043
38 -7.852647e+00 6.115073e-08
* time: 1.1295599937438965
39 -7.852647e+00 5.383324e-08
* time: 1.1558520793914795
40 -7.852647e+00 2.491115e-08
* time: 1.1948490142822266
41 -7.852647e+00 1.765566e-08
* time: 1.2218589782714844
42 -7.852647e+00 1.404257e-08
* time: 1.2488081455230713
43 -7.852647e+00 1.372267e-08
* time: 1.299713134765625
44 -7.852647e+00 1.372267e-08
* time: 1.4092340469360352Newton 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=1e-1);n Energy log10(ΔE) log10(Δρ) Diag
--- --------------- --------- --------- ----
1 -7.846791291173 -0.70 4.5
2 -7.852299197883 -2.26 -1.53 1.0Remove 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(Δρ)
--- --------------- --------- ---------
1 -7.852646686579 -2.54
2 -7.852646686730 -9.82 -5.85
3 -7.852646686730 + -Inf -11.61
4 -7.852646686730 + -Inf -15.97Comparison 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| = 1.3866702276178619e-13
|ρ_newton - ρ_scfv| = 2.045522491196992e-12
|ρ_newton - ρ_dm| = 5.411763012138437e-10