Merge pull request #288 from lcpp-org/recoil_generation

drobnyjt · web-flow · commit 5b013197d211 · 2025-08-04T15:27:36.000-07:00
Recoil generation fix for z-directed particles
diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml
@@ -28,7 +28,7 @@ E = [ 1000.0 ]
 Ec = [ 1.0 ]
 Es = [ 10.0 ]
 pos = [ [ 0.0, 0.0, 0.0,] ]
-dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ]
+dir = [ [ 1.0, 0.0, 0.0,] ]
 
 [geometry_input]
 length_unit = "ANGSTROM"
diff --git a/examples/boron_nitride_sphere.toml b/examples/boron_nitride_sphere.toml
@@ -31,7 +31,7 @@ E = [ 500.0 ]
 Ec = [ 1.0 ]
 Es = [ 10.0 ]
 pos = [ [ -100.0, 0.0, 0.0,] ]
-dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ]
+dir = [ [ 1.0, 0.0, 0.0,] ]
 interaction_index = [ 0 ]
 particle_input_filename=""
 
diff --git a/examples/boron_nitride_wire.toml b/examples/boron_nitride_wire.toml
@@ -24,7 +24,7 @@ E = [ 5000.0 ]
 Ec = [ 1.0 ]
 Es = [ 10.0 ]
 pos = [ [ 0.0, 0.0, 0.0,] ]
-dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ]
+dir = [ [ 1.0, 0.0, 0.0,] ]
 
 [geometry_input]
 length_unit = "ANGSTROM"
diff --git a/examples/boron_nitride_wire_homogeneous.toml b/examples/boron_nitride_wire_homogeneous.toml
@@ -24,7 +24,7 @@ E = [ 5000.0 ]
 Ec = [ 1.0 ]
 Es = [ 10.0 ]
 pos = [ [ 0.0, 0.0, 0.0,] ]
-dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,] ]
+dir = [ [ 1.0, 0.0, 0.0,] ]
 
 [geometry_input]
 length_unit = "ANGSTROM"
diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml
@@ -46,7 +46,7 @@ Ec = [ 1.0,]
 Es = [ 0.0,]
 interaction_index = [0]
 pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,],]
-dir = [ [ 0.999, 0.001, 0.0,],]
+dir = [ [ 1.0, 0.0, 0.0,],]
 
 [geometry_input]
 length_unit = "MICRON"
diff --git a/examples/layered_geometry_1D.toml b/examples/layered_geometry_1D.toml
@@ -46,7 +46,7 @@ Ec = [ 1.0,]
 Es = [ 0.0,]
 interaction_index = [0]
 pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,],]
-dir = [ [ 0.999, 0.001, 0.0,],]
+dir = [ [ 1.0, 0.0, 0.0,],]
 
 [geometry_input]
 length_unit = "MICRON"
diff --git a/examples/lithium_vapor_shield.toml b/examples/lithium_vapor_shield.toml
@@ -46,7 +46,7 @@ Ec = [ 1.0,]
 Es = [ 0.0,]
 interaction_index = [0]
 pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,],]
-dir = [ [ 0.999, 0.001, 0.0,],]
+dir = [ [ 1.0, 0.0, 0.0,],]
 
 [geometry_input]
 length_unit = "MICRON"
diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml
@@ -40,7 +40,7 @@ Ec = [ 0.1, 0.1]
 Es = [ 0.0, 1.0]
 interaction_index = [0, 1]
 pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,], [ -1.7453292519934434e-8, 0.0, 0.0,],]
-dir = [ [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,], [ 0.9999999999984769, 1.7453292519934434e-6, 0.0,],]
+dir = [ [ 1.0, 0.0, 0.0,], [ 1.0, 0.0, 0.0,],]
 
 [geometry_input]
 energy_barrier_thickness = 2.2E-4
diff --git a/examples/test_cube.py b/examples/test_cube.py
@@ -115,10 +115,10 @@ def main():
     run_sim = True
     track_trajectories = False
     a = 1000
-    num_samples = 100000
+    num_samples = 10000
     num_angles = 10
-    energy = 100
-    angles = np.linspace(0.01, 89.9, num_angles)
+    energy = 200
+    angles = np.linspace(0.0, 89.9, num_angles)
     
     sputtering_yields_x = np.zeros(num_angles)
     sputtering_yields_z = np.zeros(num_angles)
@@ -138,43 +138,48 @@ def main():
 
     sim_index = 0
     impl_index = num_angles*3//4
+    bins = np.linspace(0, 500, 100)
+    
     for angle_index, angle in enumerate(angles):
 
         R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, ux=np.cos(angle*np.pi/180.), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2))
         sim_index += 1
         plt.figure(1)
-        if angle_index == impl_index: plt.hist(implanted[:, 2], histtype='step', bins=100, density=False, label='x')
-
+        if angle_index == impl_index: plt.hist(implanted[:, 2], histtype='step', bins=bins, density=False, label='x')
         sputtering_yields_x[angle_index] = Y
         R_x[angle_index] = R
 
         R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a, ux=-np.cos(angle*np.pi/180.), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2))
         sim_index += 1
-        if angle_index == impl_index: plt.hist(a - implanted[:, 2], histtype='step', bins=100, density=False, label='-x')
+        if angle_index == impl_index: plt.hist(a - implanted[:, 2], histtype='step', bins=bins, density=False, label='-x')
         sputtering_yields_x_minus[angle_index] = Y
         R_x_minus[angle_index] = R
 
         R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=a/2., z0=a, ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=-np.cos(angle*np.pi/180.))
         sim_index += 1
-        if angle_index == impl_index: plt.hist(a - implanted[:, 4], histtype='step', bins=100, density=False, label='-z')
+        if angle_index == impl_index: plt.hist(a - implanted[:, 4], histtype='step', bins=bins, density=False, label='-z')
+
         sputtering_yields_z_minus[angle_index] = Y
         R_z_minus[angle_index] = R
 
         R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=a/2., z0=0.0, ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.cos(angle*np.pi/180.))
         sim_index += 1
-        if angle_index == impl_index: plt.hist(implanted[:, 4], histtype='step', bins=100, density=False, label='z')
+        if angle_index == impl_index: plt.hist(implanted[:, 4], histtype='step', bins=bins, density=False, label='z')
+
         sputtering_yields_z[angle_index] = Y
         R_z[angle_index] = R
 
         R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=0.0, z0=a/2., ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=np.cos(angle*np.pi/180.))
         sim_index += 1
-        if angle_index == impl_index: plt.hist(implanted[:, 3], histtype='step', bins=100, density=False, label='y')
+        if angle_index == impl_index: plt.hist(implanted[:, 3], histtype='step', bins=bins, density=False, label='y')
+
         sputtering_yields_y[angle_index] = Y
         R_y[angle_index] = R
 
         R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=a, z0=a/2., ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=-np.cos(angle*np.pi/180.))
         sim_index += 1
-        if angle_index == impl_index: plt.hist(a - implanted[:, 3], histtype='step', bins=100, density=False, label='-y')
+        if angle_index == impl_index: plt.hist(a - implanted[:, 3], histtype='step', bins=bins, density=False, label='-y')
+
         sputtering_yields_y_minus[angle_index] = Y
         R_y_minus[angle_index] = R
 
@@ -206,7 +211,6 @@ def main():
         plt.show()
 
     if track_trajectories and do_plots:
-
         do_trajectory_plot('cube_19', boundary=[[0.0, 0.0], [0.0, a], [a, a], [a, 0.0], [0.0, 0.0]])
 
         do_trajectory_plot('cube_20', boundary=[[0.0, 0.0], [0.0, a], [a, a], [a, 0.0], [0.0, 0.0]])
diff --git a/examples/tungsten_twist_trimesh.toml b/examples/tungsten_twist_trimesh.toml
@@ -42,7 +42,7 @@ E = [ 2000.0 ]
 Ec = [ 1.0 ]
 Es = [ 0.0 ]
 pos = [ [ 0.0, 0.0, 1.0,] ]
-dir = [ [ 0.0, 0.0001, -0.9999,] ]
+dir = [ [ 0.0, 0.0, -1.0,] ]
 interaction_index = [ 0 ]
 
 [geometry_input]
diff --git a/src/bca.rs b/src/bca.rs
@@ -127,9 +127,13 @@ pub fn single_ion_bca<T: Geometry>(particle: particle::Particle, material: &mate
                     total_energy_lost_to_recoils += binary_collision_result.recoil_energy;
                     total_asymptotic_deflection += binary_collision_result.asymptotic_deflection;
 
+                    // Rotate each particle in the appropriate plane by psi/psi_r
                     particle_1.rotate(binary_collision_result.psi,
                         binary_collision_geometry.phi_azimuthal);
 
+                    // Since psi and psi_r are absolute-valued in the BC calculation, note
+                    // that it is negated here to correctly rotate the recoil in the correct
+                    // direction (that is, away from the direction the incident atom is deflected)
                     particle_2.rotate(-binary_collision_result.psi_recoil,
                         binary_collision_geometry.phi_azimuthal);
 
@@ -360,10 +364,28 @@ pub fn choose_collision_partner<T: Geometry>(particle_1: &particle::Particle, ma
     let sinx: f64 = (1. - cosx*cosx).sqrt();
     let cosphi: f64 = phi_azimuthal.cos();
 
-    //Find recoil location
-    let x_recoil: f64 = x + mfp*cosx - impact_parameter*cosphi*sinx;
-    let y_recoil: f64 = y + mfp*cosy - impact_parameter*(sinphi*cosz - cosphi*cosy*cosx)/sinx;
-    let z_recoil: f64 = z + mfp*cosz + impact_parameter*(sinphi*cosy - cosphi*cosx*cosz)/sinx;
+    // These formulas find the recoil one mfp away at an impact parameter p at angle phi
+    // To resolve the singularity, a different set of rotations is used when cosx == -1
+    // Because of this, the recoil location is not consistent between the two formulas at a given phi
+    // Since phi is sampled uniformly from (0, 2pi), this does not matter
+    // However, if a crystalline structure is ever added, this needs to be considered
+    let x_recoil = if cosx > -1. {
+        x + mfp*cosx - impact_parameter*(cosz*sinphi + cosy*cosphi)
+    } else {
+        x + mfp*cosx - impact_parameter*((1. + cosz - cosx*cosx)*cosphi - cosx*cosy*sinphi)/(1. + cosz)
+    };
+
+    let y_recoil = if cosx > -1. {
+        y + mfp*cosy + impact_parameter*((1. + cosx - cosy*cosy)*cosphi - cosy*cosz*sinphi)/(1. + cosx)
+    } else {
+        y + mfp*cosy + impact_parameter*((1. + cosz - cosy*cosy)*sinphi - cosx*cosy*cosphi)/(1. + cosz)
+    };
+
+    let z_recoil = if cosx > -1. {
+        z + mfp*cosz + impact_parameter*((1. + cosx - cosz*cosz)*sinphi - cosy*cosz*cosphi)/(1. + cosx)
+    } else {
+        z + mfp*cosz + impact_parameter*(cosx*cosphi + cosy*sinphi)
+    };
 
     //Choose recoil Z, M
     let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, Ed_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil);
@@ -504,8 +526,9 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: &
         InteractionPotential::COULOMB{..} => 0.,
         _ => x0*a*(theta/2.).sin()
     };
-    let psi = (theta.sin().atan2(Ma/Mb + theta.cos())).abs();
-    let psi_recoil = (theta.sin().atan2(1. - theta.cos())).abs();
+
+    let psi = theta.sin().atan2(Ma/Mb + theta.cos());//.abs();
+    let psi_recoil = theta.sin().atan2(1. - theta.cos());//.abs();
     let recoil_energy = 4.*(Ma*Mb)/(Ma + Mb).powi(2)*E0*(theta/2.).sin().powi(2);
 
     Ok(BinaryCollisionResult::new(theta, psi, psi_recoil, recoil_energy, asymptotic_deflection, x0))
diff --git a/src/input.rs b/src/input.rs
@@ -581,7 +581,7 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
                             ux: match cosx {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max);  uniform.sample(&mut rand::thread_rng())*length_unit},
-                                Distributions::POINT(ux) => {assert!(ux != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); ux}
+                                Distributions::POINT(ux) => ux,
                             },
                             uy: match cosy {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
@@ -654,7 +654,7 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
                             ux: match cosx {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
                                 Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max);  uniform.sample(&mut rand::thread_rng())*length_unit},
-                                Distributions::POINT(x) => {assert!(x != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); x*length_unit},
+                                Distributions::POINT(x) => x*length_unit
                             },
                             uy: match cosy {
                                 Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
diff --git a/src/lib.rs b/src/lib.rs
diff --git a/src/particle.rs b/src/particle.rs
diff --git a/src/tests.rs b/src/tests.rs
