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.397556573161                   -0.90    5.2   85.0ms
  2   -8.400199915895       -2.58       -1.72    1.0   19.5ms
  3   -8.400396581240       -3.71       -3.02    1.5   20.2ms
  4   -8.400427729930       -4.51       -2.94    3.5   24.1ms
  5   -8.400428044172       -6.50       -3.24    1.0   19.2ms
  6   -8.400428148968       -6.98       -4.79    1.0   19.3ms
  7   -8.400428151906       -8.53       -4.56    2.5   22.5ms
  8   -8.400428152178       -9.57       -5.10    1.0   19.6ms
  9   -8.400428152207      -10.53       -6.55    1.0   54.7ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime 
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397533205989                   -0.90           4.8    1.57s
  2   -8.400387008210       -2.54       -1.79   0.80    2.2    1.05s
  3   -8.400423569382       -4.44       -2.98   0.80    1.0    200ms
  4   -8.400428110123       -5.34       -3.45   0.80    2.5   19.6ms
  5   -8.400428148582       -7.42       -4.77   0.80    1.0    108ms
  6   -8.400428152197       -8.44       -5.66   0.80    3.0   21.3ms
  7   -8.400428152209      -10.93       -6.26   0.80    2.0   19.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/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   +1.138657044274                   -1.09    3.12s
  2   -1.477336677558        0.42       -0.64    149ms
  3   -4.144120869544        0.43       -0.35    111ms
  4   -5.470665604587        0.12       -0.41   47.6ms
  5   -7.222437465467        0.24       -0.55   47.8ms
  6   -7.832355824548       -0.21       -1.01   47.1ms
  7   -8.067025250503       -0.63       -1.28   34.6ms
  8   -8.237208526217       -0.77       -1.68   71.5ms
  9   -8.318477802794       -1.09       -1.85   35.2ms
 10   -8.345881210309       -1.56       -2.16   35.0ms
 11   -8.369535784295       -1.63       -2.30   35.4ms
 12   -8.381970521675       -1.91       -2.20   35.0ms
 13   -8.392764420084       -1.97       -2.38   35.2ms
 14   -8.396754637899       -2.40       -2.61   35.2ms
 15   -8.399253313240       -2.60       -2.74   41.6ms
 16   -8.399912974859       -3.18       -3.51   35.2ms
 17   -8.400242287429       -3.48       -3.39   35.1ms
 18   -8.400331270884       -4.05       -4.22   35.1ms
 19   -8.400394165382       -4.20       -4.05   39.1ms
 20   -8.400412842072       -4.73       -4.21   35.3ms
 21   -8.400424028494       -4.95       -4.30   35.4ms
 22   -8.400426270879       -5.65       -4.17   35.2ms
 23   -8.400427565581       -5.89       -4.47   35.0ms
 24   -8.400427788826       -6.65       -4.83   41.4ms
 25   -8.400428038580       -6.60       -4.76   35.3ms
 26   -8.400428084455       -7.34       -5.59   35.3ms
 27   -8.400428127905       -7.36       -5.63   35.4ms
 28   -8.400428137897       -8.00       -6.06   35.2ms

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.397647064650                   -0.90    5.2   24.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.400427982831                   -1.79    8.92s
  2   -8.400428152209       -6.77       -4.02    3.24s
  3   -8.400428152209      -14.75       -7.85   75.8ms

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|  = 9.755870996891817e-7
|ρ_newton - ρ_scfv| = 1.6065350275384802e-7
|ρ_newton - ρ_dm|   = 3.6371298183196053e-6