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.397851538169                   -0.90    5.0   25.7ms
  2   -8.400258128203       -2.62       -1.73    1.0   18.5ms
  3   -8.400402124027       -3.84       -2.90    1.2   18.9ms
  4   -8.400427771827       -4.59       -2.91    3.2   33.2ms
  5   -8.400427947695       -6.75       -3.03    1.0   18.8ms
  6   -8.400428145750       -6.70       -4.84    1.0   18.9ms
  7   -8.400428151735       -8.22       -4.40    3.2   24.3ms
  8   -8.400428152193       -9.34       -5.78    1.0   19.1ms
  9   -8.400428152209      -10.80       -6.13    2.2   22.6ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397850193741                   -0.90           5.2   26.2ms
  2   -8.400386344144       -2.60       -1.78   0.80    2.0   18.6ms
  3   -8.400423182103       -4.43       -3.03   0.80    1.0   16.1ms
  4   -8.400428102624       -5.31       -3.41   0.80    2.5   19.8ms
  5   -8.400428147795       -7.35       -4.50   0.80    1.0   16.0ms
  6   -8.400428152145       -8.36       -5.69   0.80    2.2   25.7ms
  7   -8.400428152209      -10.20       -6.21   0.80    3.0   21.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/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.763118547381                   -1.03   57.3ms
  2   -1.632440866257        0.38       -0.60   32.5ms
  3   -4.560282562545        0.47       -0.35   43.9ms
  4   -6.122417940015        0.19       -0.43   43.7ms
  5   -7.616937714078        0.17       -0.68   48.8ms
  6   -7.948762508452       -0.48       -1.42   32.3ms
  7   -8.175708762109       -0.64       -1.56   32.4ms
  8   -8.261109397122       -1.07       -1.72   32.5ms
  9   -8.328819924426       -1.17       -2.00   34.1ms
 10   -8.359923734423       -1.51       -2.03   32.5ms
 11   -8.383162048723       -1.63       -2.61   32.2ms
 12   -8.390150130139       -2.16       -2.86   36.8ms
 13   -8.395742257156       -2.25       -2.79   32.4ms
 14   -8.397537115394       -2.75       -3.61   32.3ms
 15   -8.398732273417       -2.92       -3.11   32.2ms
 16   -8.399486878758       -3.12       -3.15   32.4ms
 17   -8.400049516250       -3.25       -3.14   37.2ms
 18   -8.400273732549       -3.65       -3.34   32.4ms
 19   -8.400368607382       -4.02       -3.56   32.3ms
 20   -8.400399643437       -4.51       -4.05   32.4ms
 21   -8.400414117089       -4.84       -3.96   35.7ms
 22   -8.400420567626       -5.19       -4.29   32.5ms
 23   -8.400424834984       -5.37       -4.45   32.3ms
 24   -8.400426531002       -5.77       -4.47   32.4ms
 25   -8.400427503396       -6.01       -5.26   35.9ms
 26   -8.400427818401       -6.50       -5.05   32.3ms
 27   -8.400427967512       -6.83       -5.14   32.5ms
 28   -8.400428069664       -6.99       -5.67   36.0ms
 29   -8.400428106623       -7.43       -5.72   32.3ms
 30   -8.400428135891       -7.53       -5.75   32.7ms
 31   -8.400428146427       -7.98       -6.02   32.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.397813804824                   -0.90    5.0   25.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.400427980959                   -1.79    590ms
  2   -8.400428152209       -6.77       -4.03    395ms
  3   -8.400428152209      -14.75       -7.83    122ms

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.845529519497723e-7
|ρ_newton - ρ_scfv| = 2.023343343927827e-7
|ρ_newton - ρ_dm|   = 3.6491326321846233e-6