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.846825135780 -0.70 4.5
2 -7.852318901218 -2.26 -1.53 1.0 22.6ms
3 -7.852615300516 -3.53 -2.55 1.5 24.2ms
4 -7.852645992516 -4.51 -2.88 2.8 31.3ms
5 -7.852646507618 -6.29 -3.18 1.0 22.3ms
6 -7.852646679643 -6.76 -4.02 1.0 23.0ms
7 -7.852646686349 -8.17 -5.12 1.5 26.2ms
8 -7.852646686721 -9.43 -5.36 2.2 29.1ms
9 -7.852646686730 -11.05 -6.59 1.0 22.5msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.846922023684 -0.70 4.8
2 -7.852526805486 -2.25 -1.64 0.80 2.0 26.1ms
3 -7.852635967438 -3.96 -2.72 0.80 1.0 21.8ms
4 -7.852646502942 -4.98 -3.27 0.80 2.2 28.0ms
5 -7.852646672802 -6.77 -4.09 0.80 1.0 21.2ms
6 -7.852646686445 -7.87 -4.75 0.80 1.8 26.2ms
7 -7.852646686725 -9.55 -5.49 0.80 1.8 25.4ms
8 -7.852646686729 -11.37 -6.67 0.80 1.5 24.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.387706e+01 3.719426e+00
* time: 0.012444019317626953
1 1.269083e+00 1.964536e+00
* time: 0.03523516654968262
2 -1.427284e+00 2.411262e+00
* time: 0.05799508094787598
3 -3.841403e+00 1.850400e+00
* time: 0.09041309356689453
4 -5.118125e+00 1.510901e+00
* time: 0.12261819839477539
5 -6.799177e+00 8.163126e-01
* time: 0.15502214431762695
6 -7.446443e+00 2.687974e-01
* time: 0.1872570514678955
7 -7.674390e+00 1.338245e-01
* time: 0.20988798141479492
8 -7.769253e+00 1.214316e-01
* time: 0.23241114616394043
9 -7.814742e+00 7.927568e-02
* time: 0.25545215606689453
10 -7.837702e+00 5.893711e-02
* time: 0.27816009521484375
11 -7.846280e+00 5.074035e-02
* time: 0.300616979598999
12 -7.850021e+00 2.384647e-02
* time: 0.3230421543121338
13 -7.851658e+00 1.389295e-02
* time: 0.34543418884277344
14 -7.852259e+00 1.072655e-02
* time: 0.3678281307220459
15 -7.852478e+00 6.281028e-03
* time: 0.3905150890350342
16 -7.852581e+00 3.552732e-03
* time: 0.41290998458862305
17 -7.852629e+00 1.890666e-03
* time: 0.4605081081390381
18 -7.852642e+00 1.009485e-03
* time: 0.48296213150024414
19 -7.852645e+00 5.664085e-04
* time: 0.5055961608886719
20 -7.852646e+00 3.366794e-04
* time: 0.5278310775756836
21 -7.852646e+00 2.018730e-04
* time: 0.5503361225128174
22 -7.852647e+00 1.489106e-04
* time: 0.5725641250610352
23 -7.852647e+00 8.637588e-05
* time: 0.5950040817260742
24 -7.852647e+00 5.915399e-05
* time: 0.617603063583374
25 -7.852647e+00 2.230207e-05
* time: 0.6397960186004639
26 -7.852647e+00 1.156498e-05
* time: 0.662121057510376
27 -7.852647e+00 8.198795e-06
* time: 0.6845271587371826
28 -7.852647e+00 5.502192e-06
* time: 0.7071561813354492
29 -7.852647e+00 2.023899e-06
* time: 0.7294771671295166
30 -7.852647e+00 1.449614e-06
* time: 0.7516160011291504
31 -7.852647e+00 7.676239e-07
* time: 0.7739729881286621
32 -7.852647e+00 4.320281e-07
* time: 0.7965030670166016
33 -7.852647e+00 2.713169e-07
* time: 0.8190250396728516
34 -7.852647e+00 1.840549e-07
* time: 0.841616153717041
35 -7.852647e+00 1.400749e-07
* time: 0.8643391132354736
36 -7.852647e+00 7.835382e-08
* time: 0.8871190547943115
37 -7.852647e+00 5.224560e-08
* time: 0.9097461700439453
38 -7.852647e+00 3.095855e-08
* time: 0.9320921897888184
39 -7.852647e+00 1.739628e-08
* time: 0.9643402099609375
40 -7.852647e+00 1.589284e-08
* time: 1.0073161125183105
41 -7.852647e+00 1.025110e-08
* time: 1.0300660133361816
42 -7.852647e+00 5.285366e-09
* time: 1.0526010990142822
43 -7.852647e+00 4.828327e-09
* time: 1.084954023361206Newton 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.846885606838 -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.852645921498 -1.64
2 -7.852646686730 -6.12 -3.71 392ms
3 -7.852646686730 -13.34 -7.25 141msComparison 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| = 5.18209311518797e-7
|ρ_newton - ρ_scfv| = 3.6672919984490966e-7
|ρ_newton - ρ_dm| = 7.897882526685183e-10