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.397820596769                   -0.90    5.2   27.6ms
  2   -8.400252045122       -2.61       -1.74    1.0   31.7ms
  3   -8.400404101378       -3.82       -2.93    1.8   21.8ms
  4   -8.400427805065       -4.63       -2.94    3.0   25.4ms
  5   -8.400427968154       -6.79       -3.05    1.0   19.7ms
  6   -8.400428145948       -6.75       -4.90    1.0   19.4ms
  7   -8.400428151910       -8.22       -4.51    3.5   26.1ms
  8   -8.400428152199       -9.54       -5.84    1.0   19.9ms
  9   -8.400428152209      -11.00       -6.40    2.0   32.9ms

Potential-based SCF

scfres_scfv = DFTK.scf_potential_mixing(basis; tol);

n     Energy            log10(ΔE)   log10(Δρ)   α      Diag   Δtime
---   ---------------   ---------   ---------   ----   ----   ------
  1   -8.397844563200                   -0.90           5.2   28.3ms
  2   -8.400388159104       -2.59       -1.79   0.80    2.0   19.8ms
  3   -8.400423429394       -4.45       -3.03   0.80    1.0   16.9ms
  4   -8.400428100402       -5.33       -3.46   0.80    2.2   20.4ms
  5   -8.400428150358       -7.30       -4.72   0.80    1.5   18.2ms
  6   -8.400428152181       -8.74       -5.83   0.80    2.2   20.9ms
  7   -8.400428152209      -10.55       -6.35   0.80    3.0   33.1ms

Direct minimization

scfres_dm = direct_minimization(basis; tol);

┌ Warning: x_tol is deprecated. Use x_abstol or x_reltol instead. The provided value (-1) will be used as x_abstol.
└ @ Optim ~/.julia/packages/Optim/7krni/src/types.jl:110
┌ Warning: f_tol is deprecated. Use f_abstol or f_reltol instead. The provided value (-1) will be used as f_reltol.
└ @ Optim ~/.julia/packages/Optim/7krni/src/types.jl:120
n     Energy            log10(ΔE)   log10(Δρ)   Δtime
---   ---------------   ---------   ---------   ------
  1   +1.482753109331                   -1.04   59.6ms
  2   -1.366289553391        0.45       -0.69   33.5ms
  3   -4.502174741030        0.50       -0.40   44.7ms
  4   -6.446302867980        0.29       -0.48   53.6ms
  5   -6.617754234882       -0.77       -1.23   34.4ms
  6   -7.439549304426       -0.09       -1.20   33.5ms
  7   -8.038323085977       -0.22       -1.15   45.3ms
  8   -8.094664491871       -1.25       -1.56   33.2ms
  9   -8.177648314968       -1.08       -1.77   33.4ms
 10   -8.196593836126       -1.72       -1.81   41.8ms
 11   -8.222288802213       -1.59       -1.70   33.6ms
 12   -8.233336270703       -1.96       -1.75   45.8ms
 13   -8.307038291456       -1.13       -1.77   56.5ms
 14   -8.356700717290       -1.30       -2.16   33.6ms
 15   -8.373312189906       -1.78       -2.22   40.5ms
 16   -8.384507149319       -1.95       -2.16   33.9ms
 17   -8.393496540461       -2.05       -2.37   33.6ms
 18   -8.396123790485       -2.58       -2.80   33.3ms
 19   -8.398482673976       -2.63       -2.79   33.0ms
 20   -8.399162181885       -3.17       -2.98   33.3ms
 21   -8.399829063218       -3.18       -3.13   33.2ms
 22   -8.400140823394       -3.51       -3.45   41.4ms
 23   -8.400303498237       -3.79       -3.66   33.6ms
 24   -8.400368463504       -4.19       -3.90   33.6ms
 25   -8.400405211098       -4.43       -3.91   34.5ms
 26   -8.400420911859       -4.80       -4.43   33.3ms
 27   -8.400424549435       -5.44       -4.35   33.2ms
 28   -8.400426426282       -5.73       -4.82   40.4ms
 29   -8.400427464493       -5.98       -4.61   33.5ms
 30   -8.400427857441       -6.41       -5.14   33.2ms
 31   -8.400427993733       -6.87       -5.09   33.5ms
 32   -8.400428073946       -7.10       -5.38   33.4ms
 33   -8.400428114490       -7.39       -5.30   33.3ms
 34   -8.400428137299       -7.64       -6.35   33.2ms

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.397736608348                   -0.90    5.2   28.2ms

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.400427978837                   -1.79    619ms
  2   -8.400428152209       -6.76       -4.04    415ms
  3   -8.400428152209      -14.27       -7.82    147ms

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.606154694900307e-7
|ρ_newton - ρ_scfv| = 1.1577508342778339e-7
|ρ_newton - ρ_dm|   = 4.474441528100554e-6