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 AtomsBuilder
using DFTK
using LinearAlgebra
using PseudoPotentialData

pseudopotentials = PseudoFamily("dojo.nc.sr.pbesol.v0_4_1.standard.upf")
model = model_DFT(bulk(:Si); functionals=PBEsol(), pseudopotentials)
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   -8.397804294746                   -0.90    5.0   28.0ms
  2   -8.400236263940       -2.61       -1.73    1.0   20.6ms
  3   -8.400403544526       -3.78       -2.93    1.8   22.0ms
  4   -8.400427804522       -4.62       -2.92    3.2   26.9ms
  5   -8.400427956381       -6.82       -3.04    1.0   20.2ms
  6   -8.400428150139       -6.71       -4.82    1.0   20.3ms
  7   -8.400428155906       -8.24       -4.46    3.2   27.2ms
  8   -8.400428156254       -9.46       -5.13    1.0   21.1ms
  9   -8.400428156275      -10.68       -6.18    1.0   21.0ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397835747985                   -0.90           5.2   28.5ms
  2   -8.400383872644       -2.59       -1.78   0.80    2.0   19.6ms
  3   -8.400422528331       -4.41       -2.99   0.80    1.0   16.5ms
  4   -8.400428103260       -5.25       -3.42   0.80    2.2   21.2ms
  5   -8.400428152263       -7.31       -4.65   0.80    1.0   17.0ms
  6   -8.400428156252       -8.40       -5.43   0.80    2.8   22.6ms
  7   -8.400428156276      -10.62       -6.24   0.80    2.2   21.6ms

Direct minimization

scfres_dm = direct_minimization(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   Δtime
---   ---------------   ---------   ---------   ------
  1   +0.875313784281                   -1.11   54.2ms
  2   -1.507011198566        0.38       -0.73   30.7ms
  3   -3.988180372354        0.39       -0.46   45.0ms
  4   -4.656650066036       -0.17       -0.56   50.1ms
  5   -6.696928969909        0.31       -0.62   50.2ms
  6   -7.452578318315       -0.12       -0.96   50.1ms
  7   -7.713772444409       -0.58       -1.31   36.4ms
  8   -8.032409515056       -0.50       -1.70   56.5ms
  9   -8.115298047215       -1.08       -1.74   30.2ms
 10   -8.225126914766       -0.96       -1.86   29.9ms
 11   -8.288751488375       -1.20       -2.07   29.9ms
 12   -8.326773040912       -1.42       -2.00   29.5ms
 13   -8.362730530823       -1.44       -1.76   29.6ms
 14   -8.377107010263       -1.84       -2.44   29.7ms
 15   -8.389983472872       -1.89       -3.23   29.7ms
 16   -8.393568369042       -2.45       -2.82   29.9ms
 17   -8.398873032217       -2.28       -3.23   30.0ms
 18   -8.399187738841       -3.50       -2.97   29.8ms
 19   -8.400026244694       -3.08       -3.21   29.7ms
 20   -8.400087606907       -4.21       -3.06   29.9ms
 21   -8.400319032604       -3.64       -3.24   29.9ms
 22   -8.400325497192       -5.19       -3.69   29.8ms
 23   -8.400394943535       -4.16       -3.68   29.9ms
 24   -8.400407153440       -4.91       -3.97   40.1ms
 25   -8.400422250992       -4.82       -4.17   30.2ms
 26   -8.400424443066       -5.66       -4.14   30.2ms
 27   -8.400427133522       -5.57       -4.60   30.1ms
 28   -8.400427348977       -6.67       -4.59   29.9ms
 29   -8.400427938729       -6.23       -4.90   30.0ms
 30   -8.400428015238       -7.12       -5.18   30.2ms
 31   -8.400428120431       -6.98       -5.53   30.6ms
 32   -8.400428132801       -7.91       -5.16   30.1ms
 33   -8.400428150305       -7.76       -5.44   30.0ms
 34   -8.400428152448       -8.67       -5.63   30.1ms
 35   -8.400428155043       -8.59       -5.85   29.9ms
 36   -8.400428155386       -9.47       -5.96   30.1ms
 37   -8.400428156043       -9.18       -6.60   46.6ms

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   -8.397802690572                   -0.90    5.0   27.6ms

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   -8.400427945539                   -1.77    538ms
  2   -8.400428156277       -6.68       -3.97    402ms
  3   -8.400428156277      -14.15       -7.73    114ms

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|  = 1.1775454389753843e-6
|ρ_newton - ρ_scfv| = 3.334668414959783e-7
|ρ_newton - ρ_dm|   = 2.3348590710568234e-6