From 3ce944556a577ab50b99f744272def79fe7aa86e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 11 Aug 2026 13:20:00 -0700 Subject: [PATCH 01/14] Use macros in remaining python functions; removes another ~20 unwrap()s. --- src/lib.rs | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index e300238..4ebb80f 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1196,19 +1196,19 @@ pub fn compound_bca_list_tracked_py<'py>(energies: Vec, ux: Vec, uy: V ///Returns: /// vx, vy, vz (float): final x, y, and z velocity in m/s. When ion implants in material, vx, vy, and vz will all be zero. #[pyfunction] -pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDict>, vx: f64, vy: f64, vz: f64) -> (f64, f64, f64){ +pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, PyDict>, vx: f64, vy: f64, vz: f64) -> PyResult<(f64, f64, f64)> { - let Z1: f64 = ion.get_item("Z").unwrap().expect("Error: Cannot get key 'Z' from ion dict.").extract().unwrap(); - let m1: f64 = ion.get_item("m").unwrap().expect("Error: Cannot get key 'm' from ion dict.").extract().unwrap(); - let Es1: f64 = ion.get_item("Es").unwrap().expect("Error: Cannot get key 'Es' from ion dict.").extract().unwrap(); - let Ec1: f64 = ion.get_item("Ec").unwrap().expect("Error: Cannot get key 'Ec' from ion dict.").extract().unwrap(); - - let Z2: f64 = target.get_item("Z").unwrap().expect("Error: Cannot get key 'Z' from target dict.").extract().unwrap(); - let m2: f64 = target.get_item("m").unwrap().expect("Error: Cannot get key 'm' from target dict.").extract().unwrap(); - let Es2: f64 = target.get_item("Es").unwrap().expect("Error: Cannot get key 'Es' from target dict.").extract().unwrap(); - let Ec2: f64 = target.get_item("Ec").unwrap().expect("Error: Cannot get key 'Ec' from target dict.").extract().unwrap(); - let Eb2: f64 = target.get_item("Eb").unwrap().expect("Error: Cannot get key 'Eb' from target dict.").extract().unwrap(); - let n2: f64 = target.get_item("n").unwrap().expect("Error: Cannot get key 'n' from target dict.").extract().unwrap(); + let Z1: f64 = get_value_from_dict!(ion, "Z")?; + let m1: f64 = get_value_from_dict!(ion, "m")?; + let Es1: f64 = get_value_from_dict!(ion, "Es")?; + let Ec1: f64 = get_value_from_dict!(ion, "Ec")?; + + let Z2: f64 = get_value_from_dict!(target, "Z")?; + let m2: f64 = get_value_from_dict!(target, "m")?; + let Es2: f64 = get_value_from_dict!(target, "Es")?; + let Ec2: f64 = get_value_from_dict!(target, "Ec")?; + let Eb2: f64 = get_value_from_dict!(target, "Eb")?; + let n2: f64 = get_value_from_dict!(target, "n")?; assert!(vx > 0.0, "Input error: vx must be greater than zero for incident particles to hit surface at x=0."); @@ -1271,9 +1271,9 @@ pub fn reflect_single_ion_py<'py>(ion: &Bound<'py, PyDict>, target: &Bound<'py, let vz2 = output[0].dir.z*reflected_velocity; if output[0].E > 0.0 && output[0].dir.x < 0.0 && output[0].left && output[0].incident { - (vx2, vy2, vz2) + Ok((vx2, vy2, vz2)) } else { - (0.0, 0.0, 0.0) + Ok((0.0, 0.0, 0.0)) } } From f801f5549e382e5c7274fbeac1dee8241e919478 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 11 Aug 2026 16:32:28 -0700 Subject: [PATCH 02/14] replace explicit loop with total_cmp ordering on minimum Ec --- src/material.rs | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/src/material.rs b/src/material.rs index 6be7b0e..802e30b 100644 --- a/src/material.rs +++ b/src/material.rs @@ -231,13 +231,7 @@ impl Material { ///The minimum cutoff energy of all species that make up the material. pub fn minimum_cutoff_energy(&self) -> f64 { - let mut min_Ec = self.Ec.iter().sum::(); - for Ec in self.Ec.iter() { - if min_Ec > *Ec { - min_Ec = *Ec; - } - } - min_Ec + self.Ec.clone().into_iter().min_by(f64::total_cmp).unwrap() } ///Choose the parameters of a target atom as a concentration-weighted random draw from the species in the triangle that contains or is nearest to (x, y). From 27e43bacdd7db48452052c688a987b3f2befd303 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 11 Aug 2026 16:33:11 -0700 Subject: [PATCH 03/14] Actually, that function was never used. Deleted. --- src/material.rs | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/material.rs b/src/material.rs index 802e30b..26610ac 100644 --- a/src/material.rs +++ b/src/material.rs @@ -229,11 +229,6 @@ impl Material { } } - ///The minimum cutoff energy of all species that make up the material. - pub fn minimum_cutoff_energy(&self) -> f64 { - self.Ec.clone().into_iter().min_by(f64::total_cmp).unwrap() - } - ///Choose the parameters of a target atom as a concentration-weighted random draw from the species in the triangle that contains or is nearest to (x, y). pub fn choose(&self, x: f64, y: f64, z: f64, rng: &mut ChaCha8Rng) -> (usize, f64, f64, f64, f64, f64, usize) { let random_number: f64 = rng.random::(); From 4a9a700b9407977ba92aba6284114d05063cc657 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 10:42:36 -0700 Subject: [PATCH 04/14] Add draft output_dir --- src/input.rs | 16 +++++-- src/output.rs | 118 +++++++++++++++++++++++++------------------------- 2 files changed, 72 insertions(+), 62 deletions(-) diff --git a/src/input.rs b/src/input.rs index 0a95d6e..1fd1e9f 100644 --- a/src/input.rs +++ b/src/input.rs @@ -212,6 +212,10 @@ fn default_rootfinder() -> Vec> { vec![vec![Rootfinder::DEFAULTNEWTON]] } +fn default_output_dir() -> String { + "./".to_string() +} + /// Rustbca's internal representation of the simulation-level options. #[cfg(not(feature = "distributions"))] #[derive(Deserialize, Clone)] @@ -250,7 +254,9 @@ pub struct Options { #[serde(default = "default_false")] pub track_energy_losses: bool, #[serde(default = "default_seed")] - pub seed: i32 + pub seed: i32, + #[serde(default = "default_output_dir")] + pub output_dir: String, } #[cfg(not(feature = "distributions"))] @@ -275,6 +281,7 @@ impl Options { track_displacements: false, track_energy_losses: false, seed: default_seed(), + output_dir: default_output_dir(), } } } @@ -331,7 +338,9 @@ pub struct Options { pub y_num: usize, pub z_num: usize, #[serde(default = "default_seed")] - pub seed: i32 + pub seed: i32, + #[serde(default = "default_output_dir")] + pub output_dir: String, } #[cfg(feature = "distributions")] @@ -370,7 +379,8 @@ impl Options { x_num: 0, y_num: 0, z_num: 0, - seed: default_seed() + seed: default_seed(), + output_dir: default_output_dir(), } } } diff --git a/src/output.rs b/src/output.rs index 463c0c5..98f94ad 100644 --- a/src/output.rs +++ b/src/output.rs @@ -1,5 +1,6 @@ use super::*; use std::fs::File; +use std::path::Path; #[derive(Clone, Debug)] pub struct OutputUnits { @@ -71,7 +72,7 @@ impl Distributions { .write(true) .create(true) .truncate(true) - .open(format!("{}{}", options.name, "distributions.toml")) + .open(format!("{}{}{}", options.output_dir, options.name, "distributions.toml")) .context("Could not open distributions output file.") .unwrap(); let mut distribution_file_stream = BufWriter::with_capacity(8000, distribution_output_file); @@ -207,7 +208,7 @@ impl SummaryPerSpecies { .write(true) .create(true) .truncate(true) - .open(format!("{}{}", options.name, "summary.output")) + .open(format!("{}{}{}", options.output_dir, options.name, "summary.output")) .context("Could not open output file.") .unwrap(); let writer = BufWriter::with_capacity(8000, summary_output_file); @@ -257,71 +258,70 @@ impl SummaryPerSpecies { } } -/// Open list output files for streaming write -pub fn open_output_lists(options: &Options) -> OutputListStreams { +macro_rules! open_file { + ($output_path:expr, $name:expr, $options:expr) => {{ //Open output files for streaming output - let reflected_file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}", options.name, "reflected.output")) - .context("Could not open output file.") - .unwrap(); - let reflected_file_stream = BufWriter::with_capacity(options.write_buffer_size, reflected_file); - - let sputtered_file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}", options.name, "sputtered.output")) - .context("Could not open output file.") - .unwrap(); - let sputtered_file_stream = BufWriter::with_capacity(options.write_buffer_size, sputtered_file); - - let deposited_file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}", options.name, "deposited.output")) - .context("Could not open output file.") - .unwrap(); - let deposited_file_stream = BufWriter::with_capacity(options.write_buffer_size, deposited_file); - - let trajectory_file = OpenOptions::new() + let file_path = $output_path.join(format!("{}{}{}", $options.name, $name, ".output")); + let file = OpenOptions::new() .write(true) .create(true) .truncate(true) - .open(format!("{}{}", options.name, "trajectories.output")) - .context("Could not open output file.") + .open(file_path) + .context(format!("Could not open {} output file.", $options.name,)) .unwrap(); - let trajectory_file_stream = BufWriter::with_capacity(options.write_buffer_size, trajectory_file); + BufWriter::with_capacity($options.write_buffer_size, file) + }} +} - let trajectory_data = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}", options.name, "trajectory_data.output")) - .context("Could not open output file.") - .unwrap(); - let trajectory_data_stream = BufWriter::with_capacity(options.write_buffer_size, trajectory_data); +/// Open list output files for streaming write +pub fn open_output_lists(options: &Options) -> OutputListStreams { - let displacements_file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}", options.name, "displacements.output")) - .context("Could not open output file.") - .unwrap(); - let displacements_file_stream = BufWriter::with_capacity(options.write_buffer_size, displacements_file); + let output_path = Path::new(&options.output_dir); + assert!(output_path.try_exists().unwrap()); + assert!(output_path.is_dir()); - let energy_loss_file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}", options.name, "energy_loss.output")) - .context("Could not open output file.") - .unwrap(); - let energy_loss_file_stream = BufWriter::with_capacity(options.write_buffer_size, energy_loss_file); + //Open output files for streaming output + let reflected_file_stream = open_file!( + output_path, + "reflected", + options + ); + + let sputtered_file_stream = open_file!( + output_path, + "sputtered", + options + ); + + let deposited_file_stream = open_file!( + output_path, + "deposited", + options + ); + + let trajectory_file_stream = open_file!( + output_path, + "trajectories", + options + ); + + let trajectory_data_stream = open_file!( + output_path, + "trajectory_data", + options + ); + + let displacements_file_stream = open_file!( + output_path, + "displacements", + options + ); + + let energy_loss_file_stream = open_file!( + output_path, + "energy_loss", + options + ); OutputListStreams { reflected_file_stream, From 9349e16f1cccffb49fd17d6bf1a43bc22d39b71f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 11:06:03 -0700 Subject: [PATCH 05/14] Fixed tests.rs by adding output_dir to Options --- src/tests.rs | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/tests.rs b/src/tests.rs index d346226..afba5de 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -230,6 +230,7 @@ fn test_distributions() { track_displacements: false, track_energy_losses: true, seed: 0, + output_dir: ".".to_string(), energy_min: 0.0, energy_max: 10.0, energy_num: 11, @@ -885,6 +886,7 @@ fn test_momentum_conservation() { track_displacements: false, track_energy_losses: false, seed: 0, + output_dir: ".".to_string(), }; #[cfg(feature = "distributions")] @@ -907,6 +909,7 @@ fn test_momentum_conservation() { track_displacements: false, track_energy_losses: false, seed: 0, + output_dir: ".".to_string(), energy_min: 0.0, energy_max: 10.0, energy_num: 11, @@ -1119,6 +1122,7 @@ fn test_quadrature() { track_displacements: false, track_energy_losses: false, seed: 0, + output_dir: ".".to_string(), }; #[cfg(feature = "distributions")] @@ -1141,6 +1145,7 @@ fn test_quadrature() { track_displacements: false, track_energy_losses: false, seed: 0, + output_dir: ".".to_string(), energy_min: 0.0, energy_max: 10.0, energy_num: 11, From 68849e5567e5f589ef9f1cf2671b04d3c69e8a60 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 11:32:09 -0700 Subject: [PATCH 06/14] Updated make_input_file_and_run.py to include tests of output_dir option --- examples/make_input_file_and_run.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py index 37e3a16..3c900f9 100644 --- a/examples/make_input_file_and_run.py +++ b/examples/make_input_file_and_run.py @@ -255,16 +255,22 @@ 'geometry_input': geometry_0D } +# Test rustbca_py and 'output_dir' option input_data['options']['name'] = 'rustbca_input_file' +if not os.path.exists(r'./outputs'): + os.makedirs('outputs') +input_data['options']['output_dir'] = 'outputs' rustbca_py(input_data, mode) -s = np.genfromtxt('rustbca_input_filesputtered.output', delimiter=',') +s = np.genfromtxt('outputs/rustbca_input_filesputtered.output', delimiter=',') arrays = rustbca_local_py(input_data, mode) sputtered = arrays['sputtered'] np.testing.assert_approx_equal(s[0, 2], np.array(arrays['energy'])[sputtered][0]) +# reset these options before running from command line input_data['options']['name'] = 'input_file' +input_data['options']['output_dir'] = '.' # Attempt to cleanup line endings input_string = dumps(input_data).replace('\r', '') From d878e1851cfb1ff32313708e8368e21b737a29bd Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 13:16:10 -0700 Subject: [PATCH 07/14] Update file opening for summaries and distributions to use std::path::Path instead of fragile string concatenation. Fix hardcoded write buffer size (although it should not matter) --- src/output.rs | 82 +++++++++++++++++++++++++++------------------------ 1 file changed, 43 insertions(+), 39 deletions(-) diff --git a/src/output.rs b/src/output.rs index 98f94ad..b94d25e 100644 --- a/src/output.rs +++ b/src/output.rs @@ -1,6 +1,25 @@ use super::*; use std::fs::File; use std::path::Path; +#[cfg(feature = "distributions")] +extern crate ndarray; +#[cfg(feature = "distributions")] +use ndarray::prelude::*; + +macro_rules! open_output_file { + ($output_path:expr, $name:expr, $options:expr) => {{ + //Open output files for streaming output + let file_path = $output_path.join(format!("{}{}{}", $options.name, $name, ".output")); + let file = OpenOptions::new() + .write(true) + .create(true) + .truncate(true) + .open(file_path) + .context(format!("Could not open {} output file.", $options.name,)) + .unwrap(); + BufWriter::with_capacity($options.write_buffer_size, file) + }} +} #[derive(Clone, Debug)] pub struct OutputUnits { @@ -22,12 +41,6 @@ pub fn energy_angle_from_particle(particle: &particle::Particle, units: &OutputU (energy, angle) } -#[cfg(feature = "distributions")] -extern crate ndarray; - -#[cfg(feature = "distributions")] -use ndarray::prelude::*; - /// Distribution tracker for tracking EADs and implantation distributions #[derive(Serialize)] #[cfg(feature = "distributions")] @@ -68,17 +81,22 @@ impl Distributions { /// Write distributions to toml pub fn print(&self, options: &Options) { + let output_path = Path::new(&options.output_dir); + assert!(output_path.try_exists().unwrap()); + assert!(output_path.is_dir()); + + let file_path = output_path.join(format!("{}{}{}", options.name, "distributions", ".toml")); let distribution_output_file = OpenOptions::new() .write(true) .create(true) .truncate(true) - .open(format!("{}{}{}", options.output_dir, options.name, "distributions.toml")) + .open(file_path) .context("Could not open distributions output file.") .unwrap(); - let mut distribution_file_stream = BufWriter::with_capacity(8000, distribution_output_file); + let mut distribution_file_stream = BufWriter::with_capacity(options.write_buffer_size, distribution_output_file); let toml = toml::to_string(&self).unwrap(); + writeln!(distribution_file_stream, "{}", toml).unwrap(); - } /// Updates distributions with a single particle @@ -204,14 +222,15 @@ pub struct SummaryPerSpecies { impl SummaryPerSpecies { pub fn new(options: &Options) -> SummaryPerSpecies { - let summary_output_file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(format!("{}{}{}", options.output_dir, options.name, "summary.output")) - .context("Could not open output file.") - .unwrap(); - let writer = BufWriter::with_capacity(8000, summary_output_file); + let output_path = Path::new(&options.output_dir); + assert!(output_path.try_exists().unwrap()); + assert!(output_path.is_dir()); + + let writer = open_output_file!( + output_path, + "summary", + options + ); SummaryPerSpecies { m: vec![], @@ -258,21 +277,6 @@ impl SummaryPerSpecies { } } -macro_rules! open_file { - ($output_path:expr, $name:expr, $options:expr) => {{ - //Open output files for streaming output - let file_path = $output_path.join(format!("{}{}{}", $options.name, $name, ".output")); - let file = OpenOptions::new() - .write(true) - .create(true) - .truncate(true) - .open(file_path) - .context(format!("Could not open {} output file.", $options.name,)) - .unwrap(); - BufWriter::with_capacity($options.write_buffer_size, file) - }} -} - /// Open list output files for streaming write pub fn open_output_lists(options: &Options) -> OutputListStreams { @@ -281,43 +285,43 @@ pub fn open_output_lists(options: &Options) -> OutputListStreams { assert!(output_path.is_dir()); //Open output files for streaming output - let reflected_file_stream = open_file!( + let reflected_file_stream = open_output_file!( output_path, "reflected", options ); - let sputtered_file_stream = open_file!( + let sputtered_file_stream = open_output_file!( output_path, "sputtered", options ); - let deposited_file_stream = open_file!( + let deposited_file_stream = open_output_file!( output_path, "deposited", options ); - let trajectory_file_stream = open_file!( + let trajectory_file_stream = open_output_file!( output_path, "trajectories", options ); - let trajectory_data_stream = open_file!( + let trajectory_data_stream = open_output_file!( output_path, "trajectory_data", options ); - let displacements_file_stream = open_file!( + let displacements_file_stream = open_output_file!( output_path, "displacements", options ); - let energy_loss_file_stream = open_file!( + let energy_loss_file_stream = open_output_file!( output_path, "energy_loss", options From 2855ee1b42c156af382471606eecdab46741a727 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 16:24:18 -0700 Subject: [PATCH 08/14] Add check on maximum number density to catch errors. --- src/consts.rs | 4 +++- src/geometry.rs | 7 ++++++- src/parry.rs | 2 ++ src/sphere.rs | 3 ++- 4 files changed, 13 insertions(+), 3 deletions(-) diff --git a/src/consts.rs b/src/consts.rs index 829fb7a..111ae14 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -44,4 +44,6 @@ pub const TABLE_SIZE: usize = Z_MAX*(Z_MAX + 1)/2; /// Gauss-Legendre Points pub const GL_X: [f64; 5] = [0./2. + 1./2., -0.538469/2. + 1./2., 0.538469/2. + 1./2., -0.90618/2. + 1./2., 0.90618/2. + 1./2.]; /// Gauss-Legendre Weights -pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.]; \ No newline at end of file +pub const GL_W: [f64; 5] = [0.568889/2., 0.478629/2., 0.478629/2., 0.236927/2., 0.236927/2.]; +/// Enforced maximum number density +pub const MAX_DENSITY: f64 = 1e32; \ No newline at end of file diff --git a/src/geometry.rs b/src/geometry.rs index fb121c6..f7f31f9 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -65,8 +65,9 @@ impl Geometry for Mesh0D { let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); assert!(!densities.is_empty(), "Input Error: density list empty."); - + let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.; @@ -154,6 +155,7 @@ impl Geometry for Mesh1D { let densities: Vec> = geometry_input.densities .iter() .map( |row| row.iter().map(|element| element/(length_unit).powi(3)).collect() ).collect(); + //Assert all layer density lists are equal length assert!( @@ -170,6 +172,7 @@ impl Geometry for Mesh1D { layer_bottom += layer_thickness*length_unit; let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); layers.push(Layer1D::new(layer_top, layer_bottom, densities, concentrations, ck)); @@ -308,6 +311,7 @@ impl Geometry for HomogeneousMesh2D { let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.; @@ -476,6 +480,7 @@ impl Geometry for Mesh2D { y3*length_unit, ); let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck)); diff --git a/src/parry.rs b/src/parry.rs index 8686683..ae1e886 100644 --- a/src/parry.rs +++ b/src/parry.rs @@ -77,6 +77,7 @@ impl Geometry for ParryBall { let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); let radius = input.radius*length_unit; @@ -201,6 +202,7 @@ impl Geometry for ParryTriMesh { let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); let points = input.vertices.iter().map(|p| Point::new(p[0]*length_unit , p[1]*length_unit , p[2]*length_unit)).collect(); diff --git a/src/sphere.rs b/src/sphere.rs index a9326fc..bb51d05 100644 --- a/src/sphere.rs +++ b/src/sphere.rs @@ -63,13 +63,14 @@ impl Geometry for Sphere { "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())), }; let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor; let densities: Vec = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect(); let total_density: f64 = densities.iter().sum(); + assert!(total_density < MAX_DENSITY, "Input Error: total density {}/m^3 exceeds realistic values; check values or units.", total_density); let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); let radius = input.radius*length_unit; From 9a19ba503d91f8aae77750bd15624af565bc7229 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 17:24:00 -0700 Subject: [PATCH 09/14] Found input bug while working on WW benchmark - a <= on interaction_index length was defaulting to 0 inappropriately. --- src/bca.rs | 2 ++ src/input.rs | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 76b80df..f3de757 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -612,6 +612,8 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame complex_threshold: f64, truncation_threshold: f64, far_from_zero: f64, interval_limit: f64, derivative_free: bool) -> Result { + println!("yes"); + //Lindhard screening length and reduced energy let a = interactions::screening_length(Za, Zb, interaction_potential); let reduced_energy = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E0; diff --git a/src/input.rs b/src/input.rs index 1fd1e9f..6f393c8 100644 --- a/src/input.rs +++ b/src/input.rs @@ -432,11 +432,11 @@ pub fn process_input_file(input: ::InputFileFormat) assert!(material.m.len() == material.Eb.len(), "Input error: material input arrays of unequal length."); assert!(material.m.len() == material.Es.len(), "Input error: material input arrays of unequal length."); - if material.interaction_index.len() <= 1 { + if material.interaction_index.is_empty() { material.interaction_index = vec![0; material.m.len()]; } - if material.Ed.len() <= 1 { + if material.Ed.is_empty() { material.Ed = vec![material.Ed[0]; material.m.len()]; } From fadbc2424b3b4a8541a4dc830e77ca9ebbaaf23e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 17:24:34 -0700 Subject: [PATCH 10/14] Relax epsilon of Kr-C Morse --- examples/test_morse.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/test_morse.py b/examples/test_morse.py index 227203c..0e8e86c 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -87,7 +87,7 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): mean_free_path_model = "LIQUID" interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=7E10, x0=0.75E-10}}}}]] scattering_integral = [["GAUSS_LEGENDRE"]] - root_finder = [[{{"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}}}}]] + root_finder = [[{{"CPR"={{n0=2, nmax=100, epsilon=1E-3, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}}}]] num_threads = 4 num_chunks = 10 @@ -172,7 +172,7 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): #Running and plotting the H-Ni simulations with the Morse potential and updated Es num_energies = 15 energies = np.logspace(-1, 4, num_energies) -run_sim = True +run_sim = False num_samples = 10000 R_N = np.zeros(num_energies) R_E = np.zeros(num_energies) @@ -180,8 +180,8 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): R_E_2 = np.zeros(num_energies) for index, energy in enumerate(energies): - R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=True) - R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=True) + R_N[index], R_E[index] = run_krc_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim) + R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim) plt.semilogx(energies, R_N, label='R_N Morse-Kr-C H-Ni, Es=1.5eV', color='purple') plt.semilogx(energies, R_N_2, label='R_N Morse H-Ni, Es=1.5eV', color='green') From 9c6021e158cc30704b37ef6c494fcc9b61f02b52 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 17:53:09 -0700 Subject: [PATCH 11/14] Previous attempt broke interaction indices; trying again by filling vec with current value if len 1 --- src/bca.rs | 7 ++----- src/input.rs | 6 +++--- 2 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index f3de757..61e747c 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -612,8 +612,6 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame complex_threshold: f64, truncation_threshold: f64, far_from_zero: f64, interval_limit: f64, derivative_free: bool) -> Result { - println!("yes"); - //Lindhard screening length and reduced energy let a = interactions::screening_length(Za, Zb, interaction_potential); let reduced_energy = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E0; @@ -624,7 +622,7 @@ 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)}; - let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential); + let upper_bound = 10.0*impact_parameter + interactions::crossing_point_doca(interaction_potential); let roots = match derivative_free { true => find_roots_with_secant_polishing(&g, &f, 1e-15, upper_bound, @@ -638,8 +636,7 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame truncation_threshold, interval_limit, far_from_zero) } }.with_context(|| format!("Numerical error: CPR Rootfinder failed to converge when calculating distance of closest approach for Er = {} eV p = {} A using {}.", - relative_energy/EV, impact_parameter/ANGSTROM, interaction_potential)) - .unwrap(); + relative_energy/EV, impact_parameter/ANGSTROM, interaction_potential))?; let max_root = roots.iter().cloned().fold(f64::NAN, f64::max)/a; diff --git a/src/input.rs b/src/input.rs index 6f393c8..b79b87b 100644 --- a/src/input.rs +++ b/src/input.rs @@ -432,11 +432,11 @@ pub fn process_input_file(input: ::InputFileFormat) assert!(material.m.len() == material.Eb.len(), "Input error: material input arrays of unequal length."); assert!(material.m.len() == material.Es.len(), "Input error: material input arrays of unequal length."); - if material.interaction_index.is_empty() { - material.interaction_index = vec![0; material.m.len()]; + if material.Ed.len() <= 1 { + material.interaction_index = vec![material.interaction_index[0]; material.m.len()]; } - if material.Ed.is_empty() { + if material.Ed.len() <= 1 { material.Ed = vec![material.Ed[0]; material.m.len()]; } From 02c6995436c46c953140c7c3f4243abe9674857f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 22:48:19 -0700 Subject: [PATCH 12/14] Minor fixes found while working on W-W benchmark --- src/interactions.rs | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index 2def3ac..ceb6117 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -176,7 +176,7 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent 1./(1. + (r*alpha).powi(2)) } InteractionPotential::WW => { - 1. + 1./(1. + (r/a).powi(2)) }, InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => { 1./(1. + (r*alpha).powi(2)) @@ -465,8 +465,8 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 { let x2 = 2.10004200084; if x <= x1 { - - let a = screening_length(74., 74., InteractionPotential::ZBL); + // + let a = screening_length(74., 74., InteractionPotential::ZBL)*1.000_250_544_359; screened_coulomb(r, a, 74., 74., InteractionPotential::ZBL) } else if x <= x2 { @@ -512,17 +512,8 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 { /// Distance of closest approach function for the W-W cublic spline potential from Bjorkas et al. pub fn doca_tungsten_tungsten_cubic_spline(r: f64, p: f64, relative_energy: f64) -> f64 { - let x = r/ANGSTROM; - let x1 = 1.10002200044; - let x2 = 2.10004200084; - - if x <= x1 { - let a = screening_length(74., 74., InteractionPotential::ZBL); - distance_of_closest_approach_function_singularity_free(r, a, 74., 74., relative_energy, p, InteractionPotential::ZBL) - } else { - (r/ANGSTROM).powi(2) - (r/ANGSTROM).powi(2)*tungsten_tungsten_cubic_spline(r)/relative_energy - p.powi(2)/ANGSTROM.powi(2) - } + x.powi(2) - x.powi(2)*tungsten_tungsten_cubic_spline(r)/relative_energy - p.powi(2)/ANGSTROM.powi(2) } fn heaviside(x: f64) -> f64 { From 916901e1e3fedf8a34aeb9811cf18142a3b4bad6 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 22:52:03 -0700 Subject: [PATCH 13/14] Typoed an interaction index check in Mesh0D; fixed. --- src/input.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/input.rs b/src/input.rs index b79b87b..9f09883 100644 --- a/src/input.rs +++ b/src/input.rs @@ -432,7 +432,7 @@ pub fn process_input_file(input: ::InputFileFormat) assert!(material.m.len() == material.Eb.len(), "Input error: material input arrays of unequal length."); assert!(material.m.len() == material.Es.len(), "Input error: material input arrays of unequal length."); - if material.Ed.len() <= 1 { + if material.interaction_index.len() <= 1 { material.interaction_index = vec![material.interaction_index[0]; material.m.len()]; } From 11cf2866b8d37451e6a8e38b00ca893cb320932c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 13 Aug 2026 23:12:43 -0700 Subject: [PATCH 14/14] Update Cargo.toml and setup.py to reflect the fact that the python library can now use attractive repulsive potentials. --- Cargo.toml | 2 +- setup.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 661bbed..ee9637a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -44,4 +44,4 @@ cpr_rootfinder = ["rcpr"] distributions = ["ndarray"] no_list_output = [] parry3d = ["parry3d-f64"] -python = ["pyo3", "pythonize"] \ No newline at end of file +python = ["pyo3", "pythonize", "rcpr"] \ No newline at end of file diff --git a/setup.py b/setup.py index 7d2962f..f6618a8 100644 --- a/setup.py +++ b/setup.py @@ -8,7 +8,7 @@ RustExtension( "libRustBCA", binding=Binding.PyO3, - features=["python", "parry3d", "pythonize"], + features=["python", "parry3d", "pythonize", "cpr_rootfinder"], ) ], # rust extensions are not zip safe, just like C-extensions.