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-61.0e-6Density-based self-consistent field
scfres_scf = self_consistent_field(basis; tol);n Energy log10(ΔE) log10(Δρ) Diag Δtime
--- --------------- --------- --------- ---- ------
1 -8.397899562728 -0.90 5.5 27.3ms
2 -8.400250018371 -2.63 -1.74 1.0 18.6ms
3 -8.400404632689 -3.81 -2.96 1.5 20.1ms
4 -8.400427816087 -4.63 -2.96 3.0 24.1ms
5 -8.400427932567 -6.93 -3.04 1.0 45.3ms
6 -8.400428149651 -6.66 -4.71 1.0 19.4ms
7 -8.400428155721 -8.22 -4.39 2.8 24.7ms
8 -8.400428156256 -9.27 -5.49 1.0 19.2ms
9 -8.400428156274 -10.72 -5.93 2.0 22.1ms
10 -8.400428156277 -11.65 -6.78 2.0 22.3msPotential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397855001034 -0.90 5.2 27.4ms
2 -8.400385108564 -2.60 -1.79 0.80 2.0 18.7ms
3 -8.400422974787 -4.42 -3.02 0.80 1.0 15.6ms
4 -8.400428122361 -5.29 -3.43 0.80 2.5 20.2ms
5 -8.400428154014 -7.50 -4.65 0.80 1.0 15.9ms
6 -8.400428156244 -8.65 -6.05 0.80 2.2 19.8msDirect minimization
scfres_dm = direct_minimization(basis; tol);n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +0.736466371039 -1.02 64.5ms
2 -2.008508186068 0.44 -0.66 29.0ms
3 -4.048867996341 0.31 -0.43 39.3ms
4 -5.129753115541 0.03 -0.51 39.3ms
5 -6.921516484543 0.25 -0.63 39.3ms
6 -7.765720205019 -0.07 -0.98 39.3ms
7 -8.025591123676 -0.59 -1.40 28.8ms
8 -8.197351838844 -0.77 -1.78 28.7ms
9 -8.246309292700 -1.31 -1.73 36.6ms
10 -8.310890985071 -1.19 -1.91 29.0ms
11 -8.341148800847 -1.52 -2.01 29.5ms
12 -8.366979149859 -1.59 -2.28 29.2ms
13 -8.381481944099 -1.84 -2.24 29.0ms
14 -8.388802789566 -2.14 -2.46 28.9ms
15 -8.394931093548 -2.21 -2.74 29.0ms
16 -8.398215400385 -2.48 -2.84 35.6ms
17 -8.399393980607 -2.93 -3.17 29.0ms
18 -8.399894263973 -3.30 -3.47 28.9ms
19 -8.400154298128 -3.58 -3.52 29.0ms
20 -8.400343914354 -3.72 -3.82 29.1ms
21 -8.400366547956 -4.65 -4.01 28.9ms
22 -8.400407990149 -4.38 -4.38 35.3ms
23 -8.400415293739 -5.14 -4.05 29.1ms
24 -8.400423850823 -5.07 -4.48 29.1ms
25 -8.400425990602 -5.67 -4.50 29.0ms
26 -8.400427200852 -5.92 -4.69 29.1ms
27 -8.400427804216 -6.22 -4.63 35.0ms
28 -8.400428037808 -6.63 -5.23 29.0ms
29 -8.400428079330 -7.38 -5.04 29.0ms
30 -8.400428118618 -7.41 -5.31 29.0ms
31 -8.400428134795 -7.79 -5.34 34.1ms
32 -8.400428147613 -7.89 -5.46 28.9ms
33 -8.400428152149 -8.34 -5.79 33.7ms
34 -8.400428154696 -8.59 -5.77 29.9ms
35 -8.400428155701 -9.00 -6.08 35.0msNewton 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.397770640602 -0.90 5.0 29.1msRemove 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.400427955926 -1.78 559ms
2 -8.400428156277 -6.70 -4.00 383ms
3 -8.400428156277 -14.45 -7.75 125msComparison 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| = 5.29491732327665e-8
|ρ_newton - ρ_scfv| = 4.7724894625198865e-6
|ρ_newton - ρ_dm| = 1.1358496482929661e-6