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.846825135780                   -0.70    4.5
  2   -7.852318901218       -2.26       -1.53    1.0   22.6ms
  3   -7.852615300516       -3.53       -2.55    1.5   24.2ms
  4   -7.852645992516       -4.51       -2.88    2.8   31.3ms
  5   -7.852646507618       -6.29       -3.18    1.0   22.3ms
  6   -7.852646679643       -6.76       -4.02    1.0   23.0ms
  7   -7.852646686349       -8.17       -5.12    1.5   26.2ms
  8   -7.852646686721       -9.43       -5.36    2.2   29.1ms
  9   -7.852646686730      -11.05       -6.59    1.0   22.5ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -7.846922023684                   -0.70           4.8
  2   -7.852526805486       -2.25       -1.64   0.80    2.0   26.1ms
  3   -7.852635967438       -3.96       -2.72   0.80    1.0   21.8ms
  4   -7.852646502942       -4.98       -3.27   0.80    2.2   28.0ms
  5   -7.852646672802       -6.77       -4.09   0.80    1.0   21.2ms
  6   -7.852646686445       -7.87       -4.75   0.80    1.8   26.2ms
  7   -7.852646686725       -9.55       -5.49   0.80    1.8   25.4ms
  8   -7.852646686729      -11.37       -6.67   0.80    1.5   24.6ms

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.387706e+01     3.719426e+00
 * time: 0.012444019317626953
     1     1.269083e+00     1.964536e+00
 * time: 0.03523516654968262
     2    -1.427284e+00     2.411262e+00
 * time: 0.05799508094787598
     3    -3.841403e+00     1.850400e+00
 * time: 0.09041309356689453
     4    -5.118125e+00     1.510901e+00
 * time: 0.12261819839477539
     5    -6.799177e+00     8.163126e-01
 * time: 0.15502214431762695
     6    -7.446443e+00     2.687974e-01
 * time: 0.1872570514678955
     7    -7.674390e+00     1.338245e-01
 * time: 0.20988798141479492
     8    -7.769253e+00     1.214316e-01
 * time: 0.23241114616394043
     9    -7.814742e+00     7.927568e-02
 * time: 0.25545215606689453
    10    -7.837702e+00     5.893711e-02
 * time: 0.27816009521484375
    11    -7.846280e+00     5.074035e-02
 * time: 0.300616979598999
    12    -7.850021e+00     2.384647e-02
 * time: 0.3230421543121338
    13    -7.851658e+00     1.389295e-02
 * time: 0.34543418884277344
    14    -7.852259e+00     1.072655e-02
 * time: 0.3678281307220459
    15    -7.852478e+00     6.281028e-03
 * time: 0.3905150890350342
    16    -7.852581e+00     3.552732e-03
 * time: 0.41290998458862305
    17    -7.852629e+00     1.890666e-03
 * time: 0.4605081081390381
    18    -7.852642e+00     1.009485e-03
 * time: 0.48296213150024414
    19    -7.852645e+00     5.664085e-04
 * time: 0.5055961608886719
    20    -7.852646e+00     3.366794e-04
 * time: 0.5278310775756836
    21    -7.852646e+00     2.018730e-04
 * time: 0.5503361225128174
    22    -7.852647e+00     1.489106e-04
 * time: 0.5725641250610352
    23    -7.852647e+00     8.637588e-05
 * time: 0.5950040817260742
    24    -7.852647e+00     5.915399e-05
 * time: 0.617603063583374
    25    -7.852647e+00     2.230207e-05
 * time: 0.6397960186004639
    26    -7.852647e+00     1.156498e-05
 * time: 0.662121057510376
    27    -7.852647e+00     8.198795e-06
 * time: 0.6845271587371826
    28    -7.852647e+00     5.502192e-06
 * time: 0.7071561813354492
    29    -7.852647e+00     2.023899e-06
 * time: 0.7294771671295166
    30    -7.852647e+00     1.449614e-06
 * time: 0.7516160011291504
    31    -7.852647e+00     7.676239e-07
 * time: 0.7739729881286621
    32    -7.852647e+00     4.320281e-07
 * time: 0.7965030670166016
    33    -7.852647e+00     2.713169e-07
 * time: 0.8190250396728516
    34    -7.852647e+00     1.840549e-07
 * time: 0.841616153717041
    35    -7.852647e+00     1.400749e-07
 * time: 0.8643391132354736
    36    -7.852647e+00     7.835382e-08
 * time: 0.8871190547943115
    37    -7.852647e+00     5.224560e-08
 * time: 0.9097461700439453
    38    -7.852647e+00     3.095855e-08
 * time: 0.9320921897888184
    39    -7.852647e+00     1.739628e-08
 * time: 0.9643402099609375
    40    -7.852647e+00     1.589284e-08
 * time: 1.0073161125183105
    41    -7.852647e+00     1.025110e-08
 * time: 1.0300660133361816
    42    -7.852647e+00     5.285366e-09
 * time: 1.0526010990142822
    43    -7.852647e+00     4.828327e-09
 * time: 1.084954023361206

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.846885606838                   -0.70    4.8

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.852645921498                   -1.64
  2   -7.852646686730       -6.12       -3.71    392ms
  3   -7.852646686730      -13.34       -7.25    141ms

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|  = 5.18209311518797e-7
|ρ_newton - ρ_scfv| = 3.6672919984490966e-7
|ρ_newton - ρ_dm|   = 7.897882526685183e-10