Refactor ThermalDiffusionLithiumDensity to remove hardcoded parameters (#578)

* Refactor ThermalDiffusionLithiumDensity to remove hardcoded parameters

* Refactor ThermalDiffusionLithiumDensity to remove hardcoded parameters

* Fix lithium diffusion model and YAML parsing

* Bug fixes

* Pass precision type to `ThermalDiffusionLithiumDensityParameters`

* Store `ThermalDiffusionLithiumDensityParameters.H` in units of `cal`

* Revert changes to `Project.toml`

* Fix order of magnitude of `lithium_density_on_contact`

The publication quotes values for calculating cm^-3, SSD expects m^-3

* cleaning the code and fixing bugs

* Remove example JSON file for inverted coaxial detectors

* Add references to `Thermal_Diffusion_Config.yaml`

* Reorganize code and add docstrings

* Remove `in_germanium` from `diffusivity` methods

The publications used for germanium also seem to provide values for silicon. I do not see why this model should not be also applicable to silicon.

* Update activation energy unit to `cal/mol`

---------

Co-authored-by: Simone Pellegrini <simonepellegrini@Mac.fritz.box>
Co-authored-by: Felix Hagemann <hagemann@berkeley.edu>
Co-authored-by: Simone Pellegrini <simonepellegrini@MacBook-Air-di-Simone-4.fritz.box>

Author: 4 people

examples/example_config_files/ThermalDiffusionLithiumDensity/Thermal_Diffusion_Config.yaml

@@ -0,0 +1,17 @@
+model: ThermalDiffusionLithiumDensity
+material: HPGe
+annealing_temperature_ranges:
+  # https://doi.org/10.1103/PhysRev.96.21
+  - T_min: 473K
+    T_max: 873K
+    D0: 25e-4cm^2/s
+    H: 11800cal/mol
+  # https://doi.org/10.1103/PhysRev.91.193
+  - T_min: 873K
+    T_max: 1273K
+    D0: 13e-4cm^2/s
+    H: 10700cal/mol
+  
+experimental_parameters:
+  a: 21.27
+  b: 2610

src/ImpurityDensities/ThermalDiffusionLithiumDensity.jl

@@ -1,56 +1,49 @@
 """
-struct ThermalDiffusionLithiumDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
+    struct ThermalDiffusionLithiumDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
 
 Lithium impurity density model. Ref: [Dai _et al._ (2023)](https://doi.org/10.1016/j.apradiso.2022.110638)
  
 ## Fields
 * `lithium_annealing_temperature::T`: lithium annealing temperature in Kelvin, when the lithium is diffused into the crystal. The default value is 623 K.
 * `lithium_annealing_time::T`: lithium annealing time in seconds. The default value is 18 minutes.
 * `distance_to_contact::Function`: the function for describing the depth to surface. 
-    1) use the built-in function `ConstructiveSolidGeometry.distance_to_surface` to calculate the distance to the surface of the doped contact 
+    1) use the built-in function `ConstructiveSolidGeometry.distance_to_surface` to calculate the distance to the surface of the doped contact.
     2) custom. Custom function might be much faster while the detector has good symmetry.
-* `lithium_density_on_contact::T`: the lithium concentration in the surface (in m^-3)
-* `lithium_diffusivity_in_germanium::T`: the diffusivity of lithium in germanium (in m^2*s^-1)
+* `lithium_density_on_contact::T`: the lithium concentration in the surface (in m^-3).
+* `lithium_diffusivity::T`: the diffusivity of lithium (in m^2/s).
 
 ## Parameters for constructing the model
 * `lithium_annealing_temperature::T`
 * `lithium_annealing_time::T`
 * `contact_with_lithium_doped::G`: the geometry of the contact with lithium doped, which is used to calculate the distance to the surface.
     - If you don't pass this parameter, the `distance_to_contact` function should be passed.
     - If you pass this parameter but don't pass the `distance_to_contact` function, the `distance_to_contact` function will use the default one.
-* `distance_to_contact::Function`: optional, default is `ConstructiveSolidGeometry.distance_to_surface(pt, contact_with_lithium_doped)`
-* `lithium_density_on_contact::T`: optional, default is the saturated lithium density at the given temperature
-* `lithium_diffusivity_in_germanium::T`: optional, default is calculated with the annealing temperature
+* `distance_to_contact::Function`: optional, default is `ConstructiveSolidGeometry.distance_to_surface(pt, contact_with_lithium_doped)`.
+* `lithium_density_on_contact::T`: optional, default is the saturated lithium density at the given temperature.
+* `lithium_diffusivity::T`: optional, default is calculated with the annealing temperature.
 """
 struct ThermalDiffusionLithiumDensity{T <: SSDFloat} <: AbstractImpurityDensity{T}
     lithium_annealing_temperature::T
     lithium_annealing_time::T
     inactive_contact_id::Int
     distance_to_contact::Function
     lithium_density_on_contact::T
-    lithium_diffusivity_in_germanium::T
+    lithium_diffusivity::T
 end
 
-function calculate_lithium_diffusivity_in_germanium(lithium_annealing_temperature::T)::T where {T <: SSDFloat}
-    # D0 [m^2*s^-1]
-    # H [cal]
-    D0::T, H::T = ifelse(lithium_annealing_temperature <= 873, T.((2.5e-7, 11800)), T.((1.3e-7, 10700)))
-    D0 * exp(-H/(R_gas*lithium_annealing_temperature))
-end
-function calculate_lithium_saturated_density(lithium_annealing_temperature::T)::T where {T <: SSDFloat}
-    exp10(27.27 - 2610.0/lithium_annealing_temperature)
-end
+include("ThermalDiffusionLithiumDensityParameters.jl")
 
-function ThermalDiffusionLithiumDensity{T}(
+function ThermalDiffusionLithiumDensity{T}( 
     lithium_annealing_temperature::T,
     lithium_annealing_time::T,
     contact_with_lithium_doped::G,
     inactive_contact_id::Int;
+    model_parameters::ThermalDiffusionLithiumDensityParameters{T} = ThermalDiffusionLithiumParameters(T=T), 
     distance_to_contact::Function = pt::AbstractCoordinatePoint{T} -> ConstructiveSolidGeometry.distance_to_surface(pt, contact_with_lithium_doped),
-    lithium_density_on_contact::T = calculate_lithium_saturated_density(lithium_annealing_temperature),
-    lithium_diffusivity_in_germanium::T = calculate_lithium_diffusivity_in_germanium(lithium_annealing_temperature),
+    lithium_density_on_contact::T = calculate_lithium_saturated_density(lithium_annealing_temperature, model_parameters.saturation),
+    lithium_diffusivity::T = calculate_lithium_diffusivity(lithium_annealing_temperature, model_parameters.diffusion),
 ) where {T <: SSDFloat, G <: Union{<:AbstractGeometry, Nothing}}
-    ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, inactive_contact_id, distance_to_contact, lithium_density_on_contact, lithium_diffusivity_in_germanium)
+    ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, inactive_contact_id, distance_to_contact, lithium_density_on_contact, lithium_diffusivity)
 end
 
 function ImpurityDensity(T::DataType, t::Val{:li_diffusion}, dict::AbstractDict, input_units::NamedTuple)
@@ -59,15 +52,16 @@ function ImpurityDensity(T::DataType, t::Val{:li_diffusion}, dict::AbstractDict,
     contact_with_lithium_doped = haskey(dict, "contact_with_lithium_doped") ? dict["contact_with_lithium_doped"] : nothing # you don't have to pass the geometry of doped contact only when the distance_to_contact is passed
     inactive_contact_id = get(dict, "doped_contact_id", -1)
     inactive_contact_id < 1 && error("Invalid doped_contact_id: missing or misspelled key")
-    ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, contact_with_lithium_doped, inactive_contact_id)
+    model_parameters = haskey(dict,"model_parameters") ? ThermalDiffusionLithiumParameters(dict["model_parameters"]; T) : ThermalDiffusionLithiumParameters(; T)
+    ThermalDiffusionLithiumDensity{T}(lithium_annealing_temperature, lithium_annealing_time, contact_with_lithium_doped, inactive_contact_id, model_parameters)
 end
 
 function get_impurity_density(li_diffusion::ThermalDiffusionLithiumDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
     pt::CartesianPoint{T} = CartesianPoint(pt)
     depth = li_diffusion.distance_to_contact(pt)
-    li_diffusion.lithium_density_on_contact * SpecialFunctions.erfc(depth/2/sqrt(li_diffusion.lithium_diffusivity_in_germanium*li_diffusion.lithium_annealing_time))
+    li_diffusion.lithium_density_on_contact * SpecialFunctions.erfc(depth/2/sqrt(li_diffusion.lithium_diffusivity*li_diffusion.lithium_annealing_time))
 end
 
-(*)(scale::Real, tidm::ThermalDiffusionLithiumDensity{T}) where {T} = ThermalDiffusionLithiumDensity{T}(tidm.lithium_annealing_temperature, tidm.lithium_annealing_time, tidm.inactive_contact_id, tidm.distance_to_contact, T(scale * tidm.lithium_density_on_contact), tidm.lithium_diffusivity_in_germanium)
+(*)(scale::Real, tidm::ThermalDiffusionLithiumDensity{T}) where {T} = ThermalDiffusionLithiumDensity{T}(tidm.lithium_annealing_temperature, tidm.lithium_annealing_time, tidm.inactive_contact_id, tidm.distance_to_contact, T(scale * tidm.lithium_density_on_contact), tidm.lithium_diffusivity)
 
-(+)(offset::Union{<:Real, <:Quantity{<:Real, Unitful.𝐋^(-3)}}, tidm::ThermalDiffusionLithiumDensity{T}) where {T} = ThermalDiffusionLithiumDensity{T}(tidm.lithium_annealing_temperature, tidm.lithium_annealing_time, tidm.inactive_contact_id, tidm.distance_to_contact, T(to_internal_units(offset) + tidm.lithium_density_on_contact), tidm.lithium_diffusivity_in_germanium)
+(+)(offset::Union{<:Real, <:Quantity{<:Real, Unitful.𝐋^(-3)}}, tidm::ThermalDiffusionLithiumDensity{T}) where {T} = ThermalDiffusionLithiumDensity{T}(tidm.lithium_annealing_temperature, tidm.lithium_annealing_time, tidm.inactive_contact_id, tidm.distance_to_contact, T(to_internal_units(offset) + tidm.lithium_density_on_contact), tidm.lithium_diffusivity)

src/ImpurityDensities/ThermalDiffusionLithiumDensityParameters.jl

@@ -0,0 +1,159 @@
+"""
+    struct LithiumDiffusionParameters{T <: SSDFloat}
+
+Set of parameters to calculate the lithium diffusivity `D` at a given annealing temperature `T_an` using
+```math
+D(T_\text{an}) = D_0 \exp(-H / RT+\text{an})
+```
+
+## Parametric types
+* `T`: Precision type.
+    
+## Fields
+* `T_min::T`: Minimum temperature (in K) for which the parameters are valid.
+* `T_max::T`: Maximum temperature (in K) for which the parameters are valid.
+* `D0::T`: Diffusivity constant in m^2/s.
+* `H::T`: Activation energy in cal/mol.
+"""
+struct LithiumDiffusionParameters{T <: SSDFloat}
+    T_min::T
+    T_max::T
+    D0::T
+    H::T
+end
+
+function calculate_lithium_diffusivity(lithium_annealing_temperature::T, parameters::NTuple{N, LithiumDiffusionParameters{T}})::T where {T <: SSDFloat, N}
+
+    if !(parameters[1].T_min ≤ lithium_annealing_temperature ≤ parameters[end].T_max)
+        throw(ArgumentError("Invalid lithium_annealing_temperature=$(lithium_annealing_temperature): expected $(parameters[1].T_min) ≤ lithium_annealing_temperature ≤ $(parameters[end].T_max)."))
+    end
+
+    for parameter in parameters
+        if lithium_annealing_temperature <= parameter.T_max
+            return parameter.D0 * exp(-parameter.H/(R_gas*lithium_annealing_temperature))
+        end
+    end
+    
+    # throw an error if thediffusivity couldn't be calculated (this should never be the case)
+    throw(ArgumentError("Could not determine lithium diffusivity"))
+end
+
+
+"""
+    struct LithiumSaturationParameters{T <: SSDFloat}
+
+Set of experimental fit parameters to calculate the lithium saturated density `N_s` in 1/cm³ at a given annealing temperature `T_an` using
+```math
+N_s(T_\text{an}) = 10^{a - b/T_\text{an})
+```
+
+## Parametric types
+* `T`: Precision type.
+    
+## Fields
+* `a::T`: Fit parameter.
+* `b::T`: Fit parameter.`
+"""
+struct LithiumSaturationParameters{T <: SSDFloat}
+    a::T
+    b::T
+end
+
+function calculate_lithium_saturated_density(lithium_annealing_temperature::T, parameters::LithiumSaturationParameters{T})::T where {T <: SSDFloat}
+    ustrip(internal_length_unit^-3, exp10(parameters.a - parameters.b/lithium_annealing_temperature) * u"cm^-3")
+end
+
+
+"""
+    ThermalDiffusionLithiumDensityParameters{T <: SSDFloat,N}
+
+Set of parameters to calculate the lithium diffusivity and the lithium saturated density.
+
+## Parametric types
+* `T`: Precision type.
+* `N`: Number of temperature ranges to calculate the lithium diffusivity.
+    
+## Fields
+* `diffusion::NTuple{N, LithiumDiffusionParameters{T}}`: Lithium diffusivity parameters.
+* `saturation::LithiumSaturationParameters{T}`: Lithium saturated density parameters.`
+
+See also [`LithiumDiffusionParameters`](@ref) and [`LithiumSaturationParameters`](@ref).
+"""
+struct ThermalDiffusionLithiumDensityParameters{T <: SSDFloat,N}
+    diffusion::NTuple{N, LithiumDiffusionParameters{T}}
+    saturation::LithiumSaturationParameters{T}
+end
+
+
+function _ThermalDiffusionLithiumParameters(
+        config::AbstractDict;
+        T::Type = Float32,
+        temperature_ranges::AbstractVector{<:AbstractDict} = config["annealing_temperature_ranges"] ,
+        a::Union{Real, String} = config["experimental_parameters"]["a"],
+        b::Union{Real, String} = config["experimental_parameters"]["b"],
+        input_units::Union{Missing, NamedTuple} = missing
+    )
+
+    temperature_unit = !ismissing(input_units) ? input_units.temperature : internal_temperature_unit
+    diffusivity_unit = !ismissing(input_units) ? input_units.length^2 / internal_time_unit : internal_length_unit^2/internal_time_unit
+
+    diffusion_params = Tuple(
+        LithiumDiffusionParameters{T}(
+            _parse_value(T, r["T_min"], temperature_unit),
+            _parse_value(T, r["T_max"], temperature_unit),
+            _parse_value(T, r["D0"], diffusivity_unit),
+            _parse_value(T, r["H"], internal_activation_energy_unit)
+        )
+        for r in temperature_ranges
+    )
+
+    for i in eachindex(diffusion_params)
+        if diffusion_params[i].T_max < diffusion_params[i].T_min
+            throw(ConfigFileError("Invalid annealing temperature range $i: T_max must be ≥ T_min."))
+        end
+        if i > 1 && diffusion_params[i-1].T_max != diffusion_params[i].T_min
+            throw(ConfigFileError("Annealing temperature ranges must be contiguous and increasing: range $(i-1) ends at a different temperature than range $i starts."))
+        end
+    end
+
+    saturation_params = LithiumSaturationParameters{T}(_parse_value(T, a, Unitful.NoUnits), _parse_value(T, b, Unitful.NoUnits))
+    
+    ThermalDiffusionLithiumDensityParameters(diffusion_params, saturation_params)
+end
+
+
+function ThermalDiffusionLithiumParameters(config::AbstractDict, input_units::Union{Missing, NamedTuple} = missing; kwargs...)
+    if !haskey(config, "annealing_temperature_ranges") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'annealing_temperature_ranges'.")) 
+    elseif isempty(config["annealing_temperature_ranges"])
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file cannot have an empty 'annealing_temperature_ranges'."))
+    end 
+
+    for (i, temperature_range) in enumerate(config["annealing_temperature_ranges"])
+        for axis in ("T_min", "T_max", "D0", "H")
+            if !haskey(temperature_range, axis)
+                throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'annealing_temperature_ranges[$i]/$(axis)'."))
+            end
+        end
+        for k in keys(temperature_range)
+            if !(k in ("T_min", "T_max", "D0", "H"))
+                @warn "Encountered unexpected key `$(k)` in 'annealing_temperature_ranges[$i]': Allowed keys are `T_min`, `T_max`, `D0` and `H`. Key `$(k)` will be ignored."
+            end
+        end
+    end
+
+    if !haskey(config, "experimental_parameters") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'experimental_parameters'.")) 
+    elseif !haskey(config["experimental_parameters"], "a") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'experimental_parameters/a'.")) 
+    elseif !haskey(config["experimental_parameters"], "b") 
+        throw(ConfigFileError("ThermalDiffusionLithiumDensity config file needs entry 'experimental_parameters/b'.")) 
+    end
+
+    _ThermalDiffusionLithiumParameters(config; input_units, kwargs...)
+end
+
+const default_ThermalDiffusionLithiumDensity_config_file = joinpath(get_path_to_example_config_files(), "ThermalDiffusionLithiumDensity/Thermal_Diffusion_Config.yaml")
+function ThermalDiffusionLithiumParameters(config_filename::AbstractString = default_ThermalDiffusionLithiumDensity_config_file, input_units::Union{Missing, NamedTuple} = missing; kwargs...)
+    ThermalDiffusionLithiumParameters(parse_config_file(config_filename), input_units; kwargs...)
+end

src/MaterialProperties/MaterialProperties.jl

@@ -5,7 +5,7 @@ const elementary_charge = Float64(ustrip(u"C", Unitful.q))
 const ϵ0 = Float64(ustrip(u"F/m", Unitful.ϵ0))
 const kB = Float64(ustrip(u"J/K", Unitful.k))
 const me = Float64(ustrip(u"kg", Unitful.me))
-const R_gas = Float64(ustrip(u"cal/K/mol", Unitful.R))
+const R_gas = Float64(ustrip(internal_activation_energy_unit / internal_temperature_unit, Unitful.R))
 
 const material_properties = Dict{Symbol, NamedTuple}()
 

src/Units.jl

@@ -17,6 +17,7 @@ const internal_diffusion_unit = internal_length_unit ^ 2 / internal_time_unit
 const internal_charge_density_unit = internal_charge_unit / internal_length_unit ^ 3
 const internal_charge_density_gradient_unit = internal_charge_unit / internal_length_unit ^ 4
 const internal_temperature_unit = u"K"
+const internal_activation_energy_unit = u"cal/mol"
 
 const external_charge_unit  = u"e_au" # elementary charge - from UnitfulAtomic.jl
 
@@ -31,6 +32,7 @@ to_internal_units(x::Quantity{<:Real, dimension(internal_diffusion_unit)})  = us
 to_internal_units(x::Quantity{<:Real, dimension(internal_charge_density_unit)})  = ustrip(internal_charge_density_unit,  x)
 to_internal_units(x::Quantity{<:Real, dimension(internal_charge_density_gradient_unit)}) = ustrip(internal_charge_density_gradient_unit, x)
 to_internal_units(x::Quantity{<:Real, dimension(internal_temperature_unit)})     = ustrip(internal_temperature_unit,  x)
+to_internal_units(x::Quantity{<:Real, dimension(internal_activation_energy_unit)}) = ustrip(internal_activation_energy_unit,  x)
 
 from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_time_unit)})    = uconvert(unit, x * internal_time_unit)
 from_internal_units(x::Real, unit::Unitful.Units{<:Any, dimension(internal_voltage_unit)}) = uconvert(unit, x * internal_voltage_unit)