Make functionals isbits types, in preparation for GPU usage

Project.toml

@@ -4,7 +4,6 @@ authors = ["Michael F. Herbst <info@michael-herbst.com>"]
 version = "0.3.1"
 
 [deps]
-ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 
@@ -16,8 +15,9 @@ Libxc_jll = "5.1.0"
 julia = "1.10"
 
 [extras]
+ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 Libxc_jll = "a56a6d9d-ad03-58af-ab61-878bf78270d6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Libxc_jll"]
+test = ["Test", "ComponentArrays", "Libxc_jll"]

src/DftFunctionals.jl

@@ -1,12 +1,10 @@
 module DftFunctionals
 using ForwardDiff
-using ComponentArrays
 
 include("interface.jl")
 include("util.jl")
 export Functional
 export family, kind, identifier
-export parameters, change_parameters
 export needs_σ, needs_τ, needs_Δρ, has_energy
 export potential_terms, kernel_terms
 export spinindex_σ

src/functionals/gga_c_pbe.jl

@@ -1,34 +1,24 @@
-struct PbeCorrelation{Tlda,CA} <:
-       Functional{:gga,:c} where {Tlda,CA<:ComponentArray{<:Number}}
-    parameters::CA
+struct PbeCorrelation{Tlda,Tβ,Tγ} <:
+       Functional{:gga,:c} where {Tlda,Tβ<:Number,Tγ<:Number}
     lda::Tlda
-    identifier::Symbol
+    β::Tβ
+    γ::Tγ
 end
-function PbeCorrelation(parameters::ComponentArray, lda=DftFunctional(:lda_c_pw))
-    PbeCorrelation(parameters, lda, :gga_c_pbe_custom)
-end
-function PbeCorrelation(parameters::ComponentArray, identifier::Symbol)
-    PbeCorrelation(parameters, DftFunctional(:lda_c_pw), identifier)
+function PbeCorrelation(; lda=DftFunctional(:lda_c_pw), β, γ)
+    PbeCorrelation(lda, β, γ)
 end
 
-identifier(pbe::PbeCorrelation) = pbe.identifier
-parameters(pbe::PbeCorrelation) = pbe.parameters
-function change_parameters(pbe::PbeCorrelation, parameters::ComponentArray;
-                           keep_identifier=false)
-    if keep_identifier
-        PbeCorrelation(parameters, pbe.lda, pbe.identifier)
-    else
-        PbeCorrelation(parameters, pbe.lda)
-    end
+function parameters_type(pbe::PbeCorrelation)
+    promote_type(parameters_type(pbe.lda), typeof(pbe.β), typeof(pbe.γ))
 end
 
 function energy(pbe::PbeCorrelation, ρ::T, σ::U) where {T<:Number,U<:Number}
     TT = arithmetic_type(pbe, T, U)
 
     # TODO This function is quite sensitive to the floating-point type ...
     #      so for now we don't bother doing this in TT, but rather convert before return
-    β = pbe.parameters.β
-    γ = pbe.parameters.γ
+    β = pbe.β
+    γ = pbe.γ
 
     # Spin-scaling factor with ζ spin polarization.
     # Yue Wang and John P. Perdew. Phys. Rev. B 43, 8911 (1991).
@@ -54,65 +44,52 @@ end
 # Concrete functionals
 #
 
-"""
-Standard PBE correlation.
-Perdew, Burke, Ernzerhof 1996 (DOI: 10.1103/PhysRevLett.77.3865)
-"""
-function DftFunctional(::Val{:gga_c_pbe})
-    β = 0.06672455060314922
-    γ = (1 - log(2)) / π^2
-    PbeCorrelation(ComponentArray(; β, γ), :gga_c_pbe)
-end
-
-"""
-XPBE correlation.
-Xu, Goddard 2004 (DOI 10.1063/1.1771632)
-"""
-function DftFunctional(::Val{:gga_c_xpbe})
-    β = 0.089809  # Fitted constants, Table I
-    α = 0.197363  # Fitted constants, Table I
-    γ = β^2 / 2α
-    PbeCorrelation(ComponentArray(; β, γ), :gga_c_xpbe)
-end
-
-"""
-PBESol correlation.
-Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
-"""
-function DftFunctional(::Val{:gga_c_pbe_sol})
-    β = 0.046  # Page 3, left column below figure 1
-    γ = (1 - log(2)) / π^2
-    PbeCorrelation(ComponentArray(; β, γ), :gga_c_pbe_sol)
-end
-
-"""
-APBE correlation.
-Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
-"""
-function DftFunctional(::Val{:gga_c_apbe})
-    μ = 0.260   # p. 1, right column, bottom
-    β = 3μ / π^2
-    γ = (1 - log(2)) / π^2  # like in PBE
-    PbeCorrelation(ComponentArray(; β, γ), :gga_c_apbe)
-end
-
-"""
-PBEmol correlation.
-del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
-"""
-function DftFunctional(::Val{:gga_c_pbe_mol})
+const KNOWN_C_PBE = [
+    # Standard PBE correlation.
+    # Perdew, Burke, Ernzerhof 1996 (DOI: 10.1103/PhysRevLett.77.3865)
+    :gga_c_pbe => (; β=0.06672455060314922, γ=(1 - log(2)) / π^2),
+    # XPBE correlation.
+    # Xu, Goddard 2004 (DOI 10.1063/1.1771632)
+    :gga_c_xpbe => let
+        β = 0.089809  # Fitted constants, Table I
+        α = 0.197363  # Fitted constants, Table I
+        γ = β^2 / 2α
+        (; β, γ)
+    end,
+    # PBESol correlation.
+    # Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
+    # Page 3, left column below figure 1
+    :gga_c_pbe_sol => (; β=0.046, γ=(1 - log(2)) / π^2),
+    # APBE correlation.
+    # Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
+    :gga_c_apbe => let
+        μ = 0.260   # p. 1, right column, bottom
+        β = 3μ / π^2
+        γ = (1 - log(2)) / π^2  # like in PBE
+        (; β, γ)
+    end,
+    # PBEmol correlation.
+    # del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
     # β made to cancel self-interaction error in hydrogen
-    β = 0.08384             # p. 4, right column, first paragraph
-    γ = (1 - log(2)) / π^2  # like in PBE
-    PbeCorrelation(ComponentArray(; β, γ), :gga_c_pbe_mol)
+    # p. 4, right column, first paragraph
+    :gga_c_pbe_mol => (; β=0.08384, γ=(1 - log(2)) / π^2),
+    # PBEfe correlation.
+    # Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)
+    # Fitted constants, Table I
+    :gga_c_pbefe => (; β=0.043, γ=0.031090690869654895034),
+]
+
+for (id, param) in KNOWN_C_PBE
+    @eval function DftFunctional(::Val{$(QuoteNode(id))})
+        PbeCorrelation(β=$(param.β), γ=$(param.γ))
+    end
 end
 
-"""
-PBEfe correlation.
-Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)
-"""
-function DftFunctional(::Val{:gga_c_pbefe})
-    β = 0.043                    # Fitted constants, Table I
-    γ = 0.031090690869654895034  # Fitted constants, Table I
-    PbeCorrelation(ComponentArray(; β, γ), :gga_c_pbefe)
+function identifier(pbe::PbeCorrelation)
+    for (id, param) in KNOWN_C_PBE
+        if pbe.β ≈ param.β && pbe.γ ≈ param.γ
+            return id
+        end
+    end
+    nothing
 end

src/functionals/gga_x_pbe.jl

@@ -1,29 +1,21 @@
-struct PbeExchange{CA} <: Functional{:gga,:x} where {CA<:ComponentArray{<:Number}}
-    parameters::CA
-    identifier::Symbol
+struct PbeExchange{Tκ, Tμ} <: Functional{:gga,:x} where {Tκ<:Number,Tμ<:Number}
+    κ::Tκ
+    μ::Tμ
 end
-function PbeExchange(parameters::ComponentArray)
-    PbeExchange(parameters, :gga_x_pbe_custom)
+function PbeExchange(; κ, μ)
+    PbeExchange(κ, μ)
 end
-
-identifier(pbe::PbeExchange) = pbe.identifier
-parameters(pbe::PbeExchange) = pbe.parameters
-function change_parameters(pbe::PbeExchange, parameters::ComponentArray;
-                           keep_identifier=false)
-    if keep_identifier
-        PbeExchange(parameters, pbe.identifier)
-    else
-        PbeExchange(parameters)
-    end
+function parameters_type(pbe::PbeExchange)
+    promote_type(typeof(pbe.κ), typeof(pbe.μ))
 end
 
 function energy(pbe::PbeExchange, ρ::T, σ::U) where {T<:Number,U<:Number}
     TT = arithmetic_type(pbe, T, U)
 
     # TODO This function is quite sensitive to the floating-point type ...
     #      so for now we don't bother doing this in TT, but rather convert before return
-    κ = pbe.parameters.κ
-    μ = pbe.parameters.μ
+    κ = pbe.κ
+    μ = pbe.μ
 
     pbe_x_f(s²) = 1 + κ - κ^2 / (κ + μ * s²)   # (14)
     # rₛ = cbrt(3 / (4π  * ρ))                 # page 2, left column, top
@@ -43,61 +35,44 @@ pbe_β_from_μ(μ) = 3μ / π^2
 # Concrete functionals
 #
 
-"""
-Standard PBE exchange.
-Perdew, Burke, Ernzerhof 1996 (DOI: 10.1103/PhysRevLett.77.3865)
-"""
-function DftFunctional(::Val{:gga_x_pbe})
-    PbeExchange(ComponentArray(κ=0.8040, μ=pbe_μ_from_β(0.06672455060314922)), :gga_x_pbe)
-end
-
-"""
-Revised PBE exchange.
-Zhang, Yang 1998 (DOI 10.1103/physrevlett.80.890)
-"""
-function DftFunctional(::Val{:gga_x_pbe_r})
-    PbeExchange(ComponentArray(κ=1.245, μ=pbe_μ_from_β(0.06672455060314922)), :gga_x_pbe_r)
-end
-
-"""
-XPBE exchange.
-Xu, Goddard 2004 (DOI 10.1063/1.1771632)
-"""
-function DftFunctional(::Val{:gga_x_xpbe})
-    PbeExchange(ComponentArray(κ=0.91954, μ=0.23214), :gga_x_xpbe)  # Table 1
-end
-
-"""
-PBESol exchange.
-Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
-"""
-function DftFunctional(::Val{:gga_x_pbe_sol})
+const KNOWN_X_PBE = [
+    # Standard PBE exchange.
+    # Perdew, Burke, Ernzerhof 1996 (DOI: 10.1103/PhysRevLett.77.3865)
+    :gga_x_pbe     => (; κ=0.8040, μ=pbe_μ_from_β(0.06672455060314922)),
+    # Revised PBE exchange.
+    # Zhang, Yang 1998 (DOI 10.1103/physrevlett.80.890)
+    :gga_x_pbe_r   => (; κ=1.245, μ=pbe_μ_from_β(0.06672455060314922)),
+    # XPBE exchange.
+    # Xu, Goddard 2004 (DOI 10.1063/1.1771632)
+    :gga_x_xpbe    => (; κ=0.91954, μ=0.23214), # Table 1
+    # PBESol exchange.
+    # Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
     # μ given below equation (2)
-    PbeExchange(ComponentArray(κ=0.8040, μ=10 / 81), :gga_x_pbe_sol)
-end
-
-"""
-APBE exchange.
-Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
-"""
-function DftFunctional(::Val{:gga_x_apbe})
+    :gga_x_pbe_sol => (; κ=0.8040, μ=10 / 81),
+    # APBE exchange.
+    # Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
     # p. 1, right column, bottom
-    PbeExchange(ComponentArray(κ=0.8040, μ=0.260), :gga_x_apbe)
-end
-
-"""
-PBEmol exchange.
-del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
-"""
-function DftFunctional(::Val{:gga_x_pbe_mol})
+    :gga_x_apbe    => (; κ=0.8040, μ=0.260),
+    # PBEmol exchange.
+    # del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
     # p. 4, left column, bottom
-    PbeExchange(ComponentArray(κ=0.8040, μ=0.27583), :gga_x_pbe_mol)
+    :gga_x_pbe_mol => (; κ=0.8040, μ=0.27583),
+    # PBEfe exchange.
+    # Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)
+    :gga_x_pbefe   => (; κ=0.437, μ=0.346), # Table 1
+]
+
+for (id, param) in KNOWN_X_PBE
+    @eval function DftFunctional(::Val{$(QuoteNode(id))})
+        PbeExchange(κ=$(param.κ), μ=$(param.μ))
+    end
 end
 
-"""
-PBEfe exchange.
-Sarmiento-Perez, Silvana, Marques 2015 (DOI 10.1021/acs.jctc.5b00529)
-"""
-function DftFunctional(::Val{:gga_x_pbefe})
-    PbeExchange(ComponentArray(κ=0.437, μ=0.346), :gga_x_pbefe)  # Table 1
+function identifier(pbe::PbeExchange)
+    for (id, param) in KNOWN_X_PBE
+        if pbe.κ ≈ param.κ && pbe.μ ≈ param.μ
+            return id
+        end
+    end
+    nothing
 end

src/interface.jl

@@ -15,9 +15,19 @@ Return the functional kind: `:x` (exchange), `:c` (correlation), `:k` (kinetic)
 """
 kind(::Functional{F,K}) where {F,K} = K
 
-"""Return the identifier corresponding to a functional"""
+"""
+Return the identifier corresponding to a functional, if available,
+and `nothing` otherwise.
+"""
 function identifier end
-Base.show(io::IO, fun::Functional) = print(io, identifier(fun))
+function Base.show(io::IO, fun::Functional)
+    id = identifier(fun)
+    if isnothing(id)
+        Base.show_default(io, fun)
+    else
+        print(io, id)
+    end
+end
 
 @doc raw"""
 True if the functional needs ``σ = 𝛁ρ ⋅ 𝛁ρ``.
@@ -42,21 +52,13 @@ i.e. `e` will be `false` (a strong zero).
 has_energy(::Functional) = true
 
 """
-Return adjustable parameters of the functional and their values.
-"""
-parameters(::Functional) = ComponentArray{Bool}()
+The type of the parameters of the functional. If the functional has multiple parameters,
+the result of `promote_type(...)` of the parameter types should be returned.
 
+This allows the functional to influence the computed type, for example if it contains
+Dual numbers.
 """
-Return a new version of the passed functional with its parameters adjusted.
-This may not be a copy in case no changes are done to its internal parameters.
-Generally the identifier of the functional will be changed to reflect the
-change in parameter values unless `keep_identifier` is true.
-To get the tuple of adjustable parameters and their current values check out
-[`parameters`](@ref). It is not checked that the correct parameters are passed.
-
-`change_parameters(f::Functional, params_new; keep_identifier=false)::Functional`
-"""
-function change_parameters end
+parameters_type(::Functional) = Bool # default: Bool, which promotes to any number type
 
 # TODO These values are read-only for now and their defaults hard-coded for Float64
 """

src/util.jl

@@ -27,7 +27,7 @@ number types always win (to ensure that `Float32` density data causes a `Float32
 functional evaluation even if the functional parameters are stored in `Float64`.
 """
 function arithmetic_type(func::Functional, T, S...)
-    arithmetic_type_(eltype(parameters(func)), T, S...)
+    arithmetic_type_(parameters_type(func), T, S...)
 end
 arithmetic_type_(::Type{<:AbstractFloat}, T::Type, S...) = promote_type(T, S...)
 arithmetic_type_(PT::Type, T, S...) = promote_type(PT, T, S...)