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 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
Density-based self-consistent field
scfres_scf = self_consistent_field(basis; tol);
n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.397863751843 -0.90 5.2 26.2ms
2 -8.400238836002 -2.62 -1.74 1.0 24.9ms
3 -8.400405164893 -3.78 -2.98 1.5 19.8ms
4 -8.400427857909 -4.64 -2.99 3.0 23.9ms
5 -8.400427921932 -7.19 -3.02 1.0 18.9ms
6 -8.400428144057 -6.65 -4.55 1.0 23.8ms
7 -8.400428151597 -8.12 -4.36 2.8 23.5ms
8 -8.400428152173 -9.24 -5.10 1.0 19.3ms
9 -8.400428152207 -10.47 -6.24 1.2 20.0ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397796881798 -0.90 5.5 28.2ms
2 -8.400384894315 -2.59 -1.79 0.80 2.2 19.1ms
3 -8.400422730030 -4.42 -3.01 0.80 1.0 16.1ms
4 -8.400428106113 -5.27 -3.43 0.80 2.2 24.0ms
5 -8.400428149660 -7.36 -4.62 0.80 1.2 16.8ms
6 -8.400428152185 -8.60 -5.70 0.80 2.5 20.0ms
7 -8.400428152209 -10.62 -6.00 0.80 2.5 20.1ms
8 -8.400428152209 -12.59 -7.10 0.80 1.0 16.3ms
Direct minimization
scfres_dm = direct_minimization(basis; tol);
┌ Warning: x_tol is deprecated. Use x_abstol or x_reltol instead. The provided value (-1) will be used as x_abstol.
└ @ Optim ~/.julia/packages/Optim/7krni/src/types.jl:110
┌ Warning: f_tol is deprecated. Use f_abstol or f_reltol instead. The provided value (-1) will be used as f_reltol.
└ @ Optim ~/.julia/packages/Optim/7krni/src/types.jl:120
n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.051648144054 -1.09 61.7ms
2 -1.859552167927 0.46 -0.71 32.7ms
3 -4.552250590873 0.43 -0.47 44.6ms
4 -6.172299567834 0.21 -0.54 47.6ms
5 -7.604966352285 0.16 -0.83 44.2ms
6 -7.968828670400 -0.44 -1.34 32.4ms
7 -8.195406774940 -0.64 -1.55 36.3ms
8 -8.257493620890 -1.21 -1.68 32.9ms
9 -8.328531831544 -1.15 -1.90 33.0ms
10 -8.359099405008 -1.51 -2.11 36.5ms
11 -8.384647230904 -1.59 -2.53 32.5ms
12 -8.393349879024 -2.06 -2.59 32.6ms
13 -8.397091360335 -2.43 -3.00 32.5ms
14 -8.398964668611 -2.73 -3.29 35.9ms
15 -8.399782260021 -3.09 -3.41 32.4ms
16 -8.400059388742 -3.56 -3.20 32.4ms
17 -8.400276814040 -3.66 -3.33 32.4ms
18 -8.400353446940 -4.12 -3.92 36.0ms
19 -8.400396969208 -4.36 -3.80 32.5ms
20 -8.400418357490 -4.67 -4.21 32.9ms
21 -8.400424594455 -5.21 -4.10 32.4ms
22 -8.400426796142 -5.66 -4.37 36.2ms
23 -8.400427691398 -6.05 -4.52 32.8ms
24 -8.400427947513 -6.59 -5.35 32.4ms
25 -8.400428070850 -6.91 -5.17 36.3ms
26 -8.400428115764 -7.35 -5.50 32.8ms
27 -8.400428136957 -7.67 -5.74 32.9ms
28 -8.400428146309 -8.03 -5.82 32.9ms
29 -8.400428149795 -8.46 -5.82 36.2ms
30 -8.400428151043 -8.90 -6.22 32.7ms
Newton 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 -8.397910971377 -0.90 5.0 25.4ms
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.400427991194 -1.79 575ms
2 -8.400428152209 -6.79 -4.05 382ms
3 -8.400428152209 + -Inf -7.86 120ms
Comparison 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| = 1.3828366765974826e-6
|ρ_newton - ρ_scfv| = 1.951945286376085e-7
|ρ_newton - ρ_dm| = 1.691272756887701e-6