Skip to main content

lindhard/electron/inelastic/
model.rs

1//! Choosing the ELF-extension model per material.
2//!
3//! [`PennAlgorithm`] names the three models (the single-pole approximation,
4//! [`super::penn::SinglePolePenn`], the full Penn algorithm,
5//! [`super::full_penn::FullPenn`], and the Mermin-ELF oscillator fit,
6//! [`super::mermin::MerminPenn`]); [`PennInelastic`] holds one of them built
7//! from one material's optical ELF and gives all the same interface. The
8//! choice is made per material (one [`PennInelastic`] per
9//! [`OpticalElf`]) and is recorded in the metadata of whatever is derived
10//! from it: [`PennInelastic::model_identity`] is the string for the `model`
11//! field of a [`crate::electron::data::CrossSectionTable`], and it includes
12//! the algorithm's label, so a table cannot be mistaken for one of the other
13//! algorithm.
14
15use super::full_penn::FullPenn;
16use super::mermin::MerminPenn;
17use super::mermin_fit::MerminFitOptions;
18use super::penn::{InelasticPoint, SinglePolePenn};
19use crate::electron::data::{ElectronDataError, OpticalElf};
20
21type Result<T> = std::result::Result<T, ElectronDataError>;
22
23/// Which ELF-extension model to use.
24#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
25pub enum PennAlgorithm {
26    /// The single-pole approximation (Shinotsuka et al. 2017 eqs. (7)-(13)).
27    SinglePole,
28    /// The full Penn algorithm over Lindhard dielectric functions (eq. (4)).
29    Full,
30    /// A fit of Mermin oscillators to the optical ELF (MELF-GOS; dV2022
31    /// eqs. (3)-(5), see [`super::mermin`]), extended to finite momentum by
32    /// the Mermin dielectric function.
33    Mermin,
34}
35
36impl PennAlgorithm {
37    /// The stable label, `"penn-single-pole"`, `"penn-full"` or
38    /// `"mermin-melf"`.
39    pub fn label(self) -> &'static str {
40        match self {
41            Self::SinglePole => "penn-single-pole",
42            Self::Full => "penn-full",
43            Self::Mermin => "mermin-melf",
44        }
45    }
46
47    /// The algorithm with this label.
48    pub fn from_label(label: &str) -> Option<Self> {
49        [Self::SinglePole, Self::Full, Self::Mermin]
50            .into_iter()
51            .find(|a| a.label() == label)
52    }
53}
54
55impl std::str::FromStr for PennAlgorithm {
56    type Err = ElectronDataError;
57
58    fn from_str(s: &str) -> Result<Self> {
59        Self::from_label(s).ok_or_else(|| ElectronDataError::Invalid {
60            what: "Penn algorithm",
61            reason: format!("unknown label {s:?}; expected \"penn-single-pole\", \"penn-full\" or \"mermin-melf\""),
62        })
63    }
64}
65
66/// One of the three models, built from one material's optical ELF.
67#[derive(Debug, Clone)]
68pub enum PennInelastic {
69    /// The single-pole model.
70    SinglePole(SinglePolePenn),
71    /// The full model.
72    Full(FullPenn),
73    /// The Mermin-ELF model (a fit of oscillators).
74    Mermin(MerminPenn),
75}
76
77impl PennInelastic {
78    /// The model of `algorithm` for `elf` (Fermi energy zero, default
79    /// tolerance). For [`PennAlgorithm::Mermin`] the oscillators are fitted
80    /// with [`MerminFitOptions::default`]; use [`Self::try_new`] to handle a
81    /// failed fit, or [`Self::mermin`] to choose the fit options.
82    ///
83    /// # Panics
84    ///
85    /// Only for [`PennAlgorithm::Mermin`], if the default fit fails (an ELF
86    /// table too short for three oscillators, or an ELF that is zero).
87    pub fn new(algorithm: PennAlgorithm, elf: OpticalElf) -> Self {
88        Self::try_new(algorithm, elf).expect("the default Mermin fit failed; use try_new")
89    }
90
91    /// As [`Self::new`], returning the error of a failed Mermin fit.
92    pub fn try_new(algorithm: PennAlgorithm, elf: OpticalElf) -> Result<Self> {
93        Ok(match algorithm {
94            PennAlgorithm::SinglePole => Self::SinglePole(SinglePolePenn::new(elf)),
95            PennAlgorithm::Full => Self::Full(FullPenn::new(elf)),
96            PennAlgorithm::Mermin => Self::mermin(elf, &MerminFitOptions::default())?,
97        })
98    }
99
100    /// The Mermin-ELF model with oscillators fitted to `elf` with `options`.
101    pub fn mermin(elf: OpticalElf, options: &MerminFitOptions) -> Result<Self> {
102        Ok(Self::Mermin(MerminPenn::fit(elf, options)?))
103    }
104
105    /// Which algorithm this is.
106    pub fn algorithm(&self) -> PennAlgorithm {
107        match self {
108            Self::SinglePole(_) => PennAlgorithm::SinglePole,
109            Self::Full(_) => PennAlgorithm::Full,
110            Self::Mermin(_) => PennAlgorithm::Mermin,
111        }
112    }
113
114    /// The optical ELF the model is built from.
115    pub fn optical_elf(&self) -> &OpticalElf {
116        match self {
117            Self::SinglePole(m) => m.optical_elf(),
118            Self::Full(m) => m.optical_elf(),
119            Self::Mermin(m) => m.optical_elf(),
120        }
121    }
122
123    /// Set the Fermi energy, eV (see the two models).
124    pub fn with_fermi_energy_ev(self, fermi_ev: f64) -> Result<Self> {
125        Ok(match self {
126            Self::SinglePole(m) => Self::SinglePole(m.with_fermi_energy_ev(fermi_ev)?),
127            Self::Full(m) => Self::Full(m.with_fermi_energy_ev(fermi_ev)?),
128            Self::Mermin(m) => Self::Mermin(m.with_fermi_energy_ev(fermi_ev)?),
129        })
130    }
131
132    /// The Fermi energy, eV (see the two models).
133    pub fn fermi_energy_ev(&self) -> f64 {
134        match self {
135            Self::SinglePole(m) => m.fermi_energy_ev(),
136            Self::Full(m) => m.fermi_energy_ev(),
137            Self::Mermin(m) => m.fermi_energy_ev(),
138        }
139    }
140
141    /// The identity string to record as the `model` of derived tables:
142    /// the algorithm label, the nonrelativistic kinematics, the Fermi energy
143    /// and the integration tolerance, e.g.
144    /// `penn-full; nonrelativistic; E_F = 0 eV; rel tol 1e-4`. For the Mermin
145    /// model the fitted oscillators follow the label (their parameters
146    /// decide the table), e.g. `mermin-melf (3 oscillators (A=...))`.
147    pub fn model_identity(&self) -> String {
148        let (tol, fermi) = match self {
149            Self::SinglePole(m) => (m.relative_tolerance(), m.fermi_energy_ev()),
150            Self::Full(m) => (m.relative_tolerance(), m.fermi_energy_ev()),
151            Self::Mermin(m) => (m.relative_tolerance(), m.fermi_energy_ev()),
152        };
153        let label = match self {
154            Self::Mermin(m) => format!(
155                "{} ({})",
156                self.algorithm().label(),
157                m.fit_report().summary()
158            ),
159            _ => self.algorithm().label().to_string(),
160        };
161        format!("{label}; nonrelativistic; E_F = {fermi} eV; rel tol {tol:e}")
162    }
163
164    /// The DIIMFP, m⁻¹ eV⁻¹ (see the two models).
165    pub fn diimfp_per_m_ev(&self, energy_ev: f64, loss_ev: f64) -> Result<f64> {
166        match self {
167            Self::SinglePole(m) => m.diimfp_per_m_ev(energy_ev, loss_ev),
168            Self::Full(m) => m.diimfp_per_m_ev(energy_ev, loss_ev),
169            Self::Mermin(m) => m.diimfp_per_m_ev(energy_ev, loss_ev),
170        }
171    }
172
173    /// The inverse IMFP and the stopping power at `energy_ev`.
174    pub fn imfp_and_stopping(&self, energy_ev: f64) -> Result<InelasticPoint> {
175        match self {
176            Self::SinglePole(m) => m.imfp_and_stopping(energy_ev),
177            Self::Full(m) => m.imfp_and_stopping(energy_ev),
178            Self::Mermin(m) => m.imfp_and_stopping(energy_ev),
179        }
180    }
181
182    /// IMFP and stopping power on a grid of energies.
183    pub fn tabulate(&self, energies_ev: &[f64]) -> Result<Vec<InelasticPoint>> {
184        match self {
185            Self::SinglePole(m) => m.tabulate(energies_ev),
186            Self::Full(m) => m.tabulate(energies_ev),
187            Self::Mermin(m) => m.tabulate(energies_ev),
188        }
189    }
190}
191
192#[cfg(test)]
193mod tests {
194    use super::*;
195
196    #[test]
197    fn labels_round_trip_and_unknown_labels_are_rejected() {
198        for a in [
199            PennAlgorithm::SinglePole,
200            PennAlgorithm::Full,
201            PennAlgorithm::Mermin,
202        ] {
203            assert_eq!(PennAlgorithm::from_label(a.label()), Some(a));
204            assert_eq!(a.label().parse::<PennAlgorithm>().unwrap(), a);
205        }
206        assert!(PennAlgorithm::from_label("penn").is_none());
207        assert!("penn".parse::<PennAlgorithm>().is_err());
208    }
209}