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 DFTK
using LinearAlgebra

a = 10.26  # Silicon lattice constant in Bohr
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms     = [Si, Si]
positions = [ones(3)/8, -ones(3)/8]

model = model_LDA(lattice, atoms, positions)
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   -7.846835863630                   -0.70    4.8
  2   -7.852323709791       -2.26       -1.53    1.0   34.6ms
  3   -7.852612797784       -3.54       -2.55    1.2   37.2ms
  4   -7.852645917000       -4.48       -2.86    2.8   49.7ms
  5   -7.852646445638       -6.28       -3.10    1.0   34.8ms
  6   -7.852646677080       -6.64       -3.95    1.0   36.2ms
  7   -7.852646686241       -8.04       -5.13    1.5   42.6ms
  8   -7.852646686715       -9.32       -5.23    2.5   48.4ms
  9   -7.852646686729      -10.85       -6.20    1.0   35.1ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -7.846839598734                   -0.70           4.5
  2   -7.852528940588       -2.24       -1.63   0.80    2.0   38.8ms
  3   -7.852637887565       -3.96       -2.70   0.80    1.0   32.2ms
  4   -7.852646479322       -5.07       -3.31   0.80    2.0   46.6ms
  5   -7.852646682881       -6.69       -4.20   0.80    1.5   36.4ms
  6   -7.852646686602       -8.43       -4.87   0.80    2.0   45.8ms
  7   -7.852646686721       -9.93       -5.77   0.80    1.5   34.5ms
  8   -7.852646686730      -11.04       -6.79   0.80    2.0   38.9ms

Direct minimization

Note: Unlike the other algorithms, tolerance for this one is in the energy, thus we square the density tolerance value to be roughly equivalent.

scfres_dm = direct_minimization(basis; tol=tol^2);
Iter     Function value   Gradient norm
     0     1.377235e+01     3.881346e+00
 * time: 0.01699209213256836
     1     1.265511e+00     1.801907e+00
 * time: 0.049504995346069336
     2    -1.543474e+00     1.960013e+00
 * time: 0.08228492736816406
     3    -3.858899e+00     1.655478e+00
 * time: 0.13008499145507812
     4    -5.175275e+00     1.663854e+00
 * time: 0.1777019500732422
     5    -6.857095e+00     1.206923e+00
 * time: 0.22472691535949707
     6    -7.446812e+00     4.656910e-01
 * time: 0.2716360092163086
     7    -7.665155e+00     1.925371e-01
 * time: 0.3676130771636963
     8    -7.752370e+00     1.642304e-01
 * time: 0.40064191818237305
     9    -7.799010e+00     1.327639e-01
 * time: 0.4340860843658447
    10    -7.824560e+00     1.045189e-01
 * time: 0.46753406524658203
    11    -7.839799e+00     6.699483e-02
 * time: 0.5019879341125488
    12    -7.845236e+00     6.118431e-02
 * time: 0.5358829498291016
    13    -7.846562e+00     4.216923e-02
 * time: 0.5742161273956299
    14    -7.848235e+00     3.144008e-02
 * time: 0.6114299297332764
    15    -7.850677e+00     2.142225e-02
 * time: 0.6461060047149658
    16    -7.851997e+00     1.781180e-02
 * time: 0.6796009540557861
    17    -7.852452e+00     6.897497e-03
 * time: 0.714224100112915
    18    -7.852589e+00     3.823951e-03
 * time: 0.7476780414581299
    19    -7.852628e+00     2.867573e-03
 * time: 0.7834999561309814
    20    -7.852641e+00     1.471738e-03
 * time: 0.818295955657959
    21    -7.852645e+00     6.774609e-04
 * time: 0.8513889312744141
    22    -7.852646e+00     5.449233e-04
 * time: 0.8840620517730713
    23    -7.852647e+00     2.691163e-04
 * time: 0.9181230068206787
    24    -7.852647e+00     1.736887e-04
 * time: 0.9518449306488037
    25    -7.852647e+00     9.470015e-05
 * time: 0.9836061000823975
    26    -7.852647e+00     5.422476e-05
 * time: 1.0163159370422363
    27    -7.852647e+00     5.237609e-05
 * time: 1.048586130142212
    28    -7.852647e+00     2.292737e-05
 * time: 1.0814039707183838
    29    -7.852647e+00     1.119824e-05
 * time: 1.1156079769134521
    30    -7.852647e+00     5.171177e-06
 * time: 1.1503491401672363
    31    -7.852647e+00     3.022153e-06
 * time: 1.183779001235962
    32    -7.852647e+00     2.204431e-06
 * time: 1.2175700664520264
    33    -7.852647e+00     1.337301e-06
 * time: 1.2507951259613037
    34    -7.852647e+00     8.326351e-07
 * time: 1.283607006072998
    35    -7.852647e+00     5.370798e-07
 * time: 1.3170990943908691
    36    -7.852647e+00     2.981644e-07
 * time: 1.351619005203247
    37    -7.852647e+00     2.108745e-07
 * time: 1.385200023651123
    38    -7.852647e+00     1.225919e-07
 * time: 1.4177899360656738
    39    -7.852647e+00     6.459460e-08
 * time: 1.4510281085968018
    40    -7.852647e+00     2.919975e-08
 * time: 1.4831700325012207
    41    -7.852647e+00     1.092091e-08
 * time: 1.516571044921875
    42    -7.852647e+00     8.208210e-09
 * time: 1.5494189262390137
    43    -7.852647e+00     8.208202e-09
 * time: 1.6864359378814697

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   -7.846746081609                   -0.70    4.5

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   -7.852645862126                   -1.64
  2   -7.852646686730       -6.08       -3.70    695ms
  3   -7.852646686730      -13.25       -7.22    261ms

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|  = 8.951759867264316e-7
|ρ_newton - ρ_scfv| = 1.5880727604563485e-7
|ρ_newton - ρ_dm|   = 1.015388878899405e-9