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 family, kind
 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).
@@ -61,7 +51,7 @@ 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)
+    PbeCorrelation(; β, γ)
 end
 
 """
@@ -72,7 +62,7 @@ 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)
+    PbeCorrelation(; β, γ)
 end
 
 """
@@ -82,7 +72,7 @@ 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)
+    PbeCorrelation(; β, γ)
 end
 
 """
@@ -93,7 +83,7 @@ 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)
+    PbeCorrelation(; β, γ)
 end
 
 """
@@ -104,7 +94,7 @@ function DftFunctional(::Val{:gga_c_pbe_mol})
     # β 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)
+    PbeCorrelation(; β, γ)
 end
 
 """
@@ -114,5 +104,5 @@ 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)
+    PbeCorrelation(; β, γ)
 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
@@ -48,23 +40,23 @@ 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)
+    PbeExchange(κ=0.8040, μ=pbe_μ_from_β(0.06672455060314922))
 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)
+    PbeExchange(κ=1.245, μ=pbe_μ_from_β(0.06672455060314922))
 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
+    PbeExchange(κ=0.91954, μ=0.23214)  # Table 1
 end
 
 """
@@ -73,7 +65,7 @@ Perdew, Ruzsinszky, Csonka and others 2008 (DOI 10.1103/physrevlett.100.136406)
 """
 function DftFunctional(::Val{:gga_x_pbe_sol})
     # μ given below equation (2)
-    PbeExchange(ComponentArray(κ=0.8040, μ=10 / 81), :gga_x_pbe_sol)
+    PbeExchange(κ=0.8040, μ=10 / 81)
 end
 
 """
@@ -82,7 +74,7 @@ Constantin, Fabiano, Laricchia 2011 (DOI 10.1103/physrevlett.106.186406)
 """
 function DftFunctional(::Val{:gga_x_apbe})
     # p. 1, right column, bottom
-    PbeExchange(ComponentArray(κ=0.8040, μ=0.260), :gga_x_apbe)
+    PbeExchange(κ=0.8040, μ=0.260)
 end
 
 """
@@ -91,13 +83,13 @@ del Campo, Gazqez, Trickey and others 2012 (DOI 10.1063/1.3691197)
 """
 function DftFunctional(::Val{:gga_x_pbe_mol})
     # p. 4, left column, bottom
-    PbeExchange(ComponentArray(κ=0.8040, μ=0.27583), :gga_x_pbe_mol)
+    PbeExchange(κ=0.8040, μ=0.27583)
 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
+    PbeExchange(κ=0.437, μ=0.346)  # Table 1
 end

src/functionals/lda.jl

@@ -1,6 +1,5 @@
 struct LdaExchange <: Functional{:lda,:x}
 end
-identifier(::LdaExchange) = :lda_x
 
 """
 LDA Slater exchange (DOI: 10.1017/S0305004100016108 and 10.1007/BF01340281)
@@ -24,7 +23,6 @@ end
 #
 struct LdaCorrelationVwn <: Functional{:lda,:c}
 end
-identifier(::LdaCorrelationVwn) = :lda_c_vwn
 
 """
 VWN5 LDA correlation according to Vosko, Wilk, and Nusair, (DOI 10.1139/p80-159).
@@ -51,7 +49,6 @@ end
 #
 struct LdaCorrelationPw{Improved} <: Functional{:lda,:c}
 end
-identifier(::LdaCorrelationPw) = :lda_c_pw
 
 """
 Perdew, Wang correlation from 1992 (10.1103/PhysRevB.45.13244)

src/interface.jl

@@ -15,10 +15,6 @@ 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"""
-function identifier end
-Base.show(io::IO, fun::Functional) = print(io, identifier(fun))
-
 @doc raw"""
 True if the functional needs ``σ = 𝛁ρ ⋅ 𝛁ρ``.
 """
@@ -42,21 +38,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...)