Python improvements and output_dir implementation

.github/workflows/rustbca_compile_check.yml

@@ -44,6 +44,7 @@ jobs:
         python3 examples/make_input_file_and_run.py
         python3 examples/test_electronic_stopping.py
         python3 examples/test_different_options.py
+        python3 examples/test_morse.py
     - name: Test Fortran and C bindings
       run : |
         cargo build --release --lib --features parry3d

Cargo.toml

@@ -44,4 +44,4 @@ cpr_rootfinder = ["rcpr"]
 distributions = ["ndarray"]
 no_list_output = []
 parry3d = ["parry3d-f64"]
-python = ["pyo3", "pythonize"]
+python = ["pyo3", "pythonize", "rcpr"]

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', '')

examples/test_morse.py

@@ -25,7 +25,7 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True):
     mean_free_path_model = "LIQUID"
     interaction_potential = [[{{"MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10}}}}]]
     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=3, nmax=100, epsilon=1E-9, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-13, derivative_free=true}}}}]]
     num_threads = 4
     num_chunks = 10
 
@@ -87,9 +87,9 @@ 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}}}}]]
-    num_threads = 4
-    num_chunks = 10
+    root_finder = [[{{"CPR"={{n0=2, nmax=200, epsilon=1E-9, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-13, derivative_free=true}}}}]]
+    num_threads = 6
+    num_chunks = 1
 
     [particle_parameters]
     length_unit = "ANGSTROM"
@@ -173,15 +173,30 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
 num_energies = 15
 energies = np.logspace(-1, 4, num_energies)
 run_sim = True
-num_samples = 10000
+num_samples = 100
 R_N = np.zeros(num_energies)
 R_E = np.zeros(num_energies)
 R_N_2 = np.zeros(num_energies)
 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)
+
+R_N_test = [
+    0.00, 0.01, 0.28, 0.60, 0.90,
+    0.95, 0.88, 0.81, 0.70, 0.49,
+    0.30, 0.24, 0.17, 0.11, 0.10
+]
+
+R_N_2_test = [
+    0.00, 0.01, 0.28, 0.60, 0.90,
+    0.95, 0.88, 0.81, 0.74, 0.60,
+    0.46, 0.23, 0.05, 0.00, 0.00
+]
+
+np.testing.assert_allclose(R_N, R_N_test)
+np.testing.assert_allclose(R_N_2, R_N_2_test)
 
 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')

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.

src/bca.rs

@@ -623,21 +623,21 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame
         interactions::scaling_function(r, impact_parameter, interaction_potential)};
 
     let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential);
+    let lower_bound = impact_parameter / 1000.0;
 
     let roots = match derivative_free {
-        true => find_roots_with_secant_polishing(&g, &f, 1e-15, upper_bound,
+        true => find_roots_with_secant_polishing(&g, &f, lower_bound, upper_bound,
             n0, epsilon, nmax, complex_threshold,
             truncation_threshold, interval_limit, far_from_zero),
 
         false => {
             let df = |r: f64| -> f64 {interactions::diff_distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)};
-            find_roots_with_newton_polishing(&g, &f, &df, 1e-15, upper_bound,
+            find_roots_with_newton_polishing(&g, &f, &df, lower_bound, upper_bound,
             n0, epsilon, nmax, complex_threshold,
             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;
 

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.];
+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;

src/geometry.rs

@@ -65,8 +65,9 @@ impl Geometry for Mesh0D {
 
         let densities: Vec<f64> = 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<Vec<f64>> = 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<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             layers.push(Layer1D::new(layer_top, layer_bottom, densities, concentrations, ck));
@@ -308,6 +311,7 @@ impl Geometry for HomogeneousMesh2D {
         let densities: Vec<f64> = 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<f64> = densities.iter().map(|&density| density/total_density).collect::<Vec<f64>>();
 
             cells.push(Cell2D::new(coordinate_set_converted, densities, concentrations, ck));

src/input.rs

@@ -212,6 +212,10 @@ fn default_rootfinder() -> Vec<Vec<Rootfinder>> {
     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(),
         }
     }
 }
@@ -423,7 +433,7 @@ pub fn process_input_file<T: Geometry>(input: <T as Geometry>::InputFileFormat)
     assert!(material.m.len() == material.Es.len(), "Input error: material input arrays of unequal length.");
 
     if material.interaction_index.len() <= 1 {
-        material.interaction_index = vec![0; material.m.len()];
+        material.interaction_index = vec![material.interaction_index[0]; material.m.len()];
     }
 
     if material.Ed.len() <= 1 {

src/interactions.rs

@@ -8,7 +8,7 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 {
     match interaction_potential {
         InteractionPotential::LENNARD_JONES_12_6{sigma, ..} | InteractionPotential::LENNARD_JONES_65_6{sigma, ..} => sigma,
         InteractionPotential::MORSE{D, alpha, r0} => (alpha*r0 - (2.0_f64).ln())/alpha,
-        InteractionPotential::WW => 50.*ANGSTROM,
+        InteractionPotential::WW => 10.*ANGSTROM,
         _ => 10.*ANGSTROM,
     }
 }
@@ -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/ANGSTROM).powi(2))
         },
         InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => {
             1./(1. + (r*alpha).powi(2))
@@ -465,8 +465,10 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 {
     let x2 = 2.10004200084;
 
     if x <= x1 {
-
-        let a = screening_length(74., 74., InteractionPotential::ZBL);
+        // This correction factor is necessary to correct for a discrepancy in the 
+        // published interaction potential - it is slightly off from the "correct"
+        // ZBL screening length. Without this, the doca function is discontinuous.
+        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 +514,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/ANGSTROM).powi(2)
 }
 
 fn heaviside(x: f64) -> f64 {

src/lib.rs

@@ -1196,19 +1196,19 @@ pub fn compound_bca_list_tracked_py<'py>(energies: Vec<f64>, ux: Vec<f64>, 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))
     }
 }
 

src/material.rs

@@ -229,17 +229,6 @@ impl <T: Geometry> Material<T> {
         }
     }
 
-    ///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::<f64>();
-        for Ec in self.Ec.iter() {
-            if min_Ec > *Ec {
-                min_Ec = *Ec;
-            }
-        }
-        min_Ec
-    }
-
     ///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::<f64>();