Minor code fixes

Cargo.toml

@@ -28,6 +28,7 @@ hdf5 = {version = "0.7.1", optional = true}
 rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
 ndarray = {version = "0.14.0", features = ["serde"], optional = true}
 parry3d-f64 = {optional = true, version="0.2.0"}
+cached = {features=["proc_macro"], version="2.0.2"}
 
 [dependencies.pyo3]
 version = "0.19.0"

src/bca.rs

@@ -82,14 +82,6 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
             //Choose impact parameters and azimuthal angles for all collisions, and determine mean free path
             let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, material, options, rng);
 
-            #[cfg(feature = "accelerated_ions")]
-            let distance_to_target = if !material.inside(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z) {
-                let (x, y, z) = material.geometry.closest_point(particle_1.pos.x, particle_1.pos.y, particle_1.pos.z);
-                ((x - particle_1.pos.x).powi(2) + (y - particle_1.pos.y).powi(2) + (z - particle_1.pos.z).powi(2)).sqrt()
-            } else {
-                0.
-            };
-
             let mut total_energy_lost_to_recoils = 0.;
             let mut total_asymptotic_deflection = 0.;
             let mut normalized_distance_of_closest_approach = 0.;
@@ -159,7 +151,7 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
 
                         let n = material.total_number_density(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z);
                         //We just need the lindhard screening length here, so the particular potential is not important
-                        let a: f64 = interactions::screening_length(Za, Zb, InteractionPotential::MOLIERE);
+                        let a: f64 = interactions::lindhard_screening_length_lookup(Za as u64, Zb as u64);
                         let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E;
                         let estimated_range_of_recoils = (reduced_energy.powf(0.3) + 0.1).powi(3)/n/a/a;
 
@@ -181,14 +173,9 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
             }
 
             //Advance particle in space and track total distance traveled
-            #[cfg(not(feature = "accelerated_ions"))]
             let distance_traveled = particle_1.advance(
                 binary_collision_geometries[0].mfp, total_asymptotic_deflection);
 
-            #[cfg(feature = "accelerated_ions")]
-            let distance_traveled = particle_1.advance(
-                binary_collision_geometries[0].mfp + distance_to_target - material.geometry.get_energy_barrier_thickness(), total_asymptotic_deflection);
-
             //Subtract total energy from all simultaneous collisions and electronic stopping
             let energy_lost_to_electronic_stopping = bca::subtract_electronic_stopping_energy(&mut particle_1, material, distance_traveled,
                 normalized_distance_of_closest_approach, strong_collision_Z,
@@ -235,7 +222,8 @@ pub fn determine_mfp_phi_impact_parameter<T: Geometry>(particle_1: &mut particle
         let E: f64  = particle_1.E;
         let Ec: f64 = particle_1.Ec;
         //We just need the Lindhard screening length here, so the particular potential is not important
-        let a: f64 = interactions::screening_length(Za, Zb, InteractionPotential::MOLIERE);
+        let a: f64 = interactions::lindhard_screening_length_lookup(Za as u64, Zb as u64);
+
         let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E;
 
         //Minimum energy transfer for generating scattering event set to cutoff energy
@@ -510,8 +498,8 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: &
         _ => x0*a*(theta/2.).sin()
     };
 
-    let psi = theta.sin().atan2(Ma/Mb + theta.cos());//.abs();
-    let psi_recoil = theta.sin().atan2(1. - theta.cos());//.abs();
+    let psi = theta.sin().atan2(Ma/Mb + theta.cos());
+    let psi_recoil = theta.sin().atan2(1. - theta.cos());
     let recoil_energy = 4.*(Ma*Mb)/(Ma + Mb).powi(2)*E0*(theta/2.).sin().powi(2);
 
     Ok(BinaryCollisionResult::new(theta, psi, psi_recoil, recoil_energy, asymptotic_deflection, x0))
@@ -525,7 +513,9 @@ fn scattering_integral_mw(x: f64, beta: f64, reduced_energy: f64, interaction_po
 }
 
 /// Gauss-Legendre scattering integrand.
-fn scattering_function_gl(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: &dyn Fn(f64) -> f64) -> Result<f64, anyhow::Error> {
+fn scattering_function_gl<F>(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: F) -> Result<f64, anyhow::Error> 
+    where F: Fn(f64) -> f64
+{
     let result = 4.*impact_parameter*u/(r0*(1. - interaction_potential(r0/(1. - u*u))/relative_energy - impact_parameter*impact_parameter*(1. - u*u).powi(2)/r0/r0).sqrt());
 
     if result.is_nan() {
@@ -537,7 +527,9 @@ fn scattering_function_gl(u: f64, impact_parameter: f64, r0: f64, relative_energ
 }
 
 /// Gauss-Mehler scattering integrand.
-fn scattering_function_gm(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: &dyn Fn(f64) -> f64) -> Result<f64, anyhow::Error> {
+fn scattering_function_gm<F>(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: F) -> Result<f64, anyhow::Error> 
+    where F: Fn(f64) -> f64    
+{
     let result = impact_parameter/r0/(1. - interaction_potential(r0/u)/relative_energy - (impact_parameter*u/r0).powi(2)).sqrt();
 
     if result.is_nan() {
@@ -549,23 +541,27 @@ fn scattering_function_gm(u: f64, impact_parameter: f64, r0: f64, relative_energ
 }
 
 /// Compute the scattering integral for a given relative energy, distance of closest approach `r0`,  and interaction potential using a Gauss-Mehler, n-point quadrature.
-fn scattering_integral_gauss_mehler(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: &dyn Fn(f64) -> f64, n_points: usize) -> f64 {
+fn scattering_integral_gauss_mehler<F>(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: F, n_points: usize) -> f64 
+    where F: Fn(f64) -> f64 + Clone
+{
     let x: Vec<f64> = (1..=n_points).map(|i| ((2.*i as f64 - 1.)/4./n_points as f64*PI).cos()).collect();
     let w: Vec<f64> = (1..=n_points).map(|i| PI/n_points as f64*((2.*i as f64 - 1.)/4./n_points as f64*PI).sin()).collect();
 
     PI - x.iter().zip(w)
-        .map(|(&x, w)| w*scattering_function_gm(x, impact_parameter, r0, relative_energy, interaction_potential)
+        .map(|(&x, w)| w*scattering_function_gm(x, impact_parameter, r0, relative_energy, interaction_potential.clone())
         .with_context(|| format!("Numerical error: NaN in Gauss-Mehler scattering integral at x = {} with Er = {} eV and p = {} A.", x, relative_energy/EV, impact_parameter/ANGSTROM))
         .unwrap()).sum::<f64>()
 }
 
 /// Compute the scattering integral for a given relative energy, distance of closest approach `r0`,  and interaction potential using a Gauss-Legendre, 5-point quadrature.
-fn scattering_integral_gauss_legendre(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: &dyn Fn(f64) -> f64) -> f64 {
+fn scattering_integral_gauss_legendre<F>(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: F) -> f64 
+    where F: Fn(f64) -> f64 + Clone
+{
     let x: Vec<f64> = [0., -0.538469, 0.538469, -0.90618, 0.90618].iter().map(|x| x/2. + 1./2.).collect();
     let w: Vec<f64> = [0.568889, 0.478629, 0.478629, 0.236927, 0.236927].iter().map(|w| w/2.).collect();
 
     PI - x.iter().zip(w)
-        .map(|(&x, w)| w*scattering_function_gl(x, impact_parameter, r0, relative_energy, interaction_potential)
+        .map(|(&x, w)| w*scattering_function_gl(x, impact_parameter, r0, relative_energy, interaction_potential.clone())
         .with_context(|| format!("Numerical error: NaN in Gauss-Legendre scattering integral at x = {} with Er = {} eV and p = {} A.", x, relative_energy/EV, impact_parameter/ANGSTROM))
         .unwrap()).sum::<f64>()
 }
@@ -583,18 +579,6 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact
     let coefficients = interactions::polynomial_coefficients(relative_energy, impact_parameter, interaction_potential);
     let roots = real_polynomial_roots(coefficients.clone(), polynom_complex_threshold).unwrap();
 
-    /*
-    println!("p={}", impact_parameter/ANGSTROM);
-
-    for coefficient in &coefficients {
-        println!("{}", {coefficient});
-    }
-
-    for root in &roots {
-        println!("{} A", {root});
-    }
-    */
-
     let max_root = roots.iter().cloned().fold(f64::NAN, f64::max);
 
     let inverse_transformed_root = interactions::inverse_transform(max_root, interaction_potential);
@@ -642,8 +626,6 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame
     let g = |r: f64| -> f64 {interactions::distance_of_closest_approach_function_singularity_free(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)*
         interactions::scaling_function(r, impact_parameter, interaction_potential)};
 
-    //Using upper bound const, ~10, construct upper bound as a plateau near 0 and linear increase away from that
-    //let upper_bound = f64::max(upper_bound_const*p, upper_bound_const*a);
     let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential);
 
     let roots = match derivative_free {
@@ -683,14 +665,16 @@ pub fn newton_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par
     let f = |r: f64| -> f64 {interactions::distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)};
     let df = |r: f64| -> f64 {interactions::diff_distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)};
 
-    //Guess for large reduced energy from Mendenhall and Weller 1991
-    //For small energies, use pure Newton-Raphson with arbitrary guess of 1
-    let mut x0 = beta;
     let mut xn: f64;
-    if reduced_energy > 5. {
+
+    //Guess for large reduced energy from Mendenhall and Weller 1991
+    //For small energies, use pure Newton-Raphson with arbitrary guess of beta
+    let mut x0 = if reduced_energy > 5. {
         let inv_er_2 = 0.5/reduced_energy;
-        x0 = inv_er_2 + (inv_er_2*inv_er_2 + beta*beta).sqrt();
-    }
+        inv_er_2 + (inv_er_2*inv_er_2 + beta*beta).sqrt()
+    } else {
+        beta
+    };
 
     //Newton-Raphson to determine distance of closest approach
     let mut err: f64 = tolerance + 1.;

src/geometry.rs

@@ -57,7 +57,7 @@ impl Geometry for Mesh0D {
             "NM" => NM,
             "M" => 1.,
             _ => input.length_unit.parse()
-                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
+                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
                         &input.length_unit.as_str())),
         };
 
@@ -295,7 +295,7 @@ impl Geometry for HomogeneousMesh2D {
             "NM" => NM,
             "M" => 1.,
             _ => input.length_unit.parse()
-                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
+                .unwrap_or_else(|_| panic!("Input errror: could nor parse length unit {}. Use a valid float or one of ANGSTROM, NM, MICRON, CM, MM, M",
                         &input.length_unit.as_str())),
         };
 
@@ -353,16 +353,14 @@ impl Geometry for HomogeneousMesh2D {
     fn inside_energy_barrier(&self, x: f64, y: f64, z: f64) -> bool {
         if self.inside(x, y, z) {
             true
-        } else {
-            if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
-                let distance = ((x - p.x()).powi(2) +  (y - p.y()).powi(2)).sqrt();
-                distance < self.energy_barrier_thickness
-            } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
-                true
-            } else {
-                panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
-            }
-        }
+        } else if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+             let distance = ((x - p.x()).powi(2) +  (y - p.y()).powi(2)).sqrt();
+             distance < self.energy_barrier_thickness
+         } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
+             true
+         } else {
+             panic!("Geometry error: closest point routine failed to find single closest point to ({}, {}, {}).", x, y, z);
+         }
     }
     fn closest_point(&self, x: f64, y: f64, z: f64) -> (f64, f64, f64) {
         if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) {
@@ -401,7 +399,7 @@ impl Mesh2D {
     /// Finds the cell that is nearest to (x, y).
     fn nearest_to(&self, x: f64, y: f64, z: f64) -> &Cell2D {
 
-        let mut min_distance: f64 = std::f64::MAX;
+        let mut min_distance: f64 = f64::MAX;
         let mut index: usize = 0;
 
         for (cell_index, cell) in self.mesh.iter().enumerate() {
@@ -760,7 +758,7 @@ impl Triangle2D {
 
     /// Calculates the shortest distance from this triangle to the point (x, y).
     pub fn distance_to(&self, x: f64, y: f64) -> f64 {
-        let mut distance_to = std::f64::MAX;
+        let mut distance_to = f64::MAX;
 
         for segment in &self.segments {
             let length_2 = (segment.2 - segment.0)*(segment.2 - segment.0) + (segment.3 - segment.1)*(segment.3 - segment.1);

src/interactions.rs

@@ -1,4 +1,5 @@
 use super::*;
+use std::sync::LazyLock;
 
 /// Analytic solutions to outermost root of the interaction potential.
 pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 {
@@ -9,7 +10,6 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 {
         InteractionPotential::WW => 50.*ANGSTROM,
         _ => 10.*ANGSTROM,
     }
-
 }
 
 fn smootherstep(x: f64, k: f64, x0: f64) -> f64 {
@@ -53,7 +53,6 @@ pub fn interaction_potential(r: f64, a: f64, Za: f64, Zb: f64, interaction_poten
         InteractionPotential::FOUR_EIGHT{alpha, beta} => {
             four_eight(r, alpha, beta)
         }
-
     }
 }
 
@@ -271,17 +270,61 @@ pub fn dphi(xi: f64, interaction_potential: InteractionPotential) -> f64 {
 pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 {
     match interaction_potential {
         //ZBL screening length, Eckstein (4.1.8)
-        InteractionPotential::ZBL => 0.88534*A0/(Za.powf(0.23) + Zb.powf(0.23)),
+        InteractionPotential::ZBL => zbl_screening_length_lookup(Za as u64, Zb as u64),
         //Lindhard/Firsov screening length, Eckstein (4.1.5)
-        InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN | InteractionPotential::WW => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.),
-        InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.),
+        InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN | InteractionPotential::WW => lindhard_screening_length_lookup(Za as u64, Zb as u64),
+        InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64),
         InteractionPotential::MORSE{D, alpha, r0} => alpha,
-        InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => 0.88534*A0/(Z1.powf(0.23) + Z2.powf(0.23)),
-        InteractionPotential::KRC_MORSE{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.),
-        InteractionPotential::FOUR_EIGHT{..} =>0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.),
+        InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => zbl_screening_length_lookup(Za as u64, Zb as u64),
+        InteractionPotential::KRC_MORSE{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64),
+        InteractionPotential::FOUR_EIGHT{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64),
     }
 }
 
+// It turns out it's faster (~10% speedup) to just generate every possible screening length as a lookup table
+// LazyLock is a thread-safe value that is initialized whenever it is first accessed
+// It will block other threads while it runs, but it should run extremely quickly and only once
+const Z_MAX: usize = 120;
+static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new(
+    ||
+    std::array::from_fn(
+        |i| {
+            // standard 2D to 1D  array indexing
+            // I always have to look this up; copied from here
+            // (e.g. https://stackoverflow.com/questions/5494974/convert-1d-array-index-to-2d-array-index)
+            let Za = i / Z_MAX;
+            let Zb = i % Z_MAX;
+            lindhard_screening_length(Za as f64, Zb as f64)
+        }
+    )
+);
+static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new(
+    ||
+    std::array::from_fn(
+        |i| {
+            let Za = i / Z_MAX;
+            let Zb = i % Z_MAX;
+            zbl_screening_length(Za as f64, Zb as f64)
+        }
+    )
+);
+
+pub fn zbl_screening_length(Za: f64, Zb: f64) -> f64{
+    0.88534*A0/(Za.powf(0.23) + Zb.powf(0.23))
+}
+
+pub fn lindhard_screening_length(Za: f64, Zb: f64) -> f64 {
+    0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.)
+}
+
+pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 {
+    LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize]
+}
+
+pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{
+    ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize]
+}
+
 /// Coefficients of inverse-polynomial interaction potentials.
 /// Note: these go in order of decreasing orders of r, e.g., [r^n, r^n-1, r^n-2, ... r, 1].
 pub fn polynomial_coefficients(relative_energy: f64, impact_parameter: f64, interaction_potential: InteractionPotential) -> Vec<f64> {
@@ -414,11 +457,11 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 {
 
     } else if x <= x2 {
 
-        let a = [1.389653276380862E4,
-            -3.596912431628216E4,
-            3.739206756369099E4,
-            -1.933748081656593E4,
-            0.495516793802426E4,
+        let a = [1.389653276380862E4,
+            -3.596912431628216E4,
+            3.739206756369099E4,
+            -1.933748081656593E4,
+            0.495516793802426E4,
             -0.050264585985867E4];
 
         (a[0] + a[1]*x + a[2]*x.powi(2) + a[3]*x.powi(3) + a[4]*x.powi(4) + a[5]*x.powi(5))*EV