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.397742764959                   -0.90    5.0   26.0ms
  2   -8.400217793318       -2.61       -1.74    1.0   19.0ms
  3   -8.400404327688       -3.73       -2.95    1.5   20.2ms
  4   -8.400427878933       -4.63       -2.99    3.0   23.9ms
  5   -8.400428018688       -6.85       -3.14    1.0   26.8ms
  6   -8.400428147300       -6.89       -4.94    1.0   19.3ms
  7   -8.400428152005       -8.33       -4.59    3.5   26.1ms
  8   -8.400428152199       -9.71       -5.42    1.0   19.5ms
  9   -8.400428152208      -11.07       -6.12    1.0   19.7ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397843265668                   -0.90           5.2   26.6ms
  2   -8.400386627204       -2.59       -1.79   0.80    2.2   19.5ms
  3   -8.400422834491       -4.44       -3.02   0.80    1.0   16.5ms
  4   -8.400428111809       -5.28       -3.39   0.80    2.5   26.3ms
  5   -8.400428149798       -7.42       -4.44   0.80    1.0   16.7ms
  6   -8.400428152176       -8.62       -5.93   0.80    2.2   19.9ms
  7   -8.400428152209      -10.49       -5.89   0.80    3.0   22.3ms
  8   -8.400428152209      -12.89       -7.21   0.80    1.0   17.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   +1.197074277617                   -0.99   59.7ms
  2   -1.821709009552        0.48       -0.63   33.4ms
  3   -4.444926343277        0.42       -0.35   49.8ms
  4   -6.114580312377        0.22       -0.40   44.5ms
  5   -7.556420010640        0.16       -0.70   44.5ms
  6   -7.900475762705       -0.46       -1.38   33.5ms
  7   -8.204033685128       -0.52       -1.54   37.9ms
  8   -8.318600047222       -0.94       -1.77   33.1ms
  9   -8.364890644987       -1.33       -2.13   33.1ms
 10   -8.382939633450       -1.74       -2.28   33.1ms
 11   -8.392331930862       -2.03       -2.49   36.7ms
 12   -8.397004371448       -2.33       -2.87   33.1ms
 13   -8.398837745888       -2.74       -2.91   33.1ms
 14   -8.399778725947       -3.03       -3.18   33.1ms
 15   -8.400068292158       -3.54       -3.79   36.9ms
 16   -8.400296599091       -3.64       -3.52   32.8ms
 17   -8.400352286013       -4.25       -3.73   33.1ms
 18   -8.400397055817       -4.35       -3.74   36.8ms
 19   -8.400410298902       -4.88       -3.81   33.0ms
 20   -8.400423340482       -4.88       -3.83   33.1ms
 21   -8.400425636821       -5.64       -4.39   33.2ms
 22   -8.400427381278       -5.76       -4.51   37.5ms
 23   -8.400427617219       -6.63       -5.06   33.3ms
 24   -8.400428010836       -6.40       -4.76   33.2ms
 25   -8.400428052684       -7.38       -5.14   33.3ms
 26   -8.400428115155       -7.20       -5.62   36.5ms
 27   -8.400428131717       -7.78       -5.52   33.0ms
 28   -8.400428143697       -7.92       -5.73   33.1ms
 29   -8.400428147151       -8.46       -5.82   36.5ms
 30   -8.400428149976       -8.55       -6.01   33.1ms

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.397840710081                   -0.90    5.5   26.4ms

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.400427982745                   -1.79    599ms
  2   -8.400428152209       -6.77       -4.03    409ms
  3   -8.400428152209   +    -Inf       -7.83    121ms

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.652650757639361e-6
|ρ_newton - ρ_scfv| = 1.9371483567762906e-8
|ρ_newton - ρ_dm|   = 3.125837292622972e-6