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.397931706513 -0.90 5.5 24.9ms
2 -8.400242222110 -2.64 -1.74 1.0 17.2ms
3 -8.400401457770 -3.80 -2.98 1.2 17.9ms
4 -8.400427815698 -4.58 -2.94 3.2 22.7ms
5 -8.400427943286 -6.89 -3.04 1.0 17.1ms
6 -8.400428150111 -6.68 -4.79 1.0 17.2ms
7 -8.400428155827 -8.24 -4.42 2.8 21.9ms
8 -8.400428156262 -9.36 -5.96 1.0 17.7ms
9 -8.400428156276 -10.85 -5.96 2.5 21.8ms
10 -8.400428156277 -12.21 -6.63 1.0 18.1ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397804189276 -0.90 5.2 31.2ms
2 -8.400384024766 -2.59 -1.78 0.80 2.2 22.4ms
3 -8.400423224439 -4.41 -3.01 0.80 1.0 19.3ms
4 -8.400428108718 -5.31 -3.48 0.80 2.5 23.6ms
5 -8.400428152170 -7.36 -4.71 0.80 1.2 20.1ms
6 -8.400428156259 -8.39 -5.44 0.80 2.8 35.5ms
7 -8.400428156276 -10.78 -6.92 0.80 1.8 17.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/8dE7C/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/8dE7C/src/types.jl:120
n Energy log10(ΔE) log10(Δρ) Δtime
--- --------------- --------- --------- ------
1 +1.100854778973 -1.06 51.0ms
2 -1.645318362252 0.44 -0.68 28.6ms
3 -4.550888789956 0.46 -0.39 38.4ms
4 -6.195456648319 0.22 -0.49 38.8ms
5 -7.601378730970 0.15 -0.76 38.6ms
6 -8.023986560934 -0.37 -1.39 28.5ms
7 -8.209256312844 -0.73 -1.65 28.3ms
8 -8.281930897926 -1.14 -1.77 28.7ms
9 -8.319365965735 -1.43 -2.36 28.1ms
10 -8.339383707243 -1.70 -2.40 28.1ms
11 -8.358547940460 -1.72 -2.16 28.0ms
12 -8.374454821829 -1.80 -2.15 28.0ms
13 -8.387923708235 -1.87 -2.35 28.0ms
14 -8.395330133532 -2.13 -2.82 40.4ms
15 -8.398010712697 -2.57 -2.92 28.3ms
16 -8.399626832535 -2.79 -3.21 28.0ms
17 -8.400052158648 -3.37 -3.61 28.2ms
18 -8.400313372302 -3.58 -3.55 28.3ms
19 -8.400369249438 -4.25 -3.80 28.2ms
20 -8.400411605930 -4.37 -3.92 28.4ms
21 -8.400418383925 -5.17 -4.17 28.5ms
22 -8.400425009704 -5.18 -4.57 28.3ms
23 -8.400426704576 -5.77 -4.56 28.3ms
24 -8.400427635503 -6.03 -4.83 28.4ms
25 -8.400427900928 -6.58 -4.74 28.3ms
26 -8.400428038447 -6.86 -5.06 28.2ms
27 -8.400428109675 -7.15 -5.09 34.3ms
28 -8.400428132706 -7.64 -6.18 28.6ms
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.397803133072 -0.90 5.2 24.7ms
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.400427976433 -1.79 547ms
2 -8.400428156277 -6.75 -4.02 369ms
3 -8.400428156277 -14.45 -7.80 124ms
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.444120027233851e-7
|ρ_newton - ρ_scfv| = 5.11477523540288e-7
|ρ_newton - ρ_dm| = 7.241798945936816e-6