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.846777262291 -0.70 4.5
2 -7.852310526801 -2.26 -1.53 1.0
3 -7.852645905985 -3.47 -2.53 3.0
4 -7.852646675666 -6.11 -3.35 2.8
5 -7.852646685780 -8.00 -4.75 1.2
6 -7.852646686726 -9.02 -5.23 4.2
7 -7.852646686729 -11.44 -5.89 2.0
8 -7.852646686730 -12.46 -7.49 2.0
9 -7.852646686730 -14.75 -8.05 2.5
10 -7.852646686730 -14.75 -8.97 1.2
11 -7.852646686730 + -15.05 -9.51 2.0
12 -7.852646686730 -14.57 -8.81 1.0
13 -7.852646686730 + -14.35 -8.28 2.0
14 -7.852646686730 -15.05 -9.06 1.2
15 -7.852646686730 + -Inf -8.55 1.0
16 -7.852646686730 -14.75 -8.24 2.0
17 -7.852646686730 + -14.57 -9.88 1.8
18 -7.852646686730 -14.57 -9.89 2.0
19 -7.852646686730 + -14.75 -10.48 1.2
20 -7.852646686730 + -14.75 -10.46 1.0
21 -7.852646686730 + -14.75 -11.89 1.0
22 -7.852646686730 -14.35 -11.28 1.0
23 -7.852646686730 + -15.05 -11.32 1.0
24 -7.852646686730 + -Inf -11.48 1.0
25 -7.852646686730 + -Inf -11.86 1.0
26 -7.852646686730 -15.05 -13.45 2.0Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; is_converged);n Energy log10(ΔE) log10(Δρ) α Diag
--- --------------- --------- --------- ---- ----
1 -7.846842889337 -0.70 4.5
2 -7.852530500139 -2.25 -1.63 0.80 2.2
3 -7.852637976478 -3.97 -2.72 0.80 1.0
4 -7.852646355774 -5.08 -3.33 0.80 2.0
5 -7.852646682776 -6.49 -4.23 0.80 1.8
6 -7.852646686592 -8.42 -4.90 0.80 1.8
7 -7.852646686722 -9.89 -5.78 0.80 1.8
8 -7.852646686730 -11.11 -6.36 0.80 2.0
9 -7.852646686730 -12.92 -7.25 0.80 1.2
10 -7.852646686730 -13.88 -7.98 0.80 2.5
11 -7.852646686730 + -14.57 -8.66 0.80 1.8
12 -7.852646686730 -15.05 -8.99 0.80 2.0
13 -7.852646686730 -14.75 -10.20 0.80 1.2
14 -7.852646686730 + -14.27 -10.70 0.80 2.2
15 -7.852646686730 -14.75 -11.19 0.80 1.2
16 -7.852646686730 -14.45 -12.22 0.80 1.5Direct minimization
scfres_dm = direct_minimization(basis; tol);Iter Function value Gradient norm
0 1.426441e+01 3.846956e+00
* time: 0.01440286636352539
1 1.292289e+00 1.863167e+00
* time: 0.041684865951538086
2 -1.880820e+00 2.121376e+00
* time: 0.06896805763244629
3 -3.609945e+00 1.843785e+00
* time: 0.10801887512207031
4 -5.177782e+00 1.499845e+00
* time: 0.14759397506713867
5 -6.721600e+00 7.658730e-01
* time: 0.1866140365600586
6 -6.795059e+00 9.442784e-01
* time: 0.2139148712158203
7 -7.524428e+00 4.985537e-01
* time: 0.24129390716552734
8 -7.692265e+00 3.946382e-01
* time: 0.26845788955688477
9 -7.775121e+00 1.142263e-01
* time: 0.29541802406311035
10 -7.807626e+00 1.177294e-01
* time: 0.32247495651245117
11 -7.825249e+00 8.087177e-02
* time: 0.3495960235595703
12 -7.838730e+00 5.654729e-02
* time: 0.37686586380004883
13 -7.847248e+00 3.967695e-02
* time: 0.40372800827026367
14 -7.849716e+00 2.783006e-02
* time: 0.43066906929016113
15 -7.850933e+00 2.038467e-02
* time: 0.4578440189361572
16 -7.851889e+00 1.997258e-02
* time: 0.4850270748138428
17 -7.852416e+00 1.278973e-02
* time: 0.5118730068206787
18 -7.852577e+00 4.805168e-03
* time: 0.5387940406799316
19 -7.852622e+00 3.132789e-03
* time: 0.5656988620758057
20 -7.852638e+00 1.650355e-03
* time: 0.5927770137786865
21 -7.852644e+00 9.781338e-04
* time: 0.6198439598083496
22 -7.852646e+00 6.153271e-04
* time: 0.6473119258880615
23 -7.852646e+00 3.496476e-04
* time: 0.6742198467254639
24 -7.852647e+00 1.850917e-04
* time: 0.7011339664459229
25 -7.852647e+00 9.859993e-05
* time: 0.728126049041748
26 -7.852647e+00 6.753300e-05
* time: 0.7553088665008545
27 -7.852647e+00 4.875409e-05
* time: 0.7821569442749023
28 -7.852647e+00 1.616915e-05
* time: 0.8090970516204834
29 -7.852647e+00 9.364682e-06
* time: 0.8359298706054688
30 -7.852647e+00 6.121000e-06
* time: 0.8628778457641602
31 -7.852647e+00 3.238488e-06
* time: 0.8896679878234863
32 -7.852647e+00 2.134854e-06
* time: 0.9165499210357666
33 -7.852647e+00 1.451816e-06
* time: 0.9436378479003906
34 -7.852647e+00 5.900540e-07
* time: 0.9710679054260254
35 -7.852647e+00 3.404332e-07
* time: 0.9980978965759277
36 -7.852647e+00 2.247658e-07
* time: 1.0251379013061523
37 -7.852647e+00 1.441273e-07
* time: 1.0521109104156494
38 -7.852647e+00 1.203774e-07
* time: 1.079028844833374
39 -7.852647e+00 6.529254e-08
* time: 1.1058149337768555
40 -7.852647e+00 2.673794e-08
* time: 1.13277006149292
41 -7.852647e+00 1.608416e-08
* time: 1.1595900058746338
42 -7.852647e+00 1.025327e-08
* time: 1.1864609718322754
43 -7.852647e+00 8.314678e-09
* time: 1.2133948802947998
44 -7.852647e+00 5.520029e-09
* time: 1.2526710033416748
45 -7.852647e+00 5.441449e-09
* time: 1.303725004196167Newton 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.846845132465 -0.70 4.8
2 -7.852322193269 -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.852646686713 -2.55
2 -7.852646686730 -10.77 -6.02
3 -7.852646686730 + -15.05 -12.68Comparison 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| = 3.843254934083239e-14
|ρ_newton - ρ_scfv| = 7.858717003474331e-13
|ρ_newton - ρ_dm| = 6.626397645110732e-10