Python compound bca; minor refactors

.github/ISSUE_TEMPLATE/benchmark_request.md

@@ -0,0 +1,14 @@
+---
+name: Benchmark request
+about: Suggest a benchmark for validation
+title: "[benchmark]"
+labels: benchmark
+assignees: ''
+
+---
+
+**Description**
+A clear and concise description of the requested benchmark.
+
+**Source**
+Source of data for benchmark; a paper, code, or database with sputtering yields, reflection coefficients, etc.

src/bca.rs

@@ -75,9 +75,6 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
 
         //Remove particle from top of vector as particle_1
         let mut particle_1 = particles.pop().unwrap();
-
-        //println!("Particle start Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM);
-
         //BCA loop
         while !particle_1.stopped & !particle_1.left {
 
@@ -114,8 +111,6 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
                             particle_1.pos.x, particle_2.pos.x, &binary_collision_geometry))
                         .unwrap();
 
-                    //println!("{}", binary_collision_result);
-
                     //Only use 0th order collision for local electronic stopping
                     if k == 0 {
                         normalized_distance_of_closest_approach = binary_collision_result.normalized_distance_of_closest_approach;
@@ -127,17 +122,14 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
                     particle_2.E = binary_collision_result.recoil_energy - material.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z);
                     particle_2.energy_origin = particle_2.E;
 
-                    //Accumulate asymptotic deflections for primary particle
+                    //Accumulate energy losses and asymptotic deflections for primary particle
                     total_energy_loss += binary_collision_result.recoil_energy;
-
-                    //total_deflection_angle += psi;
                     total_asymptotic_deflection += binary_collision_result.asymptotic_deflection;
 
-                    //Rotate particle 1, 2 by lab frame scattering angles
-                    particle::rotate_particle(&mut particle_1, binary_collision_result.psi,
+                    particle_1.rotate(binary_collision_result.psi,
                         binary_collision_geometry.phi_azimuthal);
 
-                    particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil,
+                    particle_2.rotate(-binary_collision_result.psi_recoil,
                         binary_collision_geometry.phi_azimuthal);
 
                     particle_2.dir_old.x = particle_2.dir.x;
@@ -184,28 +176,24 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
 
             //Advance particle in space and track total distance traveled
             #[cfg(not(feature = "accelerated_ions"))]
-            let distance_traveled = particle::particle_advance(&mut particle_1,
+            let distance_traveled = particle_1.advance(
                 binary_collision_geometries[0].mfp, total_asymptotic_deflection);
 
             #[cfg(feature = "accelerated_ions")]
-            let distance_traveled = particle::particle_advance(&mut particle_1,
+            let distance_traveled = particle_1.advance(
                 binary_collision_geometries[0].mfp + distance_to_target - material.geometry.get_energy_barrier_thickness(), total_asymptotic_deflection);
 
             //Subtract total energy from all simultaneous collisions and electronic stopping
             bca::update_particle_energy(&mut particle_1, &material, distance_traveled,
                 total_energy_loss, normalized_distance_of_closest_approach, strong_collision_Z,
                 strong_collision_index, &options);
-            //println!("Particle finished collision loop Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM);
-
-            //println!("{} {} {}", energy_0/EV, energy_1/EV, (energy_1 - energy_0)/EV);
 
             //Check boundary conditions on leaving and stopping
             material::boundary_condition_planar(&mut particle_1, &material);
 
             //Set particle index to topmost particle
             particle_index = particles.len();
         }
-        //println!("Particle stopped or left Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM);
         particle_output.push(particle_1);
     }
     particle_output

src/geometry.rs

@@ -97,7 +97,6 @@ impl Geometry for Mesh0D {
     }
 
     fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool {
-        //println!("x: {} energy_barrier_thickness: {}", x/ANGSTROM, self.energy_barrier_thickness/ANGSTROM);
         x > -10.*self.energy_barrier_thickness
     }
 

src/interactions.rs

@@ -8,7 +8,6 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 {
         InteractionPotential::MORSE{D, alpha, r0} => (alpha*r0 - (2.0_f64).ln())/alpha,
         InteractionPotential::WW => 50.*ANGSTROM,
         _ => 10.*ANGSTROM,
-        //_ => panic!("Input error: potential never crosses zero for r > 0. Consider using the Newton rootfinder.")
     }
 
 }

src/lib.rs

@@ -80,6 +80,7 @@ pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh,
 pub fn pybca(py: Python, m: &PyModule) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(simple_bca_py, m)?)?;
     m.add_function(wrap_pyfunction!(simple_bca_list_py, m)?)?;
+    m.add_function(wrap_pyfunction!(compound_bca_list_py, m)?)?;
     Ok(())
 }
 
@@ -980,6 +981,198 @@ pub extern "C" fn simple_bca_c(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64
     }
 }
 
+#[cfg(feature = "python")]
+///compound_tagged_bca_list_py(ux, uy,  uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2)
+/// runs a BCA simulation for a list of particles and outputs a list of sputtered, reflected, and implanted particles.
+/// Args:
+///    energies (list(f64)): initial ion energies in eV.
+///    ux (list(f64)): initial ion directions x. ux != 0.0 to avoid gimbal lock
+///    uy (list(f64)): initial ion directions y.
+///    uz (list(f64)): initial ion directions z.
+///    Z1 (list(f64)): initial ion atomic numbers.
+///    m1 (list(f64)): initial ion masses in amu.
+///    Ec1 (list(f64)): ion cutoff energies in eV. If ion energy < Ec1, it stops in the material.
+///    Es1 (list(f64)): ion surface binding energies. Assumed planar.
+///    Z2 (list(f64)): target material species atomic numbers.
+///    m2 (list(f64)): target material species masses in amu.
+///    Ec2 (list(f64)): target material species cutoff energies in eV. If recoil energy < Ec2, it stops in the material.
+///    Es2 (list(f64)): target species surface binding energies. Assumed planar.
+///    n2 (list(f64)): target material species atomic number densities in inverse cubic Angstroms.
+///    Eb2 (list(f64)): target material species bulk binding energies in eV.
+/// Returns:
+///    output (NX9 list of f64): each row in the list represents an output particle (implanted,
+///    sputtered, or reflected). Each row consists of:
+///      [Z, m (amu), E (eV), x, y, z, (angstrom), ux, uy, uz]
+///    incident (list(bool)): whether each row of output was an incident ion or originated in the target
+#[pyfunction]
+pub fn compound_bca_list_py(energies: Vec<f64>, ux: Vec<f64>, uy: Vec<f64>, uz: Vec<f64>, Z1: Vec<f64>, m1: Vec<f64>, Ec1: Vec<f64>, Es1: Vec<f64>, Z2: Vec<f64>, m2: Vec<f64>, Ec2: Vec<f64>, Es2: Vec<f64>, n2: Vec<f64>, Eb2: Vec<f64>) -> (Vec<[f64; 9]>, Vec<bool>) {
+    let mut total_output = vec![];
+    let mut incident = vec![];
+    let num_species_target = Z2.len();
+    let num_incident_ions = energies.len();
+
+    assert_eq!(ux.len(), num_incident_ions, "Input error: list of x-directions is not the same length as list of incident energies.");
+    assert_eq!(uy.len(), num_incident_ions, "Input error: list of y-directions is not the same length as list of incident energies.");
+    assert_eq!(uz.len(), num_incident_ions, "Input error: list of z-directions is not the same length as list of incident energies.");
+    assert_eq!(Z1.len(), num_incident_ions, "Input error: list of incident atomic numbers is not the same length as list of incident energies.");
+    assert_eq!(m1.len(), num_incident_ions, "Input error: list of incident atomic masses is not the same length as list of incident energies.");
+    assert_eq!(Es1.len(), num_incident_ions, "Input error: list of incident surface binding energies is not the same length as list of incident energies.");
+    assert_eq!(Ec1.len(), num_incident_ions, "Input error: list of incident cutoff energies is not the same length as list of incident energies.");
+
+    assert_eq!(m2.len(), num_species_target, "Input error: list of target atomic masses is not the same length as atomic numbers.");
+    assert_eq!(Ec2.len(), num_species_target, "Input error: list of target cutoff energies is not the same length as atomic numbers.");
+    assert_eq!(Es2.len(), num_species_target, "Input error: list of target surface binding energies is not the same length as atomic numbers.");
+    assert_eq!(Eb2.len(), num_species_target, "Input error: list of target bulk binding energies is not the same length as atomic numbers.");
+    assert_eq!(n2.len(), num_species_target, "Input error: list of target number densities is not the same length as atomic numbers.");
+
+    #[cfg(feature = "distributions")]
+    let options = Options {
+        name: "test".to_string(),
+        track_trajectories: false,
+        track_recoils: true,
+        track_recoil_trajectories: false,
+        write_buffer_size: 8000,
+        weak_collision_order: 3,
+        suppress_deep_recoils: true,
+        high_energy_free_flight_paths: true,
+        electronic_stopping_mode: ElectronicStoppingMode::LOW_ENERGY_NONLOCAL,
+        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::DEFAULTNEWTON]],
+        track_displacements: false,
+        track_energy_losses: false,
+        energy_min: 0.0,
+        energy_max: 10.0,
+        energy_num: 11,
+        angle_min: 0.0,
+        angle_max: 90.0,
+        angle_num: 11,
+        x_min: 0.0,
+        y_min: -10.0,
+        z_min: -10.0,
+        x_max: 10.0,
+        y_max: 10.0,
+        z_max: 10.0,
+        x_num: 11,
+        y_num: 11,
+        z_num: 11,
+    };
+
+    #[cfg(not(feature = "distributions"))]
+    let options = Options {
+        name: "test".to_string(),
+        track_trajectories: false,
+        track_recoils: true,
+        track_recoil_trajectories: false,
+        write_buffer_size: 8000,
+        weak_collision_order: 3,
+        suppress_deep_recoils: true,
+        high_energy_free_flight_paths: false,
+        electronic_stopping_mode: ElectronicStoppingMode::LOW_ENERGY_NONLOCAL,
+        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::DEFAULTNEWTON]],
+        track_displacements: false,
+        track_energy_losses: false,
+    };
+
+    let x = -2.*(n2.iter().sum::<f64>()*10E30).powf(-1./3.);
+    let y = 0.0;
+    let z = 0.0;
+
+    let material_parameters = material::MaterialParameters {
+        energy_unit: "EV".to_string(),
+        mass_unit: "AMU".to_string(),
+        Eb: Eb2,
+        Es: Es2,
+        Ec: Ec2,
+        Z: Z2,
+        m: m2,
+        interaction_index: vec![0; num_species_target],
+        surface_binding_model: SurfaceBindingModel::INDIVIDUAL,
+        bulk_binding_model: BulkBindingModel::INDIVIDUAL,
+    };
+
+    let geometry_input = geometry::Mesh0DInput {
+        length_unit: "ANGSTROM".to_string(),
+        densities: n2,
+        electronic_stopping_correction_factor: 1.0
+    };
+
+    let m = material::Material::<Mesh0D>::new(&material_parameters, &geometry_input);
+
+    let mut index: usize = 0;
+    for (((((((E1_, ux_), uy_), uz_), Z1_), Ec1_), Es1_), m1_) in energies.iter().zip(ux).zip(uy).zip(uz).zip(Z1).zip(Ec1).zip(Es1).zip(m1) {
+
+        let mut energy_out;
+        let p = particle::Particle {
+            m: m1_*AMU,
+            Z: Z1_,
+            E: E1_*EV,
+            Ec: Ec1_*EV,
+            Es: Es1_*EV,
+            pos: Vector::new(x, y, z),
+            dir: Vector::new(ux_, uy_, uz_),
+            pos_origin: Vector::new(x, y, z),
+            pos_old: Vector::new(x, y, z),
+            dir_old: Vector::new(ux_, uy_, uz_),
+            energy_origin: E1_*EV,
+            asymptotic_deflection: 0.0,
+            stopped: false,
+            left: false,
+            incident: true,
+            first_step: true,
+            trajectory: vec![],
+            energies: vec![],
+            track_trajectories: false,
+            number_collision_events: 0,
+            backreflected: false,
+            interaction_index : 0,
+            weight: 0.0,
+            tag: 0,
+            tracked_vector: Vector::new(0., 0., 0.),
+        };
+
+        let output = bca::single_ion_bca(p, &m, &options);
+
+        for particle in output {
+            if (particle.left) | (particle.incident) {
+
+                incident.push(particle.incident);
+
+                if particle.stopped {
+                    energy_out = 0.
+                } else {
+                    energy_out = particle.E/EV
+                }
+                total_output.push(
+                    [
+                        particle.Z,
+                        particle.m/AMU,
+                        energy_out,
+                        particle.pos.x,
+                        particle.pos.y,
+                        particle.pos.z,
+                        particle.dir.x,
+                        particle.dir.y,
+                        particle.dir.z,
+                    ]
+                );
+            }
+        }
+        index += 1;
+    }
+    (total_output, incident)
+}
+
 #[cfg(feature = "python")]
 /// simple_bca_py( x, y, z, ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2)
 /// --

src/material.rs

@@ -323,7 +323,6 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
 
     //Actual surface binding energies
     let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old, z_old);
-    //println!("Actual Es: {}", Es);
     let Ec = particle_1.Ec;
 
     let inside_now = material.inside_energy_barrier(x, y, z);
@@ -380,7 +379,6 @@ pub fn surface_binding_energy<T: Geometry>(particle_1: &mut particle::Particle,
                 particle_1.dir.x = -2.*(costheta)*dx/mag + cosx;
                 particle_1.dir.y = -2.*(costheta)*dy/mag + cosy;
                 particle_1.dir.z = -2.*(costheta)*dz/mag + cosz;
-
                 particle_1.backreflected = true;
             }
         }