# 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.]]
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.846940822172                   -0.70    4.8
2   -7.852336526202       -2.27       -1.53    1.0   38.7ms
3   -7.852616947394       -3.55       -2.56    1.2   18.0ms
4   -7.852646013551       -4.54       -2.89    2.5   23.1ms
5   -7.852646509208       -6.30       -3.19    1.0   17.2ms
6   -7.852646680349       -6.77       -4.07    1.2   18.0ms
7   -7.852646686409       -8.22       -5.22    1.8   20.5ms
8   -7.852646686724       -9.50       -5.49    2.2   22.2ms
9   -7.852646686729      -11.29       -5.91    1.2   18.2ms
10   -7.852646686730      -12.20       -6.34    1.2   18.7ms

## Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
1   -7.846876925013                   -0.70           4.5
2   -7.852528531221       -2.25       -1.63   0.80    2.0   19.7ms
3   -7.852637238439       -3.96       -2.71   0.80    1.0   15.8ms
4   -7.852646514414       -5.03       -3.29   0.80    2.2   21.1ms
5   -7.852646674111       -6.80       -4.15   0.80    1.0   16.0ms
6   -7.852646686551       -7.91       -4.80   0.80    1.8   19.2ms
7   -7.852646686724       -9.76       -5.53   0.80    1.5   18.2ms
8   -7.852646686730      -11.29       -6.59   0.80    1.8   18.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.394672e+01     3.544632e+00
* time: 0.009553909301757812
1     1.385849e+00     1.794337e+00
* time: 0.02740788459777832
2    -1.185393e+00     2.190993e+00
* time: 0.04525589942932129
3    -3.617802e+00     1.854679e+00
* time: 0.07098102569580078
4    -4.977092e+00     1.641886e+00
* time: 0.09677791595458984
5    -6.758541e+00     8.713809e-01
* time: 0.12250304222106934
6    -6.848992e+00     1.647052e+00
* time: 0.1402430534362793
7    -7.472111e+00     1.300939e+00
* time: 0.15796804428100586
8    -7.486335e+00     2.475845e+00
* time: 0.17573809623718262
9    -7.621276e+00     2.282987e+00
* time: 0.19350790977478027
10    -7.639091e+00     1.951443e+00
* time: 0.2192380428314209
11    -7.646306e+00     1.845222e+00
* time: 0.2449660301208496
12    -7.667894e+00     1.432112e+00
* time: 0.270658016204834
13    -7.683349e+00     1.207230e+00
* time: 0.29653096199035645
14    -7.715872e+00     6.314722e-01
* time: 0.32306909561157227
15    -7.729993e+00     1.486678e+00
* time: 0.34087705612182617
16    -7.760390e+00     7.062637e-01
* time: 0.35875606536865234
17    -7.772689e+00     6.727363e-01
* time: 0.3845078945159912
18    -7.793503e+00     6.439758e-01
* time: 0.40236496925354004
19    -7.817270e+00     4.881956e-01
* time: 0.42012906074523926
20    -7.830550e+00     3.453619e-01
* time: 0.4378840923309326
21    -7.840620e+00     1.092378e-01
* time: 0.4557220935821533
22    -7.845939e+00     5.563220e-02
* time: 0.4920768737792969
23    -7.851151e+00     3.876657e-02
* time: 0.5178570747375488
24    -7.852210e+00     1.980195e-02
* time: 0.5382630825042725
25    -7.852483e+00     1.456035e-02
* time: 0.5561549663543701
26    -7.852611e+00     6.933888e-03
* time: 0.5739340782165527
27    -7.852638e+00     3.404811e-03
* time: 0.591731071472168
28    -7.852643e+00     2.032716e-03
* time: 0.6094698905944824
29    -7.852645e+00     1.024032e-03
* time: 0.6271910667419434
30    -7.852646e+00     4.012938e-04
* time: 0.6449639797210693
31    -7.852646e+00     3.457974e-04
* time: 0.6627469062805176
32    -7.852647e+00     2.393512e-04
* time: 0.6805789470672607
33    -7.852647e+00     1.441701e-04
* time: 0.6984009742736816
34    -7.852647e+00     5.113364e-05
* time: 0.7161760330200195
35    -7.852647e+00     4.267315e-05
* time: 0.7339038848876953
36    -7.852647e+00     3.439727e-05
* time: 0.7516989707946777
37    -7.852647e+00     2.289903e-05
* time: 0.7695050239562988
38    -7.852647e+00     7.868150e-06
* time: 0.7872960567474365
39    -7.852647e+00     5.090627e-06
* time: 0.8050329685211182
40    -7.852647e+00     3.425136e-06
* time: 0.8227748870849609
41    -7.852647e+00     1.263145e-06
* time: 0.8405370712280273
42    -7.852647e+00     9.913267e-07
* time: 0.8583879470825195
43    -7.852647e+00     6.137473e-07
* time: 0.876147985458374
44    -7.852647e+00     3.478600e-07
* time: 0.8939189910888672
45    -7.852647e+00     2.894924e-07
* time: 0.9116950035095215
46    -7.852647e+00     9.339685e-08
* time: 0.9294009208679199
47    -7.852647e+00     7.359717e-08
* time: 0.9471628665924072
48    -7.852647e+00     4.499015e-08
* time: 0.9649350643157959
49    -7.852647e+00     1.524489e-08
* time: 0.9826650619506836
50    -7.852647e+00     1.132541e-08
* time: 1.0005860328674316
51    -7.852647e+00     9.684323e-09
* time: 1.0342068672180176
52    -7.852647e+00     7.387184e-09
* time: 1.059967041015625

## 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.846883795509                   -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.852645919875                   -1.64
2   -7.852646686730       -6.12       -3.71    290ms
3   -7.852646686730      -13.30       -7.25    121ms

## 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|  = 4.3302220192353144e-7
|ρ_newton - ρ_scfv| = 2.567017146265263e-7
|ρ_newton - ρ_dm|   = 6.109249844477483e-10