lindhard/ion/

potential.rs

//! Screened-Coulomb interatomic potentials.
//!
//! The interaction of two atoms with atomic numbers `Z1`, `Z2` at separation
//! `r` is written
//!
//! ```text
//! V(r) = Z1 Z2 e^2 / (4 pi eps0 r) * phi(r / a)
//! ```
//!
//! where `phi` is a universal screening function (`phi(0) = 1`, decreasing to 0)
//! and `a` a screening length that depends on `Z1` and `Z2`. Everything in the
//! scattering code is expressed in the reduced variables `x = r / a` and
//! `eps = a E_cm / (Z1 Z2 e^2 / (4 pi eps0))`, in which `phi` is the only
//! input; the screening length only converts to and from SI.
//!
//! Coefficients follow the cited papers. The status of each set is recorded
//! in its `docs/data-provenance.md` row. The Firsov and Lindhard screening
//! lengths have been checked against the primary papers (page references
//! below). The ZBL, Kr-C and Moliere screening sets have been cross-checked
//! against an independent MIT-licensed implementation. The Lenz-Jensen set
//! agrees with a textbook tabulation, but not with a primary source. The
//! universal-length prefactor has been seen only in secondary sources.

use crate::constants::{BOHR_RADIUS, COULOMB_E2};

/// Screening function `phi(x)`, `x = r / a`.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum Screening {
    /// ZBL "universal" screening: four exponentials fitted to Hartree-Fock-Slater
    /// solid-state pair potentials. J. F. Ziegler, J. P. Biersack, U. Littmark,
    /// *The Stopping and Range of Ions in Solids* (Pergamon, 1985), Ch. 2.
    ZblUniversal,
    /// Kr-C potential: three exponentials fitted to Hartree-Fock Kr-Kr pair
    /// interaction. W. D. Wilson, L. G. Haggmark, J. P. Biersack,
    /// Phys. Rev. B 15, 2458 (1977).
    KrC,
    /// Moliere's approximation to the Thomas-Fermi function: three exponentials.
    /// G. Moliere, Z. Naturforsch. A 2, 133 (1947).
    Moliere,
    /// Lenz-Jensen screening, a polynomial-times-exponential approximation of
    /// the Thomas-Fermi-Jensen statistical model. W. Lenz, Z. Phys. 77, 713
    /// (1932), doi:10.1007/BF01342150; H. Jensen, Z. Phys. 77, 722 (1932),
    /// doi:10.1007/BF01342151. Polynomial form and argument scale as tabulated
    /// in W. Moller, *Fundamentals of Ion-Solid Interaction*, HZDR-073
    /// (Helmholtz-Zentrum Dresden-Rossendorf, 2017), p. 11, eq. (28); see
    /// `LJ_P` (private; see its comment in the source) for the verification status.
    LenzJensen,
}

/// Prefactors and exponents `(c_i, b_i)` of the sum-of-exponentials screening
/// functions, `phi(x) = Sum_i c_i exp(-b_i x)`.
///
/// ZBL: verified against J. F. Ziegler, J. P. Biersack, M. D. Ziegler, *SRIM: The Stopping and
/// Range of Ions in Matter* (SRIM Co., 2015), Ch. 2,
/// Fig. 2-17, p. 2-32 (these 4-digit values exactly) and Fig. 2-18, p. 2-36
/// (the 5-digit set 0.18175, 0.50986, 0.28022, 0.028171 / 3.1998, 0.94229,
/// 0.4029, 0.20162). The book's Eq. 2-74, p. 2-31, misprints the fourth
/// prefactor as `.2817` and the third exponent as `.4028`; the figures and
/// the normalisation `Sum c_i = 1` settle both.
///
/// Moliere: verified against the same book, Eq. 2-62, p. 2-20, which
/// misprints the middle prefactor as `5.5` (normalisation gives 0.55).
///
/// Kr-C and Moliere: cross-checked digit for digit against
/// `include/screened_coulomb.h` of the MIT-licensed `ir2-lab/screened_coulomb`
/// (commit f84c3c8), a component of OpenTRIM. This is a secondary source, not
/// the cited papers.
const ZBL_C: [(f64, f64); 4] = [
    (0.181_8, 3.2),
    (0.509_9, 0.942_3),
    (0.280_2, 0.402_9),
    (0.028_17, 0.201_6),
];
const KRC_C: [(f64, f64); 3] = [
    (0.190_945, 0.278_544),
    (0.473_674, 0.637_174),
    (0.335_381, 1.919_249),
];
const MOLIERE_C: [(f64, f64); 3] = [(0.35, 0.3), (0.55, 1.2), (0.10, 6.0)];

/// Lenz-Jensen: `phi = P(q) exp(-q)`, `q = sqrt(LJ_Q2 x)`,
/// `P(q) = 1 + q + 0.3344 q^2 + 0.0485 q^3 + 0.002647 q^4`.
///
/// All five polynomial coefficients and the scale `9.67` are printed exactly
/// so in W. Moller, *Fundamentals of Ion-Solid Interaction*, HZDR-073 (2017),
/// p. 11, eq. (28), `phi_LJ(xi) = (1 + xi + 0.3344 xi^2 + 0.0485 xi^3 +
/// 0.002647 xi^4) exp(-xi)`, `xi = sqrt(9.67 r/a)`
/// (<https://www.hzdr.de/publications/PublDoc-10091.pdf>). There `a` is the
/// single-atom Thomas-Fermi length `0.8853 a0 Z^(-1/3)` (eq. (26)). The
/// two-atom form replaces it with the Lindhard length (eqs. (29)-(30)), which
/// is the pairing in [`Screening::default_length`]. That is a secondary
/// (textbook) source. The primary papers of Lenz and Jensen were not
/// accessible, so the set is **unverified against the primary literature**. No
/// source was found for the alternative scale `3.11126^2 = 9.680`. See
/// `docs/data-provenance.md`.
const LJ_P: [f64; 5] = [1.0, 1.0, 0.3344, 0.0485, 0.002_647];
const LJ_Q2: f64 = 9.67;

fn sum_exp(c: &[(f64, f64)], x: f64) -> f64 {
    c.iter().map(|&(a, b)| a * (-b * x).exp()).sum()
}
fn sum_exp_d(c: &[(f64, f64)], x: f64) -> f64 {
    c.iter().map(|&(a, b)| -a * b * (-b * x).exp()).sum()
}

impl Screening {
    /// All implemented screening functions.
    pub const ALL: [Screening; 4] = [
        Screening::ZblUniversal,
        Screening::KrC,
        Screening::Moliere,
        Screening::LenzJensen,
    ];

    /// Screening function value `phi(x)`.
    pub fn phi(self, x: f64) -> f64 {
        match self {
            Screening::ZblUniversal => sum_exp(&ZBL_C, x),
            Screening::KrC => sum_exp(&KRC_C, x),
            Screening::Moliere => sum_exp(&MOLIERE_C, x),
            Screening::LenzJensen => {
                let q = (LJ_Q2 * x).sqrt();
                horner(&LJ_P, q) * (-q).exp()
            }
        }
    }

    /// Derivative `dphi/dx`.
    pub fn dphi(self, x: f64) -> f64 {
        match self {
            Screening::ZblUniversal => sum_exp_d(&ZBL_C, x),
            Screening::KrC => sum_exp_d(&KRC_C, x),
            Screening::Moliere => sum_exp_d(&MOLIERE_C, x),
            Screening::LenzJensen => {
                if x <= 0.0 {
                    return f64::NEG_INFINITY;
                }
                let q = (LJ_Q2 * x).sqrt();
                // d/dq [P e^-q] = (P' - P) e^-q ;  dq/dx = q / (2x).
                let p = horner(&LJ_P, q);
                let dp = LJ_P
                    .iter()
                    .enumerate()
                    .skip(1)
                    .rev()
                    .fold(0.0, |acc, (i, &c)| acc * q + c * i as f64);
                (dp - p) * (-q).exp() * q / (2.0 * x)
            }
        }
    }

    /// The `(c_i, b_i)` terms of a sum-of-exponentials screening function
    /// `phi(x) = Sum_i c_i exp(-b_i x)`, or `None` for Lenz-Jensen. Used by the
    /// analytic stopping asymptote in the scattering tests and the validation
    /// harness (`lindhard/tests/validation`).
    pub fn exponential_terms(self) -> Option<&'static [(f64, f64)]> {
        match self {
            Screening::ZblUniversal => Some(&ZBL_C),
            Screening::KrC => Some(&KRC_C),
            Screening::Moliere => Some(&MOLIERE_C),
            Screening::LenzJensen => None,
        }
    }

    /// The screening length conventionally paired with this function (the
    /// pairing used where the function was introduced or fitted):
    /// ZBL universal with the universal length, Kr-C and Moliere with the
    /// Firsov length, Lenz-Jensen with the Lindhard (Thomas-Fermi) length.
    pub fn default_length(self) -> ScreeningLength {
        match self {
            Screening::ZblUniversal => ScreeningLength::Universal,
            Screening::KrC | Screening::Moliere => ScreeningLength::Firsov,
            Screening::LenzJensen => ScreeningLength::Lindhard,
        }
    }
}

fn horner(c: &[f64], q: f64) -> f64 {
    c.iter().rev().fold(0.0, |acc, &k| acc * q + k)
}

/// Screening length `a(Z1, Z2)`.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum ScreeningLength {
    /// Universal length, `a_U = 0.8853 a0 / (Z1^0.23 + Z2^0.23)`.
    /// Ziegler, Biersack, Littmark (1985), Ch. 2; the page numbers are those of
    /// J. F. Ziegler, J. P. Biersack and M. D. Ziegler, *SRIM: The Stopping
    /// and Range of Ions in Matter* (SRIM Co., 2015), Eq. 2-73, p. 2-30, which
    /// prints the prefactor as `0.8853`; its Figs. 2-16 and 2-17 (p. 2-32) label it
    /// `.8854`. We use the Thomas-Fermi constant (`0.88534`, the same value
    /// as the Firsov and Lindhard lengths), which agrees with Eq. 2-73 to the
    /// printed digits; see `docs/data-provenance.md`.
    Universal,
    /// Firsov length, `a_F = 0.8853 a0 / (sqrt(Z1) + sqrt(Z2))^(2/3)`.
    /// O. B. Firsov, Sov. Phys. JETP 6, 534 (1958): the scaled distance
    /// `x' = (Z1^(1/2) + Z2^(1/2))^(2/3) R / a` (p. 536), with
    /// `a = (9 pi^2/128)^(1/3) hbar^2/(m e^2) = 4.68e-9 cm` (p. 535, below
    /// eq. (4)). Verified against the scan at
    /// <http://www.jetp.ras.ru/cgi-bin/dn/e_006_03_0534.pdf>.
    Firsov,
    /// Lindhard (Thomas-Fermi) length,
    /// `a_L = 0.8853 a0 / (Z1^(2/3) + Z2^(2/3))^(1/2)`.
    /// J. Lindhard, M. Scharff, H. E. Schiott, Mat. Fys. Medd. Dan. Vid. Selsk.
    /// 33(14) (1963), p. 8, printed as `a = a0 . 0.8853 (Z1^(2/3) +
    /// Z2^(2/3))^(-1/2)`. Verified against the scan at
    /// <http://gymarkiv.sdu.dk/MFM/kdvs/mfm%2030-39/mfm-33-14.pdf>.
    Lindhard,
}

/// The Thomas-Fermi constant `(9 pi^2 / 128)^(1/3) = 0.8853`, which multiplies
/// the Bohr radius in the Firsov and Lindhard lengths. The closed form is
/// Firsov's (Sov. Phys. JETP 6, 534 (1958), p. 535), and the rounded value
/// `0.8853` is the one printed by Lindhard, Scharff and Schiott (1963), p. 8.
pub fn thomas_fermi_constant() -> f64 {
    (9.0 * std::f64::consts::PI * std::f64::consts::PI / 128.0).cbrt()
}

impl ScreeningLength {
    /// Screening length in metres for atomic numbers `z1`, `z2`.
    pub fn metres(self, z1: f64, z2: f64) -> f64 {
        let a0 = BOHR_RADIUS;
        match self {
            ScreeningLength::Universal => {
                thomas_fermi_constant() * a0 / (z1.powf(0.23) + z2.powf(0.23))
            }
            ScreeningLength::Firsov => {
                thomas_fermi_constant() * a0 / (z1.sqrt() + z2.sqrt()).powf(2.0 / 3.0)
            }
            ScreeningLength::Lindhard => {
                thomas_fermi_constant() * a0 / (z1.powf(2.0 / 3.0) + z2.powf(2.0 / 3.0)).sqrt()
            }
        }
    }
}

/// A screened-Coulomb pair potential for a specific `(Z1, Z2)`.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Potential {
    /// Screening function.
    pub screening: Screening,
    /// Screening length convention.
    pub length: ScreeningLength,
    /// Projectile atomic number.
    pub z1: f64,
    /// Target atomic number.
    pub z2: f64,
}

impl Potential {
    /// Potential with the conventional length for the chosen screening function.
    pub fn new(screening: Screening, z1: f64, z2: f64) -> Self {
        Self {
            screening,
            length: screening.default_length(),
            z1,
            z2,
        }
    }

    /// Override the screening length.
    pub fn with_length(mut self, length: ScreeningLength) -> Self {
        self.length = length;
        self
    }

    /// Screening length `a`, metres.
    pub fn screening_length(&self) -> f64 {
        self.length.metres(self.z1, self.z2)
    }

    /// Potential energy at separation `r` (metres), joules.
    pub fn energy(&self, r: f64) -> f64 {
        self.z1 * self.z2 * COULOMB_E2 / r * self.screening.phi(r / self.screening_length())
    }

    /// Reduced centre-of-mass energy `eps = a E_cm / (Z1 Z2 e^2/(4 pi eps0))`,
    /// from `E_cm` in joules. (Lindhard-Scharff-Schiott reduced energy.)
    pub fn reduced_energy(&self, e_cm_j: f64) -> f64 {
        self.screening_length() * e_cm_j / (self.z1 * self.z2 * COULOMB_E2)
    }

    /// Centre-of-mass energy in joules from lab energy `e_lab_j` for projectile
    /// mass `m1` and target mass `m2` (any common unit): `E_cm = E M2/(M1+M2)`.
    pub fn cm_energy(e_lab_j: f64, m1: f64, m2: f64) -> f64 {
        e_lab_j * m2 / (m1 + m2)
    }

    /// Impact parameter in metres from reduced `beta = b / a`.
    pub fn impact_parameter(&self, beta: f64) -> f64 {
        beta * self.screening_length()
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn phi_is_one_at_origin_and_decreasing() {
        for s in Screening::ALL {
            let p0 = s.phi(1e-12);
            assert!((p0 - 1.0).abs() < 5e-3, "{s:?} phi(0) = {p0}");
            let mut prev = p0;
            for k in 1..200 {
                let v = s.phi(k as f64 * 0.1);
                assert!(v < prev && v > 0.0, "{s:?} not decreasing at {k}");
                prev = v;
            }
        }
    }

    #[test]
    fn derivative_matches_finite_difference() {
        for s in Screening::ALL {
            for &x in &[0.05, 0.3, 1.0, 3.0, 10.0] {
                let h = 1e-6 * x;
                let fd = (s.phi(x + h) - s.phi(x - h)) / (2.0 * h);
                let an = s.dphi(x);
                assert!((fd - an).abs() < 1e-6 * an.abs().max(1e-3), "{s:?} x={x}");
            }
        }
    }

    #[test]
    fn functions_agree_with_thomas_fermi_scale() {
        // All four approximate the Thomas-Fermi function, phi_TF(1) ~ 0.40.
        for s in Screening::ALL {
            let v = s.phi(1.0);
            assert!((0.3..0.5).contains(&v), "{s:?} phi(1) = {v}");
        }
    }

    #[test]
    fn screening_lengths_match_closed_forms() {
        // Explicit closed forms for Z1 = Z2 = 14.
        let a0 = BOHR_RADIUS;
        let tf = thomas_fermi_constant();
        assert!((tf - 0.8853).abs() < 1e-4);
        let a = ScreeningLength::Lindhard.metres(14.0, 14.0);
        let expect = tf * a0 / (2.0 * 14f64.powf(2.0 / 3.0)).sqrt();
        assert!((a / expect - 1.0).abs() < 1e-12);
        let a = ScreeningLength::Firsov.metres(14.0, 14.0);
        let expect = tf * a0 / (2.0 * 14f64.sqrt()).powf(2.0 / 3.0);
        assert!((a / expect - 1.0).abs() < 1e-12);
        // Universal: ZBL (1985) Eq. 2-73 prints 0.8853; same constant as above.
        let a = ScreeningLength::Universal.metres(5.0, 14.0);
        let expect = tf * a0 / (5f64.powf(0.23) + 14f64.powf(0.23));
        assert!((a / expect - 1.0).abs() < 1e-12);
    }

    /// The Thomas-Fermi prefactor against the numbers printed in the primary
    /// sources: Firsov (1958), p. 535, `a = (9 pi^2/128)^(1/3) hbar^2/(m e^2) =
    /// 4.68e-9 cm`, and Lindhard, Scharff, Schiott (1963), p. 8, `0.8853`.
    #[test]
    fn screening_length_prefactors_match_primary_sources() {
        let tf = thomas_fermi_constant();
        // LSS (1963) prints four digits.
        assert_eq!((tf * 1e4).round(), 8853.0);
        // Firsov (1958) prints a = 4.68e-9 cm. The value is 4.6850e-9 cm, so
        // the printed figure is truncated, not rounded: require agreement to
        // one unit in its last digit.
        let a_cm = tf * BOHR_RADIUS * 100.0;
        assert!((a_cm - 4.68e-9).abs() < 0.01e-9, "a = {a_cm} cm");
        // Firsov p. 536: R = a x' / (sqrt(Z1) + sqrt(Z2))^(2/3).
        let (z1, z2) = (5.0f64, 14.0f64);
        let expect = tf * BOHR_RADIUS / (z1.sqrt() + z2.sqrt()).powf(2.0 / 3.0);
        assert!((ScreeningLength::Firsov.metres(z1, z2) / expect - 1.0).abs() < 1e-12);
    }

    #[test]
    fn exponential_sets_are_normalised() {
        // phi(0) = 1 requires the prefactors to sum to 1.
        let sum = |c: &[(f64, f64)]| c.iter().map(|&(a, _)| a).sum::<f64>();
        assert!((sum(&ZBL_C) - 1.0).abs() < 1e-4);
        assert!((sum(&KRC_C) - 1.0).abs() < 1e-5);
        assert!((sum(&MOLIERE_C) - 1.0).abs() < 1e-12);
        // The ZBL exponents are 3.2, 0.9423, 0.4029, 0.2016.
        assert_eq!(ZBL_C[2].1, 0.4029);
    }

    #[test]
    fn reduced_energy_round_trip() {
        let p = Potential::new(Screening::ZblUniversal, 5.0, 14.0);
        let e = 1.0e3 * crate::units::J_PER_EV;
        let eps = p.reduced_energy(e);
        assert!(eps > 0.0 && eps < 1.0);
        // V(r) at r = a: Z1 Z2 e^2/a * phi(1), in units of E_cm = eps^-1 ...
        let a = p.screening_length();
        let v = p.energy(a);
        let expect = p.z1 * p.z2 * COULOMB_E2 / a * p.screening.phi(1.0);
        assert!((v / expect - 1.0).abs() < 1e-12);
    }
}