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

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