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
tol = 1e-12
is_converged = DFTK.ScfConvergenceDensity(tol);

Density-based self-consistent field

scfres_scf = self_consistent_field(basis; is_converged);
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -7.846777262291                   -0.70    4.5
  2   -7.852310526801       -2.26       -1.53    1.0
  3   -7.852645905985       -3.47       -2.53    3.0
  4   -7.852646675666       -6.11       -3.35    2.8
  5   -7.852646685780       -8.00       -4.75    1.2
  6   -7.852646686726       -9.02       -5.23    4.2
  7   -7.852646686729      -11.44       -5.89    2.0
  8   -7.852646686730      -12.46       -7.49    2.0
  9   -7.852646686730      -14.75       -8.05    2.5
 10   -7.852646686730      -14.75       -8.97    1.2
 11   -7.852646686730   +  -15.05       -9.51    2.0
 12   -7.852646686730      -14.57       -8.81    1.0
 13   -7.852646686730   +  -14.35       -8.28    2.0
 14   -7.852646686730      -15.05       -9.06    1.2
 15   -7.852646686730   +    -Inf       -8.55    1.0
 16   -7.852646686730      -14.75       -8.24    2.0
 17   -7.852646686730   +  -14.57       -9.88    1.8
 18   -7.852646686730      -14.57       -9.89    2.0
 19   -7.852646686730   +  -14.75      -10.48    1.2
 20   -7.852646686730   +  -14.75      -10.46    1.0
 21   -7.852646686730   +  -14.75      -11.89    1.0
 22   -7.852646686730      -14.35      -11.28    1.0
 23   -7.852646686730   +  -15.05      -11.32    1.0
 24   -7.852646686730   +    -Inf      -11.48    1.0
 25   -7.852646686730   +    -Inf      -11.86    1.0
 26   -7.852646686730      -15.05      -13.45    2.0

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; is_converged);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag
---   ---------------   ---------   ---------   ----   ----
  1   -7.846842889337                   -0.70           4.5
  2   -7.852530500139       -2.25       -1.63   0.80    2.2
  3   -7.852637976478       -3.97       -2.72   0.80    1.0
  4   -7.852646355774       -5.08       -3.33   0.80    2.0
  5   -7.852646682776       -6.49       -4.23   0.80    1.8
  6   -7.852646686592       -8.42       -4.90   0.80    1.8
  7   -7.852646686722       -9.89       -5.78   0.80    1.8
  8   -7.852646686730      -11.11       -6.36   0.80    2.0
  9   -7.852646686730      -12.92       -7.25   0.80    1.2
 10   -7.852646686730      -13.88       -7.98   0.80    2.5
 11   -7.852646686730   +  -14.57       -8.66   0.80    1.8
 12   -7.852646686730      -15.05       -8.99   0.80    2.0
 13   -7.852646686730      -14.75      -10.20   0.80    1.2
 14   -7.852646686730   +  -14.27      -10.70   0.80    2.2
 15   -7.852646686730      -14.75      -11.19   0.80    1.2
 16   -7.852646686730      -14.45      -12.22   0.80    1.5

Direct minimization

scfres_dm = direct_minimization(basis; tol);
Iter     Function value   Gradient norm
     0     1.426441e+01     3.846956e+00
 * time: 0.01440286636352539
     1     1.292289e+00     1.863167e+00
 * time: 0.041684865951538086
     2    -1.880820e+00     2.121376e+00
 * time: 0.06896805763244629
     3    -3.609945e+00     1.843785e+00
 * time: 0.10801887512207031
     4    -5.177782e+00     1.499845e+00
 * time: 0.14759397506713867
     5    -6.721600e+00     7.658730e-01
 * time: 0.1866140365600586
     6    -6.795059e+00     9.442784e-01
 * time: 0.2139148712158203
     7    -7.524428e+00     4.985537e-01
 * time: 0.24129390716552734
     8    -7.692265e+00     3.946382e-01
 * time: 0.26845788955688477
     9    -7.775121e+00     1.142263e-01
 * time: 0.29541802406311035
    10    -7.807626e+00     1.177294e-01
 * time: 0.32247495651245117
    11    -7.825249e+00     8.087177e-02
 * time: 0.3495960235595703
    12    -7.838730e+00     5.654729e-02
 * time: 0.37686586380004883
    13    -7.847248e+00     3.967695e-02
 * time: 0.40372800827026367
    14    -7.849716e+00     2.783006e-02
 * time: 0.43066906929016113
    15    -7.850933e+00     2.038467e-02
 * time: 0.4578440189361572
    16    -7.851889e+00     1.997258e-02
 * time: 0.4850270748138428
    17    -7.852416e+00     1.278973e-02
 * time: 0.5118730068206787
    18    -7.852577e+00     4.805168e-03
 * time: 0.5387940406799316
    19    -7.852622e+00     3.132789e-03
 * time: 0.5656988620758057
    20    -7.852638e+00     1.650355e-03
 * time: 0.5927770137786865
    21    -7.852644e+00     9.781338e-04
 * time: 0.6198439598083496
    22    -7.852646e+00     6.153271e-04
 * time: 0.6473119258880615
    23    -7.852646e+00     3.496476e-04
 * time: 0.6742198467254639
    24    -7.852647e+00     1.850917e-04
 * time: 0.7011339664459229
    25    -7.852647e+00     9.859993e-05
 * time: 0.728126049041748
    26    -7.852647e+00     6.753300e-05
 * time: 0.7553088665008545
    27    -7.852647e+00     4.875409e-05
 * time: 0.7821569442749023
    28    -7.852647e+00     1.616915e-05
 * time: 0.8090970516204834
    29    -7.852647e+00     9.364682e-06
 * time: 0.8359298706054688
    30    -7.852647e+00     6.121000e-06
 * time: 0.8628778457641602
    31    -7.852647e+00     3.238488e-06
 * time: 0.8896679878234863
    32    -7.852647e+00     2.134854e-06
 * time: 0.9165499210357666
    33    -7.852647e+00     1.451816e-06
 * time: 0.9436378479003906
    34    -7.852647e+00     5.900540e-07
 * time: 0.9710679054260254
    35    -7.852647e+00     3.404332e-07
 * time: 0.9980978965759277
    36    -7.852647e+00     2.247658e-07
 * time: 1.0251379013061523
    37    -7.852647e+00     1.441273e-07
 * time: 1.0521109104156494
    38    -7.852647e+00     1.203774e-07
 * time: 1.079028844833374
    39    -7.852647e+00     6.529254e-08
 * time: 1.1058149337768555
    40    -7.852647e+00     2.673794e-08
 * time: 1.13277006149292
    41    -7.852647e+00     1.608416e-08
 * time: 1.1595900058746338
    42    -7.852647e+00     1.025327e-08
 * time: 1.1864609718322754
    43    -7.852647e+00     8.314678e-09
 * time: 1.2133948802947998
    44    -7.852647e+00     5.520029e-09
 * time: 1.2526710033416748
    45    -7.852647e+00     5.441449e-09
 * time: 1.303725004196167

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=1e-1);
n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
  1   -7.846845132465                   -0.70    4.8
  2   -7.852322193269       -2.26       -1.53    1.0

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(Δρ)
---   ---------------   ---------   ---------
  1   -7.852646686713                   -2.55
  2   -7.852646686730      -10.77       -6.02
  3   -7.852646686730   +  -15.05      -12.68

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|  = 3.843254934083239e-14
|ρ_newton - ρ_scfv| = 7.858717003474331e-13
|ρ_newton - ρ_dm|   = 6.626397645110732e-10