Add a plasma implementation

src/ElectronicStructureModels.jl

@@ -5,7 +5,21 @@ hello_world() = "Hello, World!"
 # matter model
 export AbstractMatterModel
 
-# Electron system
+### Plasma
+export AbstactPlasma, AbstractChiharaPlasma
+export electron_system, electron_temperature, electron_density
+export ion_system, ion_temperature, ion_density
+
+# OneComponentPlasma
+export OneComponentPlasma
+
+# Plasma Components
+
+export AbstractPlasmaComponent
+export temperature, number_density
+
+### Electron system
+
 export AbstractElectronSystem
 export temperature, electron_density, imag_dynamic_response, real_dynamic_response
 export fermi_wave_vector,
@@ -28,7 +42,8 @@ export AbstractResponseApproximation, NoApprox, NonDegenerated, Degenerated, Zer
 export response_approximation
 export InteractingElectronSystem
 
-# Ion System
+
+### Ion System
 
 export AbstractIonSystem
 export atomic_number, ionization_state, number_bound_electrons
@@ -95,5 +110,8 @@ include("matter/ions/impl.jl")
 include("matter/ions/subshells.jl")
 include("matter/ions/form_factors/PaulingSherman.jl")
 
+include("matter/plasma/interface.jl")
+include("matter/plasma/generic.jl")
+include("matter/plasma/impl.jl")
 
 end

src/matter/ions/impl.jl

@@ -1,7 +1,24 @@
-struct IonSystem{I <: Integer, C <: Integer} <: AbstractIonSystem
+struct IonSystem{T <: Real, I <: Integer, C <: Integer} <: AbstractIonSystem
     ZA::I # atomic number
     ZF::C # ionization state
-    function IonSystem(ZA::A, ZF::F) where {A <: Int, F <: Int}
+    density::T
+    temperature::T
+
+    function IonSystem(
+            ZA::A,
+            ZF::F,
+            density::T1,
+            temperature::T2,
+        ) where {
+            T <: Real,
+            T1 <: Union{T, Quantity{T}},
+            T2 <: Union{T, Quantity{T}},
+            A <: Int,
+            F <: Int,
+        }
+
+        n_internal = _internalize_density(density)
+        temp_internal = _internalize_temperature(temperature)
         ZA > 0 || throw(
             ArgumentError(
                 "atomic number must be striclty positive"
@@ -14,15 +31,18 @@ struct IonSystem{I <: Integer, C <: Integer} <: AbstractIonSystem
             )
         )
 
-        return new{A, F}(ZA, ZF)
+        return new{T, A, F}(ZA, ZF, n_internal, temp_internal)
     end
 end
 
 # interface with PeriodicTables.jl
-IonSystem(e::Element, ZF) = IonSystem(e.number, ZF)
+IonSystem(e::Element, ZF, temp, density) = IonSystem(e.number, ZF, temp, density)
 element(ions::IonSystem) = getindex(elements, ions.ZA)
 
-# ESM interface
+# Plasma Component Interface
 atomic_number(ions::IonSystem) = ions.ZA
 ionization_state(ions::IonSystem) = ions.ZF
 number_bound_electrons(ions::IonSystem) = ions.ZA - ions.ZF
+
+_temperature(ions::IonSystem) = ions.temperature
+_number_density(ions::IonSystem) = ions.density

src/matter/ions/interface.jl

@@ -6,6 +6,7 @@ Abstract base type representing an ion component of a Chihara plasma.
 """
 abstract type AbstractIonSystem <: AbstractPlasmaComponent end
 
+
 """
 
     atomic_number(ions::AbstractIonSystem)

src/matter/plasma/Rayleig_weight/screening_cloud/Debye_Hueckel_screeing.jl

@@ -0,0 +1,10 @@
+struct DebyeHueckelScreening <: AbstractScreeningCloud end
+
+_Debye_Hueckel_screening_length(n_e, temp_e) = sqrt(n_e / temp_e) * ELECTRON_CHARGE
+
+function screening_cloud(plasma::AbstractChiharaPlasma, ::DebyeHueckelScreening, q::Real)
+    kappa_e = _Debye_Hueckel_screening_length(electron_density(plasma), electron_temperature(plasma))
+    ZF = ionization_state(plasma)
+
+    return ZF * kappa_e^2 / (q^2 + kappa_e^2)
+end

src/matter/plasma/Rayleig_weight/screening_cloud/interface.jl

@@ -0,0 +1,8 @@
+abstract type AbstractScreeningCloud end
+
+"""
+
+    screening_cloud(::AbstractPlasma, ::AbstractScreeningCloud, q::Real)
+
+"""
+function screening_cloud end

src/matter/plasma/generic.jl

@@ -0,0 +1,6 @@
+### Delegations from IonSystem
+
+atomic_number(plasma::AbstractChiharaPlasma) = atomic_number(ion_system(plasma))
+ionization_state(plasma::AbstractChiharaPlasma) = ionization_state(ion_system(plasma))
+number_bound_electrons(plasma::AbstractChiharaPlasma) = number_bound_electrons(ion_system(plasma))
+element(plasma::AbstractChiharaPlasma) = element(ion_system(plasma))

src/matter/plasma/impl.jl

@@ -0,0 +1,11 @@
+struct OneComponentPlasma{E, I} <: AbstractChiharaPlasma
+    elec_system::E
+    ion_system::I
+end
+
+@inline electron_system(plasma::OneComponentPlasma) = plasma.elec_system
+@inline electron_temperature(plasma::OneComponentPlasma) = temperature(plasma.elec_system)
+@inline electron_density(plasma::OneComponentPlasma) = number_density(plasma.elec_system)
+@inline ion_system(plasma::OneComponentPlasma) = plasma.ion_system
+@inline ion_temperature(plasma::OneComponentPlasma) = temperature(plasma.ion_system)
+@inline ion_density(plasma::OneComponentPlasma) = number_density(plasma.ion_system)

test/ion_system/form_factor/PaulingSherman.jl

@@ -2,19 +2,36 @@ using Test
 using ElectronicStructureModels
 using XDaveWrapper
 
+using Random
+using Unitful
+
+RNG = Xoshiro(161)
+
+_transform12(x01) = x01 + one(x01)
+
 _groundtruth_ff(ions::IonSystem, q) = XDaveWrapper.groundtruth_form_factor(
     atomic_number(ions), ionization_state(ions), q
 )
 
 const QS = [0.2]
 
-@testset "ZA = $ZA" for ZA in 1:29
-    @testset "ZF = $ZF" for ZF in 1:ZA
-        ions = IonSystem(ZA, ZF)
+const NIS_ccm = 1u"cm^(-3)" .* (
+    _transform12(rand(RNG)) * 1.0e20,
+)
+const TEMPS_eV = 1u"eV" .* (
+    rand(RNG),
+)
+@testset "T = $T" for T in TEMPS_eV
+    @testset "N = $N" for N in NIS_ccm
+        @testset "ZA = $ZA" for ZA in 1:29
+            @testset "ZF = $ZF" for ZF in 1:ZA
+                ions = IonSystem(ZA, ZF, N, T)
 
-        @testset "q = $q" for q in QS
-            groundtruth = _groundtruth_ff(ions, q)
-            @test isapprox(groundtruth, form_factor(ions, PaulingSherman(), q))
+                @testset "q = $q" for q in QS
+                    groundtruth = _groundtruth_ff(ions, q)
+                    @test isapprox(groundtruth, form_factor(ions, PaulingSherman(), q))
+                end
+            end
         end
     end
 end

test/ion_system/impl.jl

@@ -1,24 +1,53 @@
+# TODO:
+# - add checks if the quantities T and N are actually temperatues and densities
+
+
 using Test
 using ElectronicStructureModels
 using PeriodicTable
+using Random
+using Unitful
 
-@testset "ZA = $ZA" for ZA in 1:100
-    @testset "ZF = $ZF" for ZF in 1:ZA
+RNG = Xoshiro(161)
 
-        @testset "Properties" begin
-            ions = IonSystem(ZA, ZF)
-            @test atomic_number(ions) == ZA
-            @test ionization_state(ions) == ZF
-            @test number_bound_electrons(ions) == ZA - ZF
-        end
+_transform12(x01) = x01 + one(x01)
+
+NIS_ccm = 1u"cm^(-3)" .* (
+    _transform12(rand(RNG)) * 1.0e20,
+)
+TEMPS_eV = 1u"eV" .* (
+    rand(RNG),
+)
+@testset "T = $T" for T in TEMPS_eV
+    @testset "N = $N" for N in NIS_ccm
+        n_internal = ElectronicStructureModels._internalize_density(N)
+        T_internal = ElectronicStructureModels._internalize_temperature(T)
+
+        @testset "ZA = $ZA" for ZA in 1:100
+            @testset "ZF = $ZF" for ZF in 1:ZA
+
+                @testset "Properties" begin
+                    ions = IonSystem(ZA, ZF, N, T)
+                    @test atomic_number(ions) == ZA
+                    @test ionization_state(ions) == ZF
+                    @test number_bound_electrons(ions) == ZA - ZF
+
+                    @test isapprox(temperature(ions), T_internal)
+                    @test isapprox(number_density(ions), n_internal)
+                end
+
+                @testset "PeriodicTable integration" begin
+                    ions = IonSystem(elements[ZA], ZF, N, T)
 
-        @testset "PeriodicTable integration" begin
-            ions = IonSystem(elements[ZA], ZF)
+                    @test atomic_number(ions) == ZA
+                    @test ionization_state(ions) == ZF
+                    @test number_bound_electrons(ions) == ZA - ZF
+                    @test element(ions) == elements[ZA]
 
-            @test atomic_number(ions) == ZA
-            @test ionization_state(ions) == ZF
-            @test number_bound_electrons(ions) == ZA - ZF
-            @test element(ions) == elements[ZA]
+                    @test isapprox(temperature(ions), T_internal)
+                    @test isapprox(number_density(ions), n_internal)
+                end
+            end
         end
     end
 end

test/plasma/impl.jl

@@ -0,0 +1,64 @@
+using Test
+using ElectronicStructureModels
+using PeriodicTable
+using Random
+using Unitful
+
+RNG = Xoshiro(161)
+
+_transform12(x01) = x01 + one(x01)
+
+N_IONS_ccm = 1u"cm^(-3)" .* (
+    _transform12(rand(RNG)) * 1.0e20,
+)
+TEMPS_IONS_eV = 1u"eV" .* (
+    rand(RNG),
+)
+
+N_ELECS_ccm = 1u"cm^(-3)" .* (
+    _transform12(rand(RNG)) * 1.0e20,
+)
+TEMPS_ELECS_eV = 1u"eV" .* (
+    rand(RNG),
+)
+
+ZA = 12
+ZF = 2
+
+@testset "Ni= $Ni" for Ni in N_IONS_ccm
+    @testset "Ti = $Ti" for Ti in TEMPS_IONS_eV
+        @testset "Ne= $Ne" for Ne in N_ELECS_ccm
+            @testset "Te = $Te" for Te in TEMPS_IONS_eV
+                Ne_internal = ElectronicStructureModels._internalize_density(Ne)
+                Te_internal = ElectronicStructureModels._internalize_temperature(Te)
+                Ni_internal = ElectronicStructureModels._internalize_density(Ni)
+                Ti_internal = ElectronicStructureModels._internalize_temperature(Ti)
+
+                electrons = IdealElectronSystem(Ne, Te)
+                ions = IonSystem(ZA, ZF, Ni, Ti)
+                plasma = OneComponentPlasma(electrons, ions)
+
+                @testset "properties" begin
+
+                    @test electron_system(plasma) == electrons
+                    @test isapprox(electron_temperature(plasma), Te_internal)
+                    @test isapprox(electron_density(plasma), Ne_internal)
+
+                    @test ion_system(plasma) == ions
+                    @test isapprox(ion_temperature(plasma), Ti_internal)
+                    @test isapprox(ion_density(plasma), Ni_internal)
+
+                end
+
+                @testset "delegations" begin
+
+                    @test atomic_number(plasma) == ZA
+                    @test ionization_state(plasma) == ZF
+                    @test number_bound_electrons(plasma) == ZA - ZF
+                    @test element(plasma) == elements[ZA]
+
+                end
+            end
+        end
+    end
+end

test/runtests.jl

@@ -8,7 +8,7 @@ Don't add your tests to runtests.jl. Instead, create files named
 
 The file will be automatically included inside a `@testset` with title "Title For My Test".
 =#
-
+#=
 for (root, dirs, files) in walkdir(@__DIR__)
     for file in files
         if isnothing(match(r"^test-.*\.jl$", file))
@@ -20,3 +20,5 @@ for (root, dirs, files) in walkdir(@__DIR__)
         end
     end
 end
+=#
+include("test-plasma.jl")

test/test-plasma.jl

@@ -0,0 +1,9 @@
+using Test
+using SafeTestsets
+
+begin
+    @safetestset "OneComponentPlasma" begin
+        include("plasma/impl.jl")
+    end
+
+end