From 2644514115434e6a9815e9884cead919ffc94599 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 11 Jun 2025 10:32:35 -0700 Subject: [PATCH 01/19] Draft of fixing numerical issues in rotate functions. --- src/lib.rs | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 3a9f53e..547661b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1502,10 +1502,10 @@ pub fn simple_compound_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1 #[cfg(feature = "parry3d")] #[no_mangle] pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut f64, uz: &mut f64) { - const DELTA: f64 = 1e-9; - let RUSTBCA_DIRECTION: Vector3:: = Vector3::::new(1.0, 0.0, 0.0); + //const DELTA: f64 = 1e-9; + //let RUSTBCA_DIRECTION: Vector3:: = Vector3::::new(1.0, 0.0, 0.0); - let into_surface = Vector3::new(-nx, -ny, -nz); + //let into_surface = Vector3::new(-nx, -ny, -nz); let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global @@ -1514,12 +1514,14 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu //That rotation is then applied to the particle direction, and can be undone later. //Algorithm is from here: //https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/180436#180436 - let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); - let c = into_surface.dot(&RUSTBCA_DIRECTION); - let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); + //let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); + //let c = into_surface.dot(&RUSTBCA_DIRECTION); + //let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let rotation_matrix = if (c + 1.0).abs() > 1e-6 { - Matrix3::identity() + vx + vx*vx/(1. + c) + let s = ny*ny + nz*nz; + + let rotation_matrix = if nx > 0.0 { + Matrix3::::new(0.0, -ny, -nz, ny, nz*nz/s, -ny*nz/s, nz, -ny*nz/s, ny*ny/s) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation //around a non-x axis; y is chosen arbitrarily @@ -1531,12 +1533,12 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu // ux must not be exactly 1.0 to avoid gimbal lock in RustBCA // simple_bca does not normalize direction before proceeding, must be done manually assert!( - incident.x + DELTA > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. c={} n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", - (c + 1.0).abs(), nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z + incident.x > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", + nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z ); - *ux = incident.x + DELTA; - *uy = incident.y - DELTA; + *ux = incident.x; + *uy = incident.y; *uz = incident.z; let mag = (ux.powf(2.) + uy.powf(2.) + uz.powf(2.)).sqrt(); From 34e31a5899199d0dc47afe1fee052491177ac99b Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 12 Jun 2025 11:21:52 -0700 Subject: [PATCH 02/19] Removed panic in rotate_given_surface_normal_py and replaced linear algebra with precomputed matrix from SymPy --- src/lib.rs | 23 ++++++----------------- 1 file changed, 6 insertions(+), 17 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 547661b..f0d7769 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1518,10 +1518,10 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu //let c = into_surface.dot(&RUSTBCA_DIRECTION); //let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let s = ny*ny + nz*nz; + //let s = ny*ny + nz*nz; - let rotation_matrix = if nx > 0.0 { - Matrix3::::new(0.0, -ny, -nz, ny, nz*nz/s, -ny*nz/s, nz, -ny*nz/s, ny*ny/s) + let rotation_matrix = if (1.0 - nx).abs() > 0.0 { + Matrix3::::new(1. + (-ny*ny - nz*nz)/(1. - nx), -ny, -nz, ny, -ny*ny/(1. - nx) + 1., -ny*nz/(1. - nx), nz, -ny*nz/(1. - nx), -nz*nz/(1. - nx) + 1.) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation //around a non-x axis; y is chosen arbitrarily @@ -1530,13 +1530,6 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu let incident = rotation_matrix*direction; - // ux must not be exactly 1.0 to avoid gimbal lock in RustBCA - // simple_bca does not normalize direction before proceeding, must be done manually - assert!( - incident.x > 0.0, "Error: RustBCA initial direction out of surface. Please check surface normals and incident direction. n = ({}, {}, {}) u = ({}, {}, {}), unew = ({}, {}, {})", - nx, ny, nz, ux, uy, uz, incident.x, incident.y, incident.z - ); - *ux = incident.x; *uy = incident.y; *uz = incident.z; @@ -1617,7 +1610,7 @@ pub fn rotate_given_surface_normal_vec_py(nx: Vec, ny: Vec, nz: Vec = Vector3::::new(1.0, 0.0, 0.0); - let into_surface = Vector3::new(-nx, -ny, -nz); + //let into_surface = Vector3::new(-nx, -ny, -nz); let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global @@ -1625,12 +1618,8 @@ pub extern "C" fn rotate_back(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut //into-the-surface vector (1.0, 0.0, 0.0) onto the local into-the-surface vector (negative normal w.r.t. ray origin). //That rotation is then applied to the particle direction, and can be undone later. //Algorithm is from here: - //https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/180436#180436 - let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); - let c = into_surface.dot(&RUSTBCA_DIRECTION); - let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - let rotation_matrix = if c != -1.0 { - Matrix3::identity() + vx + vx*vx/(1. + c) + let rotation_matrix = if (1.0 - nx).abs() > 0.0 { + Matrix3::::new(1. + (-ny*ny - nz*nz)/(1. - nx), -ny, -nz, ny, -ny*ny/(1. - nx) + 1., -ny*nz/(1. - nx), nz, -ny*nz/(1. - nx), -nz*nz/(1. - nx) + 1.) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation //around a non-x axis; y is chosen arbitrarily From 83d35d8d918428108bd981a6dca5446002138104 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Thu, 17 Jul 2025 13:04:49 -0700 Subject: [PATCH 03/19] Rederived particle rotation and transport routines. --- src/bca.rs | 6 ++++++ src/particle.rs | 14 +++++++++----- 2 files changed, 15 insertions(+), 5 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 7d3f013..07020f6 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -130,6 +130,7 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate particle_1.rotate(binary_collision_result.psi, binary_collision_geometry.phi_azimuthal); + //don't negate psi with updated recoil geometry particle_2.rotate(-binary_collision_result.psi_recoil, binary_collision_geometry.phi_azimuthal); @@ -361,9 +362,14 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma 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; + */ + let x_recoil = x + mfp*cosx - impact_parameter*(cosz*sinphi + cosy*cosphi); + let y_recoil = y + mfp*cosy + impact_parameter*((1. + cosx - cosy*cosy)*cosphi - cosy*cosz*sinphi)/(1. + cosx); + let z_recoil = z + mfp*cosz + impact_parameter*((1. + cosx - cosz*cosz)*sinphi - cosy*cosz*cosphi)/(1. + cosx); //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); diff --git a/src/particle.rs b/src/particle.rs index 9de3daa..3f15b3a 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -239,16 +239,20 @@ impl Particle { let cosx: f64 = self.dir.x; let cosy: f64 = self.dir.y; let cosz: f64 = self.dir.z; - let cphi: f64 = phi.cos(); - let sphi: f64 = phi.sin(); + let cphi: f64 = (phi + PI).cos(); + let sphi: f64 = (phi + PI).sin(); let sa = (1. - cosx*cosx).sqrt(); //Particle direction update formulas from original TRIDYN paper, see Moeller and Eckstein 1988 let cpsi: f64 = psi.cos(); let spsi: f64 = psi.sin(); - let cosx_new: f64 = cpsi*cosx + spsi*cphi*sa; - let cosy_new: f64 = cpsi*cosy - spsi/sa*(cphi*cosx*cosy - sphi*cosz); - let cosz_new: f64 = cpsi*cosz - spsi/sa*(cphi*cosx*cosz + sphi*cosy); + //let cosx_new: f64 = cpsi*cosx + spsi*cphi*sa; + //let cosy_new: f64 = cpsi*cosy - spsi/sa*(cphi*cosx*cosy - sphi*cosz); + //let cosz_new: f64 = cpsi*cosz - spsi/sa*(cphi*cosx*cosz + sphi*cosy); + + let cosx_new = cpsi*cosx - spsi*(cosz*sphi + cosy*cphi); + let cosy_new = cpsi*cosy + spsi*((1. + cosx - cosy*cosy)*cphi - cosy*cosz*sphi)/(1. + cosx); + let cosz_new = cpsi*cosz + spsi*((1. + cosx - cosz*cosz)*sphi - cosy*cosz*cphi)/(1. + cosx); let dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new}; From 210cbc03a1fce7a374a2ffd71a3804abb27ac667 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 18 Jul 2025 14:29:05 -0700 Subject: [PATCH 04/19] Update testing suite to fix momentum conservation tests. --- src/bca.rs | 5 ++++- src/tests.rs | 6 +++--- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 07020f6..0169b21 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -127,10 +127,13 @@ pub fn single_ion_bca(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); - //don't negate psi with updated recoil geometry + // Since psi is absolute-valued in the previous 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); diff --git a/src/tests.rs b/src/tests.rs index fd263f0..5bd4449 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -930,9 +930,9 @@ fn test_momentum_conservation() { println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); println!(); - assert!(approx_eq!(f64, initial_momentum.x, final_momentum.x, epsilon = 1E-12)); - assert!(approx_eq!(f64, initial_momentum.y, final_momentum.y, epsilon = 1E-12)); - assert!(approx_eq!(f64, initial_momentum.z, final_momentum.z, epsilon = 1E-12)); + assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.)); + assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.)); + assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.)); assert!(!particle_1.E.is_nan()); assert!(!particle_2.E.is_nan()); From f8312bd92b53c63282dcd166b0e2f62aff8d07cf Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 18 Jul 2025 17:44:34 -0700 Subject: [PATCH 05/19] Add cube direction test to testing suite. --- examples/test_cube.py | 180 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 180 insertions(+) create mode 100644 examples/test_cube.py diff --git a/examples/test_cube.py b/examples/test_cube.py new file mode 100644 index 0000000..4873765 --- /dev/null +++ b/examples/test_cube.py @@ -0,0 +1,180 @@ +from libRustBCA import * +import numpy as np +import matplotlib.pyplot as plt +import sys +import os +#This should allow the script to find materials and formulas from anywhere +sys.path.append(os.path.dirname(__file__)+'/../scripts') +sys.path.append('scripts') +from materials import * +from rustbca import * + +#This function simply contains an entire input file as a multi-line f-string to modify some inputs. +def run(energy, index, num_samples=10000, run_sim=True, a=1000, x0=0, y0=500, z0=500, ux=0.999, uy=0.01, uz=0.00, track_trajectories=False): + + input_file = f''' + #Use with TRIMESH geometry option only + [options] + name = "cube_{index}" + track_trajectories = {str(track_trajectories).lower()} + track_recoils = true + track_recoil_trajectories = {str(track_trajectories).lower()} + weak_collision_order = 0 + num_threads = 4 + num_chunks = 5 + + [material_parameters] + energy_unit = "EV" + mass_unit = "AMU" + Eb = [ 0.0,] + Es = [ {copper["Es"]},] + Ec = [ {copper["Es"]},] + Z = [ {copper["Z"]} ] + m = [ {copper["m"]} ] + interaction_index = [0] + surface_binding_model = {{"PLANAR"={{calculation="TARGET"}}}} + bulk_binding_model = "AVERAGE" + + [particle_parameters] + length_unit = "ANGSTROM" + energy_unit = "EV" + mass_unit = "AMU" + N = [ {num_samples} ] + m = [ 4.008 ] + Z = [ 2 ] + E = [ {energy} ] + Ec = [ 1.0 ] + Es = [ 0.0 ] + pos = [ [ {x0}, {y0}, {z0},] ] + dir = [ [ {ux}, {uy}, {uz},] ] + interaction_index = [ 0 ] + + [geometry_input] + length_unit = "ANGSTROM" + electronic_stopping_correction_factor = 1.0 + densities = [0.05] + vertices = [ + [0.0, 0.0, 0.0], + [0.0, 0.0, {a}], + [0.0, {a}, 0.0], + [{a}, 0.0, 0.0], + [0.0, {a}, {a}], + [{a}, 0.0, {a}], + [{a}, {a}, 0.0], + [{a}, {a}, {a}] + ] + indices = [ + #front + [3, 1, 0], + [3, 5, 1], + #bottom + [3, 2, 0], + [6, 2, 3], + #right + [3, 6, 7], + [3, 7, 5], + #left + [0, 1, 2], + [2, 1, 4], + #top + [5, 4, 7], + [1, 4, 5], + #back + [4, 7, 6], + [6, 2, 4] + ] + ''' + + with open(f'cube_{index}.toml', 'w') as file: + file.write(input_file) + + if run_sim: os.system(f'cargo run --release --features parry3d TRIMESH cube_{index}.toml') + + reflected_list = np.atleast_2d(np.genfromtxt(f'cube_{index}reflected.output', delimiter=',')) + if np.size(reflected_list) > 0: + num_reflected = np.shape(reflected_list)[0] + energy_reflected = np.sum(reflected_list[:, 2]) + else: + num_reflected = 0 + energy_reflected = 0. + + sputtered_list = np.atleast_2d(np.genfromtxt(f'cube_{index}sputtered.output', delimiter=',')) + if np.size(sputtered_list) > 0: + num_sputtered = np.shape(sputtered_list)[0] + energy_sputtered = np.sum(sputtered_list[:, 2]) + else: + num_sputtered = 0 + energy_sputtered = 0. + + implanted_list = np.atleast_2d(np.genfromtxt(f'cube_{index}deposited.output', delimiter=',')) + + return num_reflected/num_samples, num_sputtered/num_samples, reflected_list, sputtered_list, implanted_list + +def main(): + run_sim = True + a = 1000 + num_samples = 10 + num_angles = 5 + energy = 1500 + angles = np.linspace(0.01, 89.9, num_angles) + sputtering_yields_x = np.zeros(num_angles) + sputtering_yields_z = np.zeros(num_angles) + sputtering_yields_y = np.zeros(num_angles) + + sputtering_yields_x_minus = np.zeros(num_angles) + sputtering_yields_z_minus = np.zeros(num_angles) + sputtering_yields_y_minus = np.zeros(num_angles) + + track_trajectories = True + + sim_index = 0 + 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 + + sputtering_yields_x[angle_index] = Y + + 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 + + sputtering_yields_x_minus[angle_index] = Y + + 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 + + sputtering_yields_z_minus[angle_index] = Y + + 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 + + sputtering_yields_z[angle_index] = Y + + 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 + + sputtering_yields_y[angle_index] = Y + + 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 + + sputtering_yields_y_minus[angle_index] = Y + + plt.figure(2) + plt.plot(angles, sputtering_yields_x, label='+x directed') + plt.plot(angles, sputtering_yields_z, label='+z directed') + plt.plot(angles, sputtering_yields_y, label='+y directed') + plt.plot(angles, sputtering_yields_x_minus, label='-x directed') + plt.plot(angles, sputtering_yields_z_minus, label='-z directed') + plt.plot(angles, sputtering_yields_y_minus, label='-y directed') + plt.legend() + plt.show() + + if track_trajectories: + + 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]]) + +if __name__ == '__main__': + main() \ No newline at end of file From 725d48e23c113b83ef5d8fd47859e700007f14f3 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 18 Jul 2025 23:53:26 -0700 Subject: [PATCH 06/19] Additional comments and updated cube test. --- examples/test_cube.py | 72 +++++++++++++++++++++++++++++++------------ src/bca.rs | 6 ++-- 2 files changed, 56 insertions(+), 22 deletions(-) diff --git a/examples/test_cube.py b/examples/test_cube.py index 4873765..fba5c82 100644 --- a/examples/test_cube.py +++ b/examples/test_cube.py @@ -111,66 +111,100 @@ def run(energy, index, num_samples=10000, run_sim=True, a=1000, x0=0, y0=500, z0 return num_reflected/num_samples, num_sputtered/num_samples, reflected_list, sputtered_list, implanted_list def main(): - run_sim = True + do_plots = True + run_sim = False + track_trajectories = False a = 1000 - num_samples = 10 - num_angles = 5 - energy = 1500 + num_samples = 100000 + num_angles = 10 + energy = 500 angles = np.linspace(0.01, 89.9, num_angles) + sputtering_yields_x = np.zeros(num_angles) sputtering_yields_z = np.zeros(num_angles) sputtering_yields_y = np.zeros(num_angles) + R_x = np.zeros(num_angles) + R_y = np.zeros(num_angles) + R_z = np.zeros(num_angles) + sputtering_yields_x_minus = np.zeros(num_angles) sputtering_yields_z_minus = np.zeros(num_angles) sputtering_yields_y_minus = np.zeros(num_angles) - track_trajectories = True + R_x_minus = np.zeros(num_angles) + R_y_minus = np.zeros(num_angles) + R_z_minus = np.zeros(num_angles) sim_index = 0 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 == 0: plt.hist(implanted[:, 2], histtype='step', bins=100, density=True, 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 == 0: plt.hist(a - implanted[:, 2], histtype='step', bins=100, density=True, 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 == 0: plt.hist(a - implanted[:, 4], histtype='step', bins=100, density=True, 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 == 0: plt.hist(implanted[:, 4], histtype='step', bins=100, density=True, 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 == 0: plt.hist(implanted[:, 3], histtype='step', bins=100, density=True, 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 == 0: plt.hist(a - implanted[:, 3], histtype='step', bins=100, density=True, label='-y') sputtering_yields_y_minus[angle_index] = Y + R_y_minus[angle_index] = R - plt.figure(2) - plt.plot(angles, sputtering_yields_x, label='+x directed') - plt.plot(angles, sputtering_yields_z, label='+z directed') - plt.plot(angles, sputtering_yields_y, label='+y directed') - plt.plot(angles, sputtering_yields_x_minus, label='-x directed') - plt.plot(angles, sputtering_yields_z_minus, label='-z directed') - plt.plot(angles, sputtering_yields_y_minus, label='-y directed') + plt.figure(1) + plt.title('Implantation') plt.legend() - plt.show() - if track_trajectories: + if do_plots: + plt.figure(2) + plt.title('Sputtering') + plt.plot(angles, sputtering_yields_x, label='+x directed') + plt.plot(angles, sputtering_yields_z, label='+z directed') + plt.plot(angles, sputtering_yields_y, label='+y directed') + plt.plot(angles, sputtering_yields_x_minus, label='-x directed') + plt.plot(angles, sputtering_yields_z_minus, label='-z directed') + plt.plot(angles, sputtering_yields_y_minus, label='-y directed') + plt.legend() + + plt.figure(3) + plt.title('Reflection') + plt.plot(angles, R_x, label='+x directed') + plt.plot(angles, R_z, label='+z directed') + plt.plot(angles, R_y, label='+y directed') + plt.plot(angles, R_x_minus, label='-x directed') + plt.plot(angles, R_z_minus, label='-z directed') + plt.plot(angles, R_y_minus, label='-y directed') + plt.legend() + + 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]]) diff --git a/src/bca.rs b/src/bca.rs index 0169b21..2c84305 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -131,9 +131,9 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate particle_1.rotate(binary_collision_result.psi, binary_collision_geometry.phi_azimuthal); - // Since psi is absolute-valued in the previous 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) + // 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); From 71a16379d6aa642ed7c212f6c45fcf68460897ca Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 20 Jul 2025 11:38:17 -0700 Subject: [PATCH 07/19] Working singularity-free version. --- src/bca.rs | 26 ++++++++++++++++++-------- src/input.rs | 9 ++------- src/particle.rs | 29 ++++++++++++++++++++--------- src/tests.rs | 4 ++-- 4 files changed, 42 insertions(+), 26 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 2c84305..d81d235 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -365,14 +365,24 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma 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; - */ - let x_recoil = x + mfp*cosx - impact_parameter*(cosz*sinphi + cosy*cosphi); - let y_recoil = y + mfp*cosy + impact_parameter*((1. + cosx - cosy*cosy)*cosphi - cosy*cosz*sinphi)/(1. + cosx); - let z_recoil = z + mfp*cosz + impact_parameter*((1. + cosx - cosz*cosz)*sinphi - cosy*cosz*cosphi)/(1. + cosx); + + 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); diff --git a/src/input.rs b/src/input.rs index 8e2b81d..7027509 100644 --- a/src/input.rs +++ b/src/input.rs @@ -172,11 +172,6 @@ fn one_u64() -> u64 { 1 } -///This helper function is a workaround to issue #368 in serde -fn three() -> usize { - 3 -} - fn zero_usize() -> usize{ 0 } @@ -586,7 +581,7 @@ where ::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) => {assert!(true, "ux cannot equal exactly 1.0 to prevent gimbal lock"); ux} }, uy: match cosy { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, @@ -659,7 +654,7 @@ where ::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) => {assert!(true, "ux cannot equal exactly 1.0 to prevent gimbal lock"); 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/particle.rs b/src/particle.rs index 3f15b3a..786a0a0 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -103,7 +103,7 @@ impl Particle { let dir_mag = (dirx*dirx + diry*diry + dirz*dirz).sqrt(); - assert!((dirx/dir_mag).abs() < 1.0 - f64::EPSILON, "Input error: incident direction cannot round to exactly (1, 0, 0) due to gimbal lock. Use a non-zero y-component."); + //assert!((dirx/dir_mag).abs() < 1.0 - f64::EPSILON, "Input error: incident direction cannot round to exactly (1, 0, 0) due to gimbal lock. Use a non-zero y-component."); assert!(input.E > 0., "Input error: incident energy {}; must be greater than zero.", input.E/EV); Particle { @@ -239,20 +239,31 @@ impl Particle { let cosx: f64 = self.dir.x; let cosy: f64 = self.dir.y; let cosz: f64 = self.dir.z; - let cphi: f64 = (phi + PI).cos(); - let sphi: f64 = (phi + PI).sin(); + let cosphi: f64 = (phi + PI).cos(); + let sinphi: f64 = (phi + PI).sin(); let sa = (1. - cosx*cosx).sqrt(); //Particle direction update formulas from original TRIDYN paper, see Moeller and Eckstein 1988 let cpsi: f64 = psi.cos(); let spsi: f64 = psi.sin(); - //let cosx_new: f64 = cpsi*cosx + spsi*cphi*sa; - //let cosy_new: f64 = cpsi*cosy - spsi/sa*(cphi*cosx*cosy - sphi*cosz); - //let cosz_new: f64 = cpsi*cosz - spsi/sa*(cphi*cosx*cosz + sphi*cosy); - let cosx_new = cpsi*cosx - spsi*(cosz*sphi + cosy*cphi); - let cosy_new = cpsi*cosy + spsi*((1. + cosx - cosy*cosy)*cphi - cosy*cosz*sphi)/(1. + cosx); - let cosz_new = cpsi*cosz + spsi*((1. + cosx - cosz*cosz)*sphi - cosy*cosz*cphi)/(1. + cosx); + let cosx_new = if cosx > -1. { + cpsi*cosx - spsi*(cosz*sinphi + cosy*cosphi) + } else { + cpsi*cosx - spsi*((1. + cosz - cosx*cosx)*cosphi - cosx*cosy*sinphi)/(1. + cosz) + }; + + let cosy_new = if cosx > -1. { + cpsi*cosy + spsi*((1. + cosx - cosy*cosy)*cosphi - cosy*cosz*sinphi)/(1. + cosx) + } else { + cpsi*cosy + spsi*((1. + cosz - cosy*cosy)*sinphi - cosx*cosy*cosphi)/(1. + cosz) + }; + + let cosz_new = if cosx > -1. { + cpsi*cosz + spsi*((1. + cosx - cosz*cosz)*sinphi - cosy*cosz*cosphi)/(1. + cosx) + } else { + cpsi*cosz + spsi*(cosx*cosphi + cosy*sinphi) + }; let dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new}; diff --git a/src/tests.rs b/src/tests.rs index 5bd4449..ea54628 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -781,8 +781,8 @@ fn test_momentum_conservation() { //Arbitrary initial angle let theta = 0.974194583091052_f64; let cosx = (theta).cos(); - let cosy = (theta).sin(); - let cosz = 0.; + let cosy = (theta).sin()/(2.0_f64).sqrt(); + let cosz = (theta).sin()/(2.0_f64).sqrt(); let material_parameters = material::MaterialParameters{ energy_unit: "EV".to_string(), From 729d50693697b481ccbacf06188281fecc10f55d Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 07:54:17 -0700 Subject: [PATCH 08/19] Updated singularity-free version. --- src/input.rs | 4 +- src/lib.rs | 10 +- src/particle.rs | 3 - src/tests.rs | 358 ++++++++++++++++++++++++------------------------ 4 files changed, 187 insertions(+), 188 deletions(-) diff --git a/src/input.rs b/src/input.rs index 7027509..ec2c5fd 100644 --- a/src/input.rs +++ b/src/input.rs @@ -581,7 +581,7 @@ where ::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!(true, "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 ::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!(true, "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 index f0d7769..c6e865c 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -803,7 +803,7 @@ pub extern "C" fn simple_bca_c(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64 /// 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 +/// ux (list(f64)): initial ion directions x. /// uy (list(f64)): initial ion directions y. /// uz (list(f64)): initial ion directions z. /// Z1 (list(f64)): initial ion atomic numbers. @@ -924,7 +924,7 @@ pub fn compound_bca_list_py(energies: Vec, ux: Vec, uy: Vec, uz: /// 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 +/// ux (list(f64)): initial ion directions x. /// uy (list(f64)): initial ion directions y. /// uz (list(f64)): initial ion directions z. /// Z1 (list(f64)): initial ion atomic numbers. @@ -1143,7 +1143,7 @@ pub fn reflect_single_ion_py(ion: &PyDict, target: &PyDict, vx: f64, vy: f64, vz ///compound_bca_list_1D_py(ux, uy, uz, energies, 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: -/// ux (list(f64)): initial ion directions x. ux != 0.0 to avoid gimbal lock +/// ux (list(f64)): initial ion directions x. /// uy (list(f64)): initial ion directions y. /// uz (list(f64)): initial ion directions z. /// energies (list(f64)): initial ion energies in eV. @@ -1279,7 +1279,7 @@ pub fn compound_bca_list_1D_py(ux: Vec, uy: Vec, uz: Vec, energie /// x (f64): initial ion position x. Material target is x>0 /// y (f64): initial ion position y. /// z (f64): initial ion position z. -/// ux (f64): initial ion direction x. ux != 0.0 to avoid gimbal lock +/// ux (f64): initial ion direction x. /// uy (f64): initial ion direction y. /// uz (f64): initial ion direction z. /// energy (f64): initial ion energy in eV. @@ -1309,7 +1309,7 @@ pub fn simple_bca_py(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, /// This function runs a 0D Binary Collision Approximation simulation for the given incident ions and material. /// Args: /// energy (list(f64)): initial energies in eV. -/// ux (list(f64)): initial ion directions x. ux != 0.0 to avoid gimbal lock +/// ux (list(f64)): initial ion directions x. /// uy (list(f64)): initial ion directions y. /// uz (list(f64)): initial ion directions z. /// Z1 (f64): initial ion atomic number. diff --git a/src/particle.rs b/src/particle.rs index 786a0a0..e9bf65f 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -102,8 +102,6 @@ impl Particle { let dirz = input.uz; let dir_mag = (dirx*dirx + diry*diry + dirz*dirz).sqrt(); - - //assert!((dirx/dir_mag).abs() < 1.0 - f64::EPSILON, "Input error: incident direction cannot round to exactly (1, 0, 0) due to gimbal lock. Use a non-zero y-component."); assert!(input.E > 0., "Input error: incident energy {}; must be greater than zero.", input.E/EV); Particle { @@ -177,7 +175,6 @@ impl Particle { let y = 0.; let z = 0.; - assert!((dirx/dir_mag).abs() < 1.0 - f64::EPSILON, "Input error: incident direction cannot round to exactly (1, 0, 0) due to gimbal lock. Use a non-zero y-component."); assert!(E_eV > 0., "Input error: incident energy {}; must be greater than zero.", E_eV); Particle { diff --git a/src/tests.rs b/src/tests.rs index ea54628..e47a434 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -761,184 +761,186 @@ fn test_surface_refraction() { #[test] fn test_momentum_conservation() { - for energy_eV in vec![1., 10., 100., 1000., 1000.] { - //Aluminum - let m1 = 183.8*AMU; - let Z1 = 74.; - let E1 = energy_eV*EV; - let Ec1 = 1.*EV; - let Es1 = 1.*EV; - let x1 = 0.; - let y1 = 0.; - let z1 = 0.; - - //Aluminum - let m2 = 6.941; - let Z2 = 3.; - let Ec2 = 1.; - let Es2 = 1.; - - //Arbitrary initial angle - let theta = 0.974194583091052_f64; - let cosx = (theta).cos(); - let cosy = (theta).sin()/(2.0_f64).sqrt(); - let cosz = (theta).sin()/(2.0_f64).sqrt(); - - let material_parameters = material::MaterialParameters{ - energy_unit: "EV".to_string(), - mass_unit: "AMU".to_string(), - Eb: vec![0.0], - Es: vec![Es2], - Ec: vec![Ec2], - Ed: vec![0.0], - Z: vec![Z2], - m: vec![m2], - interaction_index: vec![0], - surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET}, - bulk_binding_model: BulkBindingModel::INDIVIDUAL, - }; - - let thickness: f64 = 1000.; - let depth: f64 = 1000.; - let geometry_input = geometry::Mesh2DInput { - length_unit: "ANGSTROM".to_string(), - triangles: vec![(0, 1, 2), (0, 2, 3)], - points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], - densities: vec![vec![0.06306], vec![0.06306]], - boundary: vec![0, 1, 2, 3], - simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], - electronic_stopping_correction_factors: vec![0.0, 0.0], - energy_barrier_thickness: 0., - }; - - let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); - - for high_energy_free_flight_paths in vec![true, false] { - for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { - for scattering_integral in vec![ScatteringIntegral::MENDENHALL_WELLER, ScatteringIntegral::GAUSS_MEHLER{n_points: 10}, ScatteringIntegral::GAUSS_LEGENDRE] { - for root_finder in vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}] { - - println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); - - let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); - - #[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: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: high_energy_free_flight_paths, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![potential]], - scattering_integral: vec![vec![scattering_integral]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![root_finder]], - track_displacements: false, - track_energy_losses: false, - }; - - #[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: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: high_energy_free_flight_paths, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![potential]], - scattering_integral: vec![vec![scattering_integral]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![root_finder]], - 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, - }; - - let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options); - - println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal, - binary_collision_geometries[0].impact_parameter/ANGSTROM, - binary_collision_geometries[0].mfp/ANGSTROM); - - let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1, - &binary_collision_geometries[0], &options); - - let mom1_0 = particle_1.get_momentum(); - let mom2_0 = particle_2.get_momentum(); - - let initial_momentum = mom1_0.add(&mom2_0); - - let binary_collision_result = bca::calculate_binary_collision(&particle_1, - &particle_2, &binary_collision_geometries[0], &options).unwrap(); - - println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV, - binary_collision_result.psi, - binary_collision_result.psi_recoil); - - println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); - - //Energy transfer to recoil - particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); - - //Rotate particle 1, 2 by lab frame scattering angles - particle_1.rotate(binary_collision_result.psi, - binary_collision_geometries[0].phi_azimuthal); - - particle_2.rotate(-binary_collision_result.psi_recoil, - binary_collision_geometries[0].phi_azimuthal); - - //Subtract total energy from all simultaneous collisions and electronic stopping - particle_1.E += -binary_collision_result.recoil_energy; - bca::subtract_electronic_stopping_energy(&mut particle_1, &material_1, 0., - 0., particle_2.Z, species_index, &options); - - let mom1_1 = particle_1.get_momentum(); - let mom2_1 = particle_2.get_momentum(); - - let final_momentum = mom1_1.add(&mom2_1); - - println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); - println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude()); - println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude()); - println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); - println!(); - - assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.)); - assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.)); - assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.)); - - assert!(!particle_1.E.is_nan()); - assert!(!particle_2.E.is_nan()); - assert!(!initial_momentum.x.is_nan()); - assert!(!initial_momentum.x.is_nan()); - assert!(!initial_momentum.x.is_nan()); + for direction in vec![(1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0), (-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, -1.0)] { + for energy_eV in vec![1., 10., 100., 1000., 1000.] { + //Aluminum + let m1 = 183.8*AMU; + let Z1 = 74.; + let E1 = energy_eV*EV; + let Ec1 = 1.*EV; + let Es1 = 1.*EV; + let x1 = 0.; + let y1 = 0.; + let z1 = 0.; + + //Aluminum + let m2 = 6.941; + let Z2 = 3.; + let Ec2 = 1.; + let Es2 = 1.; + + //Arbitrary initial angle + let theta = 0.974194583091052_f64; + let cosx = direction.0; + let cosy = direction.1; + let cosz = direction.2; + + let material_parameters = material::MaterialParameters{ + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: vec![0.0], + Es: vec![Es2], + Ec: vec![Ec2], + Ed: vec![0.0], + Z: vec![Z2], + m: vec![m2], + interaction_index: vec![0], + surface_binding_model: SurfaceBindingModel::PLANAR{calculation: SurfaceBindingCalculation::TARGET}, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let thickness: f64 = 1000.; + let depth: f64 = 1000.; + let geometry_input = geometry::Mesh2DInput { + length_unit: "ANGSTROM".to_string(), + triangles: vec![(0, 1, 2), (0, 2, 3)], + points: vec![(0., -thickness/2.), (depth, -thickness/2.), (depth, thickness/2.), (0., thickness/2.)], + densities: vec![vec![0.06306], vec![0.06306]], + boundary: vec![0, 1, 2, 3], + simulation_boundary_points: vec![(0., 1.1*thickness/2.), (depth, 1.1*thickness/2.), (depth, -1.1*thickness/2.), (0., -1.1*thickness/2.), (0., 1.1*thickness/2.)], + electronic_stopping_correction_factors: vec![0.0, 0.0], + energy_barrier_thickness: 0., + }; + + let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); + + for high_energy_free_flight_paths in vec![true, false] { + for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { + for scattering_integral in vec![ScatteringIntegral::MENDENHALL_WELLER, ScatteringIntegral::GAUSS_MEHLER{n_points: 10}, ScatteringIntegral::GAUSS_LEGENDRE] { + for root_finder in vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}] { + + println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); + + let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); + + #[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: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: high_energy_free_flight_paths, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![potential]], + scattering_integral: vec![vec![scattering_integral]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![root_finder]], + track_displacements: false, + track_energy_losses: false, + }; + + #[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: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: high_energy_free_flight_paths, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![potential]], + scattering_integral: vec![vec![scattering_integral]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![root_finder]], + 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, + }; + + let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options); + + println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal, + binary_collision_geometries[0].impact_parameter/ANGSTROM, + binary_collision_geometries[0].mfp/ANGSTROM); + + let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1, + &binary_collision_geometries[0], &options); + + let mom1_0 = particle_1.get_momentum(); + let mom2_0 = particle_2.get_momentum(); + + let initial_momentum = mom1_0.add(&mom2_0); + + let binary_collision_result = bca::calculate_binary_collision(&particle_1, + &particle_2, &binary_collision_geometries[0], &options).unwrap(); + + println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV, + binary_collision_result.psi, + binary_collision_result.psi_recoil); + + println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); + + //Energy transfer to recoil + particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); + + //Rotate particle 1, 2 by lab frame scattering angles + particle_1.rotate(binary_collision_result.psi, + binary_collision_geometries[0].phi_azimuthal); + + particle_2.rotate(-binary_collision_result.psi_recoil, + binary_collision_geometries[0].phi_azimuthal); + + //Subtract total energy from all simultaneous collisions and electronic stopping + particle_1.E += -binary_collision_result.recoil_energy; + bca::subtract_electronic_stopping_energy(&mut particle_1, &material_1, 0., + 0., particle_2.Z, species_index, &options); + + let mom1_1 = particle_1.get_momentum(); + let mom2_1 = particle_2.get_momentum(); + + let final_momentum = mom1_1.add(&mom2_1); + + println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); + println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude()); + println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude()); + println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); + println!(); + + assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.)); + assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.)); + assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.)); + + assert!(!particle_1.E.is_nan()); + assert!(!particle_2.E.is_nan()); + assert!(!initial_momentum.x.is_nan()); + assert!(!initial_momentum.x.is_nan()); + assert!(!initial_momentum.x.is_nan()); + } } } } From b2d52a05c0fe99d300fb4812b556322fb312f74a Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 08:00:09 -0700 Subject: [PATCH 09/19] update to test_cube.py script. --- examples/test_cube.py | 29 +++++++++++++++-------------- 1 file changed, 15 insertions(+), 14 deletions(-) diff --git a/examples/test_cube.py b/examples/test_cube.py index fba5c82..e0d1a64 100644 --- a/examples/test_cube.py +++ b/examples/test_cube.py @@ -20,7 +20,7 @@ def run(energy, index, num_samples=10000, run_sim=True, a=1000, x0=0, y0=500, z0 track_recoils = true track_recoil_trajectories = {str(track_trajectories).lower()} weak_collision_order = 0 - num_threads = 4 + num_threads = 5 num_chunks = 5 [material_parameters] @@ -112,13 +112,13 @@ def run(energy, index, num_samples=10000, run_sim=True, a=1000, x0=0, y0=500, z0 def main(): do_plots = True - run_sim = False + run_sim = True track_trajectories = False a = 1000 - num_samples = 100000 + num_samples = 10000 num_angles = 10 - energy = 500 - 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) @@ -137,43 +137,44 @@ def main(): R_z_minus = np.zeros(num_angles) 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 == 0: plt.hist(implanted[:, 2], histtype='step', bins=100, density=True, 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 == 0: plt.hist(a - implanted[:, 2], histtype='step', bins=100, density=True, 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 == 0: plt.hist(a - implanted[:, 4], histtype='step', bins=100, density=True, 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 == 0: plt.hist(implanted[:, 4], histtype='step', bins=100, density=True, 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 == 0: plt.hist(implanted[:, 3], histtype='step', bins=100, density=True, 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 == 0: plt.hist(a - implanted[:, 3], histtype='step', bins=100, density=True, 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,9 +207,9 @@ def main(): 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_0', 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]]) + do_trajectory_plot('cube_1', boundary=[[0.0, 0.0], [0.0, a], [a, a], [a, 0.0], [0.0, 0.0]]) if __name__ == '__main__': main() \ No newline at end of file From b44ea4b8a856f6a94d2f7b6e8d18cc91a88aad66 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 11:44:52 -0700 Subject: [PATCH 10/19] Remove now-unnecessary gimbal lock assertions, fudges, etc. --- src/lib.rs | 29 +++++++---------------------- 1 file changed, 7 insertions(+), 22 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index c6e865c..41ea254 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1502,10 +1502,7 @@ pub fn simple_compound_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1 #[cfg(feature = "parry3d")] #[no_mangle] pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut f64, uz: &mut f64) { - //const DELTA: f64 = 1e-9; - //let RUSTBCA_DIRECTION: Vector3:: = Vector3::::new(1.0, 0.0, 0.0); - //let into_surface = Vector3::new(-nx, -ny, -nz); let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global @@ -1514,11 +1511,6 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu //That rotation is then applied to the particle direction, and can be undone later. //Algorithm is from here: //https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/180436#180436 - //let v: Vector3 = into_surface.cross(&RUSTBCA_DIRECTION); - //let c = into_surface.dot(&RUSTBCA_DIRECTION); - //let vx = Matrix3::::new(0.0, -v.z, v.y, v.z, 0.0, -v.x, -v.y, v.x, 0.0); - - //let s = ny*ny + nz*nz; let rotation_matrix = if (1.0 - nx).abs() > 0.0 { Matrix3::::new(1. + (-ny*ny - nz*nz)/(1. - nx), -ny, -nz, ny, -ny*ny/(1. - nx) + 1., -ny*nz/(1. - nx), nz, -ny*nz/(1. - nx), -nz*nz/(1. - nx) + 1.) @@ -1608,9 +1600,7 @@ pub fn rotate_given_surface_normal_vec_py(nx: Vec, ny: Vec, nz: Vec = Vector3::::new(1.0, 0.0, 0.0); - //let into_surface = Vector3::new(-nx, -ny, -nz); let direction = Vector3::new(*ux, *uy, *uz); //Rotation to local RustBCA coordinates from global @@ -1626,6 +1616,7 @@ pub extern "C" fn rotate_back(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut Rotation3::from_axis_angle(&Vector3::y_axis(), PI).into() }; + // Note: transpose of R == R^-1 let u = rotation_matrix.transpose()*direction; *ux = u.x; @@ -1706,8 +1697,6 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle); - const DELTA: f64 = 1e-6; - let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); let m1 = unpack(ion.get_item("m").expect("Cannot get ion mass from dictionary. Ensure ion['m'] exists.")); let Ec1 = unpack(ion.get_item("Ec").expect("Cannot get ion cutoff energy from dictionary. Ensure ion['Ec'] exists.")); @@ -1725,8 +1714,8 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, let y = 0.0; let z = 0.0; - let ux = (angle/180.0*PI).cos() + DELTA; - let uy = (angle/180.0*PI).sin() - DELTA; + let ux = (angle/180.0*PI).cos(); + let uy = (angle/180.0*PI).sin(); let uz = 0.0; let material_parameters = material::MaterialParameters { @@ -1797,8 +1786,6 @@ pub fn sputtering_yield(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, /// R_E (f64): energy reflection coefficient (sum of reflected particle energies / total incident energy) pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: f64, num_samples: usize) -> (f64, f64) { - const DELTA: f64 = 1e-6; - assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle); let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); @@ -1818,8 +1805,8 @@ pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: let y = 0.0; let z = 0.0; - let ux = (angle/180.0*PI).cos() + DELTA; - let uy = (angle/180.0*PI).sin() - DELTA; + let ux = (angle/180.0*PI).cos(); + let uy = (angle/180.0*PI).sin(); let uz = 0.0; let mut direction = Vector::new(ux, uy, uz); @@ -1899,8 +1886,6 @@ pub fn reflection_coefficient(ion: &PyDict, target: &PyDict, energy: f64, angle: /// R_E (f64): energy reflection coefficient (sum of reflected particle energies / total incident energy) pub fn compound_reflection_coefficient(ion: &PyDict, targets: Vec<&PyDict>, target_number_densities: Vec, energy: f64, angle: f64, num_samples: usize) -> (f64, f64) { - const DELTA: f64 = 1e-6; - assert!(angle.abs() <= 90.0, "Incident angle w.r.t. surface normal, {}, cannot exceed 90 degrees.", angle); let Z1 = unpack(ion.get_item("Z").expect("Cannot get ion Z from dictionary. Ensure ion['Z'] exists.")); @@ -1921,8 +1906,8 @@ pub fn compound_reflection_coefficient(ion: &PyDict, targets: Vec<&PyDict>, targ let y = 0.0; let z = 0.0; - let ux = (angle/180.0*PI).cos() + DELTA; - let uy = (angle/180.0*PI).sin() - DELTA; + let ux = (angle/180.0*PI).cos(); + let uy = (angle/180.0*PI).sin(); let uz = 0.0; let mut direction = Vector::new(ux, uy, uz); From 6d2f2a95496056b38ee522155b93818b7548ca92 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 11:58:30 -0700 Subject: [PATCH 11/19] Update examples to take advantage of singularity-free algorithms. --- examples/boron_nitride_0D.toml | 2 +- examples/boron_nitride_sphere.toml | 2 +- examples/boron_nitride_wire.toml | 2 +- examples/boron_nitride_wire_homogeneous.toml | 2 +- examples/layered_geometry.toml | 2 +- examples/layered_geometry_1D.toml | 2 +- examples/lithium_vapor_shield.toml | 2 +- examples/multiple_interaction_potentials.toml | 2 +- examples/tungsten_twist_trimesh.toml | 2 +- 9 files changed, 9 insertions(+), 9 deletions(-) diff --git a/examples/boron_nitride_0D.toml b/examples/boron_nitride_0D.toml index e2da618..0e41b1a 100644 --- 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 index d780ce2..6bf3f81 100644 --- 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 index d172890..6b88ad6 100644 --- 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 index fbe75ad..1682d86 100644 --- 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 index 235cc5a..aaaa434 100644 --- 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 index ae0a3e0..0c6938b 100644 --- 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 index d458937..ea887b9 100644 --- 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 index f2194ab..ba1aea1 100644 --- 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/tungsten_twist_trimesh.toml b/examples/tungsten_twist_trimesh.toml index 0664f3e..e41cf9b 100644 --- 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] From 35a699fc75cfe1d9e6c971401a59f7612e99b982 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 13:26:14 -0700 Subject: [PATCH 12/19] Remove now-unnecessary absolute values on collision angles. --- src/bca.rs | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index d81d235..2b1f067 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -523,8 +523,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)) From be800b3d516b59ab1f09dad028dc6baa2dde06c2 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 13:28:32 -0700 Subject: [PATCH 13/19] Remove unnecessary theta definition. --- src/tests.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tests.rs b/src/tests.rs index e47a434..0594493 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -780,7 +780,6 @@ fn test_momentum_conservation() { let Es2 = 1.; //Arbitrary initial angle - let theta = 0.974194583091052_f64; let cosx = direction.0; let cosy = direction.1; let cosz = direction.2; From 2936a3a4a71644275297c8d66f8f3aa8e1a7d79f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 13:32:03 -0700 Subject: [PATCH 14/19] Fix comments. --- src/tests.rs | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tests.rs b/src/tests.rs index 0594493..2d84d12 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -779,7 +779,6 @@ fn test_momentum_conservation() { let Ec2 = 1.; let Es2 = 1.; - //Arbitrary initial angle let cosx = direction.0; let cosy = direction.1; let cosz = direction.2; From 47118538d5d628f7d6ebf4bd0c025f19cc66e431 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 13:51:23 -0700 Subject: [PATCH 15/19] Update tests to include morse potential when CPR is activated. --- src/tests.rs | 287 +++++++++++++++++++++++++++++---------------------- 1 file changed, 163 insertions(+), 124 deletions(-) diff --git a/src/tests.rs b/src/tests.rs index 2d84d12..6f35021 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -761,6 +761,41 @@ fn test_surface_refraction() { #[test] fn test_momentum_conservation() { + + #[cfg(not(feature = "cpr_rootfinder"))] + let potentials = vec![ + InteractionPotential::KR_C, + InteractionPotential::MOLIERE, + InteractionPotential::ZBL, + InteractionPotential::LENZ_JENSEN + ]; + + #[cfg(feature = "cpr_rootfinder")] + let potentials = vec![ + InteractionPotential::KR_C, + InteractionPotential::MOLIERE, + InteractionPotential::ZBL, + InteractionPotential::LENZ_JENSEN, + InteractionPotential::MORSE{D: 5.4971E-20, r0: 2.782E-10, alpha: 1.4198E10} + ]; + + let mut rootfinders = vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}; 4]; + + //[[{"CPR"={n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}]] + #[cfg(feature = "cpr_rootfinder")] + rootfinders.push( + Rootfinder::CPR{ + n0: 2, + nmax: 100, + epsilon: 1e-9, + complex_threshold: 1e-3, + truncation_threshold: 1e-9, + far_from_zero: 1e9, + interval_limit:1e-12, + derivative_free: true + } + ); + for direction in vec![(1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0), (-1.0, 0.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, -1.0)] { for energy_eV in vec![1., 10., 100., 1000., 1000.] { //Aluminum @@ -813,132 +848,136 @@ fn test_momentum_conservation() { let material_1: material::Material = material::Material::::new(&material_parameters, &geometry_input); for high_energy_free_flight_paths in vec![true, false] { - for potential in vec![InteractionPotential::KR_C, InteractionPotential::MOLIERE, InteractionPotential::ZBL, InteractionPotential::LENZ_JENSEN] { + for (potential, root_finder) in potentials.iter().zip(rootfinders.clone()) { for scattering_integral in vec![ScatteringIntegral::MENDENHALL_WELLER, ScatteringIntegral::GAUSS_MEHLER{n_points: 10}, ScatteringIntegral::GAUSS_LEGENDRE] { - for root_finder in vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}] { - - println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); - - let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); - - #[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: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: high_energy_free_flight_paths, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![potential]], - scattering_integral: vec![vec![scattering_integral]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![root_finder]], - track_displacements: false, - track_energy_losses: false, - }; - - #[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: 0, - suppress_deep_recoils: false, - high_energy_free_flight_paths: high_energy_free_flight_paths, - electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, - mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![potential]], - scattering_integral: vec![vec![scattering_integral]], - num_threads: 1, - num_chunks: 1, - use_hdf5: false, - root_finder: vec![vec![root_finder]], - 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, - }; - - let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options); - - println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal, - binary_collision_geometries[0].impact_parameter/ANGSTROM, - binary_collision_geometries[0].mfp/ANGSTROM); - - let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1, - &binary_collision_geometries[0], &options); - - let mom1_0 = particle_1.get_momentum(); - let mom2_0 = particle_2.get_momentum(); - - let initial_momentum = mom1_0.add(&mom2_0); - - let binary_collision_result = bca::calculate_binary_collision(&particle_1, - &particle_2, &binary_collision_geometries[0], &options).unwrap(); - - println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV, - binary_collision_result.psi, - binary_collision_result.psi_recoil); - - println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); - - //Energy transfer to recoil - particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); - - //Rotate particle 1, 2 by lab frame scattering angles - particle_1.rotate(binary_collision_result.psi, - binary_collision_geometries[0].phi_azimuthal); - - particle_2.rotate(-binary_collision_result.psi_recoil, - binary_collision_geometries[0].phi_azimuthal); - - //Subtract total energy from all simultaneous collisions and electronic stopping - particle_1.E += -binary_collision_result.recoil_energy; - bca::subtract_electronic_stopping_energy(&mut particle_1, &material_1, 0., - 0., particle_2.Z, species_index, &options); - - let mom1_1 = particle_1.get_momentum(); - let mom2_1 = particle_2.get_momentum(); - - let final_momentum = mom1_1.add(&mom2_1); - - println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); - println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude()); - println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude()); - println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); - println!(); - - assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.)); - assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.)); - assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.)); - - assert!(!particle_1.E.is_nan()); - assert!(!particle_2.E.is_nan()); - assert!(!initial_momentum.x.is_nan()); - assert!(!initial_momentum.x.is_nan()); - assert!(!initial_momentum.x.is_nan()); + + //Skip incompatible combination + match (root_finder, scattering_integral) { + (Rootfinder::CPR{..}, ScatteringIntegral::MENDENHALL_WELLER) => continue, + _ => {} } + + println!("Case: {} {} {} {}", energy_eV, high_energy_free_flight_paths, potential, scattering_integral); + + let mut particle_1 = particle::Particle::new(m1, Z1, E1, Ec1, Es1, 0.0, x1, y1, z1, cosx, cosy, cosz, false, false, 0); + + #[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: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: high_energy_free_flight_paths, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![*potential]], + scattering_integral: vec![vec![scattering_integral]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![root_finder]], + track_displacements: false, + track_energy_losses: false, + }; + + #[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: 0, + suppress_deep_recoils: false, + high_energy_free_flight_paths: high_energy_free_flight_paths, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![potential]], + scattering_integral: vec![vec![scattering_integral]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![root_finder]], + 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, + }; + + let binary_collision_geometries = bca::determine_mfp_phi_impact_parameter(&mut particle_1, &material_1, &options); + + println!("Phi: {} rad p: {} Angstrom mfp: {} Angstrom", binary_collision_geometries[0].phi_azimuthal, + binary_collision_geometries[0].impact_parameter/ANGSTROM, + binary_collision_geometries[0].mfp/ANGSTROM); + + let (species_index, mut particle_2) = bca::choose_collision_partner(&mut particle_1, &material_1, + &binary_collision_geometries[0], &options); + + let mom1_0 = particle_1.get_momentum(); + let mom2_0 = particle_2.get_momentum(); + + let initial_momentum = mom1_0.add(&mom2_0); + + let binary_collision_result = bca::calculate_binary_collision(&particle_1, + &particle_2, &binary_collision_geometries[0], &options).unwrap(); + + println!("E_recoil: {} eV Psi: {} rad Psi_recoil: {} rad", binary_collision_result.recoil_energy/EV, + binary_collision_result.psi, + binary_collision_result.psi_recoil); + + println!("Initial Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); + + //Energy transfer to recoil + particle_2.E = binary_collision_result.recoil_energy - material_1.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); + + //Rotate particle 1, 2 by lab frame scattering angles + particle_1.rotate(binary_collision_result.psi, + binary_collision_geometries[0].phi_azimuthal); + + particle_2.rotate(-binary_collision_result.psi_recoil, + binary_collision_geometries[0].phi_azimuthal); + + //Subtract total energy from all simultaneous collisions and electronic stopping + particle_1.E += -binary_collision_result.recoil_energy; + bca::subtract_electronic_stopping_energy(&mut particle_1, &material_1, 0., + 0., particle_2.Z, species_index, &options); + + let mom1_1 = particle_1.get_momentum(); + let mom2_1 = particle_2.get_momentum(); + + let final_momentum = mom1_1.add(&mom2_1); + + println!("Final Energies: {} eV {} eV", particle_1.E/EV, particle_2.E/EV); + println!("X: {} {} {}% Error", initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, 100.*(final_momentum.x - initial_momentum.x)/initial_momentum.magnitude()); + println!("Y: {} {} {}% Error", initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, 100.*(final_momentum.y - initial_momentum.y)/initial_momentum.magnitude()); + println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); + println!(); + + assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.)); + assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.)); + assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.)); + + assert!(!particle_1.E.is_nan()); + assert!(!particle_2.E.is_nan()); + assert!(!initial_momentum.x.is_nan()); + assert!(!initial_momentum.x.is_nan()); + assert!(!initial_momentum.x.is_nan()); } } } From 165d8688edf78c3126431c17339582836531475e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 14:00:13 -0700 Subject: [PATCH 16/19] Fix to mismatched types. --- src/tests.rs | 1 + 1 file changed, 1 insertion(+) diff --git a/src/tests.rs b/src/tests.rs index 6f35021..7ab7689 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -969,6 +969,7 @@ fn test_momentum_conservation() { println!("Z: {} {} {}% Error", initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, 100.*(final_momentum.z - initial_momentum.z)/initial_momentum.magnitude()); println!(); + //These values are in [angstrom amu / second], so very large. assert!(approx_eq!(f64, initial_momentum.x/ANGSTROM/AMU, final_momentum.x/ANGSTROM/AMU, epsilon = 1000.)); assert!(approx_eq!(f64, initial_momentum.y/ANGSTROM/AMU, final_momentum.y/ANGSTROM/AMU, epsilon = 1000.)); assert!(approx_eq!(f64, initial_momentum.z/ANGSTROM/AMU, final_momentum.z/ANGSTROM/AMU, epsilon = 1000.)); From 38489140852b3f45040e82812140d1eb16ea3dfe Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 14:11:30 -0700 Subject: [PATCH 17/19] Add dereferncing fix to distributions test of momentum conservation. --- src/tests.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tests.rs b/src/tests.rs index 7ab7689..b926984 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -895,7 +895,7 @@ fn test_momentum_conservation() { high_energy_free_flight_paths: high_energy_free_flight_paths, electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, mean_free_path_model: MeanFreePathModel::LIQUID, - interaction_potential: vec![vec![potential]], + interaction_potential: vec![vec![*potential]], scattering_integral: vec![vec![scattering_integral]], num_threads: 1, num_chunks: 1, From 4816723c4c4b2fbc5020fc596151d39c22c56ed6 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 17:21:06 -0700 Subject: [PATCH 18/19] Add comment elaborating on the consequences of the singularity-free algorithm --- src/bca.rs | 7 +++++-- src/particle.rs | 6 ++++-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 2b1f067..14a068b 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -364,8 +364,11 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma let sinx: f64 = (1. - cosx*cosx).sqrt(); let cosphi: f64 = phi_azimuthal.cos(); - //Find recoil location - + // 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 + // 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 { diff --git a/src/particle.rs b/src/particle.rs index e9bf65f..5d6bb69 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -238,12 +238,14 @@ impl Particle { let cosz: f64 = self.dir.z; let cosphi: f64 = (phi + PI).cos(); let sinphi: f64 = (phi + PI).sin(); - let sa = (1. - cosx*cosx).sqrt(); - //Particle direction update formulas from original TRIDYN paper, see Moeller and Eckstein 1988 let cpsi: f64 = psi.cos(); let spsi: f64 = psi.sin(); + // 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 + // 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 cosx_new = if cosx > -1. { cpsi*cosx - spsi*(cosz*sinphi + cosy*cosphi) } else { From 1b93c017647f808c07be2682323d153585bc5676 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Mon, 21 Jul 2025 17:21:40 -0700 Subject: [PATCH 19/19] Amend comment elaborating on the consequences of the singularity-free algorithm --- src/bca.rs | 2 +- src/particle.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 14a068b..353c24d 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -366,7 +366,7 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma // 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 + // 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. { diff --git a/src/particle.rs b/src/particle.rs index 5d6bb69..79cc09d 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -243,7 +243,7 @@ impl Particle { let spsi: f64 = psi.sin(); // 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 + // 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 cosx_new = if cosx > -1. {