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.397871040479 -0.90 4.8 25.7ms
2 -8.400249627308 -2.62 -1.73 1.0 18.9ms
3 -8.400403760156 -3.81 -2.95 1.5 19.6ms
4 -8.400427784361 -4.62 -2.92 3.2 33.4ms
5 -8.400427978138 -6.71 -3.07 1.0 19.3ms
6 -8.400428150548 -6.76 -4.93 1.0 19.5ms
7 -8.400428155958 -8.27 -4.50 3.5 25.8ms
8 -8.400428156266 -9.51 -5.86 1.0 19.6ms
9 -8.400428156277 -10.95 -6.25 2.0 22.4ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397784177301 -0.90 4.5 25.5ms
2 -8.400379242924 -2.59 -1.77 0.80 2.2 19.3ms
3 -8.400423258062 -4.36 -2.96 0.80 1.0 16.4ms
4 -8.400428112097 -5.31 -3.48 0.80 2.2 26.0ms
5 -8.400428152548 -7.39 -5.00 0.80 1.2 17.2ms
6 -8.400428156273 -8.43 -5.35 0.80 3.8 23.4ms
7 -8.400428156276 -11.44 -6.82 0.80 1.0 16.9ms
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 +0.811579736385 -1.05 57.7ms
2 -2.175644224432 0.48 -0.70 33.1ms
3 -4.651360318919 0.39 -0.40 44.6ms
4 -6.331133238248 0.23 -0.50 49.7ms
5 -7.761514758414 0.16 -0.76 44.8ms
6 -8.089835691822 -0.48 -1.28 33.2ms
7 -8.246842081362 -0.80 -1.71 32.8ms
8 -8.274081976884 -1.56 -1.96 32.5ms
9 -8.293563757163 -1.71 -1.82 32.8ms
10 -8.354497393615 -1.22 -1.66 49.7ms
11 -8.375125649007 -1.69 -1.90 44.8ms
12 -8.389642704499 -1.84 -2.66 33.3ms
13 -8.394791362619 -2.29 -2.82 33.3ms
14 -8.398777967727 -2.40 -2.96 38.1ms
15 -8.399642347400 -3.06 -3.12 33.4ms
16 -8.400048043555 -3.39 -3.33 33.0ms
17 -8.400263254949 -3.67 -3.65 33.0ms
18 -8.400337650572 -4.13 -3.65 36.9ms
19 -8.400396809478 -4.23 -3.91 33.3ms
20 -8.400412494755 -4.80 -3.88 33.1ms
21 -8.400423803973 -4.95 -4.08 32.9ms
22 -8.400426241651 -5.61 -4.13 36.3ms
23 -8.400427313622 -5.97 -4.56 33.1ms
24 -8.400427713565 -6.40 -4.70 32.7ms
25 -8.400427964199 -6.60 -4.88 36.6ms
26 -8.400428079449 -6.94 -4.90 33.3ms
27 -8.400428128607 -7.31 -5.10 32.9ms
28 -8.400428144307 -7.80 -5.27 33.0ms
29 -8.400428150495 -8.21 -5.56 36.5ms
30 -8.400428153936 -8.46 -5.83 33.3ms
31 -8.400428155081 -8.94 -5.90 33.2ms
32 -8.400428155760 -9.17 -6.41 37.2ms
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.397842791157 -0.90 5.0 27.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.400427994786 -1.79 601ms
2 -8.400428156277 -6.79 -4.05 405ms
3 -8.400428156277 -14.75 -7.86 105ms
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.664401693530924e-7
|ρ_newton - ρ_scfv| = 2.8974785243296963e-7
|ρ_newton - ρ_dm| = 1.3665981130416105e-6