From e0a6b55b1cec6bf88f3bcd88cd968dbd372a0ef1 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Mon, 12 May 2025 17:47:49 +0200 Subject: [PATCH 1/4] Basic LYP correlation functional --- src/DftFunctionals.jl | 3 +- src/functionals/gga_c_lyp.jl | 78 ++++++++++++++++++++++++++++++++++++ src/functionals/lda.jl | 2 +- test/runtests.jl | 13 ++++-- 4 files changed, 91 insertions(+), 5 deletions(-) create mode 100644 src/functionals/gga_c_lyp.jl 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..ab8581a --- /dev/null +++ b/src/functionals/gga_c_lyp.jl @@ -0,0 +1,78 @@ +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 + 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) * ρ^(-(9+2)/3.0) + CF = 3/10 * (3*π^2)^(2/3) + + # ζ spin polarization (== 0 for non-spin-polarised) + ζ = false + ραρβ = (1 - ζ^2) * ρ*ρ / 4 # Note that (1 - ζ^2) * ρ = 4ρ_α ρ_β / ρ + # such that ρ_α * ρ_β = (1 - ζ^2) * ρ^2 / 4 + + # 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 * ∇ρα²) + # ρ_α/ρ ρ_β/ρ + ) + + res = ( + -a * (1-ζ^2)*ρ / (1 + d*rr) # = 4ρα ρβ / ρ * (1+d * rr) + -a*b * ω(rr) * ( ραρβ * term_square_bracket + - (2/3) * ρ^2 * σ + + ((2/3) * ρ^2 - ρα^2) * ∇ρβ² + + ((2/3) * ρ^2 - ρβ^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..4bd9ef7 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-11, # Note the numerical stability issues here ) + @testset "LDA Functional construction" begin let f = DftFunctional(:lda_x) @test kind(f) == :x @@ -147,9 +154,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)) From 0292b1abb28c17d18207089ecabc9561d2ab793d Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Mon, 12 May 2025 17:58:18 +0200 Subject: [PATCH 2/4] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 20cb3ff..77bd6d4 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" From 7ae2c0fc9e2921305e4344899fecf4a71379c669 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Mon, 12 May 2025 17:59:40 +0200 Subject: [PATCH 3/4] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 77bd6d4..3cd5a33 100644 --- a/Project.toml +++ b/Project.toml @@ -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" From bc934bcea75012333f633d1ae1e1359c375d51c0 Mon Sep 17 00:00:00 2001 From: "Michael F. Herbst" Date: Mon, 12 May 2025 18:43:02 +0200 Subject: [PATCH 4/4] up --- src/functionals/gga_c_lyp.jl | 26 +++++++++++++++----------- test/runtests.jl | 5 ++--- 2 files changed, 17 insertions(+), 14 deletions(-) diff --git a/src/functionals/gga_c_lyp.jl b/src/functionals/gga_c_lyp.jl index ab8581a..3b682d9 100644 --- a/src/functionals/gga_c_lyp.jl +++ b/src/functionals/gga_c_lyp.jl @@ -22,6 +22,10 @@ function energy(lyp::LypCorrelation, ρ::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 + # + # 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 @@ -31,35 +35,35 @@ function energy(lyp::LypCorrelation, ρ::T, σ::U) where {T<:Number,U<:Number} # 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) * ρ^(-(9+2)/3.0) + ω(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 - ραρβ = (1 - ζ^2) * ρ*ρ / 4 # Note that (1 - ζ^2) * ρ = 4ρ_α ρ_β / ρ - # such that ρ_α * ρ_β = (1 - ζ^2) * ρ^2 / 4 # Again assuming no spin polarisation |∇ρ_{α}|² = |½ ∇ρ|² = ¼ |∇ρ|² = σ/4 ∇ρα² = σ/4 ∇ρβ² = σ/4 - ρα = (1 + ζ)/2 * ρ - ρβ = (1 - ζ)/2 * ρ + ρα_ρ = (1 + ζ)/2 # = ρ_α / ρ + ρβ_ρ = (1 - ζ)/2 # = ρ_β / ρ # Term in the square brackets in (2) term_square_bracket = ( - 2.0^(11/3) * CF * (ρα^(8/3) + ρβ^(8/3)) + 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) # = 4ρα ρβ / ρ * (1+d * rr) - -a*b * ω(rr) * ( ραρβ * term_square_bracket - - (2/3) * ρ^2 * σ - + ((2/3) * ρ^2 - ρα^2) * ∇ρβ² - + ((2/3) * ρ^2 - ρβ^2) * ∇ρα² ) + -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) diff --git a/test/runtests.jl b/test/runtests.jl index 4bd9ef7..35c6e7b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -12,11 +12,11 @@ include("libxc.jl") :gga_x_pbe_mol, :gga_c_pbe_mol, :gga_x_apbe, :gga_c_apbe, - :gga_c_lyp + :gga_c_lyp, ) kernel_atol = Dict( - :gga_c_lyp => 1e-11, # Note the numerical stability issues here + :gga_c_lyp => 1e-10, # Note the numerical stability issues here ) @@ -138,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)