From d80891ae05b33190838a3ae53089eac207bfda64 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 13:47:27 -0700 Subject: [PATCH 01/18] Fixed unused code warnings in main and tests. --- src/main.rs | 1 - src/tests.rs | 1 + 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main.rs b/src/main.rs index a3b0301..8b02d4b 100644 --- a/src/main.rs +++ b/src/main.rs @@ -41,7 +41,6 @@ use rand::{SeedableRng, rngs::ChaCha8Rng}; //Load internal modules pub mod material; pub mod particle; -pub mod tests; pub mod interactions; pub mod bca; pub mod geometry; diff --git a/src/tests.rs b/src/tests.rs index 0dbdcfe..67646d6 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -2,6 +2,7 @@ use super::*; #[cfg(test)] use float_cmp::*; +#[cfg(test)] use rand::RngExt; #[test] From b93772facf187581a88a8e46ece57a9f1e337d22 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 13:59:49 -0700 Subject: [PATCH 02/18] Somehow branchless duff got overwritten - fixed. Removed warning from tests about mutable vector. --- src/math.rs | 19 ++++++------------- src/tests.rs | 3 +++ 2 files changed, 9 insertions(+), 13 deletions(-) diff --git a/src/math.rs b/src/math.rs index d8c5a59..3afc54a 100644 --- a/src/math.rs +++ b/src/math.rs @@ -4,17 +4,10 @@ use super::*; /// Duff et al., JCGT 2017 /// http://jcgt.org/published/0006/01/01/ pub fn duff_orthonormal_basis(n: Vector) -> (Vector, Vector) { - if n.z < 0. { - let a = 1.0 / (1.0 - n.z); - let b = n.x * n.y * a; - let b1 = Vector::new(1.0 - n.x*n.x*a, -b, n.x); - let b2 = Vector::new(b, n.y * n.y*a- 1.0, -n.y); - (b1, b2) - } else { - let a = 1.0 / (1.0 + n.z); - let b = -n.x * n.y * a; - let b1 = Vector::new(1.0 - n.x*n.x*a, b, -n.x); - let b2 = Vector::new(b, 1.0 - n.y*n.y*a, -n.y); - (b1, b2) - } + let sign = (1.0_f64).copysign(n.z); + let a = -1.0_f64 / (sign + n.z); + let b = n.x*n.y*a; + let b1 = Vector::new(1.0 + sign*n.x*n.x*a, sign*b, -sign*n.x); + let b2 = Vector::new(b, sign + n.y*n.y*a, -n.y); + (b1, b2) } \ No newline at end of file diff --git a/src/tests.rs b/src/tests.rs index 67646d6..8d68d80 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -781,7 +781,10 @@ fn test_momentum_conservation() { InteractionPotential::MORSE{D: 5.4971E-20, r0: 2.782E-10, alpha: 1.4198E10} ]; + #[cfg(feature = "cpr_rootfinder")] let mut rootfinders = vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}; 4]; + #[cfg(not(feature = "cpr_rootfinder"))] + let rootfinders = vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}; 4]; //[[{"CPR"={n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}]] #[cfg(feature = "cpr_rootfinder")] From 996fd3ee0744c1f6332178290deef124ac0d1d1d Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 14:12:43 -0700 Subject: [PATCH 03/18] clippy suggestions. --- src/geometry.rs | 24 +++++++++++------------- src/particle.rs | 2 -- 2 files changed, 11 insertions(+), 15 deletions(-) diff --git a/src/geometry.rs b/src/geometry.rs index 616d9ad..445560b 100644 --- a/src/geometry.rs +++ b/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() { diff --git a/src/particle.rs b/src/particle.rs index cc1cf00..89039fa 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -1,7 +1,5 @@ use super::*; -/// Rustbca's internal representation of the particle_parameters input. - fn default_vec_zero() -> Vec { vec![0] } From 54d9256910acb660a56f5e05902fdb43d1bb4d3c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 14:38:49 -0700 Subject: [PATCH 04/18] Cached screening lengths (needs testing); cargo clippy suggestions. --- src/bca.rs | 32 ++++++++++---------------------- src/geometry.rs | 2 +- src/interactions.rs | 28 +++++++++++++++++++--------- src/lib.rs | 4 ++-- 4 files changed, 32 insertions(+), 34 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index e1e300e..58cbdc2 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -510,8 +510,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)) @@ -583,18 +583,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 +630,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 +669,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.; diff --git a/src/geometry.rs b/src/geometry.rs index 445560b..624d36e 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -758,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); diff --git a/src/interactions.rs b/src/interactions.rs index 4687cf2..46b6c87 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -1,4 +1,5 @@ use super::*; +use cached::macros::cached; /// Analytic solutions to outermost root of the interaction potential. pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { @@ -271,17 +272,26 @@ 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(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(Za as u64, Zb as u64), + InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => lindhard_screening_length(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::COULOMB{Za: Z1, Zb: Z2} => zbl_screening_length(Za as u64, Zb as u64), 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.), } } +#[cached] +fn zbl_screening_length(Za: u64, Zb: u64) -> f64{ + 0.88534*A0/((Za as f64).powf(0.23) + (Zb as f64).powf(0.23)) +} +#[cached] +fn lindhard_screening_length(Za: u64, Zb: u64) -> f64 { + 0.8853*A0*((Za as f64).sqrt() + (Zb as f64).sqrt()).powf(-2./3.) +} + /// 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 { @@ -414,11 +424,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 diff --git a/src/lib.rs b/src/lib.rs index 2c9c9d9..7725658 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -364,7 +364,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag } #[no_mangle] -pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut f64, uy: &mut f64, uz: &mut f64, E1: &mut f64, Z1: &mut f64, m1: &mut f64, Ec1: &mut f64, Es1: &mut f64, Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64) { +pub unsafe extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut f64, uy: &mut f64, uz: &mut f64, E1: &mut f64, Z1: &mut f64, m1: &mut f64, Ec1: &mut f64, Es1: &mut f64, Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64) { assert!(E1 > &mut 0.0); @@ -689,7 +689,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { } #[no_mangle] -pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track_recoils: &mut bool, +pub unsafe extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track_recoils: &mut bool, ux: *mut f64, uy: *mut f64, uz: *mut f64, E1: *mut f64, Z1: *mut f64, m1: *mut f64, Ec1: *mut f64, Es1: *mut f64, num_species_target: &mut c_int, From 1ed27e013c5966287fbf9a2980ccd2a2609a36c2 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 14:41:57 -0700 Subject: [PATCH 05/18] Cache in cargo.toml --- Cargo.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Cargo.toml b/Cargo.toml index 8c6213a..1bd4fa3 100644 --- a/Cargo.toml +++ b/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" From 34ecb1787257dc58ce81ecb180a9bc7778ed7983 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 15:46:56 -0700 Subject: [PATCH 06/18] Trying a static lookup table for Lindhard screening lengths. --- src/interactions.rs | 33 ++++++++++++++++++++++++--------- 1 file changed, 24 insertions(+), 9 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 46b6c87..4e0bafd 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -1,5 +1,5 @@ use super::*; -use cached::macros::cached; +use std::sync::LazyLock; /// Analytic solutions to outermost root of the interaction potential. pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { @@ -274,24 +274,39 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote //ZBL screening length, Eckstein (4.1.8) InteractionPotential::ZBL => zbl_screening_length(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 => lindhard_screening_length(Za as u64, Zb as u64), - InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => lindhard_screening_length(Za as u64, Zb as u64), + 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} => zbl_screening_length(Za as u64, Zb as u64), - 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::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), } } -#[cached] -fn zbl_screening_length(Za: u64, Zb: u64) -> f64{ +const Z_MAX: usize = 120; +static LINDHARD_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; + lindhard_screening_length(Za as u64, Zb as u64) + } + ) +); + +pub fn zbl_screening_length(Za: u64, Zb: u64) -> f64{ 0.88534*A0/((Za as f64).powf(0.23) + (Zb as f64).powf(0.23)) } -#[cached] -fn lindhard_screening_length(Za: u64, Zb: u64) -> f64 { + +pub fn lindhard_screening_length(Za: u64, Zb: u64) -> f64 { 0.8853*A0*((Za as f64).sqrt() + (Zb as f64).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] +} + /// 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 { From 5e9daf87bb09cb8381583df9a218bfb1ea6b06ec Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 15:51:49 -0700 Subject: [PATCH 07/18] LUT ZBL as well; in the TiO2 example I got a small but significant speedup (13%) --- src/interactions.rs | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/src/interactions.rs b/src/interactions.rs index 4e0bafd..49ac5d0 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -294,6 +294,16 @@ static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock: } ) ); +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 u64, Zb as u64) + } + ) +); pub fn zbl_screening_length(Za: u64, Zb: u64) -> f64{ 0.88534*A0/((Za as f64).powf(0.23) + (Zb as f64).powf(0.23)) @@ -307,6 +317,10 @@ 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 lindhard_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 { From 00e1d081c04faa111e733c85259b3f3386df5bc3 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 16:00:51 -0700 Subject: [PATCH 08/18] Added some comments and minor cleanup --- src/interactions.rs | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 49ac5d0..11070a1 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -10,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 { @@ -272,25 +271,31 @@ 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 => zbl_screening_length(Za as u64, Zb as u64), + 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 => 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} => zbl_screening_length(Za as u64, Zb as u64), + 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 u64, Zb as u64) + lindhard_screening_length(Za as f64, Zb as f64) } ) ); @@ -300,24 +305,24 @@ static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( |i| { let Za = i / Z_MAX; let Zb = i % Z_MAX; - zbl_screening_length(Za as u64, Zb as u64) + zbl_screening_length(Za as f64, Zb as f64) } ) ); -pub fn zbl_screening_length(Za: u64, Zb: u64) -> f64{ - 0.88534*A0/((Za as f64).powf(0.23) + (Zb as f64).powf(0.23)) +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: u64, Zb: u64) -> f64 { - 0.8853*A0*((Za as f64).sqrt() + (Zb as f64).sqrt()).powf(-2./3.) +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 lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64{ +pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] } From d3c87772a421dc63f618e64dc980260ad348a223 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 16:11:25 -0700 Subject: [PATCH 09/18] Further cleanup --- src/interactions.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/interactions.rs b/src/interactions.rs index 11070a1..593e4e5 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -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) } - } } From 851865263f0c1f77a05a5bf9a82dae4b96309e05 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 16:50:56 -0700 Subject: [PATCH 10/18] First steps to make a LUT for LS stopping; this *might* break seeded values - is fp mult associative? Guess we will find out --- src/material.rs | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/material.rs b/src/material.rs index 7843e87..2e41498 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,5 +1,6 @@ use super::*; use rand::RngExt; +use std::sync::LazyLock; ///This helper function is a workaround to issue #368 in serde fn default_surface_binding_model() -> SurfaceBindingModel { @@ -299,8 +300,24 @@ impl Material { } } +const Z_MAX: usize = 120; +static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( + || + std::array::from_fn( + |i| { + let Za = i / Z_MAX; + let Zb = i % Z_MAX; + lindhard_scharff_stopping_power_constant(Za as f64, Zb as f64) + } + ) +); + +fn lindhard_scharff_stopping_power_constant(Za: f64, Zb: f64) -> f64 { + LINDHARD_SCHARFF_PREFACTOR*(Za*Za.cbrt().sqrt()*Zb)/(Za.cbrt().powi(2) + Zb.cbrt().powi(2)).powi(3).sqrt()*(AMU/Q).sqrt() +} + pub fn lindhard_scharff_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { - LINDHARD_SCHARFF_PREFACTOR*(Za*Za.cbrt().sqrt()*Zb)/(Za.cbrt().powi(2) + Zb.cbrt().powi(2)).powi(3).sqrt()*(E/Q/Ma*AMU).sqrt() + lindhard_scharff_stopping_power_constant(Za, Zb)*(E/Ma).sqrt() } pub fn bethe_bloch_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { From c277a71f952ae5f953554bccaeccb1d9ee9c5611 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 16:55:05 -0700 Subject: [PATCH 11/18] LS lookup table implemented. --- src/material.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/material.rs b/src/material.rs index 2e41498..8d21628 100644 --- a/src/material.rs +++ b/src/material.rs @@ -317,7 +317,7 @@ fn lindhard_scharff_stopping_power_constant(Za: f64, Zb: f64) -> f64 { } pub fn lindhard_scharff_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { - lindhard_scharff_stopping_power_constant(Za, Zb)*(E/Ma).sqrt() + LS_STOPPING_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]*(E/Ma).sqrt() } pub fn bethe_bloch_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { From 1db6143a107a0bac3ca36ba0d78366d85ba36f1b Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 17:07:17 -0700 Subject: [PATCH 12/18] Added LUTs for LS stopping, BV mean ionization pot. Seems to gain another ~5-10% speedup... --- src/material.rs | 30 ++++++++++++++++++++---------- 1 file changed, 20 insertions(+), 10 deletions(-) diff --git a/src/material.rs b/src/material.rs index 8d21628..748e4a5 100644 --- a/src/material.rs +++ b/src/material.rs @@ -301,6 +301,7 @@ impl Material { } const Z_MAX: usize = 120; +// Generating lookup tables for all possibilities turns out to be faster than calculating on the fly static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( || std::array::from_fn( @@ -320,26 +321,35 @@ pub fn lindhard_scharff_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, M LS_STOPPING_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]*(E/Ma).sqrt() } +static BV_EMPIRICAL_MEAN_IONIZATON_POT: LazyLock<[f64; Z_MAX]> = LazyLock::new( + || + std::array::from_fn( + |Zb| { + let I0 = if (Zb as f64) < 13. { + 12. + 7./ Zb as f64 + } else { + 9.76 + 58.5*(Zb as f64).powf(-1.19) + }; + (Zb as f64)*I0*Q + } + ) +); + pub fn bethe_bloch_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { let beta = (1. - 1./(1. + E/Ma/C.powi(2)).powi(2)).sqrt(); let v = beta*C; - // This term is an empirical fit to the mean ionization potential - let I0 = match Zb < 13. { - true => 12. + 7./Zb, - false => 9.76 + 58.5*Zb.powf(-1.19), - }; - let I = Zb*I0*Q; + let I = BV_EMPIRICAL_MEAN_IONIZATON_POT[Zb as usize]; //See Biersack and Haggmark - this looks like an empirical shell correction - let B = match Zb < 3. { - true => 100.*Za/Zb, - false => 5. + let B = if Zb < 3. { + 100.*Za/Zb + } else { + 5. }; let prefactor = BETHE_BLOCH_PREFACTOR*Zb*Za*Za/beta/beta; let eb = 2.*ME*v*v/I; - prefactor*(eb + 1. + B/eb).ln() } From 2fa41ac842347b8e5504f637fa7dc28d19a4750e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 17:13:40 -0700 Subject: [PATCH 13/18] Remove reference to deprecated accelerated ions feature --- src/bca.rs | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 58cbdc2..d054975 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -82,14 +82,6 @@ pub fn single_ion_bca(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.; From 7de3c068eb8d950991220259abd0cf5622c304e2 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 17:14:37 -0700 Subject: [PATCH 14/18] Remove reference to deprecated accelerated ions feature (2) --- src/bca.rs | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index d054975..e3e1529 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -173,14 +173,9 @@ pub fn single_ion_bca(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, From f4c42915b631dba60511f3b8317ca7f54af6ef12 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 19:44:08 -0700 Subject: [PATCH 15/18] Switch from dynamic to static dispatch in quadrature --- src/bca.rs | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index e3e1529..56435bc 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -512,7 +512,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 { +fn scattering_function_gl(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: F) -> Result + 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() { @@ -524,7 +526,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 { +fn scattering_function_gm(u: f64, impact_parameter: f64, r0: f64, relative_energy: f64, interaction_potential: F) -> Result + 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() { @@ -536,23 +540,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(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: F, n_points: usize) -> f64 + where F: Fn(f64) -> f64 + Clone +{ let x: Vec = (1..=n_points).map(|i| ((2.*i as f64 - 1.)/4./n_points as f64*PI).cos()).collect(); let w: Vec = (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::() } /// 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(impact_parameter: f64, relative_energy: f64, r0: f64, interaction_potential: F) -> f64 + where F: Fn(f64) -> f64 + Clone +{ let x: Vec = [0., -0.538469, 0.538469, -0.90618, 0.90618].iter().map(|x| x/2. + 1./2.).collect(); let w: Vec = [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::() } From c869745e9f1937954b06d66565cece0caba965c1 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 19:49:11 -0700 Subject: [PATCH 16/18] Switch deep recoil suppression to screening length lookup --- src/bca.rs | 2 +- src/physics.rs | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 56435bc..206bb04 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -151,7 +151,7 @@ pub fn single_ion_bca(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; diff --git a/src/physics.rs b/src/physics.rs index c8d59a7..304fce7 100644 --- a/src/physics.rs +++ b/src/physics.rs @@ -75,5 +75,4 @@ pub fn physics_loop(particle_input_array: Vec Date: Fri, 31 Jul 2026 19:50:03 -0700 Subject: [PATCH 17/18] Switch deep high energy ffps to screening length lookup --- src/bca.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/bca.rs b/src/bca.rs index 206bb04..a277eb6 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -222,7 +222,8 @@ pub fn determine_mfp_phi_impact_parameter(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 From 5e2c426970dc139a4b462a5314c7a19b2940ae5b Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 31 Jul 2026 19:55:58 -0700 Subject: [PATCH 18/18] Replaced stray powfs in lib.rs --- src/lib.rs | 2 +- src/material.rs | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 7725658..df3955f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1152,7 +1152,7 @@ pub fn reflect_single_ion_py(ion: &PyDict, target: &PyDict, vx: f64, vy: f64, vz let options = Options::default_options(false); - let velocity2 = vx.powf(2.) + vy.powf(2.) + vz.powf(2.); //m/s + let velocity2 = vx*vx + vy*vy + vz*vz; //m^2/s^2 let energy_eV = 0.5*m1*AMU*velocity2/EV; //EV let ux = vx/velocity2.sqrt(); diff --git a/src/material.rs b/src/material.rs index 748e4a5..8904d4a 100644 --- a/src/material.rs +++ b/src/material.rs @@ -279,7 +279,7 @@ impl Material { let S_high = bethe_bloch_stopping_power_cross_section(Za, *Zb, E, Ma); - // correction applied only to LS component + // ck correction applied only to LS component 1./(1./S_high + 1./(S_low*ck)) }, //Lindhard-Scharff, Oen-Robinson, Lindhard Equipartition @@ -289,7 +289,7 @@ impl Material { ElectronicStoppingMode::INTERPOLATEDPLUS{ci} => { let S_high = bethe_bloch_stopping_power_cross_section(Za, *Zb, E, Ma); - // correction applied only to LS component + // ck correction applied only to LS component (S_high.powf(-ci) + (S_low*ck).powf(-ci)).powf(-1./ci) },