Added layered target version of Python functions.

Author: drobnyjt

examples/test_rustbca.py

@@ -81,6 +81,46 @@ def main():
     print(f'RustBCA R: {len(reflected[:, 0])/number_ions} Thomas R: {thomas}')
     print(f'Time per ion: {delta_time/number_ions*1e3} us/{ion["symbol"]}')
 
+    #layered target version
+
+    #1 keV is above the He on W sputtering threshold of ~150 eV
+    energies_eV = 1000.0*np.ones(number_ions)
+
+    #Working with angles of exactly 0 is problematic due to gimbal lock
+    angle = 0.0001
+
+    #In RustBCA's 0D geometry, +x -> into the surface
+    ux = np.cos(angle*np.pi/180.)*np.ones(number_ions)
+    uy = np.sin(angle*np.pi/180.)*np.ones(number_ions)
+    uz = np.zeros(number_ions)
+
+    print(f'Running RustBCA for {number_ions} {ion["symbol"]} ions on {target["symbol"]} with hydrogenated layer at {energies_eV[0]/1000.} keV...')
+    print(f'This may take several minutes.')
+    output, incident, stopped = compound_bca_list_1D_py(
+        ux, uy, uz, energies_eV, [ion['Z']]*number_ions,
+        [ion['m']]*number_ions, [ion['Ec']]*number_ions, [ion['Es']]*number_ions, [target['Z'], 1.0], [target['m'], 1.008],
+        [target['Ec'], 1.0], [target['Es'], 1.5], [target['Eb'], 0.0], [[target['n']/10**30, target['n']/10**30], [target['n']/10**30, 0.0]], [50.0, 1e6]
+    )
+
+    output = np.array(output)
+
+    Z = output[:, 0]
+    m = output[:, 1]
+    E = output[:, 2]
+    x = output[:, 3]
+    y = output[:, 4]
+    z = output[:, 5]
+    ux = output[:, 6]
+    uy = output[:, 7]
+    uz = output[:, 8]
+
+    plt.figure(3)
+    plt.title(f'Implanted {ion["symbol"]} Depth Distribution with 50A {target["symbol"]}-H layer')
+    plt.xlabel('x [A]')
+    plt.ylabel('f(x) [A.U.]')
+    heights, _, _ = plt.hist(x[np.logical_and(incident, stopped)], bins=100, density=True, histtype='step')
+    plt.plot([50.0, 50.0], [0.0, np.max(heights)*1.1])
+    plt.gca().set_ylim([0.0, np.max(heights)*1.1])
     plt.show()
 
 if __name__ == '__main__':

src/RustBCA.egg-info/PKG-INFO

@@ -0,0 +1,4 @@
+Metadata-Version: 2.1
+Name: RustBCA
+Version: 1.2.0
+License-File: LICENSE

src/RustBCA.egg-info/SOURCES.txt

@@ -0,0 +1,28 @@
+Cargo.toml
+LICENSE
+MANIFEST.in
+README.md
+pyproject.toml
+setup.py
+src/bca.rs
+src/consts.rs
+src/enums.rs
+src/geometry.rs
+src/gui.rs
+src/input.rs
+src/interactions.rs
+src/lib.rs
+src/main.rs
+src/material.rs
+src/output.rs
+src/parry.rs
+src/particle.rs
+src/physics.rs
+src/sphere.rs
+src/structs.rs
+src/tests.rs
+src/RustBCA.egg-info/PKG-INFO
+src/RustBCA.egg-info/SOURCES.txt
+src/RustBCA.egg-info/dependency_links.txt
+src/RustBCA.egg-info/not-zip-safe
+src/RustBCA.egg-info/top_level.txt

src/RustBCA.egg-info/dependency_links.txt

@@ -0,0 +1 @@
+

src/RustBCA.egg-info/not-zip-safe

@@ -0,0 +1 @@
+

src/RustBCA.egg-info/top_level.txt

@@ -0,0 +1 @@
+

src/lib.rs

@@ -81,6 +81,7 @@ 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)?)?;
+    m.add_function(wrap_pyfunction!(compound_bca_list_1D_py, m)?)?;
     Ok(())
 }
 
@@ -1173,6 +1174,212 @@ pub fn compound_bca_list_py(energies: Vec<f64>, ux: Vec<f64>, uy: Vec<f64>, uz:
     (total_output, incident)
 }
 
+#[cfg(feature = "python")]
+///compound_bca_list_1D_py(ux, uy,  uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, Eb2 n2, dx)
+/// 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.
+///    Eb2 (list(f64)): target material species bulk binding energies in eV.
+///    n2 (list(list(f64))): target material species atomic number densities in inverse cubic Angstroms.
+///    dx (list(f64)): target material layer thicknesses starting at surface.
+/// 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
+/// stopped (list(bool)): whether each row of output is associated with a particle that stopped in the target
+#[pyfunction]
+pub fn compound_bca_list_1D_py(ux: Vec<f64>, uy: Vec<f64>, uz: Vec<f64>, energies: 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>, Eb2: Vec<f64>, n2: Vec<Vec<f64>>,  dx: Vec<f64>) -> (Vec<[f64; 9]>, Vec<bool>, Vec<bool>) {
+    let mut total_output = vec![];
+    let mut incident = vec![];
+    let mut stopped = vec![];
+
+    let num_layers_target = n2.len();
+    let num_species = 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, "Input error: list of target atomic masses is not the same length as atomic numbers.");
+    assert_eq!(Ec2.len(), num_species, "Input error: list of target cutoff energies is not the same length as atomic numbers.");
+    assert_eq!(Es2.len(), num_species, "Input error: list of target surface binding energies is not the same length as atomic numbers.");
+    assert_eq!(Eb2.len(), num_species, "Input error: list of target bulk binding energies is not the same length as atomic numbers.");
+    assert_eq!(n2[0].len(), num_species, "Input error: first layer list of target number densities is not the same length as atomic numbers.");
+
+    assert_eq!(n2[0].len(), num_species, "Input error: first layer species list of target number densities is not the same length as atomic numbers.");
+    assert_eq!(dx.len(), num_layers_target, "Input error: number of layer thicknesses not the same as number of layers in atomic densities list.");
+
+    #[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 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],
+        surface_binding_model: SurfaceBindingModel::INDIVIDUAL,
+        bulk_binding_model: BulkBindingModel::INDIVIDUAL,
+    };
+
+    let geometry_input = geometry::Mesh1DInput {
+        length_unit: "ANGSTROM".to_string(),
+        densities: n2,
+        layer_thicknesses: dx,
+        electronic_stopping_correction_factors: vec![1.0; num_layers_target]
+    };
+
+    let m = material::Material::<Mesh1D>::new(&material_parameters, &geometry_input);
+
+    let x = -m.geometry.top_energy_barrier_thickness/2.;
+
+    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 dir = Vector::new(ux_, uy_, uz_);
+        dir.normalize();
+
+        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: dir,
+            pos_origin: Vector::new(x, y, z),
+            pos_old: Vector::new(x, y, z),
+            dir_old: dir,
+            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);
+                stopped.push(particle.stopped);
+
+                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/ANGSTROM,
+                        particle.pos.y/ANGSTROM,
+                        particle.pos.z/ANGSTROM,
+                        particle.dir.x,
+                        particle.dir.y,
+                        particle.dir.z,
+                    ]
+                );
+            }
+        }
+        index += 1;
+    }
+    (total_output, incident, stopped)
+}
+
 #[cfg(feature = "python")]
 /// simple_bca_py( x, y, z, ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2)
 /// --