added code from QEDprobing

Project.toml

@@ -1,7 +1,19 @@
 name = "MatterModels"
 uuid = "89b7eb7d-0b8d-5078-86d5-4a0f481f2984"
-authors = ["Uwe Hernandez Acosta <u.hernandez@hzdr.de>"]
 version = "0.1.0"
+authors = ["Uwe Hernandez Acosta <u.hernandez@hzdr.de>"]
+
+[deps]
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
+QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
+LogExpFunctions = "0.3.29"
+QEDcore = "0.4.0"
+QuadGK = "2.11.2"
+SpecialFunctions = "2.6.1"
+Unitful = "1.25.1"
 julia = "1.10"

src/MatterModels.jl

@@ -1,11 +1,77 @@
 module MatterModels
 
-"""
-    hi = hello_world()
-A simple function to return "Hello, World!"
-"""
-function hello_world()
-    return "Hello, World!"
-end
+hello_world() = "Hello, World!"
+
+# temperature
+export FiniteTemperature, ZeroTemperature
+
+# matter model
+export AbstractMatterModel
+
+# Electron system
+export AbstractElectronSystem
+export temperature, electron_density, imag_dynamic_response, real_dynamic_response
+export fermi_wave_vector,
+    fermi_energy, beta, betabar, dynamic_response, dynamic_structure_factor
+
+export AbstractProperElectronSystem
+export AbstractInteractingElectronSystem
+export proper_electron_system, screening
+
+# screening
+export AbstractScreening, Screening, NoScreening
+export dielectric_function,
+    pseudo_potential, local_field_correction, local_effective_potential
+export AbstractPseudoPotential, CoulombPseudoPotential
+export AbstractLocalFieldCorrection, NoLocalFieldCorrection
+
+# concrete electron systems
+export IdealElectronSystem
+export AbstractResponseApproximation, NoApprox, NonDegenerated, Degenerated
+export response_approximation
+export InteractingElectronSystem
+
+# constants
+export HBARC,
+    HBARC_eV_ANG,
+    ELECTRONMASS,
+    ALPHA,
+    ALPHA_SQUARE,
+    ELEMENTARY_CHARGE_SQUARED,
+    ELEMENTARY_CHARGE,
+    HARTREE,
+    BOHR_RADIUS_ANG
+
+using QEDcore
+using QuadGK
+using Unitful
+using LogExpFunctions
+using SpecialFunctions
+
+include("utils.jl")
+include("units.jl")
+include("constants.jl")
+include("temperature.jl")
+include("interface.jl")
+include("generic.jl")
+include("lookup.jl")
+include("data_driven/impl.jl")
+include("electron_system/utils.jl")
+include("electron_system/interface.jl")
+include("electron_system/generic.jl")
+include("electron_system/ideal/approximations/interface.jl")
+include("electron_system/ideal/approximations/no_approx.jl")
+include("electron_system/ideal/approximations/non_degenerated.jl")
+include("electron_system/ideal/approximations/degenerated.jl")
+include("electron_system/ideal/utils.jl")
+include("electron_system/ideal/interface.jl")
+include("electron_system/ideal/generic.jl")
+include("electron_system/ideal/impl.jl")
+include("electron_system/interacting/screening/interface.jl")
+include("electron_system/interacting/screening/generic.jl")
+include("electron_system/interacting/screening/impl.jl")
+include("electron_system/interacting/interface.jl")
+include("electron_system/interacting/generic.jl")
+include("electron_system/interacting/impl.jl")
 
 end

src/constants.jl

@@ -0,0 +1,9 @@
+const HBARC = 197.3269718 # MeV fm
+const HBARC_eV_ANG = HBARC * 1.0e6 * 1.0e-5 # eV Ang
+const ELECTRONMASS = 0.5109989461e6 # eV
+const ALPHA = 1 / (137.035999074)
+const ALPHA_SQUARE = ALPHA^2
+const ELEMENTARY_CHARGE_SQUARED = 4 * pi * ALPHA
+const ELEMENTARY_CHARGE = sqrt(ELEMENTARY_CHARGE_SQUARED)
+const BOHR_RADIUS_ANG = 0.529177210544 # Ang
+const HARTREE = 27.211386245981

src/data_driven/impl.jl

@@ -0,0 +1,28 @@
+# TODO
+# - rename the type (adopt MCSS)
+# - build a better interface of data informed matter models
+# - update lookup function (signature)
+struct DataDrivenMatterModel <: AbstractMatterModel
+    lt::GridLookupTable
+    method::AbstractLookupMethod
+end
+lookup_table(dsf::DataDrivenMatterModel) = dsf.lt
+lookup_method(dsf::DataDrivenMatterModel) = dsf.method
+
+function DataDrivenMatterModel(datadict::Dict, method::AbstractLookupMethod = InterpolExtrapol())
+    lt = GridLookupTable(
+        datadict["omega_me"],
+        datadict["q_me"],
+        datadict["dsf"]',
+    )
+    return DataDrivenMatterModel(lt, method)
+end
+
+function dynamic_structure_factor(
+        sys::DataDrivenMatterModel,
+        om_q::NTuple{2, T},
+    ) where {T <: Real}
+    om, q = @inbounds om_q
+
+    return lookup(lookup_method(sys), om, q, lookup_table(sys))
+end

src/electron_system/generic.jl

@@ -0,0 +1,39 @@
+function dynamic_response(esys::AbstractElectronSystem, om_q::NTuple{2, T}) where {T <: Real}
+    return real_dynamic_response(esys, om_q) + 1im * imag_dynamic_response(esys, om_q)
+end
+
+
+@inline function _dynamic_structure_factor_pos_om(
+        esys::AbstractElectronSystem,
+        om_q::NTuple{2, T},
+    ) where {T <: Real}
+    _fac = pi * electron_density(esys)
+    imag_rf = imag_dynamic_response(esys, om_q)
+
+    om = @inbounds om_q[1]
+    if iszero(temperature(esys))
+        return -imag_rf / _fac
+    end
+
+    return -imag_rf / ((one(T) - exp(-beta(esys) * om)) * _fac)
+end
+
+function dynamic_structure_factor(
+        esys::AbstractElectronSystem,
+        om_q::NTuple{2, T},
+    ) where {T <: Real}
+    om, q = @inbounds om_q
+    if om < zero(om)
+        return exp(om * beta(esys)) * _dynamic_structure_factor_pos_om(esys, (-om, q))
+    else
+        return _dynamic_structure_factor_pos_om(esys, om_q)
+    end
+end
+
+function static_response(esys::AbstractElectronSystem, q::T) where {T <: Real}
+
+end
+
+function static_structure_factor(esys::AbstractElectronSystem, q::T) where {T <: Real}
+
+end

src/electron_system/ideal/approximations/degenerated.jl

@@ -0,0 +1,68 @@
+### real part
+
+# TODO: refactor this!
+function _deg_integral(x::Real)
+    if abs(x) < 0.95
+        # Numerically stable atanh version
+        return (1 - x^2) * atanh(x)
+    elseif abs(x - 1) < 1.0e-4
+        # Taylor expansion around x = 1
+        δ = x - 1
+        return (δ^2) / 3 + (2δ^3) / 5 - (δ^4) / 7
+    elseif abs(x + 1) < 1.0e-4
+        # Taylor expansion around x = -1
+        δ = x + 1
+        return (δ^2) / 3 - (2δ^3) / 5 + (δ^4) / 7
+    else
+        # Use log1p if safe, otherwise fall back to log(abs(...))
+        log_term = if (1 + x > 0) && (1 - x > 0)
+            log1p(x) - log1p(-x)
+        else
+            log(abs(1 + x)) - log(abs(1 - x))
+        end
+        return 0.5 * (1 - x^2) * log_term
+    end
+end
+
+function _real_lindhard_nonzero_temperature(::Degenerated, ombar::T, qbar::T, bbar::T) where {T <: Real}
+
+    prefac = inv(qbar)
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+
+
+    if qbar < 1.0e-2 * ombar
+        # limiting case for large ombar
+        # source: AB84, table 1, (4,2)
+        prefac2 = -32.0 * qbar^2 * pi / (3 * ombar^2)
+        return prefac2 * (
+            1 + qbar^4 / ombar^2 + 6 * qbar^2 / (bbar * ombar^2) + 60 * qbar^4 / (bbar^2 * ombar^4)
+        )
+    end
+    if ombar < 1.0e-4 * qbar
+        # fallback on no-approx
+        # TODO: find an improved formula for this case!
+        return _real_lindhard_zero_temperature(NoApprox(), ombar, qbar)
+    end
+
+    # general degenerated case
+    # source: AB84, eq. 10a
+    return -prefac * (qbar + _deg_integral(nup) - _deg_integral(num))
+end
+
+### imaginary part
+
+# general degenerated case
+# source: AB84, eq. 26
+function _imag_lindhard_nonzero_temperature(::Degenerated, ombar::T, qbar::T, bbar::T) where {T <: Real}
+    prefac = -pi / (2 * qbar * bbar)
+
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+    mubar = _chemical_potential_normalized(bbar)
+
+    term1 = bbar * (one(num) - num^2)
+    term2 = bbar * (one(nup) - nup^2)
+
+    return prefac * (log1pexp(term1) - log1pexp(term2))
+end

src/electron_system/ideal/approximations/interface.jl

@@ -0,0 +1,19 @@
+abstract type AbstractApproximation end
+Base.broadcastable(approx::AbstractApproximation) = Ref(approx)
+
+abstract type AbstractResponseApproximation <: AbstractApproximation end
+
+# most naive implementation
+struct NoApprox <: AbstractResponseApproximation end
+
+# using integral transforms and identities to make integration more stable
+struct TransformedIntegral <: AbstractResponseApproximation end
+
+# using the limit for non-degenerated electron gases (betabar << 1)
+struct NonDegenerated <: AbstractResponseApproximation end
+
+# using the limit for degenerated electron gases (betabar >> 1)
+struct Degenerated <: AbstractResponseApproximation end
+
+const AbstractZeroTemperatureApproximation = Union{NoApprox}
+const AbstractFiniteTemperatureApproximation = Union{NoApprox, NonDegenerated, Degenerated}

src/electron_system/ideal/approximations/no_approx.jl

@@ -0,0 +1,83 @@
+# zero temperature Lindhard, without any further approximation
+
+function _imag_lindhardDSF_zeroT_minus(omb::T, qb::T) where {T <: Real}
+    num = _nu_minus(omb, qb)
+
+    if abs(num) >= one(T)
+        return zero(T)
+    end
+    return one(omb) - num^2
+
+end
+
+function _imag_lindhardDSF_zeroT_plus(omb::T, qb::T) where {T <: Real}
+    nup = _nu_plus(omb, qb)
+    if abs(nup) >= one(T)
+        return zero(T)
+    end
+    return one(omb) - nup^2
+
+end
+
+function _imag_lindhard_zero_temperature(::NoApprox, ombar::T, qbar::T) where {T <: Real}
+    return -pi / (2 * qbar) *
+        (_imag_lindhardDSF_zeroT_minus(ombar, qbar) - _imag_lindhardDSF_zeroT_plus(ombar, qbar))
+end
+
+function _real_lindhard_zero_temperature(::NoApprox, ombar::T, qbar::T) where {T <: Real}
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+
+    term1 = (1 - num^2) / 4 * _stable_log_term(num)
+    term2 = (1 - nup^2) / 4 * _stable_log_term(nup)
+
+    return -2 * (qbar / 2 - term1 + term2) / qbar
+end
+
+# finite temperature Lindhard, without any further approximation
+
+## Real part
+
+function _general_fermi(x, beta)
+    mu = _chemical_potential_normalized(beta)
+    denom = exp(beta * (x^2 - mu)) + one(beta)
+    return inv(denom)
+end
+
+function _general_integrand(x, nu, beta)
+    return x * _general_fermi(x, beta) * log(abs(x - nu) / abs(x + nu))
+end
+
+function _general_integal(nu, beta)
+    res, _ = quadgk(x -> _general_integrand(x, nu, beta), 0.0, nu, Inf)
+
+    return res
+end
+
+function _real_lindhard_nonzero_temperature(::NoApprox, ombar, qbar, bbar)
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+
+    prefac = inv(qbar)
+
+    res = -prefac * (_general_integal(num, bbar) - _general_integal(nup, bbar))
+
+    return res
+end
+
+## Imaginary part
+
+function _imag_lindhard_nonzero_temperature(::NoApprox, ombar::T, qbar::T, bbar::T) where {T <: Real}
+
+    prefac = -pi / (2 * qbar)
+
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+    mubar = _chemical_potential_normalized(bbar)
+
+    term1 = log1pexp(bbar * (mubar - num^2))
+    term2 = log1pexp(bbar * (mubar - nup^2))
+
+    return prefac * (term1 - term2) / bbar
+
+end

src/electron_system/ideal/approximations/non_degenerated.jl

@@ -0,0 +1,28 @@
+function _non_deg_integal(nu, beta)
+    mu = _chemical_potential_normalized(beta)
+    prefac = -exp(mu * beta) / beta
+
+    return prefac * dawson(nu * sqrt(beta)) * sqrt(pi)
+end
+
+function _real_lindhard_nonzero_temperature(::NonDegenerated, ombar::T, qbar::T, bbar::T) where {T <: Real}
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+
+    prefac = inv(qbar)
+    return -prefac * (_non_deg_integal(num, bbar) - _non_deg_integal(nup, bbar))
+end
+
+
+function _imag_lindhard_nonzero_temperature(::NonDegenerated, ombar::T, qbar::T, bbar::T) where {T <: Real}
+    prefac = -pi / (2 * qbar * bbar)
+
+    num = _nu_minus(ombar, qbar)
+    nup = _nu_plus(ombar, qbar)
+    mubar = _chemical_potential_normalized(bbar)
+
+    term1 = bbar * (mubar - num^2)
+    term2 = bbar * (mubar - nup^2)
+
+    return prefac * (exp(term1) - exp(term2))
+end