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.397845319135 -0.90 5.2 26.6ms
2 -8.400251875023 -2.62 -1.74 1.0 18.9ms
3 -8.400407437159 -3.81 -2.97 1.5 19.7ms
4 -8.400427836984 -4.69 -2.98 3.0 23.8ms
5 -8.400427938122 -7.00 -3.04 1.0 19.1ms
6 -8.400428145555 -6.68 -4.77 1.0 19.1ms
7 -8.400428151778 -8.21 -4.43 3.2 24.6ms
8 -8.400428152184 -9.39 -5.26 1.0 19.0ms
9 -8.400428152206 -10.64 -6.43 1.0 19.2ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397850737087 -0.90 5.0 25.8ms
2 -8.400386509718 -2.60 -1.79 0.80 2.0 21.2ms
3 -8.400423091101 -4.44 -3.04 0.80 1.0 20.4ms
4 -8.400428109497 -5.30 -3.45 0.80 2.2 24.8ms
5 -8.400428148506 -7.41 -4.81 0.80 1.0 21.3ms
6 -8.400428152193 -8.43 -5.61 0.80 3.0 26.7ms
7 -8.400428152209 -10.81 -6.24 0.80 2.2 24.2ms
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.808642976091 -1.00 72.1ms
2 -1.735241709001 0.41 -0.63 49.0ms
3 -4.468315432957 0.44 -0.36 44.5ms
4 -6.067644117166 0.20 -0.39 44.5ms
5 -7.576717567012 0.18 -0.69 44.7ms
6 -7.956103653720 -0.42 -1.30 32.7ms
7 -8.179165787509 -0.65 -1.56 33.0ms
8 -8.263346094757 -1.07 -1.85 33.4ms
9 -8.316888978739 -1.27 -2.11 32.6ms
10 -8.348294917335 -1.50 -1.90 32.7ms
11 -8.376247433220 -1.55 -2.27 32.5ms
12 -8.387639251527 -1.94 -2.46 32.8ms
13 -8.395454959556 -2.11 -3.03 32.5ms
14 -8.397896891478 -2.61 -3.07 32.4ms
15 -8.399447821643 -2.81 -3.22 32.6ms
16 -8.399972687530 -3.28 -3.17 32.6ms
17 -8.400249265716 -3.56 -3.74 39.1ms
18 -8.400362113945 -3.95 -3.49 33.1ms
19 -8.400398395550 -4.44 -3.62 32.8ms
20 -8.400420143724 -4.66 -3.96 32.8ms
21 -8.400424413403 -5.37 -4.62 32.6ms
22 -8.400426653721 -5.65 -4.81 32.6ms
23 -8.400427464677 -6.09 -5.20 32.6ms
24 -8.400427901493 -6.36 -5.05 32.6ms
25 -8.400428037740 -6.87 -5.19 32.5ms
26 -8.400428116875 -7.10 -5.13 33.2ms
27 -8.400428132364 -7.81 -5.58 32.9ms
28 -8.400428144658 -7.91 -5.64 32.9ms
29 -8.400428148735 -8.39 -5.80 32.4ms
30 -8.400428150579 -8.73 -6.50 38.5ms
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.397777615146 -0.90 5.0 25.5ms
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.400427974948 -1.79 568ms
2 -8.400428152209 -6.75 -4.03 400ms
3 -8.400428152209 -14.45 -7.81 139ms
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.7391182297066575e-6
|ρ_newton - ρ_scfv| = 2.410341469367132e-7
|ρ_newton - ρ_dm| = 1.5256117482491983e-6