# 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.847813849318                   -0.70    4.5   29.5ms
2   -7.852412540530       -2.34       -1.53    1.0   17.0ms
3   -7.852610829942       -3.70       -2.55    1.0   16.7ms
4   -7.852645464017       -4.46       -2.78    2.2   21.4ms
5   -7.852646095994       -6.20       -2.88    1.0   16.9ms
6   -7.852646666224       -6.24       -3.83    1.0   17.2ms
7   -7.852646684709       -7.73       -4.60    1.2   17.9ms
8   -7.852646686686       -8.70       -5.03    2.0   21.0ms
9   -7.852646686726      -10.41       -5.78    1.2   18.1ms
10   -7.852646686728      -11.57       -5.64    1.8   20.6ms
11   -7.852646686730      -11.86       -6.14    1.0   29.8ms

## Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
1   -7.847760240034                   -0.70           5.0   30.5ms
2   -7.852559190284       -2.32       -1.62   0.80    2.0   19.7ms
3   -7.852641752770       -4.08       -2.68   0.80    1.0   15.7ms
4   -7.852646509485       -5.32       -3.43   0.80    1.8   18.6ms
5   -7.852646681417       -6.76       -4.39   0.80    2.0   19.0ms
6   -7.852646686614       -8.28       -4.86   0.80    2.5   21.7ms
7   -7.852646686725       -9.95       -5.52   0.80    1.0   15.9ms
8   -7.852646686729      -11.38       -6.64   0.80    1.8   18.3ms

## 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.289714437899                   -1.04   44.4ms
2   -1.354926437320        0.42       -0.70   24.8ms
3   -3.656180602273        0.36       -0.47   32.8ms
4   -4.712803391798        0.02       -0.58   33.0ms
5   -6.607905445474        0.28       -0.68   32.8ms
6   -7.333718723848       -0.14       -1.12   32.8ms
7   -7.644586616196       -0.51       -1.41   24.8ms
8   -7.753806486603       -0.96       -1.86   24.7ms
9   -7.797048807441       -1.36       -2.15   24.8ms
10   -7.819300213693       -1.65       -2.14   24.7ms
11   -7.835211614059       -1.80       -2.18   24.9ms
12   -7.842470279422       -2.14       -2.41   24.7ms
13   -7.846763713019       -2.37       -2.60   24.7ms
14   -7.848553883071       -2.75       -3.10   24.7ms
15   -7.850565083726       -2.70       -3.21   24.8ms
16   -7.851866136754       -2.89       -3.36   24.6ms
17   -7.852403449189       -3.27       -3.25   24.7ms
18   -7.852570070248       -3.78       -3.31   24.7ms
19   -7.852621400661       -4.29       -4.12   24.9ms
20   -7.852637708265       -4.79       -4.13   24.7ms
21   -7.852643623282       -5.23       -4.31   24.7ms
22   -7.852645665172       -5.69       -4.89   24.7ms
23   -7.852646287425       -6.21       -5.05   24.9ms
24   -7.852646527888       -6.62       -5.12   24.8ms
25   -7.852646623050       -7.02       -5.59   24.7ms
26   -7.852646659174       -7.44       -5.58   24.7ms
27   -7.852646677709       -7.73       -5.78   25.0ms
28   -7.852646684627       -8.16       -6.09   24.7ms
29   -7.852646686014       -8.86       -6.44   24.8ms
30   -7.852646686462       -9.35       -6.64   24.8ms
31   -7.852646686633       -9.77       -6.80   24.8ms
32   -7.852646686689      -10.25       -7.08   32.2ms
33   -7.852646686717      -10.55       -7.20   25.5ms
34   -7.852646686726      -11.03       -7.69   24.9ms
35   -7.852646686729      -11.59       -7.87   24.7ms
36   -7.852646686729      -12.08       -7.84   24.7ms
37   -7.852646686730      -12.57       -8.05   24.7ms
38   -7.852646686730      -13.01       -8.34   31.4ms
39   -7.852646686730      -13.60       -8.87   24.8ms
40   -7.852646686730      -14.15       -8.98   24.9ms
41   -7.852646686730      -14.27       -8.98   24.8ms
42   -7.852646686730      -14.45       -9.54   24.8ms
┌ 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.847774083254                   -0.70    4.8   39.0ms

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.852645844027                   -1.63    383ms
2   -7.852646686730       -6.07       -3.69    289ms
3   -7.852646686730      -13.21       -7.20    119ms

## 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.246569562664095e-7
|ρ_newton - ρ_scfv| = 6.003049770832736e-7
|ρ_newton - ρ_dm|   = 2.052852744478766e-9