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.397494083412                   -0.90    5.2   26.8ms
  2   -8.400181282264       -2.57       -1.72    1.0   18.6ms
  3   -8.400394186450       -3.67       -3.03    1.5   19.4ms
  4   -8.400427804963       -4.47       -2.96    3.2   24.9ms
  5   -8.400428101041       -6.53       -3.43    1.0   70.5ms
  6   -8.400428149938       -7.31       -4.89    1.0   19.2ms
  7   -8.400428152086       -8.67       -4.80    2.8   23.8ms
  8   -8.400428152189       -9.99       -5.18    1.0   19.3ms
  9   -8.400428152208      -10.73       -6.19    1.0   19.7ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime 
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397614434192                   -0.90           5.2    1.82s
  2   -8.400385117274       -2.56       -1.79   0.80    2.0    519ms
  3   -8.400423340454       -4.42       -2.97   0.80    1.0    232ms
  4   -8.400428112150       -5.32       -3.48   0.80    2.5   21.6ms
  5   -8.400428148744       -7.44       -4.79   0.80    1.2   18.6ms
  6   -8.400428152193       -8.46       -6.11   0.80    3.0   22.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/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.016762316222                   -1.09    3.63s
  2   -1.867629350370        0.46       -0.65    157ms
  3   -4.631248220240        0.44       -0.36   44.9ms
  4   -6.221877940494        0.20       -0.43   45.1ms
  5   -7.597854008489        0.14       -0.68   86.7ms
  6   -7.942135844718       -0.46       -1.28   33.3ms
  7   -8.156467047149       -0.67       -1.58   32.8ms
  8   -8.189350420834       -1.48       -1.73   33.0ms
  9   -8.191057833656       -2.77       -1.80   32.9ms
 10   -8.256275235509       -1.19       -1.72   43.9ms
 11   -8.310929517139       -1.26       -1.66   79.4ms
 12   -8.347110265224       -1.44       -2.53   32.7ms
 13   -8.358724876381       -1.93       -2.73   32.7ms
 14   -8.375997923764       -1.76       -2.29   32.8ms
 15   -8.388726543813       -1.90       -2.20   32.8ms
 16   -8.392429851825       -2.43       -2.60   39.8ms
 17   -8.397895573741       -2.26       -3.00   33.0ms
 18   -8.398634068414       -3.13       -3.37   32.8ms
 19   -8.399836694417       -2.92       -3.33   33.7ms
 20   -8.400091384066       -3.59       -3.34   34.2ms
 21   -8.400323832786       -3.63       -3.72   39.8ms
 22   -8.400357080776       -4.48       -3.90   33.4ms
 23   -8.400404641319       -4.32       -3.88   33.4ms
 24   -8.400409775657       -5.29       -4.15   33.5ms
 25   -8.400420891757       -4.95       -4.74   33.3ms
 26   -8.400424327836       -5.46       -4.45   40.3ms
 27   -8.400426911289       -5.59       -5.06   33.3ms
 28   -8.400427524161       -6.21       -4.83   33.1ms
 29   -8.400427959990       -6.36       -5.39   33.2ms
 30   -8.400428016384       -7.25       -5.28   33.1ms
 31   -8.400428114063       -7.01       -6.04   39.5ms

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.397448065425                   -0.90    5.0   27.1ms

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.400427942983                   -1.78    12.0s
  2   -8.400428152209       -6.68       -3.98    3.85s
  3   -8.400428152209      -14.45       -7.75   91.6ms

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|  = 6.053924158844756e-7
|ρ_newton - ρ_scfv| = 3.754686729521315e-6
|ρ_newton - ρ_dm|   = 7.685142693483296e-6