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.846872559925 -0.70 5.0 2 -7.852329074596 -2.26 -1.53 1.0 22.6ms 3 -7.852613512883 -3.55 -2.56 1.2 23.9ms 4 -7.852645956388 -4.49 -2.87 2.8 31.2ms 5 -7.852646464554 -6.29 -3.12 1.0 23.2ms 6 -7.852646676639 -6.67 -3.94 1.0 23.0ms 7 -7.852646686168 -8.02 -5.06 1.5 25.3ms 8 -7.852646686721 -9.26 -5.37 2.5 31.6ms 9 -7.852646686729 -11.09 -6.00 1.0 65.5ms
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -7.846862056955 -0.70 4.8 2 -7.852524744014 -2.25 -1.64 0.80 2.2 258ms 3 -7.852635224137 -3.96 -2.73 0.80 1.0 21.3ms 4 -7.852646499059 -4.95 -3.23 0.80 2.2 29.5ms 5 -7.852646673992 -6.76 -4.05 0.80 1.2 22.4ms 6 -7.852646686371 -7.91 -4.81 0.80 1.8 26.1ms 7 -7.852646686724 -9.45 -5.57 0.80 1.8 26.2ms 8 -7.852646686730 -11.24 -7.03 0.80 2.0 28.3ms
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.376673e+01 3.446661e+00 * time: 0.05784797668457031 1 1.035668e+00 1.764910e+00 * time: 0.27242517471313477 2 -1.780741e+00 2.153443e+00 * time: 0.2955482006072998 3 -3.884170e+00 1.924643e+00 * time: 0.32851099967956543 4 -5.096888e+00 1.758705e+00 * time: 0.3610820770263672 5 -6.838055e+00 9.559709e-01 * time: 0.3939549922943115 6 -7.465364e+00 5.330211e-01 * time: 0.4918639659881592 7 -7.681976e+00 3.708038e-01 * time: 0.5148661136627197 8 -7.757692e+00 1.084613e-01 * time: 0.5376181602478027 9 -7.806085e+00 1.598147e-01 * time: 0.5601329803466797 10 -7.829592e+00 6.875602e-02 * time: 0.5826570987701416 11 -7.841884e+00 6.256555e-02 * time: 0.6049520969390869 12 -7.846960e+00 5.524165e-02 * time: 0.6271171569824219 13 -7.848877e+00 4.607307e-02 * time: 0.6494541168212891 14 -7.850422e+00 2.534032e-02 * time: 0.6719121932983398 15 -7.851783e+00 1.470556e-02 * time: 0.6945559978485107 16 -7.852414e+00 7.602897e-03 * time: 0.7171990871429443 17 -7.852586e+00 5.228225e-03 * time: 0.7398731708526611 18 -7.852628e+00 2.276126e-03 * time: 0.7624790668487549 19 -7.852641e+00 9.418454e-04 * time: 0.7851331233978271 20 -7.852645e+00 6.488569e-04 * time: 0.8080360889434814 21 -7.852646e+00 3.455321e-04 * time: 0.830517053604126 22 -7.852646e+00 2.399337e-04 * time: 0.853262186050415 23 -7.852647e+00 1.632441e-04 * time: 0.8759369850158691 24 -7.852647e+00 1.049781e-04 * time: 0.898604154586792 25 -7.852647e+00 4.875373e-05 * time: 0.9212710857391357 26 -7.852647e+00 2.425934e-05 * time: 0.9436380863189697 27 -7.852647e+00 1.430669e-05 * time: 0.9660019874572754 28 -7.852647e+00 8.628782e-06 * time: 0.9886350631713867 29 -7.852647e+00 5.037319e-06 * time: 1.0112099647521973 30 -7.852647e+00 3.506866e-06 * time: 1.0340790748596191 31 -7.852647e+00 1.753959e-06 * time: 1.0567450523376465 32 -7.852647e+00 1.271976e-06 * time: 1.0790660381317139 33 -7.852647e+00 5.539523e-07 * time: 1.1015121936798096 34 -7.852647e+00 3.544319e-07 * time: 1.1238820552825928 35 -7.852647e+00 2.011782e-07 * time: 1.14650297164917 36 -7.852647e+00 1.115273e-07 * time: 1.1691629886627197 37 -7.852647e+00 8.662922e-08 * time: 1.191802978515625 38 -7.852647e+00 3.329043e-08 * time: 1.2141931056976318 39 -7.852647e+00 2.182318e-08 * time: 1.2366251945495605 40 -7.852647e+00 1.549553e-08 * time: 1.2692999839782715 41 -7.852647e+00 6.963544e-09 * time: 1.2917990684509277 42 -7.852647e+00 6.493402e-09 * time: 1.3149960041046143 43 -7.852647e+00 2.580704e-09 * time: 1.337494134902954
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.846853629199 -0.70 4.5
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.852645914146 -1.64 2 -7.852646686730 -6.11 -3.71 1.90s 3 -7.852646686730 -13.29 -7.25 140ms
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.327201480445901e-7 |ρ_newton - ρ_scfv| = 3.5148077303630417e-7 |ρ_newton - ρ_dm| = 7.169235209148582e-10