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.]]
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 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.847754733182                   -0.70    4.5   29.1ms
  2   -7.852400871343       -2.33       -1.53    1.0   17.3ms
  3   -7.852609952983       -3.68       -2.54    1.2   17.7ms
  4   -7.852645443883       -4.45       -2.77    2.2   21.6ms
  5   -7.852646077301       -6.20       -2.87    1.2   18.0ms
  6   -7.852646666286       -6.23       -3.83    1.0   17.1ms
  7   -7.852646684960       -7.73       -4.65    1.5   19.1ms
  8   -7.852646686694       -8.76       -5.06    1.8   19.9ms
  9   -7.852646686726      -10.50       -5.72    1.0   17.8ms
 10   -7.852646686728      -11.59       -5.66    1.5   20.1ms
 11   -7.852646686730      -11.87       -6.24    1.0   34.6ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);
n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -7.847780929274                   -0.70           4.8   31.1ms
  2   -7.852560337567       -2.32       -1.62   0.80    2.0   20.5ms
  3   -7.852641769653       -4.09       -2.68   0.80    1.0   16.4ms
  4   -7.852646527308       -5.32       -3.41   0.80    1.8   19.5ms
  5   -7.852646681635       -6.81       -4.40   0.80    1.5   18.6ms
  6   -7.852646686624       -8.30       -4.84   0.80    2.5   22.5ms
  7   -7.852646686726       -9.99       -5.52   0.80    1.2   17.3ms
  8   -7.852646686730      -11.47       -6.56   0.80    1.5   19.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.362678295933                   -0.97   45.8ms
  2   -1.203990905250        0.41       -0.62   25.5ms
  3   -3.543812684930        0.37       -0.43   33.7ms
  4   -4.836743009589        0.11       -0.57   33.6ms
  5   -6.687730082863        0.27       -0.67   33.2ms
  6   -6.758871660422       -1.15       -1.57   25.2ms
  7   -7.460227505159       -0.15       -1.41   25.0ms
  8   -7.552217784959       -1.04       -1.85   25.0ms
  9   -7.710400228843       -0.80       -1.85   25.0ms
 10   -7.751415011380       -1.39       -2.18   24.9ms
 11   -7.778758673453       -1.56       -2.35   24.8ms
 12   -7.802674537301       -1.62       -2.36   24.7ms
 13   -7.830896915392       -1.55       -2.18   24.7ms
 14   -7.843221345422       -1.91       -2.53   24.8ms
 15   -7.848336473864       -2.29       -2.26   24.9ms
 16   -7.851691769535       -2.47       -2.58   24.8ms
 17   -7.852282190995       -3.23       -3.37   24.8ms
 18   -7.852450603248       -3.77       -3.80   24.8ms
 19   -7.852568707405       -3.93       -3.38   24.8ms
 20   -7.852622392684       -4.27       -3.63   25.2ms
 21   -7.852635062012       -4.90       -3.66   25.3ms
 22   -7.852641809902       -5.17       -4.13   25.2ms
 23   -7.852645429596       -5.44       -4.58   25.4ms
 24   -7.852646348164       -6.04       -4.68   25.3ms
 25   -7.852646585798       -6.62       -4.84   25.0ms
 26   -7.852646657574       -7.14       -5.31   25.1ms
 27   -7.852646677416       -7.70       -5.27   24.9ms
 28   -7.852646683442       -8.22       -5.94   24.8ms
 29   -7.852646685523       -8.68       -6.06   25.1ms
 30   -7.852646686263       -9.13       -6.18   25.1ms
 31   -7.852646686590       -9.49       -6.54   25.0ms
 32   -7.852646686686      -10.02       -6.60   24.8ms
 33   -7.852646686716      -10.52       -6.91   24.9ms
 34   -7.852646686726      -11.02       -7.29   25.3ms
 35   -7.852646686728      -11.56       -7.46   25.1ms
 36   -7.852646686729      -12.03       -7.65   25.1ms
 37   -7.852646686730      -12.51       -8.02   25.1ms
 38   -7.852646686730      -13.05       -7.97   25.0ms
 39   -7.852646686730      -13.43       -8.25   24.9ms
 40   -7.852646686730      -13.80       -9.02   25.7ms
 41   -7.852646686730      -13.91       -8.80   34.3ms
 42   -7.852646686730      -15.05       -9.18   30.3ms
 43   -7.852646686730      -15.05       -9.37   25.4ms
 44   -7.852646686730      -14.75       -9.50   25.2ms
 45   -7.852646686730      -15.05       -9.62   25.0ms
 46   -7.852646686730      -15.05      -10.20   41.0ms
 47   -7.852646686730   +    -Inf      -10.06   32.9ms
 48   -7.852646686730      -14.75      -10.08   65.1ms
 49   -7.852646686730   +    -Inf      -10.03   73.7ms
 50   -7.852646686730   +    -Inf      -10.01   73.4ms
┌ 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.847764830578                   -0.70    4.8   43.1ms

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.852645831532                   -1.63    376ms
  2   -7.852646686730       -6.07       -3.68    274ms
  3   -7.852646686730      -13.21       -7.19    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|  = 5.819742373427392e-7
|ρ_newton - ρ_scfv| = 4.557405065127045e-7
|ρ_newton - ρ_dm|   = 7.082803768398802e-10