diff --git a/Project.toml b/Project.toml index 49c3df6..8e2a086 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,19 @@ name = "MatterModels" uuid = "89b7eb7d-0b8d-5078-86d5-4a0f481f2984" -authors = ["Uwe Hernandez Acosta "] version = "0.1.0" +authors = ["Uwe Hernandez Acosta "] + +[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" diff --git a/src/MatterModels.jl b/src/MatterModels.jl index 16c3a64..97e8695 100644 --- a/src/MatterModels.jl +++ b/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 diff --git a/src/constants.jl b/src/constants.jl new file mode 100644 index 0000000..40910dc --- /dev/null +++ b/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 diff --git a/src/data_driven/impl.jl b/src/data_driven/impl.jl new file mode 100644 index 0000000..d9f121e --- /dev/null +++ b/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 diff --git a/src/electron_system/generic.jl b/src/electron_system/generic.jl new file mode 100644 index 0000000..f291c60 --- /dev/null +++ b/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 diff --git a/src/electron_system/ideal/approximations/degenerated.jl b/src/electron_system/ideal/approximations/degenerated.jl new file mode 100644 index 0000000..6a3e2ff --- /dev/null +++ b/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 diff --git a/src/electron_system/ideal/approximations/interface.jl b/src/electron_system/ideal/approximations/interface.jl new file mode 100644 index 0000000..5731a61 --- /dev/null +++ b/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} diff --git a/src/electron_system/ideal/approximations/no_approx.jl b/src/electron_system/ideal/approximations/no_approx.jl new file mode 100644 index 0000000..0040b29 --- /dev/null +++ b/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 diff --git a/src/electron_system/ideal/approximations/non_degenerated.jl b/src/electron_system/ideal/approximations/non_degenerated.jl new file mode 100644 index 0000000..4a05e30 --- /dev/null +++ b/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 diff --git a/src/electron_system/ideal/generic.jl b/src/electron_system/ideal/generic.jl new file mode 100644 index 0000000..1582c6c --- /dev/null +++ b/src/electron_system/ideal/generic.jl @@ -0,0 +1,38 @@ +function imag_dynamic_response( + elsys::AbstractIdealElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + ombar, qbar = _transform_om_q(elsys, om_q) + approx = response_approximation(elsys) + KF = fermi_wave_vector(elsys) + N0 = KF / (2 * pi^2) + + + if ombar <= zero(ombar) + if qbar < zero(qbar) + return N0 * _imag_ideal_dynamic_response(elsys, approx, -ombar, -qbar) + else + return -N0 * _imag_ideal_dynamic_response(elsys, approx, -ombar, qbar) + end + end + + return N0 * _imag_ideal_dynamic_response(elsys, approx, ombar, qbar) +end + +function real_dynamic_response( + elsys::AbstractIdealElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + ombar, qbar = _transform_om_q(elsys, om_q) + + # TODO: insert symmetries here as well + + KF = fermi_wave_vector(elsys) + N0 = KF / (2 * pi^2) + + if ombar <= zero(ombar) + return N0 * _real_ideal_dynamic_response(elsys, response_approximation(elsys), -ombar, qbar) + end + + return N0 * _real_ideal_dynamic_response(elsys, response_approximation(elsys), ombar, qbar) +end diff --git a/src/electron_system/ideal/impl.jl b/src/electron_system/ideal/impl.jl new file mode 100644 index 0000000..56ab917 --- /dev/null +++ b/src/electron_system/ideal/impl.jl @@ -0,0 +1,199 @@ +# TODO: +# - treat ZeroTemperature mode as approximation (ZeroTemperatureApproximation) + +""" + IdealElectronSystem{TM, T, A} <: AbstractIdealElectronSystem + +Represents a non-interacting (ideal) electron system, parameterized by a temperature mode (`TM`), +a numerical type (`T`), and a response approximation model (`A`). This system can describe either +a zero-temperature or finite-temperature electron gas. + +# Type Parameters +- `TM <: AbstractTemperatureMode`: Temperature mode (`ZeroTemperature` or `FiniteTemperature`). +- `T <: Real`: Internal numeric type (e.g., `Float64`). +- `A <: AbstractResponseApproximation`: Approximation used for the response function. + +# Fields +- `electron_density::T`: Electron number density (stored in atomic units). +- `temperature::T`: Temperature (stored in atomic units). +- `approx::A`: Approximation model for the response function (e.g., `NoApprox()`). + +# Temperature Modes + +- **ZeroTemperature**: Models a zero-temperature Fermi gas. Temperature is internally set to zero. +- **FiniteTemperature**: Models finite-temperature effects in the electron gas. + +# Constructors + + IdealElectronSystem{TM}(electron_density::T1, temperature::T2, approx::A) + +Construct an ideal electron system with a specified temperature mode, electron density, temperature, +and approximation model. + + IdealElectronSystem{ZeroTemperature}(electron_density::T1, approx::A) + +Construct a zero-temperature system with a custom approximation. + + IdealElectronSystem{ZeroTemperature}(electron_density::T1) + +Construct a zero-temperature system using the default approximation (`NoApprox()`). + + IdealElectronSystem(electron_density::T1, temperature::T2, approx::A) + +Construct a finite-temperature system. The temperature mode is inferred as `FiniteTemperature`. + + IdealElectronSystem{FiniteTemperature}(electron_density::T1, temperature::T2) + +Construct a finite-temperature system with the default approximation (`NoApprox()`). + + IdealElectronSystem(electron_density::T1, temperature::T2) + +Alias for `IdealElectronSystem{FiniteTemperature}(...)`. + +# Examples + +```julia +using Unitful + +# Zero-temperature system with default approximation +sys1 = IdealElectronSystem{ZeroTemperature}(1e23u"cm^(-3)") + +# Finite-temperature system with explicit approximation +sys2 = IdealElectronSystem(1e23u"cm^(-3)", 300.0u"K", NonDegenerated()) + +# Finite-temperature system with default approximation +sys3 = IdealElectronSystem(1e23u"cm^(-3)", 300.0u"K") +``` +""" +struct IdealElectronSystem{TM, T, A <: AbstractResponseApproximation} <: AbstractIdealElectronSystem + electron_density::T + temperature::T + approx::A + + function IdealElectronSystem{TM}( + electron_density::T1, + temperature::T2, + approx::A + ) where { + T <: Real, + T1 <: Union{T, Quantity{T}}, + T2 <: Union{T, Quantity{T}}, + TM <: AbstractTemperatureMode, + A <: AbstractResponseApproximation, + } + + ne_internal = _internalize_density(electron_density) + temp_internal = _internalize_temperature(temperature) + + return new{TM, T, A}(ne_internal, temp_internal) + end +end + +# TODO: think about making this a default constructor +function IdealElectronSystem{ZeroTemperature}( + electron_density::T1, + approx::A + ) where { + T <: Real, + T1 <: Union{T, Quantity{T}}, + A <: AbstractZeroTemperatureApproximation, + } + return IdealElectronSystem{ZeroTemperature}(electron_density, zero(T), approx) +end + +function IdealElectronSystem{ZeroTemperature}( + electron_density::T1, + ) where { + T <: Real, + T1 <: Union{T, Quantity{T}}, + } + return IdealElectronSystem{ZeroTemperature}(electron_density, zero(T), NoApprox()) +end + +function IdealElectronSystem( + electron_density::T1, + temperature::T2, + approx::A + ) where { + T <: Real, + T1 <: Union{T, Quantity{T}}, + T2 <: Union{T, Quantity{T}}, + A <: AbstractFiniteTemperatureApproximation, + } + return IdealElectronSystem{FiniteTemperature}(electron_density, temperature, approx) +end + +function IdealElectronSystem{FiniteTemperature}( + electron_density::T1, + temperature::T2, + ) where { + T <: Real, + T1 <: Union{T, Quantity{T}}, + T2 <: Union{T, Quantity{T}}, + } + return IdealElectronSystem{FiniteTemperature}(electron_density, temperature, NoApprox()) +end + +function IdealElectronSystem( + electron_density::T1, + temperature::T2, + ) where { + T <: Real, + T1 <: Union{T, Quantity{T}}, + T2 <: Union{T, Quantity{T}}, + } + return IdealElectronSystem{FiniteTemperature}(electron_density, temperature, NoApprox()) +end + +@inline temperature(elsys::IdealElectronSystem{ZeroTemperature, T}) where {T <: Real} = zero(T) +@inline temperature(elsys::IdealElectronSystem{FiniteTemperature, T}) where {T <: Real} = + elsys.temperature +@inline electron_density(elsys::IdealElectronSystem) = elsys.electron_density +@inline response_approximation(elsys::IdealElectronSystem) = elsys.approx + +### default implementations -> needed for separation between zero and nonzero case + +# Zero temperature case + +function _imag_ideal_dynamic_response( + elsys::IdealElectronSystem{ZeroTemperature}, + approx::AbstractZeroTemperatureApproximation, + ombar::T, + qbar::T + ) where {T <: Real} + + + return _imag_lindhard_zero_temperature(approx, ombar, qbar) +end + +function _real_ideal_dynamic_response( + elsys::IdealElectronSystem{ZeroTemperature}, + approx::AbstractZeroTemperatureApproximation, + ombar::T, + qbar::T + ) where {T <: Real} + + return _real_lindhard_zero_temperature(approx, ombar, qbar) +end + + +# finite temperature case + +function _imag_ideal_dynamic_response( + elsys::IdealElectronSystem{FiniteTemperature}, + approx::AbstractFiniteTemperatureApproximation, + ombar::T, + qbar::T + ) where {T <: Real} + + return _imag_lindhard_nonzero_temperature(approx, ombar, qbar, betabar(elsys)) +end + +function _real_ideal_dynamic_response( + elsys::IdealElectronSystem{FiniteTemperature}, + approx::AbstractFiniteTemperatureApproximation, + ombar::T, + qbar::T + ) where {T <: Real} + return _real_lindhard_nonzero_temperature(approx, ombar, qbar, betabar(elsys)) +end diff --git a/src/electron_system/ideal/interface.jl b/src/electron_system/ideal/interface.jl new file mode 100644 index 0000000..a36c693 --- /dev/null +++ b/src/electron_system/ideal/interface.jl @@ -0,0 +1,46 @@ +#TBW: +# - implement an interface for ideal electron gases +# - consider adding `approximation` to the interface +# - think about generics on this interface + +abstract type AbstractIdealElectronSystem <: AbstractProperElectronSystem end + +""" + response_approximation(sys::AbstractIdealElectronSystem) + +Return the response approximation model associated with the given ideal electron system. +""" +function response_approximation end + +""" + _imag_ideal_dynamic_response( + sys::AbstractIdealElectronSystem, + approx::AbstractResponseApproximation, + ombar::T, + qbar::T + ) where {T<:Real} + +Compute the imaginary part of the dynamic response function of an ideal electron system +using the given response approximation. The input frequencies `ombar` and wavevectors `qbar` +are given in dimensionless units (normalized to Fermi energy and Fermi wavevector, respectively). + +This is a low-level method, typically implemented by concrete approximation types. +""" +function _imag_ideal_dynamic_response end + +""" + + real_ideal_dynamic_response( + sys::AbstractIdealElectronSystem, + approx::AbstractResponseApproximation, + ombar::T, + qbar::T + ) where {T<:Real} + +Compute the real part of the dynamic response function of an ideal electron system +using the given response approximation. The input frequencies `ombar` and wavevectors `qbar` +are given in dimensionless units (normalized to Fermi energy and Fermi wavevector, respectively). + +This is a low-level interface meant for implementation by specific approximation models. +""" +function _real_ideal_dynamic_response end diff --git a/src/electron_system/ideal/utils.jl b/src/electron_system/ideal/utils.jl new file mode 100644 index 0000000..7731b95 --- /dev/null +++ b/src/electron_system/ideal/utils.jl @@ -0,0 +1 @@ +# TBW diff --git a/src/electron_system/interacting/generic.jl b/src/electron_system/interacting/generic.jl new file mode 100644 index 0000000..69a4f12 --- /dev/null +++ b/src/electron_system/interacting/generic.jl @@ -0,0 +1,126 @@ +# delegations + +""" + local_effective_potential(iesys::AbstractInteractingElectronSystem, om_q) -> Real + +Delegates to the `local_effective_potential` method of the screening model associated with `iesys`. + +# Arguments +- `iesys`: An interacting electron system. +- `om_q`: Tuple `(ω, q)` representing frequency and wavevector magnitude. + +# Returns +- Local effective potential `V_eff(q, ω)` as a real number. +""" +@inline function local_effective_potential( + iesys::AbstractInteractingElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + return local_effective_potential(screening(iesys), om_q) +end + +""" + pseudo_potential(iesys::AbstractInteractingElectronSystem, om_q) -> Real + +Delegates to the `pseudo_potential` method of the screening model associated with `iesys`. + +# Arguments +- `iesys`: An interacting electron system. +- `om_q`: Tuple `(ω, q)` representing frequency and wavevector magnitude. + +# Returns +- Pseudo potential `V(q, ω)` as a real number. +""" +@inline function pseudo_potential( + iesys::AbstractInteractingElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + return pseudo_potential(screening(iesys), om_q) +end + +""" + dielectric_function(iesys::AbstractInteractingElectronSystem, om_q) -> Real + +Delegates to the `dielectric_function` using the proper electron system and screening model associated with `iesys`. + +# Arguments +- `iesys`: An interacting electron system. +- `om_q`: Tuple `(ω, q)` representing frequency and wavevector magnitude. + +# Returns +- Dielectric function `ϵ(q, ω)` as a real number. +""" +@inline function dielectric_function( + iesys::AbstractInteractingElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + return dielectric_function(proper_electron_system(iesys), screening(iesys), om_q) +end + +# response functions +""" + dynamic_response(iesys::AbstractInteractingElectronSystem, om_q) -> Complex + +Compute the dynamically screened response function of the interacting electron system. + +The response is computed as + +```math +χ(q, ω) = \\frac{χ_0(q, ω)}{ϵ(q, ω)} +``` + +where χ₀ is the proper response and ϵ is the dielectric function derived from the screening model. + +# Arguments + +- `iesys`: An interacting electron system. +- `om_q`: Tuple (ω, q) representing frequency and wavevector magnitude. + +# Returns + +Complex-valued dynamic response χ(q, ω). +""" +function dynamic_response( + iesys::AbstractInteractingElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + + # TODO: think about implementing this directly with imag and real part + # -> this could also be used directly for imag and real part of chi_RPA + rf = dynamic_response(proper_electron_system(iesys), om_q) + lep = local_effective_potential(screening(iesys), om_q) + eps = _dielectric_function(rf, lep) + return rf / eps +end + +""" + real_dynamic_response(iesys::AbstractInteractingElectronSystem, om_q) -> Real + +Return the real part of the dynamic response function of the interacting electron system. + +# Returns +- Real part of `χ(q, ω)`. +""" +function real_dynamic_response( + iesys::AbstractInteractingElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + rf = dynamic_response(iesys, om_q) + return real(rf) +end + +""" + imag_dynamic_response(iesys::AbstractInteractingElectronSystem, om_q) -> Real + +Return the imaginary part of the dynamic response function of the interacting electron system. + +# Returns +- Imaginary part of `χ(q, ω)`. +""" +function imag_dynamic_response( + iesys::AbstractInteractingElectronSystem, + om_q::NTuple{2, T}, + ) where {T <: Real} + rf = dynamic_response(iesys, om_q) + return imag(rf) +end diff --git a/src/electron_system/interacting/impl.jl b/src/electron_system/interacting/impl.jl new file mode 100644 index 0000000..1101e46 --- /dev/null +++ b/src/electron_system/interacting/impl.jl @@ -0,0 +1,69 @@ +""" + InteractingElectronSystem{PES, S} <: AbstractInteractingElectronSystem + +Represents an interacting electron system consisting of a proper electron system (e.g., a free or ideal Fermi gas) +and a screening model that modifies its response. + +# Type Parameters +- `PES <: AbstractProperElectronSystem`: Type of the proper (non-interacting) electron system. +- `S <: AbstractScreening`: Type of the screening model. + +# Fields +- `proper_elec_system::PES`: The underlying proper electron system. +- `screening::S`: The screening model used to compute interactions. + +# Constructors + +```julia +InteractingElectronSystem(pelsys::PES, scr::S) +``` + +Construct an interacting electron system directly from a proper electron system and a screening model. + +```julia +InteractingElectronSystem{PES}(electron_density::T1, temp::T2, scr::S) +``` +Construct an interacting system by providing physical parameters, where `PES` is the desired proper electron system type. + +# Examples + +```julia +iesys = InteractingElectronSystem{IdealFermiGas}(1e23, 300.0, Screening()) +``` +""" +struct InteractingElectronSystem{PES <: AbstractProperElectronSystem, S <: AbstractScreening} <: + AbstractInteractingElectronSystem + proper_elec_system::PES + screening::S + + function InteractingElectronSystem( + pelsys::PES, + scr::S, + ) where {PES <: AbstractProperElectronSystem, S <: AbstractScreening} + + return new{PES, S}(pelsys, scr) + end + + function InteractingElectronSystem{PES}( + electron_density::T1, + temp::T2, + scr::S, + ) where { + T <: Real, + T1 <: Union{T, Quantity}, + T2 <: Union{T, Quantity}, + PES <: AbstractProperElectronSystem, + S <: AbstractScreening, + } + proper_elec_system = PES(electron_density, temp) + return new{PES, S}(proper_elec_system, scr) + end +end + +# accessors +@inline proper_electron_system(elsys::InteractingElectronSystem) = elsys.proper_elec_system +@inline screening(elsys::InteractingElectronSystem) = elsys.screening +@inline temperature(elsys::InteractingElectronSystem) = + temperature(proper_electron_system(elsys)) +@inline electron_density(elsys::InteractingElectronSystem) = + electron_density(proper_electron_system(elsys)) diff --git a/src/electron_system/interacting/interface.jl b/src/electron_system/interacting/interface.jl new file mode 100644 index 0000000..f6386b1 --- /dev/null +++ b/src/electron_system/interacting/interface.jl @@ -0,0 +1,36 @@ +""" + + AbstractInteractingElectronSystem + +Abstract base type for electronic system with interacting electrons by introducing screening potentials. +Interface functions to be implemented are: + +1. Matter model interface: + - `electron_density(::AbstractElectronSystem)::Real` + - `temperature(::AbstractElectronSystem)::Real` + +2. Electronic system interface: + - `imag_dynamic_response(::AbstractElectronSystem, om_q::NTuple{2,Real})::Real` + - `real_dynamic_response(::AbstractElectronSystem, om_q:: + +3. Interacting electronic system interface: + - `proper_electron_system(::AbstractInteractingElectronSystem)::AbstractProperElectronSystem` + - `screening(::AbstractInteractingElectronSystem)::AbstractScreening` +""" +abstract type AbstractInteractingElectronSystem <: AbstractElectronSystem end + +""" + + proper_electron_system(::AbstractInteractingElectronSystem) + +Interface function: return the proper electron system associated with the interacting system. +""" +function proper_electron_system end + +""" + + screening(::AbstractInteractingElectronSystem) + +Interface function: return the screening associated with the interacting system. +""" +function screening end diff --git a/src/electron_system/interacting/screening/generic.jl b/src/electron_system/interacting/screening/generic.jl new file mode 100644 index 0000000..79ff383 --- /dev/null +++ b/src/electron_system/interacting/screening/generic.jl @@ -0,0 +1,42 @@ +# local effective potential +@inline function _local_effective_potential(vp, lfc) + return vp * (one(lfc) - lfc) +end + +function local_effective_potential( + scr::AbstractScreening, + om_q::NTuple{2, T}, + ) where {T <: Real} + vp = pseudo_potential(scr, om_q) + lfc = local_field_correction(scr, om_q) + + return _local_effective_potential(vp, lfc) +end + +""" + + dielectric_function(::AbstractIdealElectronSystem,::AbstractScreening,om_q) + +Return the value of the dielectric function for the given proper response function and +given local effective potential. +""" +function _dielectric_function(rf::Number, lep::Number) + temp = lep * rf + return one(temp) - temp +end + +""" + + dielectric_function(::AbstractIdealElectronSystem,::AbstractScreening,om_q) + +Return the value of the dielectric function for the given electron system and given screening. +""" +function dielectric_function( + esys::AbstractProperElectronSystem, + scr::AbstractScreening, + om_q::NTuple{2, T}, + ) where {T <: Real} + rf = dynamic_response(esys, om_q) + lep = local_effective_potential(scr, om_q) + return _dielectric_function(rf, lep) +end diff --git a/src/electron_system/interacting/screening/impl.jl b/src/electron_system/interacting/screening/impl.jl new file mode 100644 index 0000000..15eeff9 --- /dev/null +++ b/src/electron_system/interacting/screening/impl.jl @@ -0,0 +1,88 @@ +# local field correction +""" + NoLocalFieldCorrection <: AbstractLocalFieldCorrection + +Represents a model with no local field correction (i.e., the correction is zero). +""" +struct NoLocalFieldCorrection <: AbstractLocalFieldCorrection end + +function _compute(::NoLocalFieldCorrection, om_q::NTuple{2, T}) where {T <: Real} + return zero(T) +end + +# pseudo potential +""" + CoulombPseudoPotential <: AbstractPseudoPotential + +Represents the bare Coulomb pseudo potential: V(q) = e² / q². +""" +struct CoulombPseudoPotential <: AbstractPseudoPotential end + +# TODO: insert 4pi! +function _compute(::CoulombPseudoPotential, om_q::NTuple{2, T}) where {T <: Real} + om, q = om_q + return ELEMENTARY_CHARGE_SQUARED / q^2 +end + +### screening +# no screening +""" + NoScreening <: AbstractScreening + +A screening model with no dielectric screening; dielectric function is unity, and corrections vanish. +""" +struct NoScreening <: AbstractScreening end + +""" + dielectric_function(esys, scr::NoScreening, om_q) -> Real + +Return the dielectric function for a `NoScreening` model, which is always 1. +""" +dielectric_function( + esys::AbstractProperElectronSystem, + scr::NoScreening, + om_q::NTuple{2, T}, +) where {T <: Real} = one(T) + +""" + local_effective_potential(scr::NoScreening, om_q) -> Real + +Return the local effective potential for the `NoScreening` model, which is zero. +""" +local_effective_potential(src::NoScreening, om_q::NTuple{2, T}) where {T <: Real} = zero(T) + +""" + pseudo_potential(scr::NoScreening, om_q) -> Real + +Return the pseudo potential for the `NoScreening` model, which is zero. +""" +@inline pseudo_potential(scr::NoScreening, om_q::NTuple{2, T}) where {T <: Real} = zero(T) + +""" + local_field_correction(scr::NoScreening, om_q) -> Real + +Return the local field correction for the `NoScreening` model, which is zero. +""" +@inline local_field_correction(scr::NoScreening, om_q::NTuple{2, T}) where {T <: Real} = + zero(T) + +# screeing with pseudo potential and local field correction +""" + Screening{PP, LFC} <: AbstractScreening + +A general screening model that combines a pseudo potential and a local field correction. + +# Fields +- `pseudo_potential`: An instance of `AbstractPseudoPotential`. +- `lfc`: An instance of `AbstractLocalFieldCorrection`. +""" +struct Screening{PP <: AbstractPseudoPotential, LFC <: AbstractLocalFieldCorrection} <: + AbstractScreening + pseudo_potential::PP + lfc::LFC +end + +Screening() = Screening(CoulombPseudoPotential(), NoLocalFieldCorrection()) + +@inline pseudo_potential(scr::Screening, om_q) = _compute(scr.pseudo_potential, om_q) +@inline local_field_correction(scr::Screening, om_q) = _compute(scr.lfc, om_q) diff --git a/src/electron_system/interacting/screening/interface.jl b/src/electron_system/interacting/screening/interface.jl new file mode 100644 index 0000000..b7e65ba --- /dev/null +++ b/src/electron_system/interacting/screening/interface.jl @@ -0,0 +1,29 @@ +""" + + AbstractScreening + +Abstract base type for screening models. Interface functions, which must be implemented: + +- `pseudo_potential(::AbstractInteractingElectronSystem)::AbstractPseudoPotential` +- `local_field_correction(::AbstractInteractingElectronSystem)::AbstractLocalFieldCorrection` + +""" +abstract type AbstractScreening end +abstract type AbstractPseudoPotential end +abstract type AbstractLocalFieldCorrection end + +""" + + pseudo_potential(::AbstractScreening, om_q::Tuple) + +Interface function: return value of the pseudo potential for interacting electron systems +""" +function pseudo_potential end + +""" + + local_field_correction(::AbstractScreening, om_q::Tuple) + +Interface function: return value of the local field correction for interacting electron systems +""" +function local_field_correction end diff --git a/src/electron_system/interface.jl b/src/electron_system/interface.jl new file mode 100644 index 0000000..36442c8 --- /dev/null +++ b/src/electron_system/interface.jl @@ -0,0 +1,42 @@ +### general electron systems + + +""" + AbstractElectronSystem <: AbstractMatterModel + +Abstract base type for electron matter models. + +## Interface requirements + +Concrete subtypes must implement the following functions: + +1. Matter model interface: +- `electron_density(sys::AbstractElectronSystem)::Real` +- `temperature(sys::AbstractElectronSystem)::Real` + +2. Electronic system interface: +- `imag_dynamic_response(sys::AbstractElectronSystem, om_q::NTuple{2, Real})::Real` +- `real_dynamic_response(sys::AbstractElectronSystem, om_q::NTuple{2, Real})::Real` +""" +abstract type AbstractElectronSystem <: AbstractMatterModel end + +""" + imag_dynamic_response(sys::AbstractElectronSystem, om_q::NTuple{2, T}) where {T<:Real} + +Return the imaginary part of the dynamic response function for the given electron system `sys` +at frequency and wavevector specified by `om_q = (ω, q)`. +""" +function imag_dynamic_response end + +""" + real_dynamic_response(sys::AbstractElectronSystem, om_q::NTuple{2, T}) where {T<:Real} + +Return the real part of the dynamic response function for the given electron system `sys` +at frequency and wavevector specified by `om_q = (ω, q)`. +""" +function real_dynamic_response end + +### proper (non-interacting) electron system +# - electron system without screening + +abstract type AbstractProperElectronSystem <: AbstractElectronSystem end diff --git a/src/electron_system/utils.jl b/src/electron_system/utils.jl new file mode 100644 index 0000000..4001ff4 --- /dev/null +++ b/src/electron_system/utils.jl @@ -0,0 +1,47 @@ +function _nu_minus(ombar, qbar) + return (ombar - qbar^2) / (2 * qbar) +end + +function _nu_plus(ombar, qbar) + return (ombar + qbar^2) / (2 * qbar) +end + +function heaviside(t) + return 0.5 * (sign(t) + 1) +end + +# Ichimaru:2018 +function _chemical_potential_normalized(bbar) + A = 0.25954 + B = 0.072 + b = 0.858 + + term1 = 3 / 2 * log(bbar) + term2 = log(4 / (3 * sqrt(pi))) + term3 = (A * bbar^(b + 1) + B * bbar^((b + 1) / 2)) / (1 + A * bbar^b) + + return (term1 + term2 + term3) / bbar +end + +# GV:2005 eq. (4.43) +function _F(x, bbar) + mubar = _chemical_potential_normalized(bbar) + denom = exp(bbar * (x^2 - mubar)) + one(bbar) + return x / denom +end + + +""" + + _stable_log_term(x::Real) + +Stable version of + +```math +\\log\\left(\\left\\vert\\frac{1+x}{1-x}\\right\\vert\\right) +``` + +""" +function _stable_log_term(x::Real) + return abs(x) < 1 ? 2 * atanh(x) : 2 * atanh(inv(x)) +end diff --git a/src/generic.jl b/src/generic.jl new file mode 100644 index 0000000..da7c475 --- /dev/null +++ b/src/generic.jl @@ -0,0 +1,68 @@ +@inline _fermi_wave_vector(ne::T) where {T <: Real} = cbrt(3 * pi^2 * ne) + +""" + fermi_wave_vector(esys::AbstractMatterModel) -> Real + +Compute the Fermi wave vector for a given electronic system `esys`, based on its electron density. + +# Arguments +- `esys`: An instance of `AbstractMatterModel` representing the electronic system. + +# Returns +- Fermi wave vector `kF` as a real number. +""" +@inline fermi_wave_vector(esys::AbstractMatterModel) = + _fermi_wave_vector(electron_density(esys)) + +@inline _fermi_energy_from_kF(kF::T) where {T <: Real} = kF^2 / 2 + +""" + fermi_energy(esys::AbstractMatterModel) -> Real + +Compute the Fermi energy of a given electronic system `esys`, using its Fermi wave vector. + +# Arguments +- `esys`: An instance of `AbstractMatterModel` representing the electronic system. + +# Returns +- Fermi energy `EF` as a real number. +""" +@inline fermi_energy(esys::AbstractMatterModel) = + _fermi_energy_from_kF(fermi_wave_vector(esys)) + +""" + beta(esys::AbstractMatterModel) -> Real + +Compute the inverse temperature β = 1 / T of the system `esys`. + +# Arguments +- `esys`: An instance of `AbstractMatterModel` representing the electronic system. + +# Returns +- Inverse temperature `β` as a real number. +""" +@inline beta(esys::AbstractMatterModel) = inv(temperature(esys)) + +""" + betabar(esys::AbstractMatterModel) -> Real + +Compute the dimensionless inverse temperature \$\\bar\\beta= E_F / T\$ of the system `esys`. + +# Arguments +- `esys`: An instance of `AbstractMatterModel` representing the electronic system. + +# Returns +- Dimensionless inverse temperature `β̄` as a real number. +""" +@inline betabar(esys::AbstractMatterModel) = fermi_energy(esys) * beta(esys) + +function _transform_om_q(elsys::AbstractMatterModel, om_q::NTuple{2, T}) where {T <: Real} + kF = fermi_wave_vector(elsys) + EF = fermi_energy(elsys) + + om, q = om_q + ombar = om / EF + qbar = q / kF + + return ombar, qbar +end diff --git a/src/interface.jl b/src/interface.jl new file mode 100644 index 0000000..cbd8097 --- /dev/null +++ b/src/interface.jl @@ -0,0 +1,34 @@ +""" + AbstractMatterModel + +Abstract base type representing a generic matter model. + +Concrete subtypes should implement interface functions such as `temperature` and `electron_density` +to provide physical properties of the matter system. +""" +abstract type AbstractMatterModel end + +""" + temperature(m::AbstractMatterModel) + +Return the temperature of the given matter model `m` in internal units. + +This is an interface function that must be implemented by all concrete matter models. +""" +function temperature end + +""" + electron_density(m::AbstractMatterModel) + +Return the electron density of the given matter model `m`. + +This is an interface function that must be implemented by all concrete matter models. +""" +function electron_density end + +# TODO: +# - consider moving `electron_density` downstream, maybe to the interface for electron +# systems +# - implement `temperature(::Unit,::AbstactMatterModel)` to return the temperature in a +# user defined unit; similar for `electron_density` +# - propagate the unitful call downstream, e.g. to `beta` diff --git a/src/lookup.jl b/src/lookup.jl new file mode 100644 index 0000000..e62a190 --- /dev/null +++ b/src/lookup.jl @@ -0,0 +1,230 @@ +""" + GridLookupTable( + cols::T1, + rows::T2, + values::T3 + ) where {T1<:AbstractVector{T},T2<:AbstractVector{T},T3<:AbstractMatrix{T}} where {T<:Real} + +TBW +""" +struct GridLookupTable{T1, T2, T3} + rows::T1 + cols::T2 + values::T3 + + function GridLookupTable( + cols::T1, rows::T2, values::T3 + ) where { + T1 <: AbstractVector{T}, T2 <: AbstractVector{T}, T3 <: AbstractMatrix{T}, + } where {T <: Real} + len_arg_rows = length(rows) + len_arg_cols = length(cols) + len_val_rows, len_val_cols = size(values) + + len_arg_cols == len_val_cols || throw( + ArgumentError( + "Length of column arguments <$len_arg_cols> does not match the the first dimension of the values matrix <$len_val_cols>.", + ), + ) + + len_arg_rows == len_val_rows || throw( + ArgumentError( + "Length of row arguments <$len_arg_rows> does not match the the second dimension of the values matrix <$len_val_rows>.", + ), + ) + + # TODO: implement check, if rows and cols are in ascending order + + return new{T1, T2, T3}(rows, cols, values) + end +end + +# TODO: +# - change signature of lookup to lookup(GLT, method, column,row) +# - check behaviour for sampling outside of the lookup +# - add some I/O for lookup tables (simple saving and loading) + +######## +# extrapolation methods +######## +abstract type AbstractLookupMethod end +struct InterpolExtrapol <: AbstractLookupMethod end +struct InterpolEndValue <: AbstractLookupMethod end +struct NearestInput <: AbstractLookupMethod end +struct BelowInput <: AbstractLookupMethod end +struct AboveInput <: AbstractLookupMethod end + +""" + lookup( + method::InterpolExtrapol, + column::Real, + row::Real, + table::GridLookupTable) + +TBW +""" +function lookup(method::InterpolExtrapol, column::Real, row::Real, table::GridLookupTable) + # Binary search in rows vector + rowlow = searchsortedfirst(table.rows, row) + columnlow = searchsortedfirst(table.cols, column) + + + # Behaviour according to method + if rowlow >= length(table.rows) + rowlow = length(table.rows) - 1 + end + if columnlow >= length(table.cols) + columnlow = length(table.cols) - 1 + end + + rowhigh = rowlow + 1 + columnhigh = columnlow + 1 + # Calculate column and row factors + rowfactor = (row - table.rows[rowlow]) / (table.rows[rowhigh] - table.rows[rowlow]) + columnfactor = + (column - table.cols[columnlow]) / (table.cols[columnhigh] - table.cols[columnlow]) + # For row low + outputcolumnΔlow = table.values[rowlow, columnhigh] - table.values[rowlow, columnlow] + # For row high + outputcolumnΔhigh = table.values[rowhigh, columnhigh] - table.values[rowhigh, columnlow] + # For this column + outputrowΔ = + table.values[rowhigh, columnlow] + outputcolumnΔhigh * columnfactor - + (table.values[rowlow, columnlow] + outputcolumnΔlow * columnfactor) + return table.values[rowlow, columnlow] + + columnfactor * outputcolumnΔlow + + rowfactor * outputrowΔ +end + +""" + lookup( + method::InterpolEndValue, + column::Real, + row::Real, + table::GridLookupTable) + +TBW +""" +function lookup(method::InterpolEndValue, column::Real, row::Real, table::GridLookupTable) + # Binary search in rows vector + rowlow = searchsortedfirst(table.rows, row) + columnlow = searchsortedfirst(table.cols, column) + + if rowlow == length(table.rows) + rowlow -= 1 + end + if columnlow == length(table.cols) + columnlow -= 1 + end + rowhigh = rowlow + 1 + columnhigh = columnlow + 1 + if column >= table.cols[columnhigh] + columnfactor = 0 + outputcolumnΔlow = 0 + columnlow = columnhigh + elseif column <= table.cols[columnlow] + columnfactor = 0 + outputcolumnΔlow = 0 + else + # For row low + outputcolumnΔlow = + table.values[rowlow, columnhigh] - table.values[rowlow, columnlow] + end + if row >= table.rows[rowhigh] + rowfactor = 0 + outputrowΔ = 0 + rowlow = rowhigh + elseif row <= table.rows[rowlow] + rowfactor = 0 + outputrowΔ = 0 + else + # For row high + outputcolumnΔhigh = + table.values[rowhigh, columnhigh] - table.values[rowhigh, columnlow] + # Calculate row factor + rowfactor = (row - table.rows[rowlow]) / (table.rows[rowhigh] - table.rows[rowlow]) + # For this column + outputrowΔ = + table.values[rowhigh, columnlow] + outputcolumnΔhigh * columnfactor - + (table.values[rowlow, columnlow] + outputcolumnΔlow * columnfactor) + end + return table.values[rowlow, columnlow] + + columnfactor * outputcolumnΔlow + + rowfactor * outputrowΔ +end + +""" + lookup( + method::NearestInput, + column::Real, + row::Real, + table::GridLookupTable) + +TBW +""" +function lookup(method::NearestInput, column::Real, row::Real, table::GridLookupTable) + # Binary search in rows vector + rowlow = searchsortedfirst(table.rows, row) + columnlow = searchsortedfirst(table.cols, column) + if rowlow == length(table.rows) + rowlow -= 1 + end + if columnlow == length(table.cols) + columnlow -= 1 + end + rowhigh = rowlow + 1 + columnhigh = columnlow + 1 + if row < (table.rows[rowlow] + table.rows[rowhigh]) / 2 + nearestrow = rowlow + else + nearestrow = rowhigh + end + if column < (table.cols[columnlow] + table.cols[columnhigh]) / 2 + nearestcolumn = columnlow + else + nearestcolumn = columnhigh + end + return table.values[nearestrow, nearestcolumn] +end + +""" + lookup( + method::BelowInput, + column::Real, + row::Real, + table::GridLookupTable) + +TBW +""" +function lookup(method::BelowInput, column::Real, row::Real, table::GridLookupTable) + # Binary search in rows vector + rowlow = searchsortedfirst(table.rows, row) + columnlow = searchsortedfirst(table.cols, column) + return table.values[rowlow, columnlow] +end + +""" + lookup( + method::AboveInput, + column::Real, + row::Real, + table::GridLookupTable) + +TBW +""" +function lookup(method::AboveInput, column::Real, row::Real, table::GridLookupTable) + # Binary search in rows vector + rowlow = searchsortedfirst(table.rows, row) + columnlow = searchsortedfirst(table.cols, column) + if table.rows[rowlow] == row + rowhigh = rowlow + else + rowhigh = rowlow + 1 + end + if table.cols[columnlow] == column + columnhigh = columnlow + else + columnhigh = columnlow + 1 + end + return table.values[rowhigh, columnhigh] +end diff --git a/src/temperature.jl b/src/temperature.jl new file mode 100644 index 0000000..92ccbb8 --- /dev/null +++ b/src/temperature.jl @@ -0,0 +1,4 @@ +# Temperature modes +abstract type AbstractTemperatureMode end +struct FiniteTemperature <: AbstractTemperatureMode end +struct ZeroTemperature <: AbstractTemperatureMode end diff --git a/src/units.jl b/src/units.jl new file mode 100644 index 0000000..7347a4f --- /dev/null +++ b/src/units.jl @@ -0,0 +1,109 @@ +# TODO: +# - refactor this! +# - add more units with energy dimention (current, time, momentum?) +# - move this to `NaturallyUnitful` (`naturally` is type unstable!) + +const speed_of_light = 1u"c" +const hbar = Unitful.ħ +const kBoltzmann = Unitful.k +const hbarc = hbar * speed_of_light + +function _get_dim_power(x::Unitful.Dimensions{T}) where {T} + return prod(_get_dim_power.(T)) +end + +_get_dim_power(x::Unitful.Dimension{:Length}) = -x.power +_get_dim_power(x::Unitful.Dimension{:Time}) = -x.power +_get_dim_power(x::Unitful.Dimension{:Mass}) = x.power +_get_dim_power(x::Unitful.Dimension{:Temperature}) = x.power +#_get_dim_power(x::Unitful.Dimension{:Current}) = x.power + +function _get_dim_factor(x::Unitful.Dimensions{T}) where {T} + return prod(@. _get_dim_factor(T, _get_dim_power(T))) +end + +function _get_dim_factor(x::Unitful.Dimension{:Length}, p::Rational) + return hbarc^p +end + +function _get_dim_factor(x::Unitful.Dimension{:Mass}, p::Rational) + return speed_of_light^(2 * p) +end + +function _get_dim_factor(x::Unitful.Dimension{:Temperature}, p::Rational) + return kBoltzmann^p +end + + +function mynatural(q; base = u"eV") + d = dimension(q) + p = _get_dim_power(d) + nunit = base^(p) + nfac = _get_dim_factor(d) + return uconvert(nunit, q * nfac) +end + +# opt out, if the quanity is already in eV +gettypeparams(::Unitful.FreeUnits{T, U, V}) where {T, U, V} = T, U, V +const eV_DIM = gettypeparams(u"GeV")[2] + +function mynatural(q::Quantity{T, eV_DIM}; base = u"eV") where {T <: Number} + return uconvert(base, q) +end + +# FIXME: generalize to powers `eV^p`, currently it errors for `p!=1 + +## conversion for structure factors + +const me = ustrip(mynatural(Unitful.me)) + +function _eV2me_dimless(t_eV, p::Int = 1) + return ustrip(t_eV) * me^p +end + +function _density2me_dimless(density) + density_eV = mynatural(density) + return _eV2me_dimless(density_eV, -3) +end + +function _energy2me_dimless(energy) + energy_eV = mynatural(energy) + return _eV2me_dimless(energy_eV, -1) +end + +# simple unit transforms +# TODO: +# - write unit tests for this +# - if other unit transforms pop up, add them here +# +_bohr2ang(x_bohr, p = 1) = x_bohr * (BOHR_RADIUS_ANG^p) +_bohr2cm(x_bohr, p = 1) = x_bohr * ((BOHR_RADIUS_ANG * 1.0e-8)^p) +_bohr2fm(x_bohr, p = 1) = x_bohr * ((BOHR_RADIUS_ANG * 1.0e5)^p) +_bohr2inv_eV(x_bohr, p = 1) = x_bohr * ((BOHR_RADIUS_ANG / HBARC_eV_ANG)^p) +_bohr2inv_me(x_bohr, p = 1) = x_bohr * ((BOHR_RADIUS_ANG / HBARC_eV_ANG * ELECTRONMASS)^p) + +_inv_bohr2ang(x_inv_bohr, p = 1) = x_inv_bohr / (BOHR_RADIUS_ANG^p) +_inv_bohr2cm(x_inv_bohr, p = 1) = x_inv_bohr / ((BOHR_RADIUS_ANG * 1.0e-8)^p) +_inv_bohr2fm(x_inv_bohr, p = 1) = x_inv_bohr / ((BOHR_RADIUS_ANG * 1.0e5)^p) +_inv_bohr2eV(x_inv_bohr, p = 1) = x_inv_bohr / ((BOHR_RADIUS_ANG / HBARC_eV_ANG)^p) +_inv_bohr2me(x_inv_bohr, p = 1) = x_inv_bohr / ((BOHR_RADIUS_ANG / HBARC_eV_ANG * ELECTRONMASS)^p) + +_ang2inv_eV(x_ang, p = 1) = x_ang / (HBARC_eV_ANG^p) +_inv_ang2eV(x_inv_ang, p = 1) = x_inv_ang * (HBARC_eV_ANG^p) + +_me2eV(x_me, p = 1) = x_me * (ELECTRONMASS)^p +_me2keV(x_me, p = 1) = x_me * (ELECTRONMASS * 1.0e-3)^p +_me2MeV(x_me, p = 1) = x_me * (ELECTRONMASS * 1.0e-6)^p +_me2hartree(x_me, p = 1) = x_me * (ELECTRONMASS / HARTREE)^p +_me2inv_ang(x_me, p = 1) = x_me * (ELECTRONMASS / HBARC_eV_ANG)^p +_me2inv_bohr(x_me, p = 1) = x_me * (ELECTRONMASS / HBARC_eV_ANG * BOHR_RADIUS_ANG)^p + +# this is inv(Ang) to eV! +function _convert_Ang_to_eV(q) + return q * HBARC_eV_ANG +end + +# this is eV to inv Ang! +function _convert_eV_to_Ang(q) + return q / HBARC_eV_ANG +end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..28ba447 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,9 @@ +@inline _internalize_density(ne::T) where {T <: Real} = ne +@inline function _internalize_density(ne::Quantity) + return _density2me_dimless(ne) +end + +@inline _internalize_temperature(temp::T) where {T <: Real} = temp +@inline function _internalize_temperature(temp::Quantity) + return _energy2me_dimless(temp) +end diff --git a/test/Project.toml b/test/Project.toml index 0c36332..505c026 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,2 +1,6 @@ [deps] +QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" diff --git a/test/electron_system/checks.jl b/test/electron_system/checks.jl new file mode 100644 index 0000000..d2ff386 --- /dev/null +++ b/test/electron_system/checks.jl @@ -0,0 +1,111 @@ +# checks for dynamic structure factors + +function _dsf_detailed_balance_check(sys, om, q) + + input_pos_om = (om, q) + input_neg_om = (-om, q) + + @test isapprox( + dynamic_structure_factor(sys, input_neg_om), + exp(-om * beta(sys)) * + dynamic_structure_factor(sys, input_pos_om), + ) + return nothing +end + +function _dsf_sanity_check(sys, om, q) + fac = pi * electron_density(sys) + + im_rf = imag_dynamic_response(sys, (om, q)) + T = temperature(sys) + groundtruth_dsf = + iszero(T) ? -im_rf / fac : + -inv(fac * (one(om) - exp(-om / T))) * im_rf + + dsf = dynamic_structure_factor(sys, (om, q)) + @test isapprox(dsf, groundtruth_dsf) + + if dsf < 0.0 + @show dsf + @show im_rf + end + @test dsf >= zero(dsf) + + return nothing +end + +# checks for dynamic response functions + +function _f_sum_rule_check(sys, q; bounded = true) + T = temperature(sys) + KF = fermi_wave_vector(sys) + EF = fermi_energy(sys) + N0 = KF / (2 * pi^2) + + qb = q / KF + lower_om_bound = iszero(T) && bounded ? EF * max(zero(qb), qb^2 - 2 * qb) : zero(T) + upper_om_bound = iszero(T) && bounded ? EF * (qb^2 + 2 * qb) : EF * (qb^2 + 2 * qb) * 300.0 + tmp, err = quadgk( + x -> x * imag_dynamic_response(sys, (x, q)), + lower_om_bound, + upper_om_bound, + ) + _first_moment = -2 * tmp / pi + @test isapprox( + _first_moment, + electron_density(sys) * q^2, + #4 * N0 * q^2 * EF / 3, + rtol = 1.0e-2, + ) + return nothing +end + +function _rf_symmetry_check(sys, om, q, rtol) + input_pos_q = (om, q) + input_pos_om = (om, q) + input_neg_om = (-om, q) + input_neg_both = (-om, -q) + + #= + # TODO: check if this is valid! + + @test isapprox( + dynamic_response(sys, input_pos_q), + dynamic_response(sys, input_neg_both), + rtol = rtol, + ) + =# + @test isapprox( + real_dynamic_response(sys, input_pos_om), + real_dynamic_response(sys, input_neg_om), + rtol = rtol, + ) + @test isapprox( + imag_dynamic_response(sys, input_pos_om), + -imag_dynamic_response(sys, input_neg_om), + rtol = rtol, + ) + + return nothing +end + +function _rf_stability_check(sys, q) + @test real_dynamic_response(sys, (0.0, q)) <= zero(q) + @test imag_dynamic_response(sys, (0.0, q)) <= zero(q) + return nothing +end + +function _rf_property_check(sys, ne_internal, T_internal) + KF_internal = cbrt(3 * pi^2 * ne_internal) + EF_internal = KF_internal^2 / 2 + BETA_internal = inv(T_internal) + BETABAR_internal = EF_internal * BETA_internal + + @test isapprox(temperature(sys), T_internal) + @test isapprox(electron_density(sys), ne_internal) + @test isapprox(beta(sys), BETA_internal) + @test isapprox(betabar(sys), BETABAR_internal) + @test isapprox(fermi_wave_vector(sys), KF_internal) + @test isapprox(fermi_energy(sys), EF_internal) + return nothing +end diff --git a/test/electron_system/ideal.jl b/test/electron_system/ideal.jl new file mode 100644 index 0000000..47670b0 --- /dev/null +++ b/test/electron_system/ideal.jl @@ -0,0 +1,188 @@ +using Test +using Random +using QuadGK +using Unitful + +using MatterModels + +include("checks.jl") + +RNG = Xoshiro(137) +ATOL = sqrt(eps()) +RTOL = 1.0e-4 + +_transform12(x01) = x01 + one(x01) + +NES_ccm = 1u"cm^(-3)" .* ( + _transform12(rand(RNG)) * 1.0e20, + _transform12(rand(RNG)) * 1.0e21, + _transform12(rand(RNG)) * 1.0e22, + _transform12(rand(RNG)) * 1.0e23, + _transform12(rand(RNG)) * 1.0e24, +) +TEMPS_eV = 1u"eV" .* ( + rand(RNG), + rand(RNG) * 10, + rand(RNG) * 100, + rand(RNG) * 1000, +) + +APPROXS = (NoApprox(), NonDegenerated(), Degenerated()) + + +@testset "ne = $ne_ccm" for ne_ccm in NES_ccm + ne_internal = MatterModels._internalize_density(ne_ccm) + + KF = MatterModels._fermi_wave_vector(ne_internal) + EF = MatterModels._fermi_energy_from_kF(KF) + N0 = KF / (2 * pi^2) + + OMS = EF .* (0.0, rand(RNG), 1 + rand(RNG), 2 + rand(RNG), 3 + rand(RNG)) + QS = KF .* (1.0e-2 * rand(RNG), rand(RNG), 2 + rand(RNG), 3 + rand(RNG)) + + ### Zero temperature + + @testset "T = 0.0 eV" begin + test_system = IdealElectronSystem{ZeroTemperature}(ne_ccm) + test_system_small_finT = IdealElectronSystem(ne_ccm, 1.0e-2 * eps()) + + @testset "properties" begin + _rf_property_check(test_system, ne_internal, zero(ne_internal)) + @test iszero(temperature(test_system)) + @test isinf(beta(test_system)) + @test isinf(betabar(test_system)) + end + + @testset "q: $q" for q in QS + + @testset "rf stability" begin + _rf_stability_check(test_system, q) + end + @testset "rf f-sum rule" begin + _f_sum_rule_check(test_system, q) + end + + @testset "om: $om" for om in OMS + + @testset "rf symmetry" begin + _rf_symmetry_check(test_system, om, q, RTOL) + end + + @testset "rf sanity check" begin + groundtruth_imag_rf = + N0 * MatterModels._imag_lindhard_zero_temperature(NoApprox(), om / EF, q / KF) + @test isapprox( + groundtruth_imag_rf, + imag_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + groundtruth_real_rf = + N0 * MatterModels._real_lindhard_zero_temperature(NoApprox(), om / EF, q / KF) + @test isapprox( + groundtruth_real_rf, + real_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + groundtruth_rf = groundtruth_real_rf + 1im * groundtruth_imag_rf + @test isapprox( + groundtruth_rf, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + end + + @testset "rf matching" begin + @test isapprox( + imag_dynamic_response(test_system, (om, q)), + imag_dynamic_response(test_system_small_finT, (om, q)), + rtol = 1.0e-2, + atol = 1.0e-6, + ) + end + + @testset "dsf sanity check" begin + if !iszero(om) + _dsf_sanity_check(test_system, om, q) + end + end + end + end + end + + ### Finite temperature + + @testset "T = $T_eV" for T_eV in TEMPS_eV + T_internal = MatterModels._internalize_temperature.(T_eV) + @testset "approx = $approx" for approx in APPROXS + + test_system = IdealElectronSystem(ne_ccm, T_eV, approx) + + @testset "properties" begin + _rf_property_check(test_system, ne_internal, T_internal) + @test response_approximation(test_system) == approx + end + + @testset "q: $q" for q in QS + + @testset "rf stability" begin + _rf_stability_check(test_system, q) + end + + if approx == NoApprox() + @testset "rf f-sum rule" begin + _f_sum_rule_check(test_system, q) + end + end + + @testset "om: $om" for om in OMS + + @testset "rf symmetry" begin + _rf_symmetry_check(test_system, om, q, RTOL) + end + + @testset "rf sanity check" begin + groundtruth_imag_rf = N0 * MatterModels._imag_lindhard_nonzero_temperature( + approx, + om / EF, + q / KF, + betabar(test_system), + ) + @test isapprox( + groundtruth_imag_rf, + imag_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + groundtruth_real_rf = N0 * MatterModels._real_lindhard_nonzero_temperature( + approx, + om / EF, + q / KF, + betabar(test_system), + ) + @test isapprox( + groundtruth_real_rf, + real_dynamic_response(test_system, (om, q)), + ) + + groundtruth_rf = groundtruth_real_rf + 1im * groundtruth_imag_rf + @test isapprox(groundtruth_rf, dynamic_response(test_system, (om, q))) + end + + @testset "dsf sanity check" begin + if !iszero(om) + _dsf_sanity_check(test_system, om, q) + end + end + + @testset "detailed balance" begin + if !iszero(om) + _dsf_detailed_balance_check(test_system, om, q) + end + end + end # om + end # q + end # approx + end # finite T + +end diff --git a/test/electron_system/interacting.jl b/test/electron_system/interacting.jl new file mode 100644 index 0000000..15228b1 --- /dev/null +++ b/test/electron_system/interacting.jl @@ -0,0 +1,228 @@ +using Test +using Random +using QuadGK +using Unitful + +using MatterModels + +include("checks.jl") + +RNG = Xoshiro(137) +ATOL = sqrt(eps()) +RTOL = 1.0e-4 + +_transform12(x01) = x01 + one(x01) + +NES_ccm = 1u"cm^(-3)" .* ( + _transform12(rand(RNG)) * 1.0e20, + _transform12(rand(RNG)) * 1.0e21, + _transform12(rand(RNG)) * 1.0e22, + _transform12(rand(RNG)) * 1.0e23, + _transform12(rand(RNG)) * 1.0e24, +) +TEMPS_eV = 1u"eV" .* ( + rand(RNG), + rand(RNG) * 10, + rand(RNG) * 100, + rand(RNG) * 1000, +) + +APPROXS = (NoApprox(), NonDegenerated(), Degenerated()) +SCREENINGS = (NoScreening(), Screening()) + + +@testset "ne = $ne_ccm" for ne_ccm in NES_ccm + ne_internal = MatterModels._internalize_density(ne_ccm) + + KF = MatterModels._fermi_wave_vector(ne_internal) + EF = MatterModels._fermi_energy_from_kF(KF) + N0 = KF / (2 * pi^2) + + OMS = EF .* (0.0, rand(RNG), 1 + rand(RNG), 2 + rand(RNG), 3 + rand(RNG)) + QS = KF .* (1.0e-2 * rand(RNG), rand(RNG), 2 + rand(RNG), 3 + rand(RNG)) + + ### Zero temperature + + @testset "scr = $scr" for scr in SCREENINGS + @testset "T = 0.0 eV" begin + test_proper_system = IdealElectronSystem{ZeroTemperature}(ne_ccm) + test_system = InteractingElectronSystem(test_proper_system, scr) + + @testset "properties" begin + _rf_property_check(test_system, ne_internal, zero(ne_internal)) + @test screening(test_system) == scr + @test iszero(temperature(test_system)) + @test isinf(beta(test_system)) + @test isinf(betabar(test_system)) + end + + @testset "q: $q" for q in QS + + @testset "rf stability" begin + _rf_stability_check(test_system, q) + end + + # TODO: find out, why f-sum rule does not work for screened rf + if scr == NoScreening() + @testset "rf f-sum rule" begin + _f_sum_rule_check(test_system, q) + end + end + + @testset "om: $om" for om in OMS + + @testset "rf symmetry" begin + _rf_symmetry_check(test_system, om, q, RTOL) + end + + @testset "rf sanity check" begin + groundtruth_imag_rf = imag(dynamic_response(test_system, (om, q))) + @test isapprox( + groundtruth_imag_rf, + imag_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + groundtruth_real_rf = real(dynamic_response(test_system, (om, q))) + + @test isapprox( + groundtruth_real_rf, + real_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + groundtruth_rf = + real_dynamic_response(test_system, (om, q)) + + 1im * imag_dynamic_response(test_system, (om, q)) + @test isapprox( + groundtruth_rf, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + proper_rf = dynamic_response(test_proper_system, (om, q)) + lep = local_effective_potential(test_system, (om, q)) + groundtruth_rf_from_proper = proper_rf / (1 - lep * proper_rf) + @test isapprox( + groundtruth_rf_from_proper, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + diel_func = dielectric_function(test_system, (om, q)) + groundtruth_rf_from_dielec_func = proper_rf / diel_func + @test isapprox( + groundtruth_rf_from_dielec_func, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + end + + @testset "dsf sanity check" begin + if !iszero(om) + _dsf_sanity_check(test_system, om, q) + end + end + end + end + end + + ### Finite temperature + + @testset "T = $T_eV" for T_eV in TEMPS_eV + T_internal = MatterModels._internalize_temperature.(T_eV) + @testset "approx = $approx" for approx in APPROXS + + test_proper_system = IdealElectronSystem(ne_ccm, T_eV, approx) + test_system = InteractingElectronSystem(test_proper_system, scr) + + @testset "properties" begin + _rf_property_check(test_system, ne_internal, T_internal) + #@test response_approximation(test_system) == approx + #@test screening(test_system) == scr + end + + @testset "q: $q" for q in QS + + @testset "rf stability" begin + _rf_stability_check(test_system, q) + end + + #= + if approx == NoApprox() + @testset "rf f-sum rule" begin + _f_sum_rule_check(test_system, q) + end + end + =# + + @testset "om: $om" for om in OMS + + @testset "rf symmetry" begin + _rf_symmetry_check(test_system, om, q, RTOL) + end + + @testset "rf sanity check" begin + groundtruth_imag_rf = imag(dynamic_response(test_system, (om, q))) + @test isapprox( + groundtruth_imag_rf, + imag_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + groundtruth_real_rf = real(dynamic_response(test_system, (om, q))) + + @test isapprox( + groundtruth_real_rf, + real_dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + groundtruth_rf = + real_dynamic_response(test_system, (om, q)) + + 1im * imag_dynamic_response(test_system, (om, q)) + @test isapprox( + groundtruth_rf, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + proper_rf = dynamic_response(test_proper_system, (om, q)) + lep = local_effective_potential(test_system, (om, q)) + groundtruth_rf_from_proper = proper_rf / (1 - lep * proper_rf) + @test isapprox( + groundtruth_rf_from_proper, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + + diel_func = dielectric_function(test_system, (om, q)) + groundtruth_rf_from_dielec_func = proper_rf / diel_func + @test isapprox( + groundtruth_rf_from_dielec_func, + dynamic_response(test_system, (om, q)), + rtol = RTOL, + ) + end + + @testset "dsf sanity check" begin + + # TODO: find out, why this is broken sometimes + #= + if !iszero(om) + _dsf_sanity_check(test_system, om, q) + end + =# + end + + @testset "detailed balance" begin + if !iszero(om) + _dsf_detailed_balance_check(test_system, om, q) + end + end + end # om + end # q + end # approx + end # finite T + end # screening +end diff --git a/test/electron_system/screening.jl b/test/electron_system/screening.jl new file mode 100644 index 0000000..1dd5583 --- /dev/null +++ b/test/electron_system/screening.jl @@ -0,0 +1,91 @@ +using Test +using Random +using QuadGK +using Unitful + +using MatterModels + +include("checks.jl") + +RNG = Xoshiro(137) +ATOL = sqrt(eps()) +RTOL = 1.0e-4 + +_transform12(x01) = x01 + one(x01) + +NES_ccm = 1u"cm^(-3)" .* ( + _transform12(rand(RNG)) * 1.0e20, + _transform12(rand(RNG)) * 1.0e21, + _transform12(rand(RNG)) * 1.0e22, + _transform12(rand(RNG)) * 1.0e23, + _transform12(rand(RNG)) * 1.0e24, +) +TEMPS_eV = 1u"eV" .* [ + rand(RNG), + rand(RNG) * 10, + rand(RNG) * 100, + rand(RNG) * 1000, +] + +APPROXS = [NoApprox(), NonDegenerated(), Degenerated()] + + +@testset "ne = $ne_ccm" for ne_ccm in NES_ccm + ne_internal = MatterModels._internalize_density(ne_ccm) + + KF = MatterModels._fermi_wave_vector(ne_internal) + EF = MatterModels._fermi_energy_from_kF(KF) + N0 = KF / (2 * pi^2) + + OMS = EF .* (0.0, rand(RNG), 1 + rand(RNG), 2 + rand(RNG), 3 + rand(RNG)) + QS = KF .* (1.0e-2 * rand(RNG), rand(RNG), 2 + rand(RNG), 3 + rand(RNG)) + + TEST_SYSTEM_zeroT = IdealElectronSystem{ZeroTemperature}(ne_ccm) + TEST_SYSTEMS_finT = [IdealElectronSystem(ne_ccm, T_eV, approx) for T_eV in TEMPS_eV for approx in APPROXS] + TEST_SYSTEMS = (TEST_SYSTEM_zeroT, TEST_SYSTEMS_finT...) + + @testset "T = $(temperature(test_system))" for test_system in TEST_SYSTEMS + @testset "no screening" begin + scr = NoScreening() + @testset "om: $om, q: $q" for (om, q) in Iterators.product(OMS, QS) + + @test isapprox(dielectric_function(test_system, scr, (om, q)), one(om)) + + @test isapprox(local_effective_potential(scr, (om, q)), zero(om)) + + @test isapprox(pseudo_potential(scr, (om, q)), zero(om)) + + @test isapprox(local_field_correction(scr, (om, q)), zero(om)) + end + end + + @testset "screening" begin + # checking only default value: coulomb+no_lfc + scr = Screening() + + @testset "default constructor" begin + @test scr == Screening(CoulombPseudoPotential(), NoLocalFieldCorrection()) + end + + @testset "om: $om, q: $q" for (om, q) in Iterators.product(OMS, QS) + + @test isapprox( + dielectric_function(test_system, scr, (om, q)), + one(om) - + pseudo_potential(scr, (om, q)) * dynamic_response(test_system, (om, q)), + ) + + @test isapprox( + local_effective_potential(scr, (om, q)), + pseudo_potential(scr, (om, q)), + ) + + @test isapprox(pseudo_potential(scr, (om, q)), ELEMENTARY_CHARGE_SQUARED / q^2) + + @test isapprox(local_field_correction(scr, (om, q)), zero(om)) + end + end + + end + +end diff --git a/test/test-electron-system.jl b/test/test-electron-system.jl new file mode 100644 index 0000000..3dcf272 --- /dev/null +++ b/test/test-electron-system.jl @@ -0,0 +1,14 @@ +using Test +using SafeTestsets + +begin + @safetestset "Ideal system" begin + include("electron_system/ideal.jl") + end + @safetestset "Screening" begin + include("electron_system/screening.jl") + end + @safetestset "Interacting system" begin + include("electron_system/interacting.jl") + end +end diff --git a/test/test-lookup.jl b/test/test-lookup.jl new file mode 100644 index 0000000..1bb06a2 --- /dev/null +++ b/test/test-lookup.jl @@ -0,0 +1,26 @@ +using MatterModels +using Random + +RNG = Xoshiro(137) + +gaussian(x, y) = exp(-x^2 / 2 - y^2) +a_arr = collect(range(-rand(RNG), rand(RNG); length = 4)) # cols +b_arr = collect(range(-rand(RNG), rand(RNG); length = 5)) # rows +vals = gaussian.(a_arr', b_arr) +lookup_methods = [ + MatterModels.InterpolExtrapol(), + MatterModels.InterpolEndValue(), + MatterModels.NearestInput(), + MatterModels.BelowInput(), + MatterModels.AboveInput(), +] +GLT = MatterModels.GridLookupTable(a_arr, b_arr, vals) +@testset "$method" for method in lookup_methods + for a in a_arr + for b in b_arr + test_val = MatterModels.lookup(method, a, b, GLT) + groundtruth = gaussian(a, b) + @test isapprox(groundtruth, test_val) + end + end +end