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-6
1.0e-6
scfres_scf = self_consistent_field(basis; tol);
n Energy log10(ΔE) log10(Δρ) Diag Δtime --- --------------- --------- --------- ---- ------ 1 -7.847758929188 -0.70 4.5 82.7ms 2 -7.852409749928 -2.33 -1.53 1.0 17.2ms 3 -7.852609850338 -3.70 -2.54 1.0 16.9ms 4 -7.852645354134 -4.45 -2.77 2.2 21.7ms 5 -7.852646044764 -6.16 -2.85 1.0 17.0ms 6 -7.852646669333 -6.20 -3.99 1.0 17.3ms 7 -7.852646686117 -7.78 -4.46 1.8 19.9ms 8 -7.852646686699 -9.24 -5.48 1.0 17.1ms 9 -7.852646686724 -10.60 -5.42 2.2 21.5ms 10 -7.852646686729 -11.28 -5.95 1.0 17.7ms 11 -7.852646686730 -12.29 -6.63 1.0 17.2ms
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -7.847764577614 -0.70 4.8 398ms 2 -7.852560871393 -2.32 -1.62 0.80 2.0 2.22s 3 -7.852640927769 -4.10 -2.71 0.80 1.0 225ms 4 -7.852646512514 -5.25 -3.39 0.80 2.0 19.5ms 5 -7.852646681566 -6.77 -4.37 0.80 1.5 17.8ms 6 -7.852646686616 -8.30 -4.86 0.80 2.0 49.2ms 7 -7.852646686726 -9.96 -5.51 0.80 1.2 16.8ms 8 -7.852646686730 -11.42 -6.66 0.80 1.8 18.2ms
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);
n Energy log10(ΔE) log10(Δρ) Δtime --- --------------- --------- --------- ------ 1 +1.586130217025 -1.03 4.17s 2 -1.354767699875 0.47 -0.64 108ms 3 -3.667485313572 0.36 -0.46 32.5ms 4 -4.867274699357 0.08 -0.61 32.7ms 5 -6.576127400939 0.23 -0.74 32.6ms 6 -7.230048187789 -0.18 -1.18 32.7ms 7 -7.520932597211 -0.54 -1.41 24.6ms 8 -7.664884680787 -0.84 -1.73 24.9ms 9 -7.751456184188 -1.06 -1.77 24.8ms 10 -7.787663944145 -1.44 -2.13 24.9ms 11 -7.818000638126 -1.52 -2.00 24.9ms 12 -7.835546794003 -1.76 -2.06 24.8ms 13 -7.843321770235 -2.11 -2.60 24.9ms 14 -7.848252092295 -2.31 -2.64 24.8ms 15 -7.850786586997 -2.60 -3.33 25.5ms 16 -7.851987825813 -2.92 -3.08 25.1ms 17 -7.852431755245 -3.35 -3.64 24.8ms 18 -7.852568146994 -3.87 -3.74 25.2ms 19 -7.852615121311 -4.33 -3.85 25.2ms 20 -7.852634215625 -4.72 -4.03 24.9ms 21 -7.852641946230 -5.11 -4.32 24.6ms 22 -7.852645046595 -5.51 -4.54 56.7ms 23 -7.852646197929 -5.94 -5.09 24.6ms 24 -7.852646555110 -6.45 -5.46 24.7ms 25 -7.852646649767 -7.02 -5.38 24.7ms 26 -7.852646673608 -7.62 -5.82 24.7ms 27 -7.852646680989 -8.13 -5.88 24.5ms 28 -7.852646684459 -8.46 -6.00 24.6ms 29 -7.852646686066 -8.79 -6.21 24.6ms 30 -7.852646686518 -9.34 -6.69 24.5ms 31 -7.852646686652 -9.88 -6.75 24.6ms 32 -7.852646686703 -10.29 -6.89 25.2ms 33 -7.852646686721 -10.74 -7.48 24.7ms 34 -7.852646686727 -11.24 -7.51 24.6ms 35 -7.852646686729 -11.76 -7.75 24.6ms 36 -7.852646686730 -12.26 -8.36 24.5ms 37 -7.852646686730 -12.71 -8.34 24.6ms 38 -7.852646686730 -13.16 -8.99 24.6ms 39 -7.852646686730 -13.60 -9.18 24.6ms 40 -7.852646686730 -14.01 -8.97 24.6ms 41 -7.852646686730 -15.05 -9.72 24.5ms 42 -7.852646686730 -14.75 -9.61 24.6ms 43 -7.852646686730 + -Inf -9.78 24.5ms 44 -7.852646686730 + -Inf -10.34 72.4ms ┌ Warning: DM not converged. └ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/scf_callbacks.jl:60
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.847767780617 -0.70 4.8 29.8ms
Remove 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.852645859423 -1.63 13.7s 2 -7.852646686730 -6.08 -3.69 3.60s 3 -7.852646686730 -13.26 -7.21 105ms
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| = 2.573390010523403e-7 |ρ_newton - ρ_scfv| = 4.861905141684168e-7 |ρ_newton - ρ_dm| = 8.137971356684988e-10