Exact exchange and Hybrid DFT

This is a quick sketch how to run a simple hybrid DFT calculation using DFTK.

Preliminary implementation

Hybrid-DFT is not yet performance-optimised and has only seen rudimentary testing so far. Further only Gamma point calculations are supported and the interface is not yet considered stable. We appreciate any issues, bug reports or PRs.

using AtomsBuilder
using DFTK
using PseudoPotentialData

pseudopotentials = PseudoFamily("dojo.nc.sr.pbe.v0_5.stringent.upf")
system = bulk(:Si)
FlexibleSystem(Si₂, periodicity = TTT):
    cell_vectors      : [       0    2.715    2.715;
                            2.715        0    2.715;
                            2.715    2.715        0]u"Å"

    Atom(Si, [       0,        0,        0]u"Å")
    Atom(Si, [  1.3575,   1.3575,   1.3575]u"Å")

First perform a PBE calculation to get a good starting point

model  = model_DFT(system; pseudopotentials, functionals=PBE())
basis  = PlaneWaveBasis(model; Ecut=15, kgrid=[1, 1, 1])
scfres = self_consistent_field(basis; tol=1e-6);
nothing  # hide
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -7.824466518862                   -0.57   11.0    100ms
  2   -7.834472412149       -2.00       -1.44    1.0   60.0ms
  3   -7.835385751248       -3.04       -2.28    3.0   75.4ms
  4   -7.835418297634       -4.49       -3.05    2.0   66.5ms
  5   -7.835419852096       -5.81       -3.84    3.0   69.1ms
  6   -7.835419910297       -7.24       -4.49    3.0   69.3ms
  7   -7.835419912416       -8.67       -4.98    3.0   71.6ms
  8   -7.835419912543       -9.90       -5.56    2.0   68.0ms
  9   -7.835419912558      -10.81       -6.33    4.0   72.0ms

Then run PBE0, see also PBE0 and HybridFunctional for more documentation.

model  = model_DFT(system; pseudopotentials, functionals=PBE0())
basis  = PlaneWaveBasis(model; basis.Ecut, basis.kgrid)
scfres = self_consistent_field(basis; tol=1e-4, scfres.ρ, scfres.ψ, scfres.occupation,
                               scfres.eigenvalues, solver=DFTK.scf_damping_solver());
nothing  # hide
n     Energy            log10(ΔE)   log10(Δρ)   Diag   Δtime 
---   ---------------   ---------   ---------   ----   ------
  1   -7.715128088935                   -2.20    1.0    4.09s
  2   -7.715171953249       -4.36       -2.77    1.0    419ms
  3   -7.715185180882       -4.88       -3.04    6.0    137ms
  4   -7.715185658949       -6.32       -3.55    5.0    128ms
  5   -7.715185690833       -7.50       -4.11    5.0    189ms