# 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.847793415581                   -0.70    4.8   39.1ms
2   -7.852411199040       -2.34       -1.53    1.0   17.1ms
3   -7.852610879115       -3.70       -2.54    1.0   16.8ms
4   -7.852645449757       -4.46       -2.77    2.2   21.5ms
5   -7.852646066618       -6.21       -2.87    1.0   16.9ms
6   -7.852646668063       -6.22       -3.84    1.0   16.9ms
7   -7.852646684876       -7.77       -4.54    1.2   17.9ms
8   -7.852646686696       -8.74       -5.09    1.8   20.0ms
9   -7.852646686727      -10.51       -5.83    1.2   18.3ms
10   -7.852646686728      -11.78       -5.66    2.0   20.8ms
11   -7.852646686730      -11.87       -6.15    1.0   17.7ms

## Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
1   -7.847731724603                   -0.70           4.5   29.7ms
2   -7.852560384301       -2.32       -1.62   0.80    2.0   20.2ms
3   -7.852640324683       -4.10       -2.71   0.80    1.0   16.0ms
4   -7.852646562329       -5.20       -3.35   0.80    2.0   19.8ms
5   -7.852646680911       -6.93       -4.39   0.80    1.2   17.1ms
6   -7.852646686602       -8.24       -4.84   0.80    2.0   20.6ms
7   -7.852646686723       -9.92       -5.48   0.80    1.2   16.9ms
8   -7.852646686730      -11.17       -6.55   0.80    1.8   18.5ms

## 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);
n     Energy            log10(ΔE)   log10(Δρ)   Δtime
---   ---------------   ---------   ---------   ------
1   +1.352591858994                   -1.01   44.5ms
2   -1.293552376780        0.42       -0.63   24.6ms
3   -3.620594405474        0.37       -0.44   32.5ms
4   -4.667042118362        0.02       -0.56   41.8ms
5   -6.477287009564        0.26       -0.68   33.8ms
6   -7.228194857748       -0.12       -1.10   33.1ms
7   -7.468689403252       -0.62       -1.48   25.1ms
8   -7.648392209908       -0.75       -1.96   24.9ms
9   -7.692041625438       -1.36       -2.19   24.9ms
10   -7.725418022551       -1.48       -2.28   24.9ms
11   -7.749622607035       -1.62       -2.22   24.9ms
12   -7.781044032048       -1.50       -2.13   25.0ms
13   -7.815746693196       -1.46       -1.90   25.2ms
14   -7.832841364537       -1.77       -2.36   25.0ms
15   -7.845681530423       -1.89       -2.38   25.0ms
16   -7.846110910922       -3.37       -2.86   33.1ms
17   -7.846189294969       -4.11       -2.36   49.3ms
18   -7.846264180461       -4.13       -2.27   41.7ms
19   -7.846380213303       -3.94       -2.51   42.2ms
20   -7.847050852519       -3.17       -2.80   25.3ms
21   -7.849588583978       -2.60       -3.46   25.1ms
22   -7.852512311485       -2.53       -3.39   25.0ms
23   -7.852580897595       -4.16       -3.85   24.9ms
24   -7.852622859810       -4.38       -4.09   25.0ms
25   -7.852638726549       -4.80       -4.34   24.8ms
26   -7.852642982334       -5.37       -4.48   24.9ms
27   -7.852644939676       -5.71       -5.21   24.9ms
28   -7.852646174423       -5.91       -4.71   24.9ms
29   -7.852646492833       -6.50       -4.91   24.8ms
30   -7.852646597775       -6.98       -5.47   24.7ms
31   -7.852646653981       -7.25       -6.07   24.8ms
32   -7.852646676008       -7.66       -5.75   24.8ms
33   -7.852646684324       -8.08       -5.86   24.8ms
34   -7.852646686102       -8.75       -6.54   24.8ms
35   -7.852646686519       -9.38       -6.69   24.8ms
36   -7.852646686661       -9.85       -6.66   24.8ms
37   -7.852646686699      -10.41       -7.25   24.8ms
38   -7.852646686717      -10.74       -7.32   24.8ms
39   -7.852646686726      -11.05       -7.37   24.8ms
40   -7.852646686729      -11.60       -7.56   24.7ms
41   -7.852646686730      -12.16       -8.03   24.7ms
42   -7.852646686730      -12.69       -8.25   24.8ms
43   -7.852646686730      -13.18       -8.22   24.8ms
44   -7.852646686730      -13.60       -8.56   24.9ms
45   -7.852646686730      -14.15       -9.04   24.7ms
46   -7.852646686730   +    -Inf       -9.52   24.7ms
47   -7.852646686730      -14.57       -9.46   24.8ms
48   -7.852646686730      -15.05       -9.68   49.0ms
49   -7.852646686730   +    -Inf       -9.65   41.1ms
50   -7.852646686730   +    -Inf       -9.73   25.0ms
┌ Warning: DM not converged.
└ @ DFTK ~/work/DFTK.jl/DFTK.jl/src/scf/scf_callbacks.jl:60

## 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.847760454769                   -0.70    4.8   30.3ms

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.852645844855                   -1.63    390ms
2   -7.852646686730       -6.07       -3.69    274ms
3   -7.852646686730      -13.22       -7.20    137ms

## 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|  = 7.251449160836048e-7
|ρ_newton - ρ_scfv| = 2.6679838064898556e-7
|ρ_newton - ρ_dm|   = 6.165357765969544e-10