diff --git a/Project.toml b/Project.toml index 20cb3ff..a88c886 100644 --- a/Project.toml +++ b/Project.toml @@ -4,7 +4,6 @@ authors = ["Michael F. Herbst "] 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"] diff --git a/src/DftFunctionals.jl b/src/DftFunctionals.jl index 5f0f81c..4ac9c8a 100644 --- a/src/DftFunctionals.jl +++ b/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_σ diff --git a/src/functionals/gga_c_pbe.jl b/src/functionals/gga_c_pbe.jl index 10c40be..d00ada3 100644 --- a/src/functionals/gga_c_pbe.jl +++ b/src/functionals/gga_c_pbe.jl @@ -1,25 +1,15 @@ -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} @@ -27,8 +17,8 @@ function energy(pbe::PbeCorrelation, ρ::T, σ::U) where {T<:Number,U<:Number} # 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 diff --git a/src/functionals/gga_x_pbe.jl b/src/functionals/gga_x_pbe.jl index 91c70df..eab0476 100644 --- a/src/functionals/gga_x_pbe.jl +++ b/src/functionals/gga_x_pbe.jl @@ -1,20 +1,12 @@ -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} @@ -22,8 +14,8 @@ function energy(pbe::PbeExchange, ρ::T, σ::U) where {T<:Number,U<:Number} # 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 diff --git a/src/interface.jl b/src/interface.jl index 684421c..13e8b26 100644 --- a/src/interface.jl +++ b/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 """ diff --git a/src/util.jl b/src/util.jl index 3b65ed3..74bd429 100644 --- a/src/util.jl +++ b/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...) diff --git a/test/runtests.jl b/test/runtests.jl index ab5c4e3..c279d6b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,6 +22,7 @@ include("libxc.jl") @test !needs_τ(f) @test !needs_Δρ(f) @test has_energy(f) + @test isbits(f) end for id in (:lda_c_vwn, :lda_c_pw) @@ -29,6 +30,7 @@ include("libxc.jl") @test kind(f) == :c @test family(f) == :lda @test identifier(f) == id + @test isbits(f) end end @@ -40,16 +42,16 @@ end @test needs_σ(f) @test !needs_τ(f) @test !needs_Δρ(f) + @test isbits(f) end -end -@testset "Parameter interface defaults" begin - struct NewExchange <: Functional{:lda,:x} - end - f = NewExchange() - @test parameters(f) == ComponentArray{Bool}() - x = ComponentArray{Bool}() - @test_throws MethodError change_parameters(f, x) + customgga = PbeExchange(κ=1.0, μ=0.5) + @test family(customgga) == :gga + @test identifier(customgga) === nothing + @test needs_σ(customgga) + @test !needs_τ(customgga) + @test !needs_Δρ(customgga) + @test isbits(customgga) end @testset "LDA potential (without spin)" begin @@ -166,18 +168,6 @@ end end @testset "PBE functionals" begin - pbe = DftFunctional(:gga_x_pbe) - @test :μ in keys(parameters(pbe)) - @test :κ in keys(parameters(pbe)) - - pbemod = change_parameters(pbe, ComponentArray(;μ=12, κ=1.2)) - @test parameters(pbemod).μ == 12 - @test parameters(pbemod).κ == 1.2 - - pbemod = change_parameters(DftFunctional(:gga_c_pbe), ComponentArray(;β=12, γ=1.2)) - @test parameters(pbemod).β == 12 - @test parameters(pbemod).γ == 1.2 - μ = rand() @test μ ≈ DftFunctionals.pbe_μ_from_β(DftFunctionals.pbe_β_from_μ(μ)) end @@ -192,17 +182,17 @@ end ρ = reshape(ρ, 1, :) σ = reshape(σ, 1, :) - θ = ComponentArray(; parameters(pbe)...) + θ = ComponentArray(; pbe.κ, pbe.μ) egrad = ForwardDiff.jacobian(θ) do θ - potential_terms(change_parameters(pbe, θ), ρ, σ).e + potential_terms(PbeExchange(; θ...), ρ, σ).e end egrad_fd = let ε=1e-5 δ = zero(θ) δ[2] = ε - ( potential_terms(change_parameters(pbe, θ + δ), ρ, σ).e - - potential_terms(change_parameters(pbe, θ - δ), ρ, σ).e) / 2ε + ( potential_terms(PbeExchange(; (θ + δ)...), ρ, σ).e + - potential_terms(PbeExchange(; (θ - δ)...), ρ, σ).e) / 2ε end @test maximum(abs, egrad[:, 2] - egrad_fd) < 1e-5