Minor code improvements

.github/workflows/rustbca_compile_check.yml

@@ -32,9 +32,6 @@ jobs:
     - name: Install Python libraries
       run: |
         python3 -m pip install numpy shapely scipy matplotlib toml tomlkit
-    - name: Install HDF5 Libraries
-      run: |
-        sudo apt install libhdf5-dev
     - name: test Python Bindings and run testing scripts
       run: |
         python3 -m pip install setuptools_rust testresources setuptools wheel
@@ -59,7 +56,7 @@ jobs:
         ./a.out
     - name: Test RustBCA
       run: |
-        cargo test --features cpr_rootfinder,hdf5_input,distributions,parry3d
+        cargo test --features cpr_rootfinder,distributions,parry3d
     - name: Run Examples
       run: |
         cargo run --release 0D examples/boron_nitride_0D.toml

.gitignore

@@ -0,0 +1,8 @@
+target/
+debug/
+Cargo.lock
+*.output
+**/*.rs.bk
+*.pdb
+src/RustBCA.egg-info
+src/RustBCA.egg-info/*

Cargo.toml

@@ -1,50 +1,48 @@
-[package]
-name = "RustBCA"
-version = "2.9.0"
-default-run = "RustBCA"
-authors = ["Jon Drobny <drobny2@illinois.edu>", "Jon Drobny <jdrobny@tae.com>"]
-edition = "2021"
-
-[[bin]]
-name = "RustBCA"
-path = "src/main.rs"
-
-[lib]
-name = "libRustBCA"
-path = "src/lib.rs"
-crate-type = ["cdylib", "lib"]
-
-[dependencies]
-rand = {version = "0.10.2", features=["chacha"]}
-rand_distr = "0.6.0"
-toml = "1.1.4"
-anyhow = "1.0.71"
-itertools = "0.15.0"
-rayon = "1.10.0"
-geo = {version = "0.25", optional = false}
-indicatif = {version = "0.18.6", features=["rayon"]}
-serde = { version = "1.0.163", features = ["derive"] }
-hdf5 = {version = "0.7.1", optional = true}
-rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
-ndarray = {version = "0.17.2", features = ["serde"], optional = true}
-parry3d-f64 = {optional = true, version="0.2.0"}
-pyo3 = {version = "0.29.0", optional=true}
-pythonize = {version = "0.29.0", optional=true}
-
-[dev-dependencies]
-float-cmp = "0.10.0"
-
-[profile.release]
-lto = "fat"
-codegen-units = 1
-opt-level = 3
-debug = false
-
-[features]
-hdf5_input = ["hdf5"]
-cpr_rootfinder = ["rcpr"]
-distributions = ["ndarray"]
-no_list_output = []
-parry3d = ["parry3d-f64"]
-python = ["pyo3", "pythonize"]
+[package]
+name = "RustBCA"
+version = "3.0.0"
+default-run = "RustBCA"
+authors = ["Jon Drobny <drobny2@illinois.edu>", "Jon Drobny <jdrobny@tae.com>"]
+edition = "2024"
+
+[[bin]]
+name = "RustBCA"
+path = "src/main.rs"
+
+[lib]
+name = "libRustBCA"
+path = "src/lib.rs"
+crate-type = ["cdylib", "lib"]
+
+[dependencies]
+rand = {version = "0.10.2", features=["chacha"]}
+rand_distr = "0.6.0"
+toml = "1.1.4"
+anyhow = "1.0.71"
+itertools = "0.15.0"
+rayon = "1.10.0"
+geo = {version = "0.25", optional = false}
+indicatif = {version = "0.18.6", features=["rayon"]}
+serde = { version = "1.0.163", features = ["derive"] }
+rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true}
+ndarray = {version = "0.17.2", features = ["serde"], optional = true}
+parry3d-f64 = {optional = true, version="0.2.0"}
+pyo3 = {version = "0.29.0", optional=true}
+pythonize = {version = "0.29.0", optional=true}
+
+[dev-dependencies]
+float-cmp = "0.10.0"
+
+[profile.release]
+lto = "fat"
+codegen-units = 1
+opt-level = 3
+debug = false
+
+[features]
+cpr_rootfinder = ["rcpr"]
+distributions = ["ndarray"]
+no_list_output = []
+parry3d = ["parry3d-f64"]
+python = ["pyo3", "pythonize"]
 extended_max_z = []

examples/layered_geometry.toml

@@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
 mean_free_path_model = "LIQUID"
 interaction_potential = [["KR_C"]]
 scattering_integral = [["MENDENHALL_WELLER"]]
-use_hdf5 = false
 root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]]
 num_threads = 4
 num_chunks = 10

examples/layered_geometry_1D.toml

@@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
 mean_free_path_model = "LIQUID"
 interaction_potential = [["KR_C"]]
 scattering_integral = [["MENDENHALL_WELLER"]]
-use_hdf5 = false
 root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]]
 num_threads = 4
 num_chunks = 10

examples/lithium_vapor_shield.toml

@@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_LOCAL"
 mean_free_path_model = {"THRESHOLD"={density=1e28}}
 interaction_potential = [["KR_C"]]
 scattering_integral = [["MENDENHALL_WELLER"]]
-use_hdf5 = false
 root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]]
 num_threads = 4
 num_chunks = 2

examples/multiple_interaction_potentials.toml

@@ -22,7 +22,6 @@ root_finder = [
   [{"NEWTON"={max_iterations=100, tolerance=1E-6}}, {"POLYNOMIAL"={complex_threshold=1E-12}}],
   [{"POLYNOMIAL"={complex_threshold=1E-12}}, {"NEWTON"={max_iterations=100, tolerance=1E-6}}]
 ]
-use_hdf5 = false
 num_threads = 4
 num_chunks = 10
 track_displacements = false

examples/titanium_dioxide_0D.toml

@@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
 mean_free_path_model = {"THRESHOLD"={density=1E28}}
 interaction_potential = [["KR_C"]]
 scattering_integral = [["MENDENHALL_WELLER"]]
-use_hdf5 = false
 root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]]
 num_threads = 4
 num_chunks = 100

examples/tungsten_tiles_trimesh.toml

@@ -10,7 +10,6 @@ suppress_deep_recoils = false
 high_energy_free_flight_paths = false
 num_threads = 4
 num_chunks = 5
-use_hdf5 = false
 electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
 mean_free_path_model = "LIQUID"
 interaction_potential = [["KR_C"]]

examples/tungsten_twist_trimesh.toml

@@ -10,7 +10,6 @@ suppress_deep_recoils = false
 high_energy_free_flight_paths = false
 num_threads = 4
 num_chunks = 5
-use_hdf5 = false
 electronic_stopping_mode = "LOW_ENERGY_NONLOCAL"
 mean_free_path_model = "LIQUID"
 interaction_potential = [["KR_C"]]

setup.py

@@ -3,15 +3,12 @@
 
 setup(
     name="RustBCA",
-    version="2.8.4",
+    version="3.0.0",
     rust_extensions=[
         RustExtension(
             "libRustBCA",
             binding=Binding.PyO3,
-            features=["python", "parry3d"],
-            #args=["+nightly", "--edition 2021", "-Z unstable-options"],
-            #optional=True,
-            #rust_version="1.57.0"
+            features=["python", "parry3d", "pythonize"],
         )
     ],
     # rust extensions are not zip safe, just like C-extensions.

src/bca.rs

@@ -425,8 +425,6 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa
 pub fn subtract_electronic_stopping_energy<T: Geometry>(particle_1: &mut particle::Particle, material: &material::Material<T>, distance_traveled: f64,
     x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) -> f64 {
 
-    //assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following collision.");
-
     let x = particle_1.pos.x;
     let y = particle_1.pos.y;
     let z = particle_1.pos.z;
@@ -452,9 +450,7 @@ pub fn subtract_electronic_stopping_energy<T: Geometry>(particle_1: &mut particl
 
         particle_1.E += -delta_energy_electronic;
         //Make sure particle energy doesn't become negative again
-        /*
-        assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following electronic stopping.");
-        */
+
         if particle_1.E < 0. {
             particle_1.E = 0.;
         }
@@ -479,17 +475,12 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: &
     let x0 = distance_of_closest_approach(particle_1, particle_2, binary_collision_geometry, options);
 
     let theta: f64 = match  scattering_integral {
-        ScatteringIntegral::MENDENHALL_WELLER => mendenhall_weller(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential),
-        ScatteringIntegral::GAUSS_MEHLER{n_points} => gauss_mehler(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential, n_points),
-        ScatteringIntegral::GAUSS_LEGENDRE => gauss_legendre(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential),
-        ScatteringIntegral::MAGIC => magic(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential),
+        ScatteringIntegral::MENDENHALL_WELLER => mendenhall_weller(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential),
+        ScatteringIntegral::GAUSS_MEHLER{n_points} => gauss_mehler(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential, n_points),
+        ScatteringIntegral::GAUSS_LEGENDRE => gauss_legendre(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential),
+        ScatteringIntegral::MAGIC => magic(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential),
     };
 
-    /*
-    if theta.is_nan() {
-        return Err(anyhow!("Numerical error: CoM deflection angle is NaN for {}. Check input parameters.", binary_collision_geometry));
-    }
-    */
     let (sin_theta, cos_theta) = theta.sin_cos();
     let sin_2_theta = (theta/2.).sin();
 
@@ -695,25 +686,23 @@ pub fn newton_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par
 }
 
 /// Gauss-Mehler quadrature.
-pub fn gauss_mehler(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential, n_points: usize) -> f64 {
-    let a: f64 = interactions::screening_length(Za, Zb, interaction_potential);
+pub fn gauss_mehler(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential, n_points: usize) -> f64 {
     let r0 = x0*a;
     let V = |r| {interactions::interaction_potential(r, a, Za, Zb, interaction_potential)};
     let relative_energy = E0*Mb/(Ma + Mb);
     scattering_integral_gauss_mehler(impact_parameter, relative_energy, r0, &V, n_points)
 }
 
 /// Gauss-Legendre quadrature.
-pub fn gauss_legendre(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential) -> f64 {
-    let a: f64 = interactions::screening_length(Za, Zb, interaction_potential);
+pub fn gauss_legendre(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential) -> f64 {
     let r0 = x0*a;
     let V = |r| {interactions::interaction_potential(r, a, Za, Zb, interaction_potential)};
     let relative_energy = E0*Mb/(Ma + Mb);
     scattering_integral_gauss_legendre(impact_parameter, relative_energy, r0, &V)
 }
 
 /// Ziegler's MAGIC algorithm for approximating the scattering integral.
-pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential) -> f64 {
+pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential) -> f64 {
     //MAGIC algorithm
     //Since this is legacy code I don't think I will clean this up
     let C_ = match  interaction_potential {
@@ -723,7 +712,6 @@ pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64,
         InteractionPotential::TRIDYN => [1.0144, 0.235809, 0.126, 69350., 83550.], //Undocumented Tridyn constants
         _ => panic!("Input error: unimplemented interaction potential {} for MAGIC algorithm. Use a screened Coulomb potential.",  interaction_potential)
     };
-    let a: f64 = interactions::screening_length(Za, Zb, interaction_potential);
     let beta: f64 = impact_parameter/a;
     let V0 = Za*Zb*Q*Q/4.0/PI/EPS0/a;
     let relative_energy = E0*Mb/(Ma + Mb);
@@ -740,9 +728,8 @@ pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64,
 }
 
 /// Mendenhall-Weller quadrature.
-pub fn mendenhall_weller(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential) -> f64 {
+pub fn mendenhall_weller(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential) -> f64 {
     //Lindhard screening length and reduced energy
-    let a: f64 = interactions::screening_length(Za, Zb, interaction_potential);
     let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E0;
     let beta: f64 = impact_parameter/a;