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.846835863630 -0.70 4.8
2 -7.852323709791 -2.26 -1.53 1.0 34.6ms
3 -7.852612797784 -3.54 -2.55 1.2 37.2ms
4 -7.852645917000 -4.48 -2.86 2.8 49.7ms
5 -7.852646445638 -6.28 -3.10 1.0 34.8ms
6 -7.852646677080 -6.64 -3.95 1.0 36.2ms
7 -7.852646686241 -8.04 -5.13 1.5 42.6ms
8 -7.852646686715 -9.32 -5.23 2.5 48.4ms
9 -7.852646686729 -10.85 -6.20 1.0 35.1msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -7.846839598734 -0.70 4.5
2 -7.852528940588 -2.24 -1.63 0.80 2.0 38.8ms
3 -7.852637887565 -3.96 -2.70 0.80 1.0 32.2ms
4 -7.852646479322 -5.07 -3.31 0.80 2.0 46.6ms
5 -7.852646682881 -6.69 -4.20 0.80 1.5 36.4ms
6 -7.852646686602 -8.43 -4.87 0.80 2.0 45.8ms
7 -7.852646686721 -9.93 -5.77 0.80 1.5 34.5ms
8 -7.852646686730 -11.04 -6.79 0.80 2.0 38.9msDirect 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.377235e+01 3.881346e+00
* time: 0.01699209213256836
1 1.265511e+00 1.801907e+00
* time: 0.049504995346069336
2 -1.543474e+00 1.960013e+00
* time: 0.08228492736816406
3 -3.858899e+00 1.655478e+00
* time: 0.13008499145507812
4 -5.175275e+00 1.663854e+00
* time: 0.1777019500732422
5 -6.857095e+00 1.206923e+00
* time: 0.22472691535949707
6 -7.446812e+00 4.656910e-01
* time: 0.2716360092163086
7 -7.665155e+00 1.925371e-01
* time: 0.3676130771636963
8 -7.752370e+00 1.642304e-01
* time: 0.40064191818237305
9 -7.799010e+00 1.327639e-01
* time: 0.4340860843658447
10 -7.824560e+00 1.045189e-01
* time: 0.46753406524658203
11 -7.839799e+00 6.699483e-02
* time: 0.5019879341125488
12 -7.845236e+00 6.118431e-02
* time: 0.5358829498291016
13 -7.846562e+00 4.216923e-02
* time: 0.5742161273956299
14 -7.848235e+00 3.144008e-02
* time: 0.6114299297332764
15 -7.850677e+00 2.142225e-02
* time: 0.6461060047149658
16 -7.851997e+00 1.781180e-02
* time: 0.6796009540557861
17 -7.852452e+00 6.897497e-03
* time: 0.714224100112915
18 -7.852589e+00 3.823951e-03
* time: 0.7476780414581299
19 -7.852628e+00 2.867573e-03
* time: 0.7834999561309814
20 -7.852641e+00 1.471738e-03
* time: 0.818295955657959
21 -7.852645e+00 6.774609e-04
* time: 0.8513889312744141
22 -7.852646e+00 5.449233e-04
* time: 0.8840620517730713
23 -7.852647e+00 2.691163e-04
* time: 0.9181230068206787
24 -7.852647e+00 1.736887e-04
* time: 0.9518449306488037
25 -7.852647e+00 9.470015e-05
* time: 0.9836061000823975
26 -7.852647e+00 5.422476e-05
* time: 1.0163159370422363
27 -7.852647e+00 5.237609e-05
* time: 1.048586130142212
28 -7.852647e+00 2.292737e-05
* time: 1.0814039707183838
29 -7.852647e+00 1.119824e-05
* time: 1.1156079769134521
30 -7.852647e+00 5.171177e-06
* time: 1.1503491401672363
31 -7.852647e+00 3.022153e-06
* time: 1.183779001235962
32 -7.852647e+00 2.204431e-06
* time: 1.2175700664520264
33 -7.852647e+00 1.337301e-06
* time: 1.2507951259613037
34 -7.852647e+00 8.326351e-07
* time: 1.283607006072998
35 -7.852647e+00 5.370798e-07
* time: 1.3170990943908691
36 -7.852647e+00 2.981644e-07
* time: 1.351619005203247
37 -7.852647e+00 2.108745e-07
* time: 1.385200023651123
38 -7.852647e+00 1.225919e-07
* time: 1.4177899360656738
39 -7.852647e+00 6.459460e-08
* time: 1.4510281085968018
40 -7.852647e+00 2.919975e-08
* time: 1.4831700325012207
41 -7.852647e+00 1.092091e-08
* time: 1.516571044921875
42 -7.852647e+00 8.208210e-09
* time: 1.5494189262390137
43 -7.852647e+00 8.208202e-09
* time: 1.6864359378814697Newton 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.846746081609 -0.70 4.5Remove 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.852645862126 -1.64
2 -7.852646686730 -6.08 -3.70 695ms
3 -7.852646686730 -13.25 -7.22 261msComparison 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| = 8.951759867264316e-7
|ρ_newton - ρ_scfv| = 1.5880727604563485e-7
|ρ_newton - ρ_dm| = 1.015388878899405e-9