lindhard/ion/stopping/
oen_robinson.rs

1//! Oen-Robinson local (impact-parameter dependent) electronic energy loss.
2//!
3//! O. S. Oen and M. T. Robinson, Nucl. Instrum. Methods 132, 647 (1976).
4//!
5//! The nonlocal part is Lindhard-Scharff with the same `k_L`. The local model
6//! assigns a collision with closest approach `r_min` the loss
7//!
8//! `ΔE_e(E, r_min) = S_LS(E) · c² / (2π a²) · exp(-c r_min / a)`, `c = 0.3`,
9//!
10//! with `a` the Firsov screening length `0.8853 a0 (Z1^(1/2) + Z2^(1/2))^(-2/3)`.
11//! The prefactor normalises `∫ ΔE_e 2π p dp = S_LS` when `r_min` is identified
12//! with the impact parameter `p`; that identification is the approximation made
13//! in the original paper, so the impact-averaged stopping equals the LS value.
14//!
15//! The constants `c = 0.3` and the choice of length were entered from the
16//! literature as remembered by the contributor and have **not** been verified
17//! against the 1976 paper; see `docs/data-provenance.md` and
18//! [`OR_CONSTANTS_UNVERIFIED`]. A secondary source, J. F. Ziegler,
19//! J. P. Biersack, M. D. Ziegler, *SRIM: The Stopping and Range of Ions in
20//! Matter* (SRIM Co., 2015), Eq. 7-45, p. 7-17, prints
21//! `ΔE = 0.045 k √E / (π a²) exp(-0.3 r0 / a)`, which confirms `c = 0.3` and
22//! the prefactor `c² / (2π)`, but not which screening length `a` is.
23
24use super::lindhard_scharff::LindhardScharff;
25use super::{ElectronicStopping, Ion, StoppingError, ValidityRange};
26use crate::constants::BOHR_RADIUS;
27use std::f64::consts::PI;
28
29/// Whether the constants of the local model ([`OR_EXPONENT_COEFFICIENT`] and
30/// the length of [`local_length`]) are still **not verified** against Oen and
31/// Robinson, Nucl. Instrum. Methods 132, 647 (1976). `true` while the
32/// Oen-Robinson row of `docs/data-provenance.md` says "not verified".
33///
34/// This is the one place to change when the constants are checked against the
35/// paper: run metadata that reports the flag
36/// ([`crate::ion::bca::CrystalMetadata::electronic_constants_unverified`])
37/// follows it.
38pub const OR_CONSTANTS_UNVERIFIED: bool = true;
39
40/// Exponent coefficient `c` in `exp(-c r_min / a)`.
41pub const OR_EXPONENT_COEFFICIENT: f64 = 0.3;
42
43/// Firsov-type screening length used by the local model, m.
44pub fn local_length(z1: u8, z2: u8) -> f64 {
45    0.8853 * BOHR_RADIUS / (f64::from(z1).sqrt() + f64::from(z2).sqrt()).powf(2.0 / 3.0)
46}
47
48/// Oen-Robinson model. `stopping` is the impact-averaged value; use
49/// [`OenRobinson::local_loss`] for the per-collision loss.
50#[derive(Debug, Clone, Default)]
51pub struct OenRobinson {
52    ls: LindhardScharff,
53}
54
55impl OenRobinson {
56    /// New model with plain LS nonlocal coefficient.
57    pub fn new() -> Self {
58        Self::default()
59    }
60
61    /// Local electronic energy loss, J, for a collision with closest approach
62    /// `r_min_m` (metres), which must be finite and non-negative.
63    pub fn local_loss(
64        &self,
65        ion: &Ion,
66        target_z: u8,
67        energy_ev: f64,
68        r_min_m: f64,
69    ) -> Result<f64, StoppingError> {
70        let s = self.ls.stopping(ion, target_z, energy_ev)?;
71        if !(r_min_m.is_finite() && r_min_m >= 0.0) {
72            return Err(StoppingError::InvalidParameter {
73                name: "r_min_m",
74                value: r_min_m,
75            });
76        }
77        let a = local_length(ion.z(), target_z);
78        let c = OR_EXPONENT_COEFFICIENT;
79        Ok(s * c * c / (2.0 * PI * a * a) * (-c * r_min_m / a).exp())
80    }
81}
82
83impl ElectronicStopping for OenRobinson {
84    fn name(&self) -> &'static str {
85        "oen-robinson"
86    }
87
88    fn stopping(&self, ion: &Ion, target_z: u8, energy_ev: f64) -> Result<f64, StoppingError> {
89        self.ls.stopping(ion, target_z, energy_ev)
90    }
91
92    fn validity(&self, ion: &Ion) -> ValidityRange {
93        self.ls.validity(ion)
94    }
95}
96
97#[cfg(test)]
98mod tests {
99    use super::*;
100
101    #[test]
102    fn impact_average_recovers_stopping() {
103        // Integrate local_loss * 2 pi p dp over p numerically (r_min = p).
104        let ion = Ion::new(15).unwrap();
105        let m = OenRobinson::new();
106        let s = m.stopping(&ion, 14, 1.0e4).unwrap();
107        let a = local_length(15, 14);
108        let (n, pmax) = (200_000, 60.0 * a);
109        let dp = pmax / f64::from(n);
110        let mut sum = 0.0;
111        for i in 0..n {
112            let p = (f64::from(i) + 0.5) * dp;
113            sum += m.local_loss(&ion, 14, 1.0e4, p).unwrap() * 2.0 * PI * p * dp;
114        }
115        assert!((sum / s - 1.0).abs() < 1e-4, "{}", sum / s);
116    }
117
118    #[test]
119    fn loss_decreases_with_r_min() {
120        let ion = Ion::new(15).unwrap();
121        let m = OenRobinson::new();
122        let a = m.local_loss(&ion, 14, 1e4, 1e-11).unwrap();
123        let b = m.local_loss(&ion, 14, 1e4, 1e-10).unwrap();
124        assert!(a > b && b > 0.0);
125    }
126
127    #[test]
128    fn rejects_bad_closest_approach() {
129        let ion = Ion::proton();
130        let m = OenRobinson::new();
131        for bad in [-1e-10, f64::NAN, f64::INFINITY, f64::NEG_INFINITY] {
132            assert!(matches!(
133                m.local_loss(&ion, 14, 1e4, bad),
134                Err(StoppingError::InvalidParameter {
135                    name: "r_min_m",
136                    ..
137                })
138            ));
139        }
140    }
141
142    #[test]
143    fn zero_closest_approach_is_the_maximum_loss() {
144        let ion = Ion::proton();
145        let m = OenRobinson::new();
146        let s = m.stopping(&ion, 14, 1e4).unwrap();
147        let a = local_length(1, 14);
148        let c = OR_EXPONENT_COEFFICIENT;
149        let l0 = m.local_loss(&ion, 14, 1e4, 0.0).unwrap();
150        assert!((l0 / (s * c * c / (2.0 * PI * a * a)) - 1.0).abs() < 1e-14);
151    }
152}
lindhard/ion/stopping/oen_robinson.rs

lindhard/ion/stopping/
oen_robinson.rs
