From 723cb5b2684a880ab14f4416ca0df7143f747e08 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 12:04:07 -0700 Subject: [PATCH 01/20] Experimental triangular lookup tables for Lindhard screening length and LS stopping constant --- src/interactions.rs | 44 ++++++++++++++++++++++++-------------------- src/material.rs | 39 +++++++++++++++++++++++++++------------ src/math.rs | 7 +++++++ 3 files changed, 58 insertions(+), 32 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 36ca4c83..d634ada3 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -1,7 +1,9 @@ use super::*; use std::sync::LazyLock; +use crate::math::triangular_index; -const Z_MAX: usize = 120; +const Z_MAX: usize = 88; +const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2; /// Analytic solutions to outermost root of the interaction potential. pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { @@ -304,20 +306,31 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote // 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 -static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( +static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = 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) + { + let mut array = [0.0; TABLE_SIZE]; + for i in 0..Z_MAX { + for j in 0..=i { + let index = (i * (i + 1))/2 + j; + array[index] = lindhard_screening_length(i as f64, j as f64); + } } - ) + array + } ); +pub fn lindhard_screening_length(Za: f64, Zb: f64) -> f64 { + 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.) +} + +#[inline] +pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { + let mut i = Za as usize; + let mut j = Zb as usize; + LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] +} + static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( || std::array::from_fn( @@ -333,15 +346,6 @@ 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.) -} - -#[inline] -pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { - LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] -} - #[inline] pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] diff --git a/src/material.rs b/src/material.rs index b55bb5cb..9c030c8e 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,6 +1,8 @@ use super::*; use rand::RngExt; use std::sync::LazyLock; +use crate::math::triangular_index; +const Z_MAX: usize = 88; ///This helper function is a workaround to issue #368 in serde fn default_surface_binding_model() -> SurfaceBindingModel { @@ -300,26 +302,39 @@ impl Material { } } -const Z_MAX: usize = 120; +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() +} + +//https://math.stackexchange.com/questions/2388887/ +//num elements in a triangular NxN matrix (including diag) +const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2; + // 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( +// Tables for Za, Zb are upper-triangular +static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; TABLE_SIZE]> = 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) + { + let mut array = [0.0; TABLE_SIZE]; + for i in 0..Z_MAX { + for j in 0..=i { + let index = (i * (i + 1))/2 + j; + array[index] = lindhard_scharff_stopping_power_constant(i as f64, j as f64); + } } - ) + array + } ); -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() -} #[inline] pub fn lindhard_scharff_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { - LS_STOPPING_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]*(E/Ma).sqrt() + + let mut i = Za as usize; + let mut j = Zb as usize; + + LS_STOPPING_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]*(E/Ma).sqrt() } + static BV_EMPIRICAL_MEAN_IONIZATON_POT: LazyLock<[f64; Z_MAX]> = LazyLock::new( || std::array::from_fn( diff --git a/src/math.rs b/src/math.rs index 3afc54a5..2b99ea69 100644 --- a/src/math.rs +++ b/src/math.rs @@ -10,4 +10,11 @@ pub fn duff_orthonormal_basis(n: Vector) -> (Vector, Vector) { 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) +} + +pub fn triangular_index(i: &mut usize, j: &mut usize) -> usize { + if i < j { + std::mem::swap(i, j); + } + (*i*(*i + 1)/2) + *j } \ No newline at end of file From f5bc2bfd4797013b52e7bcf48fbaab2dfcc3561c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 12:17:22 -0700 Subject: [PATCH 02/20] Draft working version of LUTs --- src/consts.rs | 6 +++++ src/interactions.rs | 54 +++++++++++++++++++++++++-------------------- src/material.rs | 5 ----- 3 files changed, 36 insertions(+), 29 deletions(-) diff --git a/src/consts.rs b/src/consts.rs index 052627db..b8383232 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -35,3 +35,9 @@ pub const BETHE_BLOCH_PREFACTOR: f64 = 4.*PI*(Q*Q/(4.*PI*EPS0))*(Q*Q/(4.*PI*EPS0 pub const LINDHARD_SCHARFF_PREFACTOR: f64 = 1.212*ANGSTROM*ANGSTROM*Q; /// Lindhard reduced energy prefactor, in SI units. pub const LINDHARD_REDUCED_ENERGY_PREFACTOR: f64 = 4.*PI*EPS0/Q/Q; +/// Maximum atomic number in RustBCA. +pub const Z_MAX: usize = 120; +//https://math.stackexchange.com/questions/2388887/ +//num elements in a triangular NxN matrix (including diag) +/// LUT size; triangular matrix of size Z_MAX * Z_MAX +pub const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2; \ No newline at end of file diff --git a/src/interactions.rs b/src/interactions.rs index d634ada3..13c88b62 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -2,9 +2,6 @@ use super::*; use std::sync::LazyLock; use crate::math::triangular_index; -const Z_MAX: usize = 88; -const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2; - /// Analytic solutions to outermost root of the interaction potential. pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { @@ -236,18 +233,18 @@ pub fn diff_distance_of_closest_approach_function_singularity_free(r: f64, a: f6 } } -static COULOMB_CONSTANT_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( +static COULOMB_CONSTANT_TABLE: LazyLock<[f64; TABLE_SIZE]> = 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; - coulomb_constant(Za as f64, Zb as f64) + { + let mut array = [0.0; TABLE_SIZE]; + for i in 0..Z_MAX { + for j in 0..=i { + let index = (i * (i + 1))/2 + j; + array[index] = coulomb_constant(i as f64, j as f64); + } } - ) + array + } ); fn coulomb_constant(Za: f64, Zb: f64) -> f64 { @@ -256,12 +253,16 @@ fn coulomb_constant(Za: f64, Zb: f64) -> f64 { /// Screened coulomb interaction potential. pub fn screened_coulomb(r: f64, a: f64, Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 { - COULOMB_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]/r*phi(r/a, interaction_potential) + let mut i = Za as usize; + let mut j = Zb as usize; + COULOMB_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]/r*phi(r/a, interaction_potential) } /// Coulombic interaction potential. pub fn coulomb(r: f64, Za: f64, Zb: f64) -> f64 { - COULOMB_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]/r + let mut i = Za as usize; + let mut j = Zb as usize; + COULOMB_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]/r } /// Screening functions for screened-coulomb interaction potentials. @@ -303,7 +304,7 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote } } -// It turns out it's faster (~10% speedup) to just generate every possible screening length as a lookup table +// It turns out it's faster 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 static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( @@ -331,15 +332,18 @@ pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] } -static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( +static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = 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) + { + let mut array = [0.0; TABLE_SIZE]; + for i in 0..Z_MAX { + for j in 0..=i { + let index = (i * (i + 1))/2 + j; + array[index] = zbl_screening_length(i as f64, j as f64); + } } - ) + array + } ); pub fn zbl_screening_length(Za: f64, Zb: f64) -> f64{ @@ -348,7 +352,9 @@ pub fn zbl_screening_length(Za: f64, Zb: f64) -> f64{ #[inline] pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ - ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] + let mut i = Za as usize; + let mut j = Zb as usize; + ZBL_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] } /// Coefficients of inverse-polynomial interaction potentials. diff --git a/src/material.rs b/src/material.rs index 9c030c8e..468e045f 100644 --- a/src/material.rs +++ b/src/material.rs @@ -2,7 +2,6 @@ use super::*; use rand::RngExt; use std::sync::LazyLock; use crate::math::triangular_index; -const Z_MAX: usize = 88; ///This helper function is a workaround to issue #368 in serde fn default_surface_binding_model() -> SurfaceBindingModel { @@ -306,10 +305,6 @@ 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() } -//https://math.stackexchange.com/questions/2388887/ -//num elements in a triangular NxN matrix (including diag) -const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2; - // Generating lookup tables for all possibilities turns out to be faster than calculating on the fly // Tables for Za, Zb are upper-triangular static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( From 43d33405106e87e270c15fe810f6ac21a0b1fb63 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 12:26:37 -0700 Subject: [PATCH 03/20] Actually, can't make LS stopping power constant triangular - not symmetric in Za Zb. --- src/material.rs | 46 +++++++++++++++++++++++++--------------------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/src/material.rs b/src/material.rs index 468e045f..d9de74c9 100644 --- a/src/material.rs +++ b/src/material.rs @@ -3,6 +3,7 @@ use rand::RngExt; use std::sync::LazyLock; use crate::math::triangular_index; + ///This helper function is a workaround to issue #368 in serde fn default_surface_binding_model() -> SurfaceBindingModel { SurfaceBindingModel::TARGET @@ -301,34 +302,37 @@ impl Material { } } -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() -} - -// Generating lookup tables for all possibilities turns out to be faster than calculating on the fly -// Tables for Za, Zb are upper-triangular -static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( +static LS_STOPPING_CONSTANT_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( || - { - let mut array = [0.0; TABLE_SIZE]; - for i in 0..Z_MAX { - for j in 0..=i { - let index = (i * (i + 1))/2 + j; - array[index] = lindhard_scharff_stopping_power_constant(i as f64, j as f64); - } + 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) } - array - } + ) ); +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() +} #[inline] pub fn lindhard_scharff_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f64) -> f64 { - - let mut i = Za as usize; - let mut j = Zb as usize; - - LS_STOPPING_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]*(E/Ma).sqrt() + 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 + } + ) +); static BV_EMPIRICAL_MEAN_IONIZATON_POT: LazyLock<[f64; Z_MAX]> = LazyLock::new( || From 801e30f06daeb0ea13700bc62bd0f939d5cc1e5c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 12:27:39 -0700 Subject: [PATCH 04/20] Remove accidentally pasted BV MIP table. --- src/material.rs | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/src/material.rs b/src/material.rs index d9de74c9..723eb1bf 100644 --- a/src/material.rs +++ b/src/material.rs @@ -320,21 +320,8 @@ 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 { 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 - } - ) -); -static BV_EMPIRICAL_MEAN_IONIZATON_POT: LazyLock<[f64; Z_MAX]> = LazyLock::new( +static BV_EMPIRICAL_MEAN_IONIZATON_POT_TABLE: LazyLock<[f64; Z_MAX]> = LazyLock::new( || std::array::from_fn( |Zb| { @@ -352,7 +339,7 @@ pub fn bethe_bloch_stopping_power_cross_section(Za: f64, Zb: f64, E: f64, Ma: f6 let beta = (1. - 1./(1. + E/Ma/C.powi(2)).powi(2)).sqrt(); let v = beta*C; - let I = BV_EMPIRICAL_MEAN_IONIZATON_POT[Zb as usize]; + let I = BV_EMPIRICAL_MEAN_IONIZATON_POT_TABLE[Zb as usize]; //See Biersack and Haggmark - this looks like an empirical shell correction let B = if Zb < 3. { From 51b16f7f363a5649fab5868816496cd8c9472af6 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 12:34:10 -0700 Subject: [PATCH 05/20] removed triangular LUT --- src/material.rs | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/material.rs b/src/material.rs index 723eb1bf..5fac2100 100644 --- a/src/material.rs +++ b/src/material.rs @@ -1,8 +1,6 @@ use super::*; use rand::RngExt; use std::sync::LazyLock; -use crate::math::triangular_index; - ///This helper function is a workaround to issue #368 in serde fn default_surface_binding_model() -> SurfaceBindingModel { From 51677cef4648e821dffd2bb3d2ce1c96e71cf6f9 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 12:42:27 -0700 Subject: [PATCH 06/20] Added cfg to swap between 'reasonable Z and full periodic table --- Cargo.toml | 1 + src/consts.rs | 3 +++ 2 files changed, 4 insertions(+) diff --git a/Cargo.toml b/Cargo.toml index d5a0a2f5..6b2c3744 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -46,3 +46,4 @@ distributions = ["ndarray"] no_list_output = [] parry3d = ["parry3d-f64"] python = ["pyo3"] +extended_max_z = [] \ No newline at end of file diff --git a/src/consts.rs b/src/consts.rs index b8383232..c260084a 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -36,6 +36,9 @@ pub const LINDHARD_SCHARFF_PREFACTOR: f64 = 1.212*ANGSTROM*ANGSTROM*Q; /// Lindhard reduced energy prefactor, in SI units. pub const LINDHARD_REDUCED_ENERGY_PREFACTOR: f64 = 4.*PI*EPS0/Q/Q; /// Maximum atomic number in RustBCA. +#[cfg(not(feature = "extended_max_z"))] +pub const Z_MAX: usize = 92; +#[cfg(feature(="extended_max_z"))] pub const Z_MAX: usize = 120; //https://math.stackexchange.com/questions/2388887/ //num elements in a triangular NxN matrix (including diag) From 251dbb6078449266d89d6cc17034e7c360e31b4d Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 14:12:55 -0700 Subject: [PATCH 07/20] Fix to cfg statement re Z_MAX --- src/consts.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/consts.rs b/src/consts.rs index c260084a..dc87f273 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -38,7 +38,7 @@ pub const LINDHARD_REDUCED_ENERGY_PREFACTOR: f64 = 4.*PI*EPS0/Q/Q; /// Maximum atomic number in RustBCA. #[cfg(not(feature = "extended_max_z"))] pub const Z_MAX: usize = 92; -#[cfg(feature(="extended_max_z"))] +#[cfg(feature="extended_max_z")] pub const Z_MAX: usize = 120; //https://math.stackexchange.com/questions/2388887/ //num elements in a triangular NxN matrix (including diag) From b5629990d78812260aded01e9cef90f926c2c41a Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 2 Aug 2026 16:47:15 -0700 Subject: [PATCH 08/20] Add tests for Z1 > Z2 to test_rustbca.py --- examples/test_rustbca.py | 28 +++++++++++++++++++++++++++- src/interactions.rs | 6 ++++++ 2 files changed, 33 insertions(+), 1 deletion(-) diff --git a/examples/test_rustbca.py b/examples/test_rustbca.py index 5b748d8b..e374135f 100644 --- a/examples/test_rustbca.py +++ b/examples/test_rustbca.py @@ -110,7 +110,6 @@ def main(): print(f'Sputtering yield for {ion["symbol"]} on {target["symbol"]} at {energy} eV is {Y} at/ion. Yamamura predicts { np.round(yamamura(ion, target, energy),3)} at/ion.') np.testing.assert_approx_equal(Y, 0.044) - R_N, R_E = reflection_coefficient(ion, target, energy, angle, num_samples) print(f'Particle reflection coefficient for {ion["symbol"]} on {target["symbol"]} at {energy} eV is {R_N}. Thomas predicts {np.round(thomas_reflection(ion, target, energy), 3)}.') print(f'Energy reflection coefficient for {ion["symbol"]} on {target["symbol"]} at {energy} eV is {R_E}') @@ -123,6 +122,33 @@ def main(): np.testing.assert_approx_equal(R_N, 0.424) np.testing.assert_approx_equal(R_E, 0.22840032456593984) + # test of triangular LUTs correctly handling Za > Zb + ion = neon + target = boron + angle = 60.0 + num_samples = 10000 + energy = 2500.0 + + Y = sputtering_yield(ion, target, energy, angle, num_samples) + R_N, R_E = reflection_coefficient(ion, target, energy, angle, num_samples) + + np.testing.assert_approx_equal(Y, 3.3481) + np.testing.assert_approx_equal(R_N, 0.0878) + np.testing.assert_approx_equal(R_E, 0.013734709021659743) + + ion = copper # testing with Es > 0 + Y = sputtering_yield(ion, target, energy, angle, num_samples) + R_N, R_E = reflection_coefficient(ion, target, energy, angle, num_samples) + + np.testing.assert_approx_equal(Y, 4.9431) + np.testing.assert_approx_equal(R_N, 0.0066) + np.testing.assert_approx_equal(R_E, 0.000298720196409247) + + # reset species + ion = helium + ion['Eb'] = 0.0 + target = tungsten + vx0 = 1e5 vy0 = 1e5 vz0 = 0.0 diff --git a/src/interactions.rs b/src/interactions.rs index 13c88b62..fec1a0b6 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -239,6 +239,8 @@ static COULOMB_CONSTANT_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( let mut array = [0.0; TABLE_SIZE]; for i in 0..Z_MAX { for j in 0..=i { + //going from 1D to linear triangular upper array + //https://stackoverflow.com/questions/27086195 let index = (i * (i + 1))/2 + j; array[index] = coulomb_constant(i as f64, j as f64); } @@ -313,6 +315,8 @@ static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock:: let mut array = [0.0; TABLE_SIZE]; for i in 0..Z_MAX { for j in 0..=i { + //going from 1D to linear triangular upper array + //https://stackoverflow.com/questions/27086195 let index = (i * (i + 1))/2 + j; array[index] = lindhard_screening_length(i as f64, j as f64); } @@ -338,6 +342,8 @@ static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( let mut array = [0.0; TABLE_SIZE]; for i in 0..Z_MAX { for j in 0..=i { + //going from 1D to linear triangular upper array + //https://stackoverflow.com/questions/27086195 let index = (i * (i + 1))/2 + j; array[index] = zbl_screening_length(i as f64, j as f64); } From 31f0beb0b70f6f1924db688bddec237142500827 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 10:32:22 -0700 Subject: [PATCH 09/20] Add Python function to compute screened coulomb scattering integrals and test. --- examples/test_scattering_integrals.py | 46 +++++++++++++++++++++++++++ src/lib.rs | 41 ++++++++++++++++++++++++ 2 files changed, 87 insertions(+) create mode 100644 examples/test_scattering_integrals.py diff --git a/examples/test_scattering_integrals.py b/examples/test_scattering_integrals.py new file mode 100644 index 00000000..7e5d9c4e --- /dev/null +++ b/examples/test_scattering_integrals.py @@ -0,0 +1,46 @@ +import numpy as np +import matplotlib.pyplot as plt +import sys, os + +from libRustBCA import * +#This should allow the script to find materials and formulas from anywhere +sys.path.append(os.path.dirname(__file__)+'/../scripts') +sys.path.append('scripts') +from materials import * +from formulas import * + +energies = np.logspace(0, 4, 4) +impact_parameters = np.logspace(-3, 3, 100) + +ion = helium +target = boron + +Za = ion['Z'] +Zb = target['Z'] +Ma = ion['m'] +Mb = target['m'] + +show_plots = True + +linestyles = ['-', '--', ':', '-.'] + +for linestyle, energy in zip(linestyles, energies): + gm = np.zeros(100) + gl = np.zeros(100) + mw = np.zeros(100) + magic = np.zeros(100) + for index, p in enumerate(impact_parameters): + gm[index], gl[index], mw[index], magic[index] = scattering_integrals(Za, Zb, Ma, Mb, energy, p) + plt.semilogx(impact_parameters, gm, label=f'Gauss-Mehler, E={np.round(energy/1000)} keV', linestyle=linestyle) + plt.semilogx(impact_parameters, gl, label=f'Gauss-Legendre, E={np.round(energy/1000)} keV', linestyle=linestyle) + plt.semilogx(impact_parameters, mw, label=f'Mendenhall-Weller, E={np.round(energy/1000)} keV', linestyle=linestyle) + plt.semilogx(impact_parameters, magic, label=f'MAGIC, E={np.round(energy/1000)} keV', linestyle=linestyle) + plt.gca().set_prop_cycle(None) + + np.testing.assert_allclose(gm, gl, atol=5e-3) # 0.5% seems reasonable? max is ~0.3% + np.testing.assert_allclose(gm, mw, atol=5e-3) + np.testing.assert_allclose(mw, gl, atol=5e-3) + plt.legend() + plt.xlabel('p [A]') + plt.ylabel('theta [rad]') +if show_plots: plt.show() \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index bc03c53f..72276330 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -143,6 +143,9 @@ mod libRustBCA { #[pymodule_export] use super::electronic_stopping_cross_sections; + + #[pymodule_export] + use super::scattering_integrals; } #[derive(Debug)] @@ -2164,4 +2167,42 @@ fn moller_knuth_two_sum(a: f64, b: f64) -> (f64, f64) { let delta_a = a - a_prime; let r = delta_a + delta_b; (s, r) +} + +#[pyfunction] +fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64) -> (f64, f64, f64, f64) { + let E0 = E0*EV; + let p = p*ANGSTROM; + + let options = Options { + name: "test".to_string(), + track_trajectories: false, + track_recoils: true, + track_recoil_trajectories: false, + write_buffer_size: 8000, + weak_collision_order: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::KR_C]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-14}]], + track_displacements: false, + track_energy_losses: false, + seed: 0, + }; + + let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 1000, 1E-12).unwrap(); + + //Compute center of mass deflection angle with each algorithm + let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C, 100); + let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); + let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); + let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); + + (theta_gm, theta_gl, theta_mw, theta_magic) } \ No newline at end of file From 3533c75591173816894507bf83ff872c3fa5d16c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 10:37:07 -0700 Subject: [PATCH 10/20] Missing cfg flag for scat int func --- src/lib.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib.rs b/src/lib.rs index 72276330..67c50d20 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -2168,7 +2168,7 @@ fn moller_knuth_two_sum(a: f64, b: f64) -> (f64, f64) { let r = delta_a + delta_b; (s, r) } - +#[cfg(feature = "python")] #[pyfunction] fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64) -> (f64, f64, f64, f64) { let E0 = E0*EV; From 149a836dc12dc964e92e3bf052d101d094450219 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 10:46:57 -0700 Subject: [PATCH 11/20] cleanup of lib.rs --- src/lib.rs | 34 +++------------------------------- 1 file changed, 3 insertions(+), 31 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 67c50d20..9d24cae9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -52,13 +52,7 @@ use std::f64::consts::SQRT_2; #[cfg(feature = "python")] use pyo3::prelude::*; #[cfg(feature = "python")] -use pyo3::wrap_pyfunction; -#[cfg(feature = "python")] use pyo3::types::*; -#[cfg(feature = "python")] -use pyo3::exceptions::PyTypeError; -#[cfg(feature = "python")] -use pyo3::*; //Load internal modules pub mod material; @@ -96,7 +90,6 @@ pub use parry3d_f64::na::{Point3, Vector3, Matrix3, Rotation3}; #[cfg(feature = "python")] #[pymodule] mod libRustBCA { - use pyo3::prelude::*; #[pymodule_export] use super::simple_bca_py; @@ -2170,36 +2163,15 @@ fn moller_knuth_two_sum(a: f64, b: f64) -> (f64, f64) { } #[cfg(feature = "python")] #[pyfunction] -fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64) -> (f64, f64, f64, f64) { +#[pyo3(signature = (Za, Zb, Ma, Mb, E0, p, n_gl_points=100))] +fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_gl_points: usize) -> (f64, f64, f64, f64) { let E0 = E0*EV; let p = p*ANGSTROM; - let options = Options { - name: "test".to_string(), - track_trajectories: false, - track_recoils: true, - track_recoil_trajectories: false, - write_buffer_size: 8000, - weak_collision_order: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: false, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![InteractionPotential::KR_C]], - scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-14}]], - track_displacements: false, - track_energy_losses: false, - seed: 0, - }; - let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 1000, 1E-12).unwrap(); //Compute center of mass deflection angle with each algorithm - let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C, 100); + let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C, n_gl_points); let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); From 53c2b58c69c207d9c9ff0c96150da5870fe5d677 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 10:47:47 -0700 Subject: [PATCH 12/20] Added scattering integral tests to workflow --- .github/workflows/rustbca_compile_check.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index fed9b85e..167f1200 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -42,6 +42,7 @@ jobs: python3 -m pip install . python3 -c "from libRustBCA import *;" python3 examples/test_rustbca.py + python3 examples/test_scattering_integrals.py python3 examples/test_cube.py python3 examples/make_input_file_and_run.py python3 examples/test_electronic_stopping.py From 7c108a94b03e554f407623f708c4e13ca0154dd6 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 11:03:08 -0700 Subject: [PATCH 13/20] fixed accidental overwrite of triangular lut --- src/interactions.rs | 78 +++++++++++++++++++-------------------------- 1 file changed, 32 insertions(+), 46 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index fec1a0b6..2720a4b6 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -1,6 +1,6 @@ use super::*; use std::sync::LazyLock; -use crate::math::triangular_index; +use math::triangular_index; /// Analytic solutions to outermost root of the interaction potential. pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { @@ -239,8 +239,6 @@ static COULOMB_CONSTANT_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( let mut array = [0.0; TABLE_SIZE]; for i in 0..Z_MAX { for j in 0..=i { - //going from 1D to linear triangular upper array - //https://stackoverflow.com/questions/27086195 let index = (i * (i + 1))/2 + j; array[index] = coulomb_constant(i as f64, j as f64); } @@ -257,14 +255,13 @@ fn coulomb_constant(Za: f64, Zb: f64) -> f64 { pub fn screened_coulomb(r: f64, a: f64, Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 { let mut i = Za as usize; let mut j = Zb as usize; - COULOMB_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]/r*phi(r/a, interaction_potential) + let index = triangular_index(&mut i, &mut j); + COULOMB_CONSTANT_TABLE[index]/r*phi(r/a, interaction_potential) } /// Coulombic interaction potential. pub fn coulomb(r: f64, Za: f64, Zb: f64) -> f64 { - let mut i = Za as usize; - let mut j = Zb as usize; - COULOMB_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]/r + COULOMB_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]/r } /// Screening functions for screened-coulomb interaction potentials. @@ -306,61 +303,50 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote } } -// It turns out it's faster to just generate every possible screening length as a lookup table +// 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 -static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( +static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( || - { - let mut array = [0.0; TABLE_SIZE]; - for i in 0..Z_MAX { - for j in 0..=i { - //going from 1D to linear triangular upper array - //https://stackoverflow.com/questions/27086195 - let index = (i * (i + 1))/2 + j; - array[index] = lindhard_screening_length(i as f64, j as f64); - } + 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) } - array - } + ) ); -pub fn lindhard_screening_length(Za: f64, Zb: f64) -> f64 { - 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.) -} - -#[inline] -pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { - let mut i = Za as usize; - let mut j = Zb as usize; - LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] -} - -static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( +static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( || - { - let mut array = [0.0; TABLE_SIZE]; - for i in 0..Z_MAX { - for j in 0..=i { - //going from 1D to linear triangular upper array - //https://stackoverflow.com/questions/27086195 - let index = (i * (i + 1))/2 + j; - array[index] = zbl_screening_length(i as f64, j as f64); - } + std::array::from_fn( + |i| { + let Za = i / Z_MAX; + let Zb = i % Z_MAX; + zbl_screening_length(Za as f64, Zb as f64) } - array - } + ) ); 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.) +} + +#[inline] +pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { + LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] +} + #[inline] pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ - let mut i = Za as usize; - let mut j = Zb as usize; - ZBL_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] + ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] } /// Coefficients of inverse-polynomial interaction potentials. From 315dc7014f9af8a32efe3351dffbfb4389c8d53e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 11:28:25 -0700 Subject: [PATCH 14/20] Add testing of various input options. --- examples/test_different_options.py | 434 +++++++++++++++++++++++++++++ 1 file changed, 434 insertions(+) create mode 100644 examples/test_different_options.py diff --git a/examples/test_different_options.py b/examples/test_different_options.py new file mode 100644 index 00000000..6762ed0d --- /dev/null +++ b/examples/test_different_options.py @@ -0,0 +1,434 @@ +from libRustBCA import * +import numpy as np +import matplotlib.pyplot as plt +import sys +import os +#This should allow the script to find materials and formulas from anywhere +sys.path.append(os.path.dirname(__file__)+'/../scripts') +sys.path.append('scripts') +import time +from tomlkit import parse, dumps + +''' +This script is a first draft of a comprehensive RustBCA input +file creation script using tomlkit. + +It includes two geometry modes, 1D and 0D. + +It simulates the following situations: + +if mode == '1D': + H+ (1 keV) + | + V +__________ +| | +| B | dx = 100 A +_________| +| | +| TiB2 | dx = 100 A +_________| +| | +| Ti | dx = 1000 A + + +if mode == '0D': + H+ (1 keV) + | + V +__________ +| | +| TiB2 | +| | +| | + +And calculates implantation profiles, reflection coefficients, +and sputtering yields. It also uses the ergonomic python functions +to compare the result of using the default values for H on B with +the custom values of this input file. + +It creates an input file as a nested dictionary which is written to +a TOML file using tomlkit. + +It runs the input file with cargo run --release and reads the output files. +''' + +def run_test( + interaction_potential='KR_C', + high_energy_free_flight_path=False, + electronic_stopping_mode='LOW_ENERGY_NONLOCAL', + mfp='LIQUID', + weak_collision_order=0, + num_threads=6, + index=0, + scattering_integral={'GAUSS_MEHLER': {'n_points': 6}}, + run_sim=True, + num_samples=100000 + ): + mode = '1D' + incident_energy = 1000.0 # eV + angle = 45.0 # degrees; measured from surface normal + + ''' + For organizational purposes, species are commonly defined in dictionaries. + Additional examples can be found in scripts/materials.py, but values + should be checked for correctness before use. Values are explained + in the relevant sections below. + ''' + hydrogen = { + 'symbol': 'H', + 'name': 'hydrogen', + 'Z': 1.0, + 'm': 1.008, # AMU + 'Ec': 0.95, # eV + 'Es': 1.5, # eV + } + + titanium = { + 'symbol': 'Ti', + 'name': 'titanium', + 'Z': 22.0, + 'm': 47.867, # AMU + 'Es': 4.84, # eV + 'Ec': 3.5, # eV + 'Eb': 0., # eV + 'Ed': 19.0, # eV + 'n': 5.67e28, # 1/m^3 + } + + boron = { + 'symbol': 'B', + 'name': 'boron', + 'Z': 5.0, + 'm': 10.811, # AMU + 'n': 1.309E29, # 1/m^3 + 'Es': 5.77, # eV + 'Eb': 0., # eV + 'Ec': 5., # eV + 'Ed': 25.0 # eV + } + + # species definitions + ion = hydrogen + target1 = boron + target2 = titanium + + # geometry definitions + n_i = 0.0328 # 1 / A^3 from n_i = rho_TiB2 / (mB * 2 + mTi) + layer_thicknesses = [100.0, 100.0, 1000.0] # A + layer_1_densities = [boron["n"]/10**30, 0.0] # 1/A^3 + layer_2_densities = [n_i * 2, n_i] # 1/A^3 + layer_3_densities = [0.0, titanium["n"]/10**30] # 1/A^3 + + options = { + 'name': f'input_file_{index}', + 'track_trajectories': False, # whether to track trajectories for plotting; memory intensive + 'track_recoils': True, # whether to track recoils; must enable for sputtering + 'track_recoil_trajectories': False, # whether to track recoil trajectories for plotting + 'track_displacements': False, # whether to track collisions with T > Ed for each species + 'track_energy_losses': False, # whether to track detailed collision energies; memory intensive + 'write_buffer_size': 8192, # how big the buffer is for file writing + 'weak_collision_order': weak_collision_order, # weak collisions at radii (k + 1)*r; enable only when required + 'suppress_deep_recoils': False, # suppress recoils too deep to ever sputter + 'high_energy_free_flight_paths': high_energy_free_flight_path, # SRIM-style high energy free flight distances; use with caution + 'num_threads': num_threads, # number of threads to run in parallel + 'num_chunks': 10, # code will write to file every nth chunk; for very large simulations, increase num_chunks + 'electronic_stopping_mode': electronic_stopping_mode, + 'mean_free_path_model': mfp, # liquid is amorphous (constant mean free path); gas is exponentially-distributed mean free paths + 'interaction_potential': [[interaction_potential]], + 'scattering_integral': [ + [ + scattering_integral + ] + ], + + 'root_finder': [ + [ + { + 'NEWTON': { + 'max_iterations': 100, + 'tolerance': 1e-6 + } + } + ] + ], + } + + # material parameters are per-species + material_parameters = { + 'energy_unit': 'EV', + 'mass_unit': 'AMU', + # bulk binding energy; typically zero as a model choice + 'Eb': [ + target1["Eb"], + target2["Eb"] + ], + # surface binding energy + 'Es': [ + target1["Es"], + target2["Es"] + ], + # cutoff energy - particles with E < Ec stop + 'Ec': [ + target1["Ec"], + target2["Ec"] + ], + # displacement energy - only used to track displacements + 'Ed': [ + target1["Ed"], + target2["Ed"] + ], + # atomic number + 'Z': [ + target1["Z"], + target2["Z"] + ], + # atomic mass + 'm': [ + target1["m"], + target2["m"] + ], + # used to pick interaction potential from matrix in [options] + 'interaction_index': [0, 0], + 'surface_binding_model': { + "PLANAR": {'calculation': "INDIVIDUAL"} + }, + 'bulk_binding_model': 'INDIVIDUAL' + } + + particle_parameters = { + 'length_unit': 'ANGSTROM', + 'energy_unit': 'EV', + 'mass_unit': 'AMU', + # number of computational ions of this species to run at this energy + 'N': [num_samples], + # atomic mass + 'm': [ion["m"]], + # atomic number + 'Z': [ion["Z"]], + # incidenet energy + 'E': [incident_energy], + # cutoff energy - if E < Ec, particle stops + 'Ec': [ion["Ec"]], + # surface binding energy + 'Es': [ion["Es"]], + # initial position - if Es significant and E low, start (n)^(-1/3) above surface + # otherwise 0, 0, 0 is fine; most geometry modes have surface at x=0 with target x>0 + 'pos': [[0.0, 0.0, 0.0]], + # initial direction unit vector; most geometry modes have x-axis into the surface + 'dir': [ + [ + np.cos(angle*np.pi/180.0), + np.sin(angle*np.pi/180.0), + 0.0 + ] + ], + } + + geometry_0D = { + 'length_unit': 'ANGSTROM', + # used to correct nonlocal stopping for known compound discrpancies + 'electronic_stopping_correction_factor': 1.0, + # number densities of each species + 'densities': [2 * n_i, n_i] + } + + geometry_1D = { + 'length_unit': 'ANGSTROM', + # used to correct nonlocal stopping for known compound discrpancies + 'electronic_stopping_correction_factors': [1.0, 1.0, 1.0], + # thickness of each layer in order from top (x=0) to bottom + 'layer_thicknesses': layer_thicknesses, + # number densitiy of each layer in order from top to bottom + 'densities': [ + layer_1_densities, + layer_2_densities, + layer_3_densities, + ] + } + + if mode == '1D': + input_data = { + 'options': options, + 'material_parameters': material_parameters, + 'particle_parameters': particle_parameters, + 'geometry_input': geometry_1D + } + elif mode == '0D': + input_data = { + 'options': options, + 'material_parameters': material_parameters, + 'particle_parameters': particle_parameters, + 'geometry_input': geometry_0D + } + + # Attempt to cleanup line endings + input_string = dumps(input_data).replace('\r', '') + with open(f'examples/input_file_{index}.toml', 'w') as input_file: + input_file.write(input_string) + + if run_sim: + os.system(f'cargo run --release {mode} examples/input_file_{index}.toml') + + # Read output files - ensure arrays are at least 2D for indexing + sputtered = np.atleast_2d(np.genfromtxt(f'input_file_{index}sputtered.output', delimiter=',')) + reflected = np.atleast_2d(np.genfromtxt(f'input_file_{index}reflected.output', delimiter=',')) + implanted = np.atleast_2d(np.genfromtxt(f'input_file_{index}deposited.output', delimiter=',')) + + return sputtered, reflected, implanted + +num_bins = 75 +num_samples = 100000 +run_sim = True +show_plots = True + +# interaction potentials +interaction_potentials = ['KR_C', 'ZBL', 'MOLIERE', 'LENZ_JENSEN'] + +sim_index = 0 + +if not run_sim: + sim_times = np.genfromtxt('sim_times.txt') +else: + sim_times = [] + +Y_test = np.array([ + 0.02593, 0.02104, 0.02467, 0.02731, 0.02593, 0.0218, 0.02211, 0.02288, + 0.03577, 0.02593, 0.02587, 0.0296, 0.0255, 0.02568, 0.02587, 0.02593, + 0.02905, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593, + 0.02593 +]) + +R_N_test = np.array([ + 0.17442, 0.17032, 0.16802, 0.18216, 0.17442, 0.17437, 0.17449, 0.17777, + 0.28801, 0.17442, 0.17452, 0.22074, 0.17506, 0.17446, 0.17445, 0.17442, + 0.17722, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442, + 0.17442 +]) + +Y = np.zeros_like(Y_test) +R_N = np.zeros_like(R_N_test) + +for interaction_potential in interaction_potentials: + start = time.time() + s, r, i = run_test(interaction_potential=interaction_potential, index=sim_index, num_samples=num_samples, run_sim=run_sim) + stop = time.time() + Y[sim_index] = np.shape(s)[0]/num_samples + R_N[sim_index] = np.shape(r)[0]/num_samples + sim_time = (stop - start)/1e-3 + if run_sim: sim_times.append(sim_time) + x = i[:, 2] + plt.figure(1) + plt.title('Interaction Potentials') + plt.hist(x, bins=num_bins, histtype='step', label=f'{interaction_potential} [{np.round(sim_time)} ms]') + plt.legend() + plt.xlabel('x [A]') + plt.ylabel(f'f(x) [counts]') + print(f'{sim_index} Y: {np.shape(s)[0]/num_samples} R_N: {np.shape(r)[0]/num_samples}') + sim_index += 1 + +for weak_collision_order in [0, 1, 2, 3]: + + start = time.time() + s, r, i = run_test(weak_collision_order=weak_collision_order, index=sim_index, num_samples=num_samples, run_sim=run_sim) + stop = time.time() + Y[sim_index] = np.shape(s)[0]/num_samples + R_N[sim_index] = np.shape(r)[0]/num_samples + sim_time = (stop - start)/1e-3 + if run_sim: sim_times.append(sim_time) + x = i[:, 2] + plt.figure(2) + plt.title('Weak Collision Orders') + plt.hist(x, bins=num_bins, histtype='step', label=f'k={weak_collision_order} [{np.round(sim_time)} ms]') + plt.legend() + plt.xlabel('x [A]') + plt.ylabel(f'f(x) [counts]') + print(f'{sim_index} Y: {np.shape(s)[0]/num_samples} R_N: {np.shape(r)[0]/num_samples}') + sim_index += 1 + +for electronic_stopping_mode in ['LOW_ENERGY_LOCAL', 'LOW_ENERGY_NONLOCAL', 'INTERPOLATED', 'LOW_ENERGY_EQUIPARTITION']: + start = time.time() + s, r, i = run_test(electronic_stopping_mode=electronic_stopping_mode, index=sim_index, num_samples=num_samples, run_sim=run_sim) + stop = time.time() + Y[sim_index] = np.shape(s)[0]/num_samples + R_N[sim_index] = np.shape(r)[0]/num_samples + sim_time = (stop - start)/1e-3 + if run_sim: sim_times.append(sim_time) + x = i[:, 2] + plt.figure(3) + plt.title('Electronic Stopping Modes') + plt.hist(x, bins=num_bins, histtype='step', label=f'{electronic_stopping_mode} [{np.round(sim_time)} ms]') + plt.legend() + plt.xlabel('x [A]') + plt.ylabel(f'f(x) [counts]') + print(f'{sim_index} Y: {np.shape(s)[0]/num_samples} R_N: {np.shape(r)[0]/num_samples}') + sim_index += 1 + +for scattering_integral in ['MENDENHALL_WELLER', {'GAUSS_MEHLER': {'n_points': 6}}, 'GAUSS_LEGENDRE']: + start = time.time() + s, r, i = run_test(scattering_integral=scattering_integral, index=sim_index, num_samples=num_samples, run_sim=run_sim) + stop = time.time() + Y[sim_index] = np.shape(s)[0]/num_samples + R_N[sim_index] = np.shape(r)[0]/num_samples + sim_time = (stop - start)/1e-3 + if run_sim: sim_times.append(sim_time) + x = i[:, 2] + plt.figure(4) + plt.title('Scattering Integrals') + plt.hist(x, bins=num_bins, histtype='step', label=f'{scattering_integral} [{np.round(sim_time)} ms]') + plt.legend() + plt.xlabel('x [A]') + plt.ylabel(f'f(x) [counts]') + print(f'{sim_index} Y: {np.shape(s)[0]/num_samples} R_N: {np.shape(r)[0]/num_samples}') + sim_index += 1 + +for mfp in ['LIQUID', 'GASEOUS']: + start = time.time() + s, r, i = run_test(mfp=mfp, index=sim_index, num_samples=num_samples, run_sim=run_sim) + stop = time.time() + Y[sim_index] = np.shape(s)[0]/num_samples + R_N[sim_index] = np.shape(r)[0]/num_samples + sim_time = (stop - start)/1e-3 + if run_sim: sim_times.append(sim_time) + x = i[:, 2] + plt.figure(5) + plt.title('MFP Distribution') + plt.hist(x, bins=num_bins, histtype='step', label=f'{mfp} [{np.round(sim_time)} ms]') + plt.legend() + plt.xlabel('x [A]') + plt.ylabel(f'f(x) [counts]') + print(f'{sim_index} Y: {np.shape(s)[0]/num_samples} R_N: {np.shape(r)[0]/num_samples}') + sim_index += 1 + +num_threads = [1, 2, 3, 4, 5, 6, 7, 8] +sim_index_threads_start = sim_index +for n in num_threads: + start = time.time() + s, r, i = run_test(num_threads=n, index=sim_index, num_samples=num_samples, run_sim=run_sim) + stop = time.time() + Y[sim_index] = np.shape(s)[0]/num_samples + R_N[sim_index] = np.shape(r)[0]/num_samples + sim_time = (stop - start)/1e-3 + if run_sim: sim_times.append(sim_time) + print(f'{sim_index} Y: {np.shape(s)[0]/num_samples} R_N: {np.shape(r)[0]/num_samples}') + sim_index += 1 + +sim_index_threads_stop = sim_index +plt.figure(6) +plt.plot(num_threads, sim_times[sim_index_threads_start]/np.array(sim_times[sim_index_threads_start:sim_index_threads_stop]), label="Amdahl's law; s=0.06") +s = 0.06 +p = 1 - s +plt.plot(num_threads, 1/(s + p/np.array(num_threads)), label='RustBCA (i5-8600k, 6 cores)') +plt.xlabel('n threads') +plt.ylabel('t [ms]') +plt.legend() +plt.plot([6, 6], [0, 10], linestyle='--', color='gray') +plt.gca().set_ylim([1, 5]) + +if run_sim: np.savetxt('sim_times.txt', sim_times) +if show_plots: plt.show() + +np.testing.assert_allclose(Y_test, Y) +np.testing.assert_allclose(R_N_test, R_N) + + From df8785f4310f7294070996cef56bee609361eade Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 11:29:25 -0700 Subject: [PATCH 15/20] Add different option test script to tests --- .github/workflows/rustbca_compile_check.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 167f1200..a50af72a 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -46,6 +46,7 @@ jobs: python3 examples/test_cube.py python3 examples/make_input_file_and_run.py python3 examples/test_electronic_stopping.py + python3 examples/test_different_options.py - name: Test Fortran and C bindings run : | cargo build --release --lib --features parry3d From e16189c4cb0dceb68b5b5f054307e3e11230691e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 11:49:24 -0700 Subject: [PATCH 16/20] Fix test data in different options test --- examples/test_different_options.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/test_different_options.py b/examples/test_different_options.py index 6762ed0d..2f9ce812 100644 --- a/examples/test_different_options.py +++ b/examples/test_different_options.py @@ -279,7 +279,7 @@ def run_test( num_bins = 75 num_samples = 100000 -run_sim = True +run_sim = False show_plots = True # interaction potentials @@ -294,14 +294,14 @@ def run_test( Y_test = np.array([ 0.02593, 0.02104, 0.02467, 0.02731, 0.02593, 0.0218, 0.02211, 0.02288, - 0.03577, 0.02593, 0.02587, 0.0296, 0.0255, 0.02568, 0.02587, 0.02593, + 0.03577, 0.02593, 0.02587, 0.0296, 0.0255, 0.02593, 0.02587, 0.02593, 0.02905, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593, 0.02593 ]) R_N_test = np.array([ 0.17442, 0.17032, 0.16802, 0.18216, 0.17442, 0.17437, 0.17449, 0.17777, - 0.28801, 0.17442, 0.17452, 0.22074, 0.17506, 0.17446, 0.17445, 0.17442, + 0.28801, 0.17442, 0.17452, 0.22074, 0.17506, 0.17442, 0.17445, 0.17442, 0.17722, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442, 0.17442 ]) @@ -415,10 +415,10 @@ def run_test( sim_index_threads_stop = sim_index plt.figure(6) -plt.plot(num_threads, sim_times[sim_index_threads_start]/np.array(sim_times[sim_index_threads_start:sim_index_threads_stop]), label="Amdahl's law; s=0.06") +plt.plot(num_threads, sim_times[sim_index_threads_start]/np.array(sim_times[sim_index_threads_start:sim_index_threads_stop]), label='RustBCA (i5-8600k, 6 cores)') s = 0.06 p = 1 - s -plt.plot(num_threads, 1/(s + p/np.array(num_threads)), label='RustBCA (i5-8600k, 6 cores)') +plt.plot(num_threads, 1/(s + p/np.array(num_threads)), label="Amdahl's law; s=0.06" ) plt.xlabel('n threads') plt.ylabel('t [ms]') plt.legend() From 600b1dc03ab2a868651da9f00b7788036bd1c4d5 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 11:59:44 -0700 Subject: [PATCH 17/20] Fix run_sim flag to True --- examples/test_different_options.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/test_different_options.py b/examples/test_different_options.py index 2f9ce812..587d6ce1 100644 --- a/examples/test_different_options.py +++ b/examples/test_different_options.py @@ -279,7 +279,7 @@ def run_test( num_bins = 75 num_samples = 100000 -run_sim = False +run_sim = True show_plots = True # interaction potentials From 8b2e34cc83ee61faaebd70d0edb57cdc2746e4da Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 15:19:28 -0700 Subject: [PATCH 18/20] Undo trilut --- src/interactions.rs | 51 +++++++++++++++++++++++++++------------------ 1 file changed, 31 insertions(+), 20 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 2720a4b6..50f46d15 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -261,7 +261,9 @@ pub fn screened_coulomb(r: f64, a: f64, Za: f64, Zb: f64, interaction_potential: /// Coulombic interaction potential. pub fn coulomb(r: f64, Za: f64, Zb: f64) -> f64 { - COULOMB_CONSTANT_TABLE[Za as usize * Z_MAX + Zb as usize]/r + let mut i = Za as usize; + let mut j = Zb as usize; + COULOMB_CONSTANT_TABLE[triangular_index(&mut i, &mut j)]/r } /// Screening functions for screened-coulomb interaction potentials. @@ -306,29 +308,32 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote // 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 -static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( +static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = 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) + { + let mut array = [0.0; TABLE_SIZE]; + for i in 0..Z_MAX { + for j in 0..=i { + let index = (i * (i + 1))/2 + j; + array[index] = lindhard_screening_length(i as f64, j as f64); + } } - ) + array + } ); -static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; Z_MAX*Z_MAX]> = LazyLock::new( +static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = 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) + { + let mut array = [0.0; TABLE_SIZE]; + for i in 0..Z_MAX { + for j in 0..=i { + let index = (i * (i + 1))/2 + j; + array[index] = zbl_screening_length(i as f64, j as f64); + } } - ) + array + } ); pub fn zbl_screening_length(Za: f64, Zb: f64) -> f64{ @@ -341,12 +346,18 @@ pub fn lindhard_screening_length(Za: f64, Zb: f64) -> f64 { #[inline] pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { - LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] + let mut i = Za as usize; + let mut j = Zb as usize; + + LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] } #[inline] pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ - ZBL_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] + let mut i = Za as usize; + let mut j = Zb as usize; + + ZBL_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] } /// Coefficients of inverse-polynomial interaction potentials. From 54c8ce642e6ffccd97c690eb2e24ddbcdc4776de Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 15:19:40 -0700 Subject: [PATCH 19/20] Undo trilut --- src/interactions.rs | 33 +++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 50f46d15..1c338017 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -308,6 +308,7 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote // 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 +/* static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( || { @@ -321,6 +322,30 @@ static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock:: array } ); +#[inline] +pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { + let mut i = Za as usize; + let mut j = Zb as usize; + + LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] +} +*/ + +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 f64, Zb as f64) + } + ) +); + +pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { + LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] +} + static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( || @@ -344,14 +369,6 @@ pub fn lindhard_screening_length(Za: f64, Zb: f64) -> f64 { 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.) } -#[inline] -pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { - let mut i = Za as usize; - let mut j = Zb as usize; - - LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] -} - #[inline] pub fn zbl_screening_length_lookup(Za: u64, Zb: u64) -> f64{ let mut i = Za as usize; From 9849bf99e889d45d00ae9e8cfb5f678437c5e7a7 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 3 Aug 2026 15:26:52 -0700 Subject: [PATCH 20/20] Final version of triangular LUTs. --- src/interactions.rs | 19 +------------------ 1 file changed, 1 insertion(+), 18 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 1c338017..da619479 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -308,7 +308,7 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote // 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 -/* + static LINDHARD_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( || { @@ -329,23 +329,6 @@ pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { LINDHARD_SCREENING_LENGTH_TABLE[triangular_index(&mut i, &mut j)] } -*/ - -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 f64, Zb as f64) - } - ) -); - -pub fn lindhard_screening_length_lookup(Za: u64, Zb: u64) -> f64 { - LINDHARD_SCREENING_LENGTH_TABLE[Za as usize * Z_MAX + Zb as usize] -} - static ZBL_SCREENING_LENGTH_TABLE: LazyLock<[f64; TABLE_SIZE]> = LazyLock::new( ||