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.397880059610 -0.90 5.2 24.5ms
2 -8.400244575049 -2.63 -1.74 1.0 17.6ms
3 -8.400407389741 -3.79 -2.97 1.5 28.2ms
4 -8.400427865074 -4.69 -3.00 2.8 22.3ms
5 -8.400427969566 -6.98 -3.08 1.0 17.7ms
6 -8.400428149762 -6.74 -4.40 1.0 17.7ms
7 -8.400428155841 -8.22 -4.47 2.0 20.4ms
8 -8.400428156236 -9.40 -5.02 1.0 18.0ms
9 -8.400428156274 -10.42 -6.26 1.5 18.8ms
Potential-based SCF
scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n Energy log10(ΔE) log10(Δρ) α Diag Δtime
--- --------------- --------- --------- ---- ---- ------
1 -8.397840930287 -0.90 5.2 24.8ms
2 -8.400386054131 -2.59 -1.79 0.80 2.0 17.8ms
3 -8.400423050608 -4.43 -3.03 0.80 1.0 15.1ms
4 -8.400428111756 -5.30 -3.42 0.80 2.8 19.3ms
5 -8.400428154618 -7.37 -4.49 0.80 1.5 23.3ms
6 -8.400428156245 -8.79 -5.65 0.80 2.0 17.8ms
7 -8.400428156276 -10.50 -6.42 0.80 2.8 19.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 +0.865547824673 -1.01 50.7ms
2 -1.288313265980 0.33 -0.66 28.7ms
3 -4.097123720406 0.45 -0.34 38.1ms
4 -5.360682000074 0.10 -0.42 37.9ms
5 -7.367113135079 0.30 -0.53 44.3ms
6 -7.962102776699 -0.23 -1.11 38.9ms
7 -8.145764786062 -0.74 -1.27 28.4ms
8 -8.264683713413 -0.92 -1.68 28.6ms
9 -8.314504990127 -1.30 -1.84 28.3ms
10 -8.344764904650 -1.52 -1.99 33.2ms
11 -8.367016866406 -1.65 -2.09 28.9ms
12 -8.382028788590 -1.82 -2.15 28.6ms
13 -8.392178482550 -1.99 -2.41 28.7ms
14 -8.396965448672 -2.32 -2.73 28.6ms
15 -8.398962434081 -2.70 -2.72 28.5ms
16 -8.399723479189 -3.12 -3.06 33.0ms
17 -8.400137936897 -3.38 -3.20 28.1ms
18 -8.400311930274 -3.76 -3.80 28.1ms
19 -8.400387458276 -4.12 -3.70 28.3ms
20 -8.400408624430 -4.67 -4.26 33.3ms
21 -8.400421025257 -4.91 -3.93 28.5ms
22 -8.400425118159 -5.39 -4.43 28.8ms
23 -8.400427447175 -5.63 -4.34 32.3ms
24 -8.400427802121 -6.45 -4.69 28.9ms
25 -8.400428031797 -6.64 -5.37 29.0ms
26 -8.400428087830 -7.25 -5.02 28.6ms
27 -8.400428127487 -7.40 -5.36 32.9ms
28 -8.400428140728 -7.88 -5.50 28.2ms
29 -8.400428149684 -8.05 -5.56 28.3ms
30 -8.400428152596 -8.54 -5.82 28.4ms
31 -8.400428154766 -8.66 -5.97 32.6ms
32 -8.400428155568 -9.10 -6.32 28.7ms
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.397840924757 -0.90 5.2 29.3ms
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.400427982103 -1.78 530ms
2 -8.400428156277 -6.76 -4.03 355ms
3 -8.400428156277 -14.45 -7.82 108ms
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| = 7.264133849910721e-7
|ρ_newton - ρ_scfv| = 4.5754765066279303e-7
|ρ_newton - ρ_dm| = 1.3915523196257542e-6