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
tol = 1e-12
is_converged = DFTK.ScfConvergenceDensity(tol);
scfres_scf = self_consistent_field(basis; is_converged);
n Energy log10(ΔE) log10(Δρ) Diag --- --------------- --------- --------- ---- 1 -7.846832249821 -0.70 4.8 2 -7.852321060476 -2.26 -1.53 1.0 3 -7.852646061824 -3.49 -2.52 3.2 4 -7.852646678176 -6.21 -3.37 2.5 5 -7.852646686299 -8.09 -4.73 1.5 6 -7.852646686727 -9.37 -5.44 3.5 7 -7.852646686730 -11.54 -6.27 1.2 8 -7.852646686730 -13.01 -7.71 2.2 9 -7.852646686730 -14.57 -7.74 2.2 10 -7.852646686730 -15.05 -8.71 2.5 11 -7.852646686730 + -Inf -8.94 1.0 12 -7.852646686730 + -15.05 -10.67 2.2 13 -7.852646686730 + -Inf -9.47 1.0 14 -7.852646686730 -15.05 -9.24 1.0 15 -7.852646686730 -14.75 -9.10 1.0 16 -7.852646686730 + -14.35 -8.50 2.0 17 -7.852646686730 + -15.05 -9.80 1.5 18 -7.852646686730 -14.75 -10.57 1.8 19 -7.852646686730 -14.75 -10.95 1.8 20 -7.852646686730 + -15.05 -10.82 1.0 21 -7.852646686730 + -15.05 -11.71 1.0 22 -7.852646686730 -14.75 -12.45 2.0
scfres_scfv = DFTK.scf_potential_mixing(basis; is_converged);
n Energy log10(ΔE) log10(Δρ) α Diag --- --------------- --------- --------- ---- ---- 1 -7.846826529307 -0.70 4.8 2 -7.852525112695 -2.24 -1.63 0.80 2.0 3 -7.852635705254 -3.96 -2.72 0.80 1.0 4 -7.852646446816 -4.97 -3.27 0.80 2.2 5 -7.852646672886 -6.65 -4.10 0.80 1.2 6 -7.852646686341 -7.87 -4.72 0.80 1.5 7 -7.852646686721 -9.42 -5.57 0.80 1.8 8 -7.852646686729 -11.08 -6.98 0.80 1.8 9 -7.852646686730 -12.41 -7.57 0.80 3.2 10 -7.852646686730 -14.75 -7.99 0.80 1.8 11 -7.852646686730 -15.05 -9.12 0.80 1.0 12 -7.852646686730 -14.75 -9.60 0.80 3.2 13 -7.852646686730 + -15.05 -10.10 0.80 1.0 14 -7.852646686730 + -15.05 -10.82 0.80 1.2 15 -7.852646686730 + -Inf -11.97 0.80 2.2 16 -7.852646686730 + -Inf -12.64 0.80 2.5
scfres_dm = direct_minimization(basis; tol);
Iter Function value Gradient norm 0 1.403903e+01 3.503994e+00 * time: 0.44902586936950684 1 1.342235e+00 1.778317e+00 * time: 0.651695966720581 2 -1.828507e+00 1.896481e+00 * time: 0.6792130470275879 3 -3.754821e+00 1.820452e+00 * time: 0.7940008640289307 4 -5.188438e+00 1.584808e+00 * time: 0.8335838317871094 5 -6.918980e+00 8.684053e-01 * time: 0.8725728988647461 6 -6.992047e+00 9.572685e-01 * time: 0.899763822555542 7 -7.689852e+00 7.397375e-01 * time: 0.9269559383392334 8 -7.759766e+00 8.800055e-01 * time: 0.9543449878692627 9 -7.778897e+00 7.374077e-01 * time: 0.9812648296356201 10 -7.785095e+00 1.073822e+00 * time: 1.0081868171691895 11 -7.798981e+00 8.160733e-01 * time: 1.0350289344787598 12 -7.819909e+00 4.539929e-01 * time: 1.0741620063781738 13 -7.824708e+00 6.909547e-01 * time: 1.1010818481445312 14 -7.843710e+00 2.596031e-01 * time: 1.1279468536376953 15 -7.846046e+00 2.959001e-01 * time: 1.15488600730896 16 -7.851371e+00 2.329684e-02 * time: 1.1817529201507568 17 -7.852230e+00 8.761840e-03 * time: 1.2086799144744873 18 -7.852569e+00 6.798428e-03 * time: 1.23555588722229 19 -7.852604e+00 4.930048e-03 * time: 1.262833833694458 20 -7.852625e+00 2.813465e-03 * time: 1.289719820022583 21 -7.852641e+00 1.634270e-03 * time: 1.3174378871917725 22 -7.852645e+00 9.289906e-04 * time: 1.3442578315734863 23 -7.852646e+00 4.364427e-04 * time: 1.3711888790130615 24 -7.852646e+00 3.232755e-04 * time: 1.3981859683990479 25 -7.852647e+00 1.153466e-04 * time: 1.4249730110168457 26 -7.852647e+00 6.613080e-05 * time: 1.451746940612793 27 -7.852647e+00 4.115242e-05 * time: 1.4787218570709229 28 -7.852647e+00 2.545155e-05 * time: 1.505707025527954 29 -7.852647e+00 1.371649e-05 * time: 1.5330719947814941 30 -7.852647e+00 8.255275e-06 * time: 1.5598359107971191 31 -7.852647e+00 3.840623e-06 * time: 1.586789846420288 32 -7.852647e+00 2.238605e-06 * time: 1.6139118671417236 33 -7.852647e+00 1.239133e-06 * time: 1.6414918899536133 34 -7.852647e+00 6.781613e-07 * time: 1.6685888767242432 35 -7.852647e+00 3.939354e-07 * time: 1.6955578327178955 36 -7.852647e+00 2.162569e-07 * time: 1.7229039669036865 37 -7.852647e+00 1.289076e-07 * time: 1.7498118877410889 38 -7.852647e+00 6.530468e-08 * time: 1.7768030166625977 39 -7.852647e+00 3.634828e-08 * time: 1.8040690422058105 40 -7.852647e+00 2.074358e-08 * time: 1.831233024597168 41 -7.852647e+00 1.169768e-08 * time: 1.858158826828003 42 -7.852647e+00 7.731757e-09 * time: 1.8854949474334717 43 -7.852647e+00 4.208548e-09 * time: 1.913322925567627 44 -7.852647e+00 4.207983e-09 * time: 1.988473892211914
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=1e-1);
n Energy log10(ΔE) log10(Δρ) Diag --- --------------- --------- --------- ---- 1 -7.846843511052 -0.70 4.5 2 -7.852322679588 -2.26 -1.53 1.0
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(Δρ) --- --------------- --------- --------- 1 -7.852646686713 -2.54 2 -7.852646686730 -10.77 -6.03 3 -7.852646686730 + -15.05 -12.64
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| = 9.899854352494933e-14 |ρ_newton - ρ_scfv| = 2.1712353697319453e-13 |ρ_newton - ρ_dm| = 2.420188489527844e-10