lindhard/electron/inelastic/
drude.rs

1//! An analytic Drude-Lorentz energy-loss function with closed-form sum rules.
2//!
3//! This is a **synthetic** model ELF for tests, examples and benchmarks: its
4//! parameters are chosen by the caller, not fitted to any material, and its
5//! sum rules are known exactly, so a numerical sum-rule evaluation can be
6//! checked against them.
7//!
8//! # Model
9//!
10//! Each oscillator `i` has a strength `A_i` (eV²), a resonance energy `E_i`
11//! (eV) and a width `γ_i` (eV). The inverse dielectric function is
12//!
13//! ```text
14//! 1/ε(W) = 1 + Σ_i A_i / (W² - E_i² + i γ_i W),
15//! ```
16//!
17//! so the energy-loss function is
18//!
19//! ```text
20//! ELF(W) = Im[-1/ε(W)] = Σ_i A_i γ_i W / ((W² - E_i²)² + γ_i² W²),   W > 0.
21//! ```
22//!
23//! One term with `A = E_p²` and `E_i = E_p` is the Mermin-type optical ELF of
24//! Shinotsuka et al., Surf. Interface Anal. 49, 238 (2017),
25//! doi:10.1002/sia.6123, eq. (16) (`a_i = 1`); there it describes a free-electron
26//! plasmon of energy `E_p` and width `γ`.
27//!
28//! # Closed-form sum rules
29//!
30//! The two sum rules of the same paper, eqs. (18) and (19), are, for the full
31//! range `0 < W < ∞` and in energy units,
32//!
33//! ```text
34//! ∫_0^∞ W ELF(W) dW         = (π/2) Σ_i A_i           (f-sum)
35//! (2/π) ∫_0^∞ ELF(W) / W dW = Σ_i A_i / E_i²          (P_eff, without the n(0)⁻² term)
36//! ```
37//!
38//! Both follow from the two elementary integrals (by residues: the integrands
39//! are even in `W`, and `(W² - E²)² + γ² W² = |W² - E² + iγW|²` vanishes at
40//! `W = ±sqrt(E² - γ²/4) ± iγ/2`)
41//!
42//! ```text
43//! ∫_0^∞ γ W² / ((W² - E²)² + γ² W²) dW = π/2,
44//! ∫_0^∞ γ    / ((W² - E²)² + γ² W²) dW = π / (2 E²).
45//! ```
46//!
47//! The derivation is ours; no external source for these two integrals was
48//! consulted. The unit tests check both against an independent numerical
49//! quadrature. With `A = E_p²` and `E_i = E_p` (a free-electron plasmon), the
50//! f-sum gives the plasma energy `E_p` back and `P_eff = 1`, the values the
51//! paper states for a complete ELF (`P_eff -> 1` as `ΔE_max -> ∞`, text below
52//! eq. (19)).
53
54use crate::electron::data::{ElectronDataError, OpticalElf};
55
56/// One Drude-Lorentz oscillator (see the module docs).
57#[derive(Debug, Clone, Copy, PartialEq)]
58pub struct DrudeLorentzOscillator {
59    /// Strength `A`, eV². `A = E²` is a free-electron plasmon of energy `E`.
60    pub strength_ev2: f64,
61    /// Resonance energy `E`, eV.
62    pub energy_ev: f64,
63    /// Width `γ`, eV.
64    pub width_ev: f64,
65}
66
67impl DrudeLorentzOscillator {
68    /// A free-electron plasmon of energy `plasmon_ev` and width `width_ev`:
69    /// strength `plasmon_ev²`, resonance at `plasmon_ev`.
70    pub fn plasmon(plasmon_ev: f64, width_ev: f64) -> Self {
71        Self {
72            strength_ev2: plasmon_ev * plasmon_ev,
73            energy_ev: plasmon_ev,
74            width_ev,
75        }
76    }
77
78    fn elf(&self, w: f64) -> f64 {
79        let d = w * w - self.energy_ev * self.energy_ev;
80        self.strength_ev2 * self.width_ev * w / (d * d + self.width_ev * self.width_ev * w * w)
81    }
82}
83
84/// A sum of Drude-Lorentz oscillators: an analytic, synthetic ELF with
85/// closed-form sum rules (see the module docs).
86#[derive(Debug, Clone, PartialEq)]
87pub struct DrudeLorentz {
88    oscillators: Vec<DrudeLorentzOscillator>,
89}
90
91impl DrudeLorentz {
92    /// At least one oscillator; every strength, energy and width finite and
93    /// positive.
94    pub fn new(oscillators: Vec<DrudeLorentzOscillator>) -> Result<Self, ElectronDataError> {
95        if oscillators.is_empty() {
96            return Err(ElectronDataError::Invalid {
97                what: "Drude-Lorentz model",
98                reason: "needs at least one oscillator".into(),
99            });
100        }
101        for o in &oscillators {
102            for (name, v) in [
103                ("strength", o.strength_ev2),
104                ("energy", o.energy_ev),
105                ("width", o.width_ev),
106            ] {
107                if !(v.is_finite() && v > 0.0) {
108                    return Err(ElectronDataError::Invalid {
109                        what: "Drude-Lorentz oscillator",
110                        reason: format!("{name} must be finite and positive, got {v}"),
111                    });
112                }
113            }
114        }
115        Ok(Self { oscillators })
116    }
117
118    /// The oscillators.
119    pub fn oscillators(&self) -> &[DrudeLorentzOscillator] {
120        &self.oscillators
121    }
122
123    /// `ELF(W)` at energy loss `W` (eV); zero for `W <= 0`.
124    pub fn elf(&self, energy_ev: f64) -> f64 {
125        if energy_ev <= 0.0 {
126            return 0.0;
127        }
128        self.oscillators.iter().map(|o| o.elf(energy_ev)).sum()
129    }
130
131    /// `∫_0^∞ W ELF(W) dW = (π/2) Σ A_i`, eV² (closed form).
132    pub fn f_sum_ev2(&self) -> f64 {
133        std::f64::consts::FRAC_PI_2 * self.oscillators.iter().map(|o| o.strength_ev2).sum::<f64>()
134    }
135
136    /// The plasma energy `ħΩ_p = sqrt(Σ A_i)`, eV, defined by the f-sum
137    /// `∫ W ELF dW = (π/2) (ħΩ_p)²` (closed form).
138    pub fn plasma_energy_ev(&self) -> f64 {
139        self.oscillators
140            .iter()
141            .map(|o| o.strength_ev2)
142            .sum::<f64>()
143            .sqrt()
144    }
145
146    /// `P_eff = (2/π) ∫_0^∞ ELF(W)/W dW = Σ A_i / E_i²` (closed form).
147    pub fn p_eff(&self) -> f64 {
148        self.oscillators
149            .iter()
150            .map(|o| o.strength_ev2 / (o.energy_ev * o.energy_ev))
151            .sum()
152    }
153
154    /// Sample the model on `n` log-spaced energies from `lo_ev` to `hi_ev`
155    /// into an [`OpticalElf`] (interpolated linearly between the samples).
156    /// The provenance says that the table is synthetic and lists the
157    /// parameters.
158    pub fn to_optical_elf(
159        &self,
160        material: &str,
161        lo_ev: f64,
162        hi_ev: f64,
163        n: usize,
164    ) -> Result<OpticalElf, ElectronDataError> {
165        if !(lo_ev.is_finite() && hi_ev.is_finite() && lo_ev > 0.0 && hi_ev > lo_ev && n >= 2) {
166            return Err(ElectronDataError::Invalid {
167                what: "Drude-Lorentz sampling grid",
168                reason: format!("need 0 < lo < hi and n >= 2, got {lo_ev}, {hi_ev}, {n}"),
169            });
170        }
171        let ratio = hi_ev / lo_ev;
172        let energy: Vec<f64> = (0..n)
173            .map(|i| {
174                if i == n - 1 {
175                    hi_ev
176                } else {
177                    lo_ev * ratio.powf(i as f64 / (n - 1) as f64)
178                }
179            })
180            .collect();
181        let elf = energy.iter().map(|&w| self.elf(w)).collect();
182        let params: Vec<String> = self
183            .oscillators
184            .iter()
185            .map(|o| {
186                format!(
187                    "(A = {} eV^2, E = {} eV, gamma = {} eV)",
188                    o.strength_ev2, o.energy_ev, o.width_ev
189                )
190            })
191            .collect();
192        let provenance = format!(
193            "synthetic Drude-Lorentz ELF, not physical data: {}; {n} log-spaced samples from {lo_ev} to {hi_ev} eV",
194            params.join(", ")
195        );
196        OpticalElf::new(material, provenance, energy, elf)
197    }
198}
199
200#[cfg(test)]
201mod tests {
202    use super::*;
203    use crate::electron::inelastic::quadrature::{integrate_segments, GaussLegendre};
204
205    /// `∫_0^∞ f(W) dW` by the substitution `W = e^u`, `u` in [-30, 30] (the
206    /// integrands below fall off as `W` and `W^-2` or faster at the ends).
207    fn numeric(f: impl Fn(f64) -> f64, peak: f64) -> f64 {
208        let gl = GaussLegendre::new(10);
209        let breaks: Vec<f64> = [-30.0, peak.ln() - 1.0, peak.ln(), peak.ln() + 1.0, 30.0].to_vec();
210        integrate_segments(&gl, &mut |u: f64| [f(u.exp()) * u.exp()], &breaks, 1e-12)[0]
211    }
212
213    #[test]
214    fn closed_form_integrals_match_quadrature() {
215        for (e, g) in [(20.0, 5.0), (16.7, 0.5), (5.0, 12.0), (100.0, 300.0)] {
216            let den = |w: f64| (w * w - e * e).powi(2) + g * g * w * w;
217            let first = numeric(|w| g * w * w / den(w), e);
218            let second = numeric(|w| g / den(w), e);
219            assert!(
220                (first / std::f64::consts::FRAC_PI_2 - 1.0).abs() < 1e-10,
221                "{e} {g}: {first}"
222            );
223            assert!(
224                (second / (std::f64::consts::PI / (2.0 * e * e)) - 1.0).abs() < 1e-10,
225                "{e} {g}: {second}"
226            );
227        }
228    }
229
230    #[test]
231    fn sum_rules_of_a_two_oscillator_model() {
232        let m = DrudeLorentz::new(vec![
233            DrudeLorentzOscillator::plasmon(16.0, 4.0),
234            DrudeLorentzOscillator {
235                strength_ev2: 300.0,
236                energy_ev: 110.0,
237                width_ev: 60.0,
238            },
239        ])
240        .unwrap();
241        let f = numeric(|w| w * m.elf(w), 16.0);
242        let p = numeric(|w| m.elf(w) / w, 16.0) * 2.0 / std::f64::consts::PI;
243        assert!((f / m.f_sum_ev2() - 1.0).abs() < 1e-10);
244        assert!((p / m.p_eff() - 1.0).abs() < 1e-10);
245        assert!((m.plasma_energy_ev() - (256.0_f64 + 300.0_f64).sqrt()).abs() < 1e-12);
246    }
247
248    #[test]
249    fn rejects_bad_parameters() {
250        assert!(DrudeLorentz::new(vec![]).is_err());
251        assert!(DrudeLorentz::new(vec![DrudeLorentzOscillator::plasmon(10.0, 0.0)]).is_err());
252        assert!(DrudeLorentz::new(vec![DrudeLorentzOscillator::plasmon(f64::NAN, 1.0)]).is_err());
253        let m = DrudeLorentz::new(vec![DrudeLorentzOscillator::plasmon(10.0, 1.0)]).unwrap();
254        assert!(m.to_optical_elf("x", 1.0, 1.0, 10).is_err());
255        assert!(m.to_optical_elf("x", 1.0, 10.0, 1).is_err());
256        let t = m.to_optical_elf("x", 0.1, 1e3, 50).unwrap();
257        assert!(t.provenance().contains("synthetic"));
258        assert_eq!(t.energy_range_ev(), (0.1, 1e3));
259    }
260}
lindhard/electron/inelastic/drude.rs

lindhard/electron/inelastic/
drude.rs
