[WIP] Hybrid functionals

.github/workflows/CompatHelper.yaml

@@ -6,6 +6,7 @@ on:
 
 jobs:
   CompatHelper:
+    if: github.repository_owner = "JuliaMolSim"
     runs-on: ubuntu-latest
     steps:
       - uses: julia-actions/setup-julia@latest

src/DFTK.jl

@@ -89,6 +89,7 @@ export Hamiltonian
 export HamiltonianBlock
 export energy_hamiltonian
 export Kinetic
+export ExactExchange
 export ExternalFromFourier
 export ExternalFromReal
 export AtomicLocal
@@ -132,7 +133,7 @@ include("eigen/preconditioners.jl")
 include("eigen/diag.jl")
 
 export model_atomic
-export model_DFT, model_PBE, model_LDA, model_SCAN
+export model_DFT, model_PBE, model_LDA, model_SCAN, model_PBE0, model_HF
 include("standard_models.jl")
 
 export KerkerMixing, KerkerDosMixing, SimpleMixing, DielectricMixing

src/DispatchFunctional.jl

@@ -12,12 +12,19 @@ function LibxcFunctional(identifier::Symbol)
     @assert fun.kind   in (:exchange, :correlation, :exchange_correlation)
     kind = Dict(:exchange => :x, :correlation => :c, :exchange_correlation => :xc)[fun.kind]
 
-    @assert fun.family in (:lda, :gga, :mgga)  # Hybrids not supported yet.
-    if fun.family == :mgga && Libxc.needs_laplacian(fun)
-        family = :mggal
+    @assert fun.family in (:lda, :gga, :mgga, :hyb_lda, :hyb_gga, :hyb_mgga)
+
+    # Libxc maintains the distinction between hybrid and non-hybrid equivalents,
+    # but this distinction is not relevant for the functional interface
+    if startswith(string(fun.family), "hyb")
+        family = Symbol(string(fun.family)[5:end])
     else
         family = fun.family
     end
+
+    if family == :mgga && Libxc.needs_laplacian(fun)
+        family = :mggal
+    end
     LibxcFunctional{family,kind}(identifier)
 end
 
@@ -154,3 +161,10 @@ for fun in (:potential_terms, :kernel_terms)
         end
     end
 end
+
+# TODO This is hackish for now until Libxc has fully picked up the DftFunctionals.jl interface
+exx_coefficient(::Functional{:lda})      = nothing
+exx_coefficient(::Functional{:gga})      = nothing
+exx_coefficient(::Functional{:mgga})     = nothing
+exx_coefficient(fun::DispatchFunctional) = exx_coefficient(fun.inner)
+exx_coefficient(fun::LibxcFunctional)    = Libxc.Functional(fun.identifier).exx_coefficient

src/fft.jl

@@ -25,7 +25,7 @@ function ifft!(f_real::AbstractArray3, basis::PlaneWaveBasis, f_fourier::Abstrac
     f_real .*= basis.ifft_normalization
 end
 function ifft!(f_real::AbstractArray3, basis::PlaneWaveBasis,
-                 kpt::Kpoint, f_fourier::AbstractVector; normalize=true)
+               kpt::Kpoint, f_fourier::AbstractVector; normalize=true)
     @assert length(f_fourier) == length(kpt.mapping)
     @assert size(f_real) == basis.fft_size
 

src/standard_models.jl

@@ -33,8 +33,15 @@ function model_DFT(lattice::AbstractMatrix,
                    xc::Xc;
                    extra_terms=[], kwargs...)
     model_name = isempty(xc.functionals) ? "rHF" : join(string.(xc.functionals), "+")
+    exx = [ExactExchange(; scaling_factor=exx_coefficient(f))
+           for f in xc.functionals if !isnothing(exx_coefficient(f))]
+    if !isempty(exx)
+        @assert length(exx) == 1
+        @warn "Exact exchange in DFTK is hardly optimised and not yet production-ready."
+    end
+
     model_atomic(lattice, atoms, positions;
-                 extra_terms=[Hartree(), xc, extra_terms...], model_name, kwargs...)
+                 extra_terms=[Hartree(), xc, exx..., extra_terms...], model_name, kwargs...)
 end
 function model_DFT(lattice::AbstractMatrix,
                    atoms::Vector{<:Element},
@@ -75,8 +82,31 @@ function model_SCAN(lattice::AbstractMatrix, atoms::Vector{<:Element},
 end
 
 
+"""
+Build an PBE0 model from the specified atoms.
+DOI:10.1103/PhysRevLett.77.3865
+"""
+function model_PBE0(lattice::AbstractMatrix, atoms::Vector{<:Element},
+                   positions::Vector{<:AbstractVector}; kwargs...)
+    model_DFT(lattice, atoms, positions, [:hyb_gga_xc_pbeh]; kwargs...)
+end
+
+
+"""
+Build an Hartree-Fock model from the specified atoms.
+"""
+function model_HF(lattice::AbstractMatrix, atoms::Vector{<:Element},
+                  positions::Vector{<:AbstractVector}; extra_terms=[], kwargs...)
+    @warn "Exact exchange in DFTK is hardly optimised and not yet production-ready."
+    model_atomic(lattice, atoms, positions;
+                 extra_terms=[Hartree(), ExactExchange(), extra_terms...],
+                 model_name="HF", kwargs...)
+end
+
+
 # Generate equivalent functions for AtomsBase
-for fun in (:model_atomic, :model_DFT, :model_LDA, :model_PBE, :model_SCAN)
+for fun in (:model_atomic, :model_DFT, :model_LDA, :model_PBE, :model_SCAN,
+            :model_PBE0, :model_HF)
     @eval function $fun(system::AbstractSystem, args...; kwargs...)
         @assert !(:magnetic_moments in keys(kwargs))
         parsed = parse_system(system)

src/terms/Hamiltonian.jl

@@ -67,8 +67,8 @@ function Base.size(block::HamiltonianBlock)
 end
 
 import Base: Matrix, Array
-Array(block::HamiltonianBlock)  = Matrix(block)
-Matrix(block::HamiltonianBlock) = sum(Matrix, block.operators)
+Array(block::HamiltonianBlock) = Matrix(block)
+Matrix(block::HamiltonianBlock)        = sum(Matrix, block.operators)
 Matrix(block::GenericHamiltonianBlock) = sum(Matrix, block.optimized_operators)
 
 struct Hamiltonian

src/terms/exact_exchange.jl

@@ -0,0 +1,83 @@
+"""
+Exact exchange term: the Hartree-Exact exchange energy of the orbitals
+
+-1/2 ∑ ∫∫ ϕ_i^*(r)ϕ_j^*(r')ϕ_i(r')ϕ_j(r) / |r - r'| dr dr'
+
+"""
+struct ExactExchange
+    scaling_factor::Real  # to scale the term (e.g. for hybrid models)
+end
+ExactExchange(; scaling_factor=1) = ExactExchange(scaling_factor)
+(exchange::ExactExchange)(basis) = TermExactExchange(basis, exchange.scaling_factor)
+function Base.show(io::IO, exchange::ExactExchange)
+    fac = isone(exchange.scaling_factor) ? "" : "scaling_factor=$(exchange.scaling_factor)"
+    print(io, "ExactExchange($fac)")
+end
+struct TermExactExchange <: Term
+    scaling_factor::Real  # scaling factor, absorbed into poisson_green_coeffs
+    poisson_green_coeffs::AbstractArray
+end
+function TermExactExchange(basis::PlaneWaveBasis{T}, scaling_factor) where T
+    scale = T(scaling_factor)
+
+    # q = T[0.5, 0.5, 0.5]
+    # Gqs = map(G_vectors(basis)) do G
+    #     recip_vector_red_to_cart(basis.model, G + q)
+    # end
+    # poisson_green_coeffs     = 4T(π) * scale ./ norm2.(Gqs)
+
+    poisson_green_coeffs    = 4T(π) * scale ./ norm2.(G_vectors_cart(basis))
+    poisson_green_coeffs[1] = 0  # Compensating charge background => Zero DC
+
+    @assert iszero(G_vectors(basis, basis.kpoints[1])[1])
+
+    TermExactExchange(scale, poisson_green_coeffs)
+end
+
+# Note: Implementing exact exchange in a scalable and numerically stable way, such that it
+# rapidly converges with k-points is tricky. This implementation here is far too simple and
+# slow to be useful.
+#
+# For further information (in particular on regularising the Coulomb), consider the following
+#      https://www.vasp.at/wiki/index.php/Coulomb_singularity
+#      https://journals.aps.org/prb/pdf/10.1103/PhysRevB.34.4405   (QE default)
+#      https://journals.aps.org/prb/pdf/10.1103/PhysRevB.73.205119
+#      https://docs.abinit.org/topics/documents/hybrids-2017.pdf (Abinit apparently
+#           uses a short-ranged Coulomb)
+
+@timing "ene_ops: ExactExchange" function ene_ops(term::TermExactExchange,
+                                                  basis::PlaneWaveBasis{T}, ψ, occupation;
+                                                  kwargs...) where {T}
+    if isnothing(ψ) || isnothing(occupation)
+        return (; E=T(0), ops=NoopOperator.(basis, basis.kpoints))
+    end
+
+    @assert iszero(basis.model.temperature)  # ground state
+    ψ, occupation = select_occupied_orbitals(basis, ψ, occupation; threshold=0.1)
+    E = zero(T)
+
+    @assert length(ψ) == 1  # TODO: make it work for more kpoints
+    ik   = 1
+    kpt  = basis.kpoints[ik]
+    occk = occupation[ik]
+    ψk   = ψ[ik]
+
+    for (n, ψn) in enumerate(eachcol(ψk))
+        for (m, ψm) in enumerate(eachcol(ψk))
+            ψn_real = ifft(basis, kpt, ψn)
+            ψm_real = ifft(basis, kpt, ψm)
+
+            ρnm_real = conj(ψn_real) .* ψm_real
+            ρmn_real = conj(ψm_real) .* ψn_real
+
+            Vx_nm_four = fft(basis, ρnm_real) .* term.poisson_green_coeffs
+            Vx_nm_real = ifft(basis, Vx_nm_four)
+
+            occ_mn = occk[n] * occk[m]
+            E -= real(dot(Vx_nm_real, ρmn_real)) * basis.dvol * occ_mn / 2
+        end
+    end
+
+    ops = [ExchangeOperator(basis, kpt, term.poisson_green_coeffs, occk, ψk)]
+    (E=E, ops=ops)
+end

src/terms/hartree.jl

@@ -17,7 +17,7 @@ end
 Hartree(; scaling_factor=1) = Hartree(scaling_factor)
 (hartree::Hartree)(basis)   = TermHartree(basis, hartree.scaling_factor)
 function Base.show(io::IO, hartree::Hartree)
-    fac = isone(hartree.scaling_factor) ? "" : ", scaling_factor=$scaling_factor"
+    fac = isone(hartree.scaling_factor) ? "" : ", scaling_factor=$(hartree.scaling_factor)"
     print(io, "Hartree($fac)")
 end
 

src/terms/operators.jl

@@ -9,6 +9,22 @@ They also implement mul! and Matrix(op) for exploratory use.
 """
 abstract type RealFourierOperator end
 # RealFourierOperator contain fields `basis` and `kpoint`
+#
+Array(op::RealFourierOperator) = Matrix(op)
+
+# Slow fallback for getting operator as a dense matrix
+function Matrix(op::RealFourierOperator)
+    T = complex(eltype(op.basis))
+    n_G = length(G_vectors(op.basis, op.kpoint))
+    H = zeros(T, n_G, n_G)
+    ψ = zeros(T, n_G)
+    @views for i in 1:n_G
+        ψ[i] = one(T)
+        mul!(H[:, i], op, ψ)
+        ψ[i] = zero(T)
+    end
+    H
+end
 
 # Unoptimized fallback, intended for exploratory use only.
 # For performance, call through Hamiltonian which saves FFTs.
@@ -162,3 +178,27 @@ function optimize_operators_(ops)
                                                sum([op.potential for op in RSmults]))
     [nonRSmults..., combined_RSmults]
 end
+
+struct ExchangeOperator{T <: Real} <: RealFourierOperator
+    basis::PlaneWaveBasis{T}
+    kpoint::Kpoint{T}
+    poisson_green_coeffs  # TODO This needs to be typed
+    occk  # TODO This needs to be typed
+    ψk    # TODO This needs to be typed
+end
+function apply!(Hψ, op::ExchangeOperator, ψ)
+    # Hψ = - ∑_n f_n ψ_n(r) ∫ (ψ_n)†(r') * ψ(r') / |r-r'|dr'
+    for (n, ψnk) in enumerate(eachcol(op.ψk))
+        ψnk_real = ifft(op.basis, op.kpoint, ψnk)
+        x_real   = conj(ψnk_real) .* ψ.real
+
+        # Compute integral by Poisson solve
+        x_four   = fft(op.basis, x_real)
+        Vx_four  = x_four .* op.poisson_green_coeffs
+        Vx_real  = ifft(op.basis, Vx_four)
+
+        # Real-space multiply and accumulate
+        Hψ.real .-= op.occk[n] .* ψnk_real .* Vx_real
+    end
+end
+# TODO Implement  Matrix(op::ExchangeOperator)

src/terms/terms.jl

@@ -51,6 +51,7 @@ include("ewald.jl")
 include("psp_correction.jl")
 include("entropy.jl")
 include("pairwise.jl")
+include("exact_exchange.jl")
 
 include("magnetic.jl")
 breaks_symmetries(::Magnetic) = true