Merge pull request #189 from lcpp-org/dev

Merge dev into main

.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.

examples/test_rustbca.py

@@ -81,6 +81,49 @@ 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"]}')
 
+    #Next up is the layered target version. I'll add a 50 Angstrom layer of W-H to the top of the target.
+
+    #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.')
+    #Not the different argument order; when a breaking change is due, this will
+    #be back-ported to the other bindings as well for consistency.
+    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/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/consts.rs

@@ -15,6 +15,8 @@ pub const MICRON: f64 = 1E-6;
 pub const NM: f64 = 1E-9;
 /// One centimeter in meters.
 pub const CM: f64 = 1E-2;
+/// One milimter in meters.
+pub const MM: f64 = 1E-3;
 /// Vacuum permitivity in Farads/meter.
 pub const EPS0: f64 = 8.8541878128E-12;
 /// Bohr radius in meters.

src/geometry.rs

@@ -52,6 +52,7 @@ impl Geometry for Mesh0D {
         let length_unit: f64 = match input.length_unit.as_str() {
             "MICRON" => MICRON,
             "CM" => CM,
+            "MM" => MM,
             "ANGSTROM" => ANGSTROM,
             "NM" => NM,
             "M" => 1.,
@@ -97,7 +98,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
     }
 
@@ -143,6 +143,7 @@ impl Geometry for Mesh1D {
         let length_unit: f64 = match geometry_input.length_unit.as_str() {
             "MICRON" => MICRON,
             "CM" => CM,
+            "MM" => MM,
             "ANGSTROM" => ANGSTROM,
             "NM" => NM,
             "M" => 1.,
@@ -324,6 +325,7 @@ impl Geometry for Mesh2D {
         let length_unit: f64 = match geometry_input.length_unit.as_str() {
             "MICRON" => MICRON,
             "CM" => CM,
+            "MM" => MM,
             "ANGSTROM" => ANGSTROM,
             "NM" => NM,
             "M" => 1.,

src/input.rs

@@ -381,6 +381,7 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
     let length_unit: f64 = match particle_parameters.length_unit.as_str() {
         "MICRON" => MICRON,
         "CM" => CM,
+        "MM" => MM,
         "ANGSTROM" => ANGSTROM,
         "NM" => NM,
         "M" => 1.,

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.")
     }
 
 }