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.397539502634                   -0.90    5.0   27.1ms
  2   -8.400195871695       -2.58       -1.72    1.0   18.9ms
  3   -8.400398110394       -3.69       -3.02    1.5   19.9ms
  4   -8.400427766454       -4.53       -2.94    3.2   98.0ms
  5   -8.400428048210       -6.55       -3.24    1.0   19.3ms
  6   -8.400428149314       -7.00       -4.80    1.0   19.3ms
  7   -8.400428152017       -8.57       -4.65    2.8   23.8ms
  8   -8.400428152192       -9.76       -5.24    1.0   19.8ms
  9   -8.400428152208      -10.80       -6.44    1.0   20.0ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime 
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397521905396                   -0.90           5.8    1.82s
  2   -8.400387290136       -2.54       -1.78   0.80    2.2    553ms
  3   -8.400423721157       -4.44       -3.00   0.80    1.0    261ms
  4   -8.400428117489       -5.36       -3.52   0.80    2.5   21.5ms
  5   -8.400428149100       -7.50       -5.01   0.80    1.0   17.6ms
  6   -8.400428152201       -8.51       -5.38   0.80    3.5   24.5ms
  7   -8.400428152208      -11.12       -6.13   0.80    1.0   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/gmigl/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/gmigl/src/types.jl:120
n     Energy            log10(ΔE)   log10(Δρ)   Δtime 
---   ---------------   ---------   ---------   ------
  1   +0.964874060889                   -1.07    4.98s
  2   -1.402789038380        0.37       -0.66    158ms
  3   -4.021552573625        0.42       -0.38   45.7ms
  4   -5.359226123539        0.13       -0.45   46.0ms
  5   -7.114977115250        0.24       -0.57   46.4ms
  6   -7.737812476674       -0.21       -1.03   45.6ms
  7   -8.001476161793       -0.58       -1.37   33.8ms
  8   -8.213132832119       -0.67       -1.54   34.2ms
  9   -8.309687717281       -1.02       -1.65   34.1ms
 10   -8.355478584503       -1.34       -2.09   33.7ms
 11   -8.378987017113       -1.63       -2.83   33.9ms
 12   -8.388264796329       -2.03       -2.56   34.1ms
 13   -8.394146655577       -2.23       -2.75   34.0ms
 14   -8.396235210067       -2.68       -3.04   34.3ms
 15   -8.397832460939       -2.80       -3.21   33.9ms
 16   -8.398857114642       -2.99       -3.13   33.7ms
 17   -8.399653855209       -3.10       -3.09   33.8ms
 18   -8.400159188710       -3.30       -3.18   33.6ms
 19   -8.400325656600       -3.78       -3.67   34.4ms
 20   -8.400393043016       -4.17       -3.75   33.7ms
 21   -8.400413442919       -4.69       -3.81   92.5ms
 22   -8.400422332062       -5.05       -4.45   34.0ms
 23   -8.400426068171       -5.43       -4.18   34.1ms
 24   -8.400427214899       -5.94       -4.58   34.1ms
 25   -8.400427774090       -6.25       -5.05   34.2ms
 26   -8.400428026003       -6.60       -5.33   34.0ms
 27   -8.400428086486       -7.22       -5.15   34.0ms
 28   -8.400428122199       -7.45       -5.27   33.9ms
 29   -8.400428136371       -7.85       -5.76   33.6ms
 30   -8.400428144976       -8.07       -5.51   34.5ms
 31   -8.400428149373       -8.36       -6.23   34.3ms

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.397523099214                   -0.90    5.5   28.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.400427982952                   -1.79    12.0s
  2   -8.400428152209       -6.77       -4.04    3.84s
  3   -8.400428152209      -14.45       -7.85    102ms

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|  = 4.7551806089688964e-7
|ρ_newton - ρ_scfv| = 4.7313696123730925e-7
|ρ_newton - ρ_dm|   = 5.409855931616314e-6