Merge dev into main

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "RustBCA"
-version = "2.8.4"
+version = "2.9.0"
 default-run = "RustBCA"
 authors = ["Jon Drobny <drobny2@illinois.edu>", "Jon Drobny <jdrobny@tae.com>"]
 edition = "2021"
@@ -20,7 +20,7 @@ rand_distr = "0.4.3"
 toml = "0.7.4"
 anyhow = "1.0.71"
 itertools = "0.10.5"
-rayon = "1.7.0"
+rayon = "1.10.0"
 geo = {version = "0.25", optional = false}
 indicatif = {version = "0.15.0", features=["rayon"]}
 serde = { version = "1.0.163", features = ["derive"] }

examples/boron_nitride_0D.toml

@@ -1,3 +1,4 @@
+#Use with 0D geometry option only
 [options]
 name = "boron_nitride_"
 weak_collision_order = 0

examples/boron_nitride_sphere.toml

@@ -1,3 +1,4 @@
+#Use with SPHERE geometry option only
 [options]
 name = "boron_nitride_"
 track_trajectories = true

examples/boron_nitride_wire.toml

@@ -1,3 +1,4 @@
+#Use with 2D geometry option only
 [options]
 name = "boron_nitride_"
 weak_collision_order = 0

examples/boron_nitride_wire_homogeneous.toml

@@ -1,3 +1,4 @@
+#Use with 0D geometry option only
 [options]
 name = "boron_nitride_h_"
 weak_collision_order = 0

examples/layered_geometry.toml

@@ -1,3 +1,4 @@
+#Use with 2D geometry option only
 [options]
 name = "2000.0eV_0.0001deg_He_TiO2_Al_Si"
 track_trajectories = false

examples/layered_geometry_1D.toml

@@ -1,3 +1,4 @@
+#Use with 2D geometry option only
 [options]
 name = "2000.0eV_0.0001deg_He_TiO2_Al_Si"
 track_trajectories = false

examples/lithium_vapor_shield.toml

@@ -1,3 +1,4 @@
+#Use with 1D geometry option only
 [options]
 name = "lithium_vapor_shield_"
 track_trajectories = true

examples/multiple_interaction_potentials.toml

@@ -1,3 +1,4 @@
+#Use with 2D geometry input only
 [options]
 name = "2000.0eV_0.0001deg_He_H_TiO2_Al_Si"
 track_trajectories = true

examples/test_morse.py

@@ -180,8 +180,8 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True):
 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=False)
-    R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=False)
+    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)
 
 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')

examples/test_rustbca.py

@@ -12,6 +12,7 @@
 
 def main():
 
+
     #test rotation to and from RustBCA coordinates
 
     #nx, ny, nz is the normal vector (out of surface)
@@ -35,6 +36,31 @@ def main():
     assert(abs(uy) < 1e-6)
     assert(abs(uz) < 1e-6)
 
+    #test vectorized rotation to and from RustBCA coordinates
+
+    num_rot_test = 1000000
+
+    #nx, ny, nz is the normal vector (out of surface)
+    nx = [-np.sqrt(2)/2]*num_rot_test
+    ny = [-np.sqrt(2)/2]*num_rot_test
+    nz = [0.0]*num_rot_test
+
+    #ux, uy, uz is the particle direction (simulation coordinates)
+    ux = [1.0]*num_rot_test
+    uy = [0.0]*num_rot_test
+    uz = [0.0]*num_rot_test
+
+    start = time.time()
+    ux, uy, uz = rotate_given_surface_normal_vec_py(nx, ny, nz, ux, uy, uz)
+    ux, uy, uz = rotate_back_vec_py(nx, ny, nz, ux, uy, uz)
+    stop = time.time()
+    print(f'Time to rotate: {(stop - start)/num_rot_test} sec/vector')
+
+    #After rotating and rotating back, effect should be where you started (minus fp error)
+    assert(abs(ux[0] - 1.0) < 1e-6)
+    assert(abs(uy[0]) < 1e-6)
+    assert(abs(uz[0]) < 1e-6)
+
     #scripts/materials.py has a number of potential ions and targets
     ion = helium
     ion['Eb'] = 0.0
@@ -128,7 +154,7 @@ def main():
     print('Data processing complete.')
     print(f'RustBCA Y: {len(sputtered[:, 0])/number_ions} Yamamura Y: {yamamura_yield}')
     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"]}')
+    print(f'Time per ion: {delta_time/number_ions} s/{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.
 
@@ -145,14 +171,15 @@ def main():
 
     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
+    #Note 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]
@@ -173,6 +200,80 @@ def main():
     plt.plot([50.0, 50.0], [0.0, np.max(heights)*1.1])
     plt.gca().set_ylim([0.0, np.max(heights)*1.1])
 
+    number_ions = 10000
+
+    #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)
+
+    Z1 = np.array([ion['Z']]*number_ions)
+    m1 = np.array([ion['m']]*number_ions)
+    Ec1 = np.array([ion['Ec']]*number_ions)
+    Es1 = np.array([ion['Es']]*number_ions)
+
+    print(f'Running tracked RustBCA for {number_ions} {ion["symbol"]} ions on {target["symbol"]} with 1% {ion["symbol"]} at {energies_eV[0]/1000.} keV...')
+    print(f'This may take several minutes.')
+
+    start = time.time()
+    # Note that simple_bca_list_py expects number densities in 1/Angstrom^3
+    output, incident, incident_index = compound_bca_list_tracked_py(
+        energies_eV, ux, uy, uz, Z1, m1, Ec1, Es1,
+        [target['Z'], ion['Z']], [target['m'], ion['m']],
+        [target['Ec'], ion['Ec']], [target['Es'], ion['Es']], [target['n']/10**30, 0.01*target['n']/10**30],
+        [target['Eb'], 0.0])
+    stop = time.time()
+    delta_time = stop - start
+
+    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]
+
+    #For the python bindings, these conditionals can be used to distinguish
+    #between sputtered, reflected, and implanted particles in the output list
+    sputtered = output[np.logical_and(Z == target['Z'], E > 0), :]
+    reflected = output[np.logical_and(Z == ion['Z'], x < 0), :]
+    implanted = output[np.logical_and(Z == ion['Z'], x > 0), :]
+
+    plt.figure(1)
+    plt.title(f'Sputtered {target["symbol"]} Energy Distribution')
+    plt.xlabel('E [eV]')
+    plt.ylabel('f(E) [A.U.]')
+    plt.hist(sputtered[:, 2], bins=100, density=True, histtype='step')
+
+    plt.figure(2)
+    plt.title(f'Implanted {ion["symbol"]} Depth Distribution')
+    plt.xlabel('x [A]')
+    plt.ylabel('f(x) [A.U.]')
+    plt.hist(implanted[:, 3], bins=100, density=True, histtype='step')
+
+    thomas = thomas_reflection(ion, target, energies_eV[0])
+    yamamura_yield = yamamura(ion, target, energies_eV[0])
+
+    print('Data processing complete.')
+    print(f'RustBCA Y: {len(sputtered[:, 0])/number_ions} Yamamura Y: {yamamura_yield}')
+    print(f'RustBCA R: {len(reflected[:, 0])/number_ions} Thomas R: {thomas}')
+    print(f'Time per ion: {delta_time/number_ions} s/{ion["symbol"]}')
+
+    plt.figure()
+    plt.plot(incident_index)
+    plt.xlabel('Particle number')
+    plt.ylabel('Particle index')
+    plt.legend(['Incident', 'Indicies'])
+
     plt.show()
 
 if __name__ == '__main__':

examples/titanium_dioxide_0D.toml

@@ -1,3 +1,4 @@
+#Use with 0D geometry option only
 [options]
 name = "10.0"
 track_trajectories = false

examples/tungsten_tiles_trimesh.toml

@@ -1,3 +1,4 @@
+#Use with TRIMESH geometry option only
 [options]
 name = "tungsten_tiles_"
 track_trajectories = true

examples/tungsten_twist_trimesh.toml

@@ -1,3 +1,4 @@
+#Use with TRIMESH geometry option only
 [options]
 name = "tungsten_twist_"
 track_trajectories = true

src/geometry.rs

@@ -66,6 +66,7 @@ impl Geometry for Mesh0D {
         let electronic_stopping_correction_factor = input.electronic_stopping_correction_factor;
 
         let densities: Vec<f64> = input.densities.iter().map(|element| element/(length_unit).powi(3)).collect();
+        assert!(densities.len() > 0, "Input Error: density list empty.");
 
         let total_density: f64 = densities.iter().sum();
 
@@ -135,6 +136,7 @@ impl Geometry for Mesh1D {
 
         let layer_thicknesses = geometry_input.layer_thicknesses.clone();
         let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone();
+        assert!(electronic_stopping_correction_factors.len() > 0, "Input Error: Electronic stopping correction factor list empty.");
         let n = layer_thicknesses.len();
 
         let mut layers: Vec<Layer1D> =  Vec::with_capacity(n);
@@ -433,7 +435,7 @@ impl Geometry for Mesh2D {
 
         let simulation_boundary_points = geometry_input.simulation_boundary_points.clone();
         let electronic_stopping_correction_factors = geometry_input.electronic_stopping_correction_factors.clone();
-
+        assert!(electronic_stopping_correction_factors.len() > 0, "Input Error: Electronic stopping correction factor list empty.");
         let n = triangles.len();
 
         let mut cells: Vec<Cell2D> =  Vec::with_capacity(n);