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.397805109440                   -0.90    5.2   27.5ms
  2   -8.400249679235       -2.61       -1.74    1.0   19.5ms
  3   -8.400406854707       -3.80       -2.97    1.5   20.5ms
  4   -8.400427886620       -4.68       -3.00    3.0   24.8ms
  5   -8.400427979893       -7.03       -3.08    1.0   19.8ms
  6   -8.400428146381       -6.78       -4.80    1.0   19.5ms
  7   -8.400428151899       -8.26       -4.50    3.2   54.6ms
  8   -8.400428152198       -9.52       -6.04    1.0   20.3ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397873638813                   -0.90           5.2   28.6ms
  2   -8.400383079695       -2.60       -1.78   0.80    2.0   20.0ms
  3   -8.400422336437       -4.41       -3.00   0.80    1.0   16.9ms
  4   -8.400428110397       -5.24       -3.38   0.80    2.5   21.3ms
  5   -8.400428150000       -7.40       -4.42   0.80    1.2   17.8ms
  6   -8.400428152143       -8.67       -5.92   0.80    2.2   20.8ms
  7   -8.400428152209      -10.18       -6.07   0.80    3.2   24.0ms

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   +0.466030040445                   -1.11   58.8ms
  2   -1.732663266432        0.34       -0.65   32.9ms
  3   -4.777219994828        0.48       -0.42   44.1ms
  4   -6.349398786930        0.20       -0.52   44.3ms
  5   -7.732889752533        0.14       -0.82   44.2ms
  6   -8.023348192913       -0.54       -1.34   32.8ms
  7   -8.197180676327       -0.76       -1.68   32.7ms
  8   -8.276781692779       -1.10       -2.21   32.7ms
  9   -8.327935980610       -1.29       -2.12   44.5ms
 10   -8.351173826465       -1.63       -2.35   32.7ms
 11   -8.372394313965       -1.67       -2.28   33.4ms
 12   -8.385506059602       -1.88       -2.41   33.4ms
 13   -8.391931772392       -2.19       -2.41   33.0ms
 14   -8.396913417064       -2.30       -2.44   33.0ms
 15   -8.398457680490       -2.81       -3.27   33.0ms
 16   -8.399649861876       -2.92       -3.11   32.8ms
 17   -8.400031782925       -3.42       -3.69   32.9ms
 18   -8.400261847139       -3.64       -3.37   32.8ms
 19   -8.400345883327       -4.08       -3.94   32.7ms
 20   -8.400389686583       -4.36       -3.79   32.6ms
 21   -8.400410415942       -4.68       -3.83   32.7ms
 22   -8.400420274207       -5.01       -4.20   32.6ms
 23   -8.400424395399       -5.38       -4.25   42.1ms
 24   -8.400426702476       -5.64       -4.35   33.2ms
 25   -8.400427545717       -6.07       -4.55   32.8ms
 26   -8.400427857415       -6.51       -5.07   33.0ms
 27   -8.400428038459       -6.74       -4.95   32.9ms
 28   -8.400428095392       -7.24       -5.80   32.8ms
 29   -8.400428127360       -7.50       -5.41   32.8ms
 30   -8.400428140609       -7.88       -5.82   32.7ms
 31   -8.400428147467       -8.16       -5.80   33.5ms
 32   -8.400428150292       -8.55       -5.91   32.7ms
 33   -8.400428151455       -8.93       -6.20   32.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.397809710138                   -0.90    5.5   36.5ms

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.400427990853                   -1.80    603ms
  2   -8.400428152209       -6.79       -4.04    393ms
  3   -8.400428152209      -14.75       -7.87    107ms

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|  = 3.5230247129346427e-6
|ρ_newton - ρ_scfv| = 1.4917457135556732e-7
|ρ_newton - ρ_dm|   = 1.7797462751467772e-6