Basic LYP correlation functional

Project.toml

@@ -1,7 +1,7 @@
 name = "DftFunctionals"
 uuid = "6bd331d2-b28d-4fd3-880e-1a1c7f37947f"
 authors = ["Michael F. Herbst <info@michael-herbst.com>"]
-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"
 

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

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

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
 
 #

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))