Make apply_kernel error when not implemented, instead of assuming zero (#1141)

Author: niklasschmitz

src/terms/entropy.jl

@@ -37,3 +37,5 @@ function ene_ops(term::TermEntropy, basis::PlaneWaveBasis{T}, ψ, occupation;
 
     (; E, ops)
 end
+compute_kernel(term::TermEntropy, basis::AbstractBasis; kwargs...) = nothing
+apply_kernel(term::TermEntropy, basis::AbstractBasis, δρ; kwargs...) = nothing

src/terms/ewald.jl

@@ -11,7 +11,7 @@ Base.@kwdef struct Ewald
 end
 (ewald::Ewald)(basis) = TermEwald(basis; η=something(ewald.η, default_η(basis.model.lattice)))
 
-struct TermEwald{T} <: Term
+struct TermEwald{T} <: TermLinear
     energy::T                # precomputed energy
     forces::Vector{Vec3{T}}  # and forces
     η::T                     # Parameter used for the splitting

src/terms/kinetic.jl

@@ -16,7 +16,7 @@ function Base.show(io::IO, kin::Kinetic)
     print(io, "Kinetic($bup$fac)")
 end
 
-struct TermKinetic <: Term
+struct TermKinetic <: TermLinear
     scaling_factor::Real  # scaling factor, absorbed into kinetic_energies
     # kinetic energies 1/2(k+G)^2 *blowup(|k+G|, Ecut) for each k-point.
     kinetic_energies::Vector{<:AbstractVector}

src/terms/local.jl

@@ -4,7 +4,7 @@
 # potential in real space on the grid. If the potential is different in the α and β
 # components then it should be a 4d-array with the last axis running over the
 # two spin components.
-abstract type TermLocalPotential <: Term end
+abstract type TermLocalPotential <: TermLinear end
 
 @timing "ene_ops: local" function ene_ops(term::TermLocalPotential,
                                           basis::PlaneWaveBasis{T}, ψ, occupation;

src/terms/magnetic.jl

@@ -22,7 +22,7 @@ function (M::MagneticFromValues)(Apotential)
     TermMagnetic(basis, Apotential)
 end
 
-struct TermMagnetic <: Term
+struct TermMagnetic <: TermLinear
     # Apotential[α] is an array of size fft_size for α=1:3
     Apotential::Vector{<:AbstractArray}
 end

src/terms/nonlocal.jl

@@ -24,7 +24,7 @@ function (::AtomicNonlocal)(basis::PlaneWaveBasis{T}) where {T}
     TermAtomicNonlocal(ops)
 end
 
-struct TermAtomicNonlocal <: Term
+struct TermAtomicNonlocal <: TermLinear
     ops::Vector{NonlocalOperator}
 end
 

src/terms/pairwise.jl

@@ -25,7 +25,7 @@ end
     TermPairwisePotential(P.V, P.params, T(P.max_radius), energy, forces)
 end
 
-struct TermPairwisePotential{TV, Tparams, T} <:Term
+struct TermPairwisePotential{TV, Tparams, T} <: TermLinear
     V::TV
     params::Tparams
     max_radius::T

src/terms/psp_correction.jl

@@ -4,7 +4,7 @@ Pseudopotential correction energy. TODO discuss the need for this.
 struct PspCorrection end
 (::PspCorrection)(basis) = TermPspCorrection(basis)
 
-struct TermPspCorrection{T} <: Term
+struct TermPspCorrection{T} <: TermLinear
     energy::T  # precomputed energy
 end
 function TermPspCorrection(basis::PlaneWaveBasis)

src/terms/terms.jl

@@ -16,6 +16,11 @@ include("operators.jl")
 # In particular, dE/dψn = 2 fn |Hψn> (plus weighting for k-point sampling)
 abstract type Term end
 
+# Terms that are linear in the density matrix, i.e. have zero second derivative
+abstract type TermLinear <: Term end
+compute_kernel(term::TermLinear, basis::AbstractBasis; kwargs...) = nothing
+apply_kernel(term::TermLinear, basis::AbstractBasis, δρ; kwargs...) = nothing
+
 # Terms that are non-linear in the density (i.e. which give rise to a Hamiltonian
 # contribution that is density-dependent or orbital-dependent as well)
 abstract type TermNonlinear <: Term end
@@ -33,7 +38,7 @@ DftFunctionals.needs_τ(t::Term) = false
 """
 A term with a constant zero energy.
 """
-struct TermNoop <: Term end
+struct TermNoop <: TermLinear end
 function ene_ops(term::TermNoop, basis::PlaneWaveBasis{T}, ψ, occupation; kwargs...) where {T}
     (; E=zero(eltype(T)), ops=[NoopOperator(basis, kpt) for kpt in basis.kpoints])
 end
@@ -101,7 +106,6 @@ In this case the matrix has effectively 4 blocks
     end
     kernel
 end
-compute_kernel(::Term, ::AbstractBasis{T}; kwargs...) where {T} = nothing  # By default no kernel
 
 
 """
@@ -128,4 +132,3 @@ as a 4D `(i,j,k,σ)` array.
     end
     δV
 end
-apply_kernel(::Term, ::AbstractBasis{T}, δρ; kwargs...) where {T} = nothing  # by default, no kernel