# 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

In [1]:
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
tol = 1e-12
is_converged = DFTK.ScfConvergenceDensity(tol);


## Density-based self-consistent field¶

In [2]:
scfres_scf = self_consistent_field(basis; is_converged);

n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
1   -7.846832249821                   -0.70    4.8
2   -7.852321060476       -2.26       -1.53    1.0
3   -7.852646061824       -3.49       -2.52    3.2
4   -7.852646678176       -6.21       -3.37    2.5
5   -7.852646686299       -8.09       -4.73    1.5
6   -7.852646686727       -9.37       -5.44    3.5
7   -7.852646686730      -11.54       -6.27    1.2
8   -7.852646686730      -13.01       -7.71    2.2
9   -7.852646686730      -14.57       -7.74    2.2
10   -7.852646686730      -15.05       -8.71    2.5
11   -7.852646686730   +    -Inf       -8.94    1.0
12   -7.852646686730   +  -15.05      -10.67    2.2
13   -7.852646686730   +    -Inf       -9.47    1.0
14   -7.852646686730      -15.05       -9.24    1.0
15   -7.852646686730      -14.75       -9.10    1.0
16   -7.852646686730   +  -14.35       -8.50    2.0
17   -7.852646686730   +  -15.05       -9.80    1.5
18   -7.852646686730      -14.75      -10.57    1.8
19   -7.852646686730      -14.75      -10.95    1.8
20   -7.852646686730   +  -15.05      -10.82    1.0
21   -7.852646686730   +  -15.05      -11.71    1.0
22   -7.852646686730      -14.75      -12.45    2.0


## Potential-based SCF¶

In [3]:
scfres_scfv = DFTK.scf_potential_mixing(basis; is_converged);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag
---   ---------------   ---------   ---------   ----   ----
1   -7.846826529307                   -0.70           4.8
2   -7.852525112695       -2.24       -1.63   0.80    2.0
3   -7.852635705254       -3.96       -2.72   0.80    1.0
4   -7.852646446816       -4.97       -3.27   0.80    2.2
5   -7.852646672886       -6.65       -4.10   0.80    1.2
6   -7.852646686341       -7.87       -4.72   0.80    1.5
7   -7.852646686721       -9.42       -5.57   0.80    1.8
8   -7.852646686729      -11.08       -6.98   0.80    1.8
9   -7.852646686730      -12.41       -7.57   0.80    3.2
10   -7.852646686730      -14.75       -7.99   0.80    1.8
11   -7.852646686730      -15.05       -9.12   0.80    1.0
12   -7.852646686730      -14.75       -9.60   0.80    3.2
13   -7.852646686730   +  -15.05      -10.10   0.80    1.0
14   -7.852646686730   +  -15.05      -10.82   0.80    1.2
15   -7.852646686730   +    -Inf      -11.97   0.80    2.2
16   -7.852646686730   +    -Inf      -12.64   0.80    2.5


## Direct minimization¶

In [4]:
scfres_dm = direct_minimization(basis; tol);

Iter     Function value   Gradient norm
0     1.403903e+01     3.503994e+00
* time: 0.44902586936950684
1     1.342235e+00     1.778317e+00
* time: 0.651695966720581
2    -1.828507e+00     1.896481e+00
* time: 0.6792130470275879
3    -3.754821e+00     1.820452e+00
* time: 0.7940008640289307
4    -5.188438e+00     1.584808e+00
* time: 0.8335838317871094
5    -6.918980e+00     8.684053e-01
* time: 0.8725728988647461
6    -6.992047e+00     9.572685e-01
* time: 0.899763822555542
7    -7.689852e+00     7.397375e-01
* time: 0.9269559383392334
8    -7.759766e+00     8.800055e-01
* time: 0.9543449878692627
9    -7.778897e+00     7.374077e-01
* time: 0.9812648296356201
10    -7.785095e+00     1.073822e+00
* time: 1.0081868171691895
11    -7.798981e+00     8.160733e-01
* time: 1.0350289344787598
12    -7.819909e+00     4.539929e-01
* time: 1.0741620063781738
13    -7.824708e+00     6.909547e-01
* time: 1.1010818481445312
14    -7.843710e+00     2.596031e-01
* time: 1.1279468536376953
15    -7.846046e+00     2.959001e-01
* time: 1.15488600730896
16    -7.851371e+00     2.329684e-02
* time: 1.1817529201507568
17    -7.852230e+00     8.761840e-03
* time: 1.2086799144744873
18    -7.852569e+00     6.798428e-03
* time: 1.23555588722229
19    -7.852604e+00     4.930048e-03
* time: 1.262833833694458
20    -7.852625e+00     2.813465e-03
* time: 1.289719820022583
21    -7.852641e+00     1.634270e-03
* time: 1.3174378871917725
22    -7.852645e+00     9.289906e-04
* time: 1.3442578315734863
23    -7.852646e+00     4.364427e-04
* time: 1.3711888790130615
24    -7.852646e+00     3.232755e-04
* time: 1.3981859683990479
25    -7.852647e+00     1.153466e-04
* time: 1.4249730110168457
26    -7.852647e+00     6.613080e-05
* time: 1.451746940612793
27    -7.852647e+00     4.115242e-05
* time: 1.4787218570709229
28    -7.852647e+00     2.545155e-05
* time: 1.505707025527954
29    -7.852647e+00     1.371649e-05
* time: 1.5330719947814941
30    -7.852647e+00     8.255275e-06
* time: 1.5598359107971191
31    -7.852647e+00     3.840623e-06
* time: 1.586789846420288
32    -7.852647e+00     2.238605e-06
* time: 1.6139118671417236
33    -7.852647e+00     1.239133e-06
* time: 1.6414918899536133
34    -7.852647e+00     6.781613e-07
* time: 1.6685888767242432
35    -7.852647e+00     3.939354e-07
* time: 1.6955578327178955
36    -7.852647e+00     2.162569e-07
* time: 1.7229039669036865
37    -7.852647e+00     1.289076e-07
* time: 1.7498118877410889
38    -7.852647e+00     6.530468e-08
* time: 1.7768030166625977
39    -7.852647e+00     3.634828e-08
* time: 1.8040690422058105
40    -7.852647e+00     2.074358e-08
* time: 1.831233024597168
41    -7.852647e+00     1.169768e-08
* time: 1.858158826828003
42    -7.852647e+00     7.731757e-09
* time: 1.8854949474334717
43    -7.852647e+00     4.208548e-09
* time: 1.913322925567627
44    -7.852647e+00     4.207983e-09
* time: 1.988473892211914


## 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.

In [5]:
scfres_start = self_consistent_field(basis; tol=1e-1);

n     Energy            log10(ΔE)   log10(Δρ)   Diag
---   ---------------   ---------   ---------   ----
1   -7.846843511052                   -0.70    4.5
2   -7.852322679588       -2.26       -1.53    1.0


Remove the virtual orbitals (which Newton cannot treat yet)

In [6]:
ψ = DFTK.select_occupied_orbitals(basis, scfres_start.ψ, scfres_start.occupation).ψ
scfres_newton = newton(basis, ψ; tol);

n     Energy            log10(ΔE)   log10(Δρ)
---   ---------------   ---------   ---------
1   -7.852646686713                   -2.54
2   -7.852646686730      -10.77       -6.03
3   -7.852646686730   +  -15.05      -12.64


## Comparison of results¶

In [7]:
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|  = 9.899854352494933e-14
|ρ_newton - ρ_scfv| = 2.1712353697319453e-13
|ρ_newton - ρ_dm|   = 2.420188489527844e-10