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 AtomsBuilder
using DFTK
using LinearAlgebra
using PseudoPotentialData
pseudopotentials = PseudoFamily("dojo.nc.sr.pbesol.v0_4_1.standard.upf")
model = model_DFT(bulk(:Si); functionals=PBEsol(), pseudopotentials)
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 -8.397850398111 -0.90 5.2 27.0ms 2 -8.400256071796 -2.62 -1.74 1.0 18.7ms 3 -8.400406663600 -3.82 -2.96 1.5 19.8ms 4 -8.400427833193 -4.67 -2.98 3.0 57.9ms 5 -8.400427949655 -6.93 -3.05 1.0 18.9ms 6 -8.400428149316 -6.70 -4.70 1.0 19.0ms 7 -8.400428155862 -8.18 -4.45 3.2 26.0ms 8 -8.400428156256 -9.40 -5.21 1.0 19.6ms 9 -8.400428156275 -10.72 -6.13 1.0 19.7ms
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime --- --------------- --------- --------- ---- ---- ------ 1 -8.397726925764 -0.90 4.8 553ms 2 -8.400374639656 -2.58 -1.77 0.80 2.0 495ms 3 -8.400422170514 -4.32 -3.00 0.80 1.0 180ms 4 -8.400428104597 -5.23 -3.46 0.80 2.2 19.9ms 5 -8.400428151685 -7.33 -4.97 0.80 1.2 17.8ms 6 -8.400428156269 -8.34 -5.35 0.80 3.0 23.0ms 7 -8.400428156276 -11.12 -6.86 0.80 1.0 16.4ms
scfres_dm = direct_minimization(basis; tol);
n Energy log10(ΔE) log10(Δρ) Δtime --- --------------- --------- --------- ------ 1 +0.926488639176 -1.09 2.77s 2 -1.848359491024 0.44 -0.65 106ms 3 -4.582656986494 0.44 -0.39 39.3ms 4 -6.169239948250 0.20 -0.47 73.3ms 5 -7.680944399394 0.18 -0.73 39.6ms 6 -7.933004251426 -0.60 -1.45 29.2ms 7 -8.218204309711 -0.54 -1.62 29.0ms 8 -8.227933182068 -2.01 -1.94 28.8ms 9 -8.275254161630 -1.32 -1.79 49.0ms 10 -8.319166652806 -1.36 -1.67 39.4ms 11 -8.353780606343 -1.46 -1.75 39.5ms 12 -8.374376330935 -1.69 -2.43 29.0ms 13 -8.389239941286 -1.83 -2.40 33.4ms 14 -8.396952636531 -2.11 -2.50 29.1ms 15 -8.398593985646 -2.78 -2.59 29.0ms 16 -8.399631212256 -2.98 -3.11 29.0ms 17 -8.400022884144 -3.41 -3.11 32.3ms 18 -8.400262490188 -3.62 -3.37 29.0ms 19 -8.400353215508 -4.04 -3.59 29.1ms 20 -8.400407239636 -4.27 -4.05 32.2ms 21 -8.400417995679 -4.97 -3.82 28.9ms 22 -8.400425286419 -5.14 -4.36 29.0ms 23 -8.400426555823 -5.90 -4.17 32.1ms 24 -8.400427343278 -6.10 -4.71 29.1ms 25 -8.400427745067 -6.40 -4.57 32.0ms 26 -8.400427971511 -6.65 -5.01 29.0ms 27 -8.400428076202 -6.98 -5.06 29.5ms 28 -8.400428128024 -7.29 -5.21 33.7ms 29 -8.400428141728 -7.86 -5.46 29.0ms 30 -8.400428151034 -8.03 -5.36 29.0ms 31 -8.400428154120 -8.51 -5.67 32.2ms 32 -8.400428155458 -8.87 -5.85 29.0ms 33 -8.400428155871 -9.38 -6.61 32.2ms
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 -8.397856585042 -0.90 5.0 26.7ms
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 -8.400427983085 -1.78 24.2s 2 -8.400428156277 -6.76 -4.02 2.19s 3 -8.400428156277 -14.75 -7.82 115ms
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| = 1.318651071213495e-6 |ρ_newton - ρ_scfv| = 1.2366889710809778e-7 |ρ_newton - ρ_dm| = 1.5628380220615619e-6