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.397871684260                   -0.89    5.0   27.2ms
  2   -8.400258931030       -2.62       -1.73    1.0   19.6ms
  3   -8.400400904333       -3.85       -2.90    1.2   31.0ms
  4   -8.400427707761       -4.57       -2.89    3.2   25.1ms
  5   -8.400427903515       -6.71       -2.99    1.0   19.7ms
  6   -8.400428143778       -6.62       -4.78    1.0   19.7ms
  7   -8.400428151576       -8.11       -4.33    3.2   25.7ms
  8   -8.400428152186       -9.21       -5.85    1.0   20.5ms
  9   -8.400428152209      -10.64       -5.91    2.2   32.9ms
 10   -8.400428152209   +  -13.50       -5.92    1.0   20.0ms
 11   -8.400428152209      -12.60       -6.32    1.0   20.5ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397731805557                   -0.90           5.0   26.9ms
  2   -8.400381210642       -2.58       -1.80   0.80    2.2   20.0ms
  3   -8.400422528823       -4.38       -3.03   0.80    1.0   24.0ms
  4   -8.400428095043       -5.25       -3.44   0.80    2.5   21.1ms
  5   -8.400428148384       -7.27       -4.57   0.80    1.2   17.5ms
  6   -8.400428152184       -8.42       -5.64   0.80    2.5   20.7ms
  7   -8.400428152209      -10.61       -6.12   0.80    2.5   21.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   +0.899934478728                   -1.03   58.0ms
  2   -1.691384719218        0.41       -0.68   32.4ms
  3   -4.468274557807        0.44       -0.38   50.9ms
  4   -6.081411528329        0.21       -0.47   43.8ms
  5   -7.615560542601        0.19       -0.69   43.7ms
  6   -8.030811593196       -0.38       -1.41   32.2ms
  7   -8.210587557168       -0.75       -1.68   39.2ms
  8   -8.243422912498       -1.48       -2.00   32.8ms
  9   -8.292211535855       -1.31       -2.06   32.9ms
 10   -8.310926579201       -1.73       -2.23   32.5ms
 11   -8.344017106425       -1.48       -2.04   32.4ms
 12   -8.373056081001       -1.54       -1.85   32.0ms
 13   -8.387466044340       -1.84       -2.56   38.3ms
 14   -8.395368217210       -2.10       -2.72   32.4ms
 15   -8.398149198979       -2.56       -2.74   32.5ms
 16   -8.399526321548       -2.86       -2.99   32.3ms
 17   -8.400032316983       -3.30       -2.99   32.0ms
 18   -8.400277846626       -3.61       -3.23   31.7ms
 19   -8.400360387826       -4.08       -3.27   38.7ms
 20   -8.400402653517       -4.37       -3.59   33.0ms
 21   -8.400413874892       -4.95       -3.71   32.6ms
 22   -8.400422787880       -5.05       -3.98   32.4ms
 23   -8.400425148563       -5.63       -4.02   32.3ms
 24   -8.400427036235       -5.72       -4.53   38.0ms
 25   -8.400427725488       -6.16       -4.28   32.3ms
 26   -8.400427990509       -6.58       -5.13   32.4ms
 27   -8.400428076523       -7.07       -4.71   32.4ms
 28   -8.400428124223       -7.32       -5.22   32.5ms
 29   -8.400428139541       -7.81       -5.10   32.7ms
 30   -8.400428147065       -8.12       -5.74   38.7ms
 31   -8.400428149770       -8.57       -5.44   32.5ms
 32   -8.400428151276       -8.82       -6.23   32.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.397864708749                   -0.90    5.8   27.9ms

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.400427987447                   -1.79    582ms
  2   -8.400428152209       -6.78       -4.04    406ms
  3   -8.400428152209      -14.45       -7.85    106ms

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|  = 5.068046991347762e-7
|ρ_newton - ρ_scfv| = 1.8503596117007803e-7
|ρ_newton - ρ_dm|   = 2.4855398025396908e-6