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.846879465649                   -0.70    4.5
  2   -7.852323288389       -2.26       -1.53    1.0
  3   -7.852646227872       -3.49       -2.52    3.0
  4   -7.852646678225       -6.35       -3.37    2.5
  5   -7.852646686016       -8.11       -4.70    1.2
  6   -7.852646686726       -9.15       -5.31    3.5
  7   -7.852646686730      -11.44       -6.22    1.2
  8   -7.852646686730      -12.74       -7.63    2.5
  9   -7.852646686730   +  -15.05       -7.75    3.0
 10   -7.852646686730      -15.05       -8.70    2.0
 11   -7.852646686730   +    -Inf       -8.88    1.0
 12   -7.852646686730      -15.05      -10.61    2.0
 13   -7.852646686730   +  -15.05      -10.50    2.0
 14   -7.852646686730      -15.05      -11.65    1.5
 15   -7.852646686730   +  -15.05      -12.72    2.8

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; is_converged);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag
---   ---------------   ---------   ---------   ----   ----
  1   -7.846855597255                   -0.70           4.8
  2   -7.852525576107       -2.25       -1.64   0.80    2.0
  3   -7.852634641447       -3.96       -2.74   0.80    1.0
  4   -7.852646498626       -4.93       -3.22   0.80    2.2
  5   -7.852646672112       -6.76       -4.02   0.80    1.0
  6   -7.852646686240       -7.85       -4.69   0.80    1.5
  7   -7.852646686718       -9.32       -5.60   0.80    1.8
  8   -7.852646686729      -10.95       -6.76   0.80    1.8
  9   -7.852646686730      -12.43       -7.49   0.80    2.8
 10   -7.852646686730      -14.27       -7.96   0.80    2.0
 11   -7.852646686730      -15.05       -9.13   0.80    1.0
 12   -7.852646686730   +  -15.05       -9.50   0.80    2.8
 13   -7.852646686730   +    -Inf      -10.63   0.80    1.0
 14   -7.852646686730      -14.75      -10.62   0.80    2.8
 15   -7.852646686730   +  -14.45      -12.05   0.80    1.0

Direct minimization

scfres_dm = direct_minimization(basis; tol);
Iter     Function value   Gradient norm
     0     1.362834e+01     3.671831e+00
 * time: 0.013882160186767578
     1     1.041927e+00     2.087498e+00
 * time: 0.04104304313659668
     2    -1.838201e+00     2.461604e+00
 * time: 0.06888198852539062
     3    -3.889363e+00     2.130812e+00
 * time: 0.11160898208618164
     4    -5.223937e+00     1.673251e+00
 * time: 0.1499011516571045
     5    -6.797038e+00     1.250613e+00
 * time: 0.18706011772155762
     6    -7.347661e+00     4.849194e-01
 * time: 0.22365903854370117
     7    -7.532242e+00     1.615324e-01
 * time: 0.249406099319458
     8    -7.625804e+00     2.244724e-01
 * time: 0.27483415603637695
     9    -7.690670e+00     1.126843e-01
 * time: 0.30014514923095703
    10    -7.737975e+00     1.402996e-01
 * time: 0.32523202896118164
    11    -7.787923e+00     1.372769e-01
 * time: 0.3504910469055176
    12    -7.820711e+00     8.621622e-02
 * time: 0.3774580955505371
    13    -7.844840e+00     3.536448e-02
 * time: 0.4040100574493408
    14    -7.849804e+00     3.987997e-02
 * time: 0.4292171001434326
    15    -7.851889e+00     2.367369e-02
 * time: 0.4546520709991455
    16    -7.852439e+00     7.460525e-03
 * time: 0.48195505142211914
    17    -7.852581e+00     8.377179e-03
 * time: 0.5072400569915771
    18    -7.852610e+00     4.217875e-03
 * time: 0.5326650142669678
    19    -7.852631e+00     1.868054e-03
 * time: 0.5580241680145264
    20    -7.852642e+00     1.823901e-03
 * time: 0.5841591358184814
    21    -7.852645e+00     1.021985e-03
 * time: 0.6095001697540283
    22    -7.852646e+00     3.768700e-04
 * time: 0.6392440795898438
    23    -7.852647e+00     3.033489e-04
 * time: 0.6668670177459717
    24    -7.852647e+00     1.611514e-04
 * time: 0.6950380802154541
    25    -7.852647e+00     1.158842e-04
 * time: 0.7224111557006836
    26    -7.852647e+00     3.677708e-05
 * time: 0.7546300888061523
    27    -7.852647e+00     2.004913e-05
 * time: 0.7818961143493652
    28    -7.852647e+00     1.334430e-05
 * time: 0.8101541996002197
    29    -7.852647e+00     6.533916e-06
 * time: 0.8373889923095703
    30    -7.852647e+00     7.269142e-06
 * time: 0.8644170761108398
    31    -7.852647e+00     2.925066e-06
 * time: 0.9419059753417969
    32    -7.852647e+00     1.671245e-06
 * time: 0.969520092010498
    33    -7.852647e+00     1.091535e-06
 * time: 0.9981181621551514
    34    -7.852647e+00     6.211008e-07
 * time: 1.0241641998291016
    35    -7.852647e+00     5.391051e-07
 * time: 1.0497171878814697
    36    -7.852647e+00     3.728997e-07
 * time: 1.0750751495361328
    37    -7.852647e+00     9.735543e-08
 * time: 1.102931022644043
    38    -7.852647e+00     6.115073e-08
 * time: 1.1295599937438965
    39    -7.852647e+00     5.383324e-08
 * time: 1.1558520793914795
    40    -7.852647e+00     2.491115e-08
 * time: 1.1948490142822266
    41    -7.852647e+00     1.765566e-08
 * time: 1.2218589782714844
    42    -7.852647e+00     1.404257e-08
 * time: 1.2488081455230713
    43    -7.852647e+00     1.372267e-08
 * time: 1.299713134765625
    44    -7.852647e+00     1.372267e-08
 * time: 1.4092340469360352

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.846791291173                   -0.70    4.5
  2   -7.852299197883       -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.852646686579                   -2.54
  2   -7.852646686730       -9.82       -5.85
  3   -7.852646686730   +    -Inf      -11.61
  4   -7.852646686730   +    -Inf      -15.97

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|  = 1.3866702276178619e-13
|ρ_newton - ρ_scfv| = 2.045522491196992e-12
|ρ_newton - ρ_dm|   = 5.411763012138437e-10