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.397811886237 -0.90 5.5 26.9ms
2 -8.400253533221 -2.61 -1.74 1.0 20.7ms
3 -8.400397482785 -3.84 -2.91 1.2 19.6ms
4 -8.400427779469 -4.52 -2.92 3.0 31.8ms
5 -8.400427674702 + -6.98 -2.84 1.0 19.1ms
6 -8.400428142288 -6.33 -4.31 1.0 19.1ms
7 -8.400428155176 -7.89 -4.24 2.2 22.5ms
8 -8.400428156231 -8.98 -5.07 1.0 24.0ms
9 -8.400428156273 -10.37 -6.24 1.2 19.7ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397799737108 -0.90 5.2 30.8ms
2 -8.400385248102 -2.59 -1.79 0.80 2.2 19.4ms
3 -8.400422774312 -4.43 -3.05 0.80 1.0 16.4ms
4 -8.400428110320 -5.27 -3.40 0.80 2.5 20.2ms
5 -8.400428152968 -7.37 -4.45 0.80 1.0 16.5ms
6 -8.400428156189 -8.49 -5.61 0.80 2.2 19.7ms
7 -8.400428156277 -10.06 -6.13 0.80 3.2 28.1ms
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.015142122654 -1.09 58.4ms
2 -1.496222509640 0.40 -0.69 32.8ms
3 -4.291963730684 0.45 -0.42 47.9ms
4 -5.669940207943 0.14 -0.53 44.1ms
5 -7.338045367684 0.22 -0.68 44.5ms
6 -7.497436872008 -0.80 -1.41 36.4ms
7 -8.021504091459 -0.28 -1.55 32.8ms
8 -8.122686219232 -0.99 -1.79 32.6ms
9 -8.215991592012 -1.03 -1.68 32.6ms
10 -8.237914579347 -1.66 -1.89 36.1ms
11 -8.255931323429 -1.74 -1.83 32.7ms
12 -8.293781212957 -1.42 -1.64 44.4ms
13 -8.310198098145 -1.78 -1.63 47.8ms
14 -8.352927358479 -1.37 -1.85 44.3ms
15 -8.385001005696 -1.49 -2.22 33.3ms
16 -8.391529762152 -2.19 -2.30 33.7ms
17 -8.396870965654 -2.27 -2.91 36.8ms
18 -8.398407530843 -2.81 -3.08 33.1ms
19 -8.399541194369 -2.95 -2.99 33.5ms
20 -8.400091474958 -3.26 -3.37 32.9ms
21 -8.400248621688 -3.80 -3.47 36.7ms
22 -8.400358543643 -3.96 -3.74 33.0ms
23 -8.400395850582 -4.43 -3.60 33.2ms
24 -8.400417587224 -4.66 -3.92 36.1ms
25 -8.400423759293 -5.21 -4.16 33.2ms
26 -8.400426225927 -5.61 -4.31 32.7ms
27 -8.400427319521 -5.96 -4.62 33.0ms
28 -8.400427745308 -6.37 -5.02 36.2ms
29 -8.400427997733 -6.60 -4.83 33.3ms
30 -8.400428091321 -7.03 -5.27 32.8ms
31 -8.400428131338 -7.40 -5.16 32.9ms
32 -8.400428146928 -7.81 -5.76 36.3ms
33 -8.400428152179 -8.28 -5.51 33.4ms
34 -8.400428154666 -8.60 -6.50 32.9ms
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.397823981485 -0.90 5.0 25.6ms
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.400427995172 -1.79 575ms
2 -8.400428156277 -6.79 -4.04 380ms
3 -8.400428156277 -14.45 -7.87 103ms
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| = 9.079478852921349e-7
|ρ_newton - ρ_scfv| = 9.047653528040841e-8
|ρ_newton - ρ_dm| = 2.061922894339901e-6