diff --git a/Project.toml b/Project.toml index 20cb3ff..3cd5a33 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "DftFunctionals" uuid = "6bd331d2-b28d-4fd3-880e-1a1c7f37947f" authors = ["Michael F. Herbst "] -version = "0.3.1" +version = "0.3.2" [deps] ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" @@ -11,7 +11,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" [compat] ComponentArrays = "0.15" DiffResults = "1" -ForwardDiff = "1" +ForwardDiff = "0.10, 1" Libxc_jll = "5.1.0" julia = "1.10" diff --git a/src/DftFunctionals.jl b/src/DftFunctionals.jl index 5f0f81c..cbf122e 100644 --- a/src/DftFunctionals.jl +++ b/src/DftFunctionals.jl @@ -16,7 +16,8 @@ include("DftFunctional.jl") include("functionals/lda.jl") include("functionals/gga_x_pbe.jl") include("functionals/gga_c_pbe.jl") +include("functionals/gga_c_lyp.jl") export DftFunctional -export LdaExchange, LdaCorrelationVwn, LdaCorrelationPw, PbeExchange, PbeCorrelation +export LdaExchange, LdaCorrelationVwn, LdaCorrelationPw, PbeExchange, PbeCorrelation, LypCorrelation end diff --git a/src/functionals/gga_c_lyp.jl b/src/functionals/gga_c_lyp.jl new file mode 100644 index 0000000..3b682d9 --- /dev/null +++ b/src/functionals/gga_c_lyp.jl @@ -0,0 +1,82 @@ +struct LypCorrelation{CA} <: Functional{:gga,:c} where {CA<:ComponentArray{<:Number}} + parameters::CA + identifier::Symbol +end +function LypCorrelation(parameters::ComponentArray) + LypCorrelation(parameters, :gga_c_lyp_custom) +end + +identifier(pbe::LypCorrelation) = pbe.identifier +parameters(pbe::LypCorrelation) = pbe.parameters +function change_parameters(pbe::LypCorrelation, parameters::ComponentArray; + keep_identifier=false) + if keep_identifier + PbeCorrelation(parameters, pbe.identifier) + else + PbeCorrelation(parameters) + end +end + +function energy(lyp::LypCorrelation, ρ::T, σ::U) where {T<:Number,U<:Number} + TT = arithmetic_type(lyp, 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 + # + # TODO Comparing to libxc, this implementation has quite some numerical stability issues + # in the kernel + # + a = lyp.parameters.a + b = lyp.parameters.b + c = lyp.parameters.c + d = lyp.parameters.d + + rr = cbrt(1 / ρ) # (Wigner radius rₛ without prefactors) + + # Follows the partial integration trick in DOI 10.1016/0009-2614(89)87234-3, equation (2) + δ(rr) = c * rr + d * rr / (1 + d * rr) + ω(rr) = exp(-c * rr) / (1 + d * rr) # Factor ρ^(-11/3) absorbed into other terms + # for stability reasons + CF = 3/10 * (3*π^2)^(2/3) + + # ζ spin polarization (== 0 for non-spin-polarised) + ζ = false + + # Again assuming no spin polarisation |∇ρ_{α}|² = |½ ∇ρ|² = ¼ |∇ρ|² = σ/4 + ∇ρα² = σ/4 + ∇ρβ² = σ/4 + ρα_ρ = (1 + ζ)/2 # = ρ_α / ρ + ρβ_ρ = (1 - ζ)/2 # = ρ_β / ρ + + # Term in the square brackets in (2) + term_square_bracket = ( + 2.0^(11/3) * CF * ((ρα_ρ*ρ)^(8/3) + (ρβ_ρ*ρ)^(8/3)) + + (47/18 - 7/18 * δ(rr)) * σ + - (5/2 - δ(rr)/18) * (∇ρα² + ∇ρβ²) + - (δ(rr) - 11)/9 * ((1+ζ)/2 * ∇ρβ² + (1-ζ)/2 * ∇ρα²) + # ρ_α/ρ ρ_β/ρ + ) + + # Note, that (1 - ζ^2) = 4ρ_α ρ_β / ρ^2 + res = ( + -a * (1-ζ^2)*ρ / (1 + d*rr) + -a*b * ω(rr) * rr^5 * ( term_square_bracket * (1 - ζ^2)/4 + - (2/3) * σ + + ((2/3) - ρα_ρ^2) * ∇ρβ² + + ((2/3) - ρβ_ρ^2) * ∇ρα² ) + ) + + TT(res) +end + +# +# Concrete functionals +# + +""" +Standard LYP correlation. +Lee, Yang Parr 1988 (DOI: 10.1103/PhysRevB.37.785) +""" +function DftFunctional(::Val{:gga_c_lyp}) + LypCorrelation(ComponentArray(; a=0.04918, b=0.132, c=0.2533, d=0.349), :gga_c_lyp) +end diff --git a/src/functionals/lda.jl b/src/functionals/lda.jl index 04a4632..80944af 100644 --- a/src/functionals/lda.jl +++ b/src/functionals/lda.jl @@ -43,7 +43,7 @@ function energy_per_particle(::LdaCorrelationVwn, ρ::T) where {T<:Number} Xx0 = x0^2 + b * x0 + c Q = sqrt(4c - b^2) A * (log(x^2 / Xx) + 2b / Q * atan(Q / (2x + b)) - - b * x0 / Xx0 * (log((x - x0)^2 / Xx) + 2 * (b + 2x0) / Q * atan(Q / (2x + b)))) + b * x0 / Xx0 * (log((x - x0)^2 / Xx) + 2 * (b + 2x0) / Q * atan(Q / (2x + b)))) end # diff --git a/test/runtests.jl b/test/runtests.jl index ab5c4e3..35c6e7b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,8 +11,15 @@ include("libxc.jl") :gga_x_xpbe, :gga_c_xpbe, :gga_x_pbe_mol, :gga_c_pbe_mol, :gga_x_apbe, :gga_c_apbe, + + :gga_c_lyp, + ) + + kernel_atol = Dict( + :gga_c_lyp => 1e-10, # Note the numerical stability issues here ) + @testset "LDA Functional construction" begin let f = DftFunctional(:lda_x) @test kind(f) == :x @@ -131,7 +138,6 @@ end Vρσref = similar(ρ) Vσσref = similar(ρ) - ptr = xc_functional_alloc(identifier(func)) xc_gga(ptr, n_p, ρ, σ, εref, Vρref, Vσref, Vρρref, Vρσref, Vσσref, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL) @@ -147,9 +153,9 @@ end @test result.e ≈ eref atol=5e-13 @test result.Vρ ≈ Vρref atol=5e-13 @test result.Vσ ≈ Vσref atol=5e-13 - @test result.Vρρ ≈ Vρρref atol=5e-13 - @test result.Vρσ ≈ Vρσref atol=5e-13 - @test result.Vσσ ≈ Vσσref atol=5e-13 + @test result.Vρρ ≈ Vρρref atol=get(kernel_atol, func_name, 5e-13) + @test result.Vρσ ≈ Vρσref atol=get(kernel_atol, func_name, 5e-13) + @test result.Vσσ ≈ Vσσref atol=get(kernel_atol, func_name, 5e-13) # Ensure Float32 evaluation works e = DftFunctionals.energy(func, rand(Float32), rand(Float32))