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.397492041558                   -0.90    5.0   26.9ms
  2   -8.400200471588       -2.57       -1.72    1.0   18.7ms
  3   -8.400393874502       -3.71       -2.98    1.5   20.0ms
  4   -8.400427747836       -4.47       -2.91    3.2   25.1ms
  5   -8.400428047941       -6.52       -3.21    1.0   19.2ms
  6   -8.400428148555       -7.00       -4.75    1.0   19.0ms
  7   -8.400428152024       -8.46       -4.66    2.8   64.0ms
  8   -8.400428152180       -9.80       -5.04    1.0   20.4ms
  9   -8.400428152208      -10.57       -6.11    1.0   19.8ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime 
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397531305810                   -0.90           5.2    1.77s
  2   -8.400387857574       -2.54       -1.79   0.80    2.0    530ms
  3   -8.400423486440       -4.45       -3.00   0.80    1.0    244ms
  4   -8.400428099229       -5.34       -3.45   0.80    2.5   21.0ms
  5   -8.400428149422       -7.30       -4.76   0.80    1.2   17.7ms
  6   -8.400428152196       -8.56       -5.62   0.80    3.0   22.6ms
  7   -8.400428152209      -10.90       -6.28   0.80    2.5   20.6ms

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.203376302240                   -1.08    3.60s
  2   -1.253902026321        0.39       -0.65    149ms
  3   -4.156767131602        0.46       -0.36   44.0ms
  4   -5.662816819378        0.18       -0.41   43.6ms
  5   -7.338349865905        0.22       -0.67    1.18s
  6   -7.591626136808       -0.60       -1.29   32.6ms
  7   -8.088616607607       -0.30       -1.58   32.6ms
  8   -8.179152058705       -1.04       -1.87   32.5ms
  9   -8.306200486114       -0.90       -1.74   32.5ms
 10   -8.352997856725       -1.33       -2.30   32.3ms
 11   -8.378381586124       -1.60       -2.34   32.4ms
 12   -8.389014059735       -1.97       -2.47   32.3ms
 13   -8.395837580205       -2.17       -2.69   32.4ms
 14   -8.398194655417       -2.63       -2.66   32.2ms
 15   -8.399492038612       -2.89       -3.04   32.2ms
 16   -8.400038250430       -3.26       -3.23   32.4ms
 17   -8.400277272243       -3.62       -3.89   32.2ms
 18   -8.400345741617       -4.16       -3.73   32.0ms
 19   -8.400392834280       -4.33       -4.44   32.3ms
 20   -8.400414484554       -4.66       -4.09   32.1ms
 21   -8.400423267330       -5.06       -4.29   32.2ms
 22   -8.400425300473       -5.69       -4.27   32.0ms
 23   -8.400427155301       -5.73       -5.21   31.9ms
 24   -8.400427513474       -6.45       -4.58   32.1ms
 25   -8.400427880212       -6.44       -5.41   32.1ms
 26   -8.400427987661       -6.97       -4.91   32.1ms
 27   -8.400428096515       -6.96       -5.39   32.3ms
 28   -8.400428120808       -7.61       -5.32   59.1ms
 29   -8.400428138332       -7.76       -5.96   32.4ms
 30   -8.400428145940       -8.12       -5.50   32.8ms
 31   -8.400428149965       -8.40       -5.97   33.0ms
 32   -8.400428151124       -8.94       -6.64   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.397540456440                   -0.90    5.0   27.2ms

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.400427942139                   -1.78    11.7s
  2   -8.400428152209       -6.68       -3.98    3.70s
  3   -8.400428152209      -14.27       -7.74   95.9ms

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.86275260185875e-7
|ρ_newton - ρ_scfv| = 2.8896459288936294e-7
|ρ_newton - ρ_dm|   = 2.2555026197706125e-6