From db015aa406b3c86ca30b3e768e63842311b825f0 Mon Sep 17 00:00:00 2001 From: Jon Drobny <37962344+drobnyjt@users.noreply.github.com> Date: Mon, 9 Jun 2025 17:33:02 -0700 Subject: [PATCH 01/31] Update lib.rs --- src/lib.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lib.rs b/src/lib.rs index 3a9f53e..4e52076 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1518,7 +1518,7 @@ 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 rotation_matrix = if (c + 1.0).abs() > 1e-6 { + let rotation_matrix = if (c + 1.0).abs() > DELTA { Matrix3::identity() + vx + vx*vx/(1. + c) } else { //If c == -1.0, the correct rotation should simply be a 180 degree rotation From 2644514115434e6a9815e9884cead919ffc94599 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 11 Jun 2025 10:32:35 -0700 Subject: [PATCH 02/31] 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 03/31] 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 04/31] 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 05/31] 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 06/31] 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 07/31] 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 8e26b98b73ef5cae56bac7f75572b0402b314a60 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sat, 19 Jul 2025 14:07:27 -0700 Subject: [PATCH 08/31] Updates to test_cube.py --- examples/test_cube.py | 215 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 215 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..13e6f7b --- /dev/null +++ b/examples/test_cube.py @@ -0,0 +1,215 @@ +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(): + do_plots = True + run_sim = True + track_trajectories = False + a = 1000 + num_samples = 100000 + num_angles = 10 + energy = 100 + 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) + + R_x_minus = np.zeros(num_angles) + R_y_minus = np.zeros(num_angles) + R_z_minus = np.zeros(num_angles) + + sim_index = 0 + impl_index = num_angles*3//4 + 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') + + 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') + 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') + 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') + 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') + 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') + sputtering_yields_y_minus[angle_index] = Y + R_y_minus[angle_index] = R + + plt.figure(1) + plt.title('Implantation') + plt.legend() + + 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]]) + + 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 71a16379d6aa642ed7c212f6c45fcf68460897ca Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Sun, 20 Jul 2025 11:38:17 -0700 Subject: [PATCH 09/31] 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 10/31] 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 11/31] 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 12/31] 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 13/31] 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 14/31] 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 15/31] 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 16/31] 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 17/31] 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 18/31] 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 19/31] 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 20/31] 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 21/31] 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. { From 1c3f4ac0fc767d0fad16d88e9c8cbdd3cb8d0090 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Fri, 22 Aug 2025 14:20:03 -0700 Subject: [PATCH 22/31] Replace sigmoid with smootherstep in Kr-C-Morse. --- src/interactions.rs | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/interactions.rs b/src/interactions.rs index b050294..921aa7e 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -12,8 +12,18 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { } -fn sigmoid(x: f64, k: f64, x0: f64) -> f64 { - 1./(1. + (-k*(x - x0)).exp()) +fn smootherstep(x: f64, k: f64, x0: f64) -> f64 { + let x_transformed = x*k/8.0; //Constant to match sigmoid, determined by eye + let x0_transformed = x0*k/8.0; + + if x_transformed <= x0_transformed - 0.5 { + 0. + } else if x_transformed >= 0.5 + x0_transformed { + 1. + } else { + let x1 = x_transformed - x0_transformed + 0.5; + return x1 * x1 * x1 * (x1 * (6. * x1 - 15.) + 10.); + } } /// Interaction potential between two particles a and b at a distance `r`. @@ -346,9 +356,9 @@ pub fn morse(r: f64, D: f64, alpha: f64, r0: f64) -> f64 { D*((-2.*alpha*(r - r0)).exp() - 2.*(-alpha*(r - r0)).exp()) } -/// Kr-C + Morse potential with sigmoid interpolation at crossing point of two potentials +/// Kr-C + Morse potential with smootherstep interpolation at crossing point of two potentials pub fn krc_morse(r: f64, a: f64, Za: f64, Zb: f64, D: f64, alpha: f64, r0: f64, k: f64, x0: f64) -> f64 { - sigmoid(r, k, x0)*morse(r, D, alpha, r0) + sigmoid(r, -k, x0)*screened_coulomb(r, a, Za, Zb, InteractionPotential::KR_C) + smootherstep(r, k, x0)*morse(r, D, alpha, r0) + smootherstep(r, -k, x0)*screened_coulomb(r, a, Za, Zb, InteractionPotential::KR_C) } /// Distance of closest approach function for four-eight potential (Born repulsion) From 3fc46939039852a68846782d526581eabfc36322 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 26 Aug 2025 16:51:45 -0700 Subject: [PATCH 23/31] Add comparison to 2024 V Shulga paper to examples. --- examples/test_sputtered_atom_spectra.py | 482 ++++++++++++++++++++++++ 1 file changed, 482 insertions(+) create mode 100644 examples/test_sputtered_atom_spectra.py diff --git a/examples/test_sputtered_atom_spectra.py b/examples/test_sputtered_atom_spectra.py new file mode 100644 index 0000000..9dccda5 --- /dev/null +++ b/examples/test_sputtered_atom_spectra.py @@ -0,0 +1,482 @@ +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 formulas import * + +data_exp_0_10_deg = np.array([ + [1.795219872324359, 0.9533726147006264], + [2.431545232813113, 1.0000000000000002], + [3.552942960337357, 0.968668361722603], + [5.135565058063388, 0.839388729523552], + [6.518049391281579, 0.727362910870967], + [8.73326162382844, 0.5549318240872879], + [10.964781961431852, 0.4101127070551304], + [14.532921950830353, 0.27113134189526106], + [18.246369720400814, 0.1880154542858142], + [21.466671275519875, 0.12832017135233023], + [25.808615404180767, 0.09186152281210946], + [31.70855642178006, 0.05882815844474357], + [36.90275785618217, 0.04144875792022527], + [40.68288702863619, 0.03063198320387291], + [46.33213925076186, 0.021928747399631854], + [52.19718220435652, 0.014966326223296179], + [57.54399373371572, 0.009284145445194745], +]) +data_oksana_0_10_deg = np.array([ + [1.0108946133097567, 0.7045734392511628], + [1.3254182346670706, 0.8130893054759956], + [1.795219872324359, 0.9383183950023508], + [2.4580359778697947, 1.0000000000000002], + [3.591650900001854, 1.0000000000000002], + [5.135565058063388, 0.908919342461109], + [6.660846290809163, 0.7629366521078121], + [8.73326162382844, 0.6203356325845223], + [10.964781961431852, 0.48086996251685976], + [14.532921950830353, 0.35537857716782895], + [18.246369720400814, 0.2584891963499745], + [21.466671275519875, 0.1850465983402524], + [25.808615404180767, 0.13675546811805692], + [31.70855642178006, 0.08898335080774952], + [36.90275785618217, 0.0668166646296556], + [40.68288702863619, 0.05432791543794923], + [45.83280852498209, 0.04015010043002134], + [51.63464273833496, 0.02451325153735838], + [58.17091329374364, 0.018406743419819866], + [69.1830970918937, 0.011418358823154931], + [82.27996866848511, 0.006137880562106176], + [97.85617482684427, 0.002998876115879274], + [115.1268098518791, 0.001374829517774794], +]) + +data_oksana = np.array([ + [1.1708281368909705, 0.8931414506078489], + [1.5308460109164432, 0.9721427940738623], + [1.970246424053571, 0.9721427940738623], + [2.7008897040734974, 0.9721427940738623], + [3.8211297149093113, 0.8682610252970971], + [6.037010313703754, 0.6363279236793021], + [8.954777399948457, 0.4407291533212], + [13.282739984823985, 0.25047997204261246], + [19.09070284923147, 0.13078695156316839], + [27.438234558085025, 0.06638743505522715], + [42.671409285185625, 0.02844380745224324], + [61.3297554384229, 0.011196434778829558], + [90.97123897530997, 0.003720078867528625], + [130.7489940395067, 0.0011039382600225845], + [190.90702849231468, 0.00022057601289773862], +]) +data_trim_sp = np.array([ + [1.1708281368909705, 0.8682610252970971], + [1.5308460109164432, 0.9721427940738623], + [1.970246424053571, 0.9721427940738623], + [2.7008897040734974, 0.9721427940738623], + [3.8211297149093113, 0.8205601642645971], + [6.037010313703754, 0.6013690933573494], + [8.954777399948457, 0.4165162041238635], + [13.282739984823985, 0.2650408913489223], + [18.791982647375733, 0.16395472095739783], + [27.438234558085025, 0.09058477599821053], + [42.671409285185625, 0.039923351629481865], + [61.3297554384229, 0.01809958893761981], + [90.97123897530997, 0.006926132531994934], + [130.7489940395067, 0.0023671900616977384], + [187.91982647375713, 0.0005929337004735993], + [253.57684207388002, 0.00010000000000000021], +]) + +oksana_500 = np.array([ + [3.2154340836012807, 0.9298674652605262], + [12.861736334405123, 0.257351270001691], + [24.115755627009662, 0.1], + [40.192926045016065, 0.032008340465997674], + [62.700964630225144, 0.010496696290308798], + [86.81672025723475, 0.0036132227567962547], + [114.14790996784569, 0.001000000000000002], + [138.2636655948553, 0.00024517358879792915], + [167.20257234726688, 0.000058670670659931154], + [191.31832797427649, 0.000010245338593872244], +]) +oksana_1000 = np.array([ + [1.6077170418006403, 0.9298674652605262], + [12.861736334405123, 0.25118864315095796], + [30.54662379421228, 0.06622967617148331], + [40.192926045016065, 0.039810717055349734], + [62.700964630225144, 0.016636142493842227], + [86.81672025723475, 0.008040131611167865], + [114.14790996784569, 0.003792690190732254], + [144.69453376205792, 0.0015469407652462018], + [184.88745980707398, 0.0006309573444801943], + [234.7266881028939, 0.0001924024182760241], + [294.21221864951764, 0.000037018690558462135], + [337.6205787781351, 0.00001049669629030881], +]) +oksana_2000 = np.array([ + [1.6077170418006403, 0.9298674652605262], + [12.861736334405123, 0.25118864315095796], + [36.977491961414785, 0.06309573444801936], + [54.662379421221885, 0.031241857136026677], + [88.42443729903539, 0.011849068510006685], + [135.04823151125402, 0.004493984590721679], + [186.49517684887456, 0.0018779508018514951], + [245.98070739549843, 0.0008439481965654006], + [303.8585209003215, 0.00039810717055349773], + [366.5594855305466, 0.0001924024182760241], + [435.6913183279742, 0.00009076005216818159], + [499.99999999999994, 0.000035266992141746665], + [559.4855305466237, 0.000017044291274532015], +]) +oksana_5000 = np.array([ + [12.861736334405123, 0.25118864315095796], + [30.54662379421228, 0.11018063301098226], + [56.270096463022526, 0.0379269019073225], + [94.85530546623801, 0.016237767391887224], + [146.30225080385856, 0.006785454573393582], + [226.6881028938907, 0.002701338121133005], + [327.97427652733126, 0.0011565326417902512], + [406.7524115755627, 0.0006622967617148338], + [464.63022508038586, 0.00044939845907216705], + [553.0546623794212, 0.00026366508987303664], + [638.2636655948554, 0.00016636142493842227], + [729.903536977492, 0.00010754220761125644], + [834.4051446945336, 0.0000615848211066028], + [1000, 0.000023930257311805045], +]) + +srim_500 = np.array([ + [3.2154340836013944, 1.024599382282019], + [12.861736334405123, 0.20605547994561435], + [35.369774919613974, 0.040444564028863234], + [77.17041800643074, 0.008964098010170358], + [144.69453376205774, 0.00253334263054688], + [244.3729903536978, 0.0009353545286845042], + [344.0514469453374, 0.0006039522986305644], + [424.4372990353695, 0.0005753001127001087], + [450.16077170417975, 0.00022298978070390968], + [456.59163987138254, 0.000036921753789087356], + [456.59163987138254, 0.000010183933793683], +]) +srim_1000 = np.array([ + [3.2154340836013944, 1.024599382282019], + [12.861736334405123, 0.20605547994561435], + [38.58520900321537, 0.040444564028863234], + [80.38585209003213, 0.008964098010170358], + [247.58842443729918, 0.000958363672304831], + [340.83601286173644, 0.000522009131604022], + [446.9453376205788, 0.00034534929613596964], + [578.7781350482314, 0.00022298978070390968], + [697.7491961414789, 0.00016658499644765562], + [803.8585209003218, 0.00013715232249995635], + [897.1061093247586, 0.0001625854937337872], + [926.0450160771702, 0.00006150610223796307], + [938.9067524115758, 0.000019156753004056052], + [938.9067524115758, 0.00001043445227420857], +]) +srim_2000 = np.array([ + [3.2154340836013944, 1.0000000000000004], + [16.077170418006062, 0.20605547994561435], + [38.58520900321537, 0.04143947532063884], + [77.17041800643074, 0.008964098010170358], + [144.69453376205774, 0.002595661294367036], + [250.80385852090012, 0.0009353545286845042], + [347.2668810289392, 0.0005753001127001087], + [446.9453376205788, 0.00034534929613596964], + [578.7781350482314, 0.00022298978070390968], + [700.9646302250803, 0.00015487225246188873], + [797.427652733119, 0.00012444768075857784], + [900.32154340836, 0.00010245993822820222], + [922.8295819935688, 0.00009296890822370853], + [1083.601286173633, 0.00007842582625510765], + [1286.1736334405145, 0.000058588182592362524], + [1466.2379421221872, 0.00005188473200240593], + [1665.5948553054664, 0.00005580869238300307], + [1803.8585209003218, 0.000057181551741592105], + [1871.3826366559483, 0.00002758247577257015], + [1897.1061093247586, 0.000010183933793683], +]) +srim_5000 = np.array([ + [3.2154340836013944, 1.0000000000000004], + [12.861736334405123, 0.20110833903362435], + [41.80064308681631, 0.042458860815617526], + [77.17041800643074, 0.008748880943306109], + [144.69453376205774, 0.0032302402811878513], + [244.3729903536978, 0.001192661363416099], + [347.2668810289392, 0.0006987586173067845], + [450.16077170417975, 0.00041946071507487087], + [578.7781350482314, 0.00027750561041465377], + [691.3183279742766, 0.00019273492632738962], + [807.0739549839222, 0.00014752493615384886], + [922.8295819935688, 0.00012444768075857784], + [1086.8167202572345, 0.00009525588593744557], + [1292.6045016077164, 0.00007291163443479598], + [1466.2379421221872, 0.000058588182592362524], + [1665.5948553054664, 0.0000482366688774338], + [1803.8585209003218, 0.00004376842812354139], + [1999.9999999999995, 0.000034325751394760175], +]) + +yield_trim_sp = np.array([ + [99.44868379777449, 0.14099733492521718], + [148.8912402825362, 0.3306501492695028], + [199.58139231677055, 0.5250849312986958], + [297.15921039359404, 0.8449466108283331], + [496.91088635950564, 1.289654737002661], + [696.2049308523546, 1.6251889984207488], + [1002.7680317396586, 1.9297867743942156], + [2012.433270130691, 2.496994027573269], + [2996.3368557559875, 2.738985101913617], + [5010.487142071969, 2.8876542595284915], + [10055.437254790288, 3.2096393886804577], +]) + +yield_oksana_zbl = np.array([ + [98.90040709109736, 0.21806902115934437], + [137.80205758914963, 0.3798629815834076], + [198.48106776430075, 0.5875029383534014], + [396.1312785291842, 1.0647795078569717], + [700.0644998661448, 1.433456456230312], + [1008.3271024267672, 1.657723885589316], + [2012.433270130691, 2.034533519429473], + [4016.442339972745, 2.3528407832211435], + [7020.023006215825, 2.496994027573269], + [10111.181838502403, 2.5638656594088123], + +]) +yield_oksana_krc = np.array([ + [100, 0.2329637662380699], + [140.1065965676314, 0.41120677198225625], + [199.58139231677055, 0.6660846290809157], + [400.53553891270326, 1.1835067295194683], + [696.2049308523546, 1.6359623861500767], + [1008.3271024267672, 1.9044538138367328], + [2023.5896477251576, 2.3066633648505985], + [4016.442339972745, 2.632528170635725], + [6981.3204819826, 2.7938172299514883], + [10111.181838502403, 2.906796537132128], + +]) +yield_srim = np.array([ + [98.90040709109736, 0.4947714040850636], + [150.54664046568803, 0.7306368599685215], + [199.58139231677055, 0.9268328852273553], + [398.3273215910063, 1.6251889984207488], + [703.9454652710157, 2.4158518681039256], + [997.239609109939, 2.964988105897036], + [2001.3383994534818, 4.125674667864929], + [4016.442339972745, 5.199066764846401], + [6981.3204819826, 5.894470021435907], + [10055.437254790288, 6.132836773591143], +]) + +yield_exp = np.array([ + [98.90040709109736, 0.24560878824173277], + [98.90040709109736, 0.28030620499491193], + [149.71665244489452, 0.4600878303875116], + [149.71665244489452, 0.5013528314599843], + [198.48106776430075, 0.5875029383534014], + [198.48106776430075, 0.6660846290809157], + [250.35681838914445, 0.7162972054537278], + [248.97656068598883, 0.867575002428819], + [298.80657947958076, 0.9516542611656196], + [350.7664175200366, 1.0932952121194448], + [496.91088635950564, 1.3068096673810672], + [599.6684007122172, 1.350701993655997], + [599.6684007122172, 1.5013107289081726], + [647.9248097041271, 1.5932926487627252], + [803.8260284178258, 1.8547811900706128], + [1002.7680317396586, 2.1308404505536047], + [3029.6506797955367, 2.597970087225646], + [9835.515312273716, 3.1466462350737068], +]) + +#This function simply contains an entire input file as a multi-line f-string to modify some inputs. +def run_rustbca(ion, target, energy, index=0, num_samples=10000, run_sim=True, estop='LOW_ENERGY_NONLOCAL', weak_collisions=0, Eb=3.0, correction=1.0, interaction='KR_C'): + + input_file = f''' + [options] + name = "{ion["symbol"]}_{target["symbol"]}_{index}" + weak_collision_order = {weak_collisions} + electronic_stopping_mode = "{estop}" + mean_free_path_model = "LIQUID" + scattering_integral = [["MENDENHALL_WELLER"]] + interaction_potential = [["{interaction}"]] + num_threads = 4 + num_chunks = {int(min(num_samples//10, 100))} + + [particle_parameters] + length_unit = "ANGSTROM" + energy_unit = "EV" + mass_unit = "AMU" + N = [ {num_samples} ] + m = [ {ion["m"]}] + Z = [ {ion["Z"]} ] + E = [ {energy} ] + Ec = [ {ion["Ec"]} ] + Es = [ {ion["Es"]} ] + interaction_index = [ 0 ] + pos = [ [ -4.4, 0.0, 0.0,] ] + dir = [ [ 1.0, 0.0, 0.0] ] + + [geometry_input] + length_unit = "ANGSTROM" + electronic_stopping_correction_factor = {correction} + densities = [ {target["n"]/10**30} ] + + [material_parameters] + energy_unit = "EV" + mass_unit = "AMU" + Eb = [ {Eb} ] + Es = [ {target["Es"]} ] + Ec = [ {target["Ec"]} ] + Z = [ {target["Z"]} ] + m = [ {target["m"]} ] + ''' + + with open(f'{ion["symbol"]}_{target["symbol"]}_{index}.toml', 'w') as file: + file.write(input_file) + + if run_sim: os.system(f'cargo run --release --features cpr_rootfinder 0D {ion["symbol"]}_{target["symbol"]}_{index}.toml') + + sputtered_list = np.atleast_2d(np.genfromtxt(f'{ion["symbol"]}_{target["symbol"]}_{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_reflected = 0. + + + return sputtered_list, num_sputtered/num_samples + +# Figure 4(a), V Shulga 2024 + +run_sim = True +sim_index = 0 +num_samples = 100000 +energy = 1000.0 +nickel["Es"] = 4.46 # From Fig. 4(a) in V Shulga 2024 +sputtered, Y = run_rustbca(argon, nickel, energy, sim_index, num_samples, run_sim) + +hist, bin_edges = np.histogram(sputtered[:, 2], bins=1000, range=[0.0, 1000.0]) +bin_centers = (bin_edges[1:] + bin_edges[:-1])/2. + +plt.figure() +plt.plot(bin_centers, hist/np.max(hist), label='RustBCA', linewidth=1) +plt.plot(data_oksana[:, 0], data_oksana[:, 1], label='TRIM.SP', linestyle='', marker='*', markersize=7) +plt.plot(data_trim_sp[:, 0], data_trim_sp[:, 1], label='OKSANA', linestyle='--', linewidth=3) +plt.xlim([1.0, 1000.0]) +plt.ylim([1e-4, 2.0]) +plt.legend() +plt.xscale('log') +plt.yscale('log') +plt.title('Comparison to Fig. 4(a) in V Shulga 2024; Ar on Ni 1 keV, Es=4.46 eV, sputtered energies') +plt.xlabel('E [eV]') +plt.ylabel('N/N_max') + +argon['Es'] = 5.0 +vanadium = { + 'Z': 23.0, + 'm': 50.9415, + 'n': 0.07223*10**30, + 'Es': 5.33, #From Fig. 5 in V Shulga 2024 + 'Ec': 5.33, + 'Eb': 0.0, + 'symbol': 'V', + 'name': 'vanadium' +} + +# Figure 5, V Shulga 2024 +run_sim = True +num_samples = 1000000 +energy = 1000.0 +sim_index = 0 +sputtered, Y = run_rustbca(argon, vanadium, energy, sim_index, num_samples, run_sim) +cutoff = np.cos(10.0*np.pi/180.0) +sputtered = sputtered[sputtered[:, 6] < -cutoff] + +plt.figure() +plt.plot(bin_centers, hist/np.max(hist), label='RustBCA', linewidth=1) +plt.plot(data_exp_0_10_deg[:, 0], data_exp_0_10_deg[:, 1], label='Exp. (Dembowski86)', marker='^', linestyle='', markersize=7) +plt.plot(data_oksana_0_10_deg[:, 0], data_oksana_0_10_deg[:, 1], label='OKSANA', linestyle='--', linewidth=3) +plt.xlim([1.0, 1000.0]) +plt.ylim([1e-3, 2.0]) +plt.legend() +plt.xscale('log') +plt.yscale('log') +plt.title('Comparison to Fig. 5 in V Shulga 2024; Ar on V 1 keV, Es=5.33 eV, sputtered energies theta=0-10') +plt.xlabel('E [eV]') +plt.ylabel('N/N_max') + +# Figure 7, V Shulga 2024 +run_sim = True +energies = [500, 1000, 2000, 5000] +oksana_results = [oksana_500, oksana_1000, oksana_2000, oksana_5000] +srim_results = [srim_500, srim_1000, srim_2000, srim_5000] +nickel['Es'] = 4.46 +num_samples = 100000 + +plt.figure() +index = 1 +for energy, oksana, srim in zip(energies, oksana_results, srim_results): + sputtered, Y = run_rustbca(argon, nickel, energy, index, num_samples, run_sim) + hist, bin_edges = np.histogram(sputtered[:, 2], bins=100, range=[0.0, 1000.0]) + bin_centers = (bin_edges[1:] + bin_edges[:-1])/2. + color = plt.semilogy(bin_centers, hist/np.max(hist), label=f'RustBCA {energy/1000.0} keV')[0].get_color() + plt.semilogy(oksana[:, 0], oksana[:, 1], color=color, linestyle='--', linewidth=3, label=f'Oksana {energy/1000.0} keV') + plt.semilogy(srim[:, 0], srim[:, 1], color=color, linestyle=':', linewidth=3, label=f'SRIM {energy/1000.0} keV') + index += 1 + + + +plt.title('Comparison to Fig. 7 in V Shulga 2024; 0.5-5 keV Ar on Si, Es=4.46 eV, sputtered energies') +plt.xlabel('E [eV]') +plt.ylabel('N/N_max') +plt.legend() +plt.xlim([0.0, 1000.0]) +plt.ylim([1e-5, 2.0]) + + +# Figure 8 A V Shulga 2024: Sputtering yields +index = 5 +run_sim = True +num_energies = 5 +energies = np.logspace(2, 4, num_energies) +sputtering_yield_0_eV = np.zeros(num_energies) +sputtering_yield_2_eV = np.zeros(num_energies) +sputterinng_yield_zbl = np.zeros(num_energies) +num_samples = 1000 +nickel['Es'] = 4.46 + +plt.figure() + +for energy_index, energy in enumerate(energies): + _, Y = run_rustbca(argon, nickel, energy, index, num_samples, run_sim, Eb=2.0, estop='LOW_ENERGY_NONLOCAL', weak_collisions=3, correction=1.09) + + sputtering_yield_2_eV[energy_index] = Y + + index += 1 + + _, Y = run_rustbca(argon, nickel, energy, index, num_samples, run_sim, Eb=0.0, estop='LOW_ENERGY_NONLOCAL', weak_collisions=3, correction=1.09) + + sputtering_yield_0_eV[energy_index] = Y + + index += 1 + +color = plt.loglog(energies, sputtering_yield_2_eV, label='RustBCA Kr-C Eb=2 eV')[0].get_color() +plt.loglog(energies, sputtering_yield_0_eV, label='RustBCA Kr-C Eb=0 eV', color=color, linestyle='--') +plt.loglog(yield_oksana_krc[:, 0], yield_oksana_krc[:, 1], label='OKSANA Kr-C') +plt.loglog(yield_trim_sp[:, 0], yield_trim_sp[:, 1], label='TRIM.SP Kr-C n=3') +plt.loglog(yield_srim[:, 0], yield_srim[:, 1], label='SRIM') +plt.scatter(yield_exp[:, 0], yield_exp[:, 1], label='Experiments compiled by V Shulga 2024', color='black') +plt.xlabel('E [eV]') +plt.ylabel('Y [at/ion]') +plt.title('Comparison to Fig. 8 V Shulga 2024 Ar->Ni Sputtering Yields') +plt.legend() + +plt.show() \ No newline at end of file From 259bf795598a9eb5707207076f00320bfe631e55 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 26 Aug 2025 16:52:51 -0700 Subject: [PATCH 24/31] Add test_cube.py to tests. --- .github/workflows/rustbca_compile_check.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 584476c..99de991 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -41,6 +41,7 @@ jobs: python3 -m pip install . python3 -c "from libRustBCA import *;" python3 examples/test_rustbca.py + python3 examples/test_cube.py - name: Test Fortran and C bindings run : | cargo build --release --lib --features parry3d From f18e381c29e661e4857570a0cb3ebd437934e64d Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 26 Aug 2025 18:48:31 -0700 Subject: [PATCH 25/31] Fix default Eb --- examples/test_sputtered_atom_spectra.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/test_sputtered_atom_spectra.py b/examples/test_sputtered_atom_spectra.py index 9dccda5..171bbd3 100644 --- a/examples/test_sputtered_atom_spectra.py +++ b/examples/test_sputtered_atom_spectra.py @@ -295,7 +295,7 @@ ]) #This function simply contains an entire input file as a multi-line f-string to modify some inputs. -def run_rustbca(ion, target, energy, index=0, num_samples=10000, run_sim=True, estop='LOW_ENERGY_NONLOCAL', weak_collisions=0, Eb=3.0, correction=1.0, interaction='KR_C'): +def run_rustbca(ion, target, energy, index=0, num_samples=10000, run_sim=True, estop='LOW_ENERGY_NONLOCAL', weak_collisions=0, Eb=0.0, correction=1.0, interaction='KR_C'): input_file = f''' [options] From b82d0a0ddba838411f2e0f8a8ae6807e45be19cc Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 23 Sep 2025 10:49:55 -0700 Subject: [PATCH 26/31] Add compound_bca_list_1d_py example. --- examples/compound_bca_list_1d_py.py | 116 ++++++++++++++++++++++++++++ 1 file changed, 116 insertions(+) create mode 100644 examples/compound_bca_list_1d_py.py diff --git a/examples/compound_bca_list_1d_py.py b/examples/compound_bca_list_1d_py.py new file mode 100644 index 0000000..e830335 --- /dev/null +++ b/examples/compound_bca_list_1d_py.py @@ -0,0 +1,116 @@ +from libRustBCA import * +import numpy as np +import matplotlib.pyplot as plt +import sys, os +sys.path.append(os.path.dirname(__file__)+'/../scripts') +sys.path.append('scripts') +from materials import * +from formulas import * + + +''' +This script is intended to serve as an example of using the Python +interface for a multi-layer, multi-species simulation. + +In this simulation, 50 He and 50 H ions are incident on a target +of two layers, one titanium boride and one aluminum: + +_____________ +| | +| TiB2 | dx = 100 A +_____________ +| | +| Al | dx = 1 um +| | + +''' +number_ions = 10000 +angle = 45.0 # angles in BCA codes are typically measured from the surface normal. +energy = 1000.0 # eV + +ux = [np.cos(angle*np.pi/180.0)]*number_ions +uy = [np.sin(angle*np.pi/180.0)]*number_ions +uz = [0.0]*number_ions + +energies = [energy]*number_ions +ion1 = hydrogen +ion2 = helium + +# Ion properties are per-ion; so we have a list of number_ions/2 of each: +Z1 = [ion1["Z"]]*(number_ions//2) + [ion2["Z"]]*(number_ions//2) +m1 = [ion1["m"]]*(number_ions//2) + [ion2["m"]]*(number_ions//2) +Ec1 = [ion1["Ec"]]*(number_ions//2) + [ion2["Ec"]]*(number_ions//2) +Es1 = [ion1["Es"]]*(number_ions//2) + [ion2["Es"]]*(number_ions//2) + +# Material properties are per species; so we have a list of 3 species: +Z2 = [titanium["Z"], boron["Z"], aluminum["Z"]] +m2 = [titanium["m"], boron["m"], aluminum["m"]] +Ec2 = [titanium["Ec"], boron["Ec"], aluminum["Ec"]] +Es2 = [titanium["Es"], boron["Es"], aluminum["Es"]] +Eb2 = [titanium["Eb"], boron["Eb"], aluminum["Eb"]] + +# Densities (n2) are specified as a list of layers, +# each of which has a list of the densities per-species: + +#top layer, titanium diboride: +# n_i calculated from ni = rho_TiB2 / (m_Ti + 2 * m_B) +ni = 3.8e28 # 1/m3 + +nTi1 = ni/10**30 # 1/A^3 +nB1 = 2 * ni/10**30 # 1/A^3 +nAl1 = 0.0 + +# bottom layer, pure aluminum: +nTi2 = 0.0 +nB2 = 0.0 +nAl2 = aluminum["n"]/10**30 # 1/A^3 + +n2 = [ + [nTi1, nB1, nAl1], # top layer + [nTi2, nB2, nAl2] # bottom layer +] + +dx = [ + 100.0, # Angstrom, top layer + 1.0*1e-6/1e-10 # Angstrom; bottom layer +] + +# compound_bca_list_py provides three return values +output, incident, stopped = compound_bca_list_1D_py( + ux, + uy, + uz, + energies, + Z1, + m1, + Ec1, + Es1, + Z2, + m2, + Ec2, + Es2, + Eb2, + n2, + dx +) + +# output columns = [Z, m (amu), E (eV), x, y, z, (angstrom), ux, uy, uz] +output = np.array(output) +Z = output[:, 0] +E = output[:, 0] + +# implanted ions can be selected using the stopped value +x_He = output[np.logical_and(Z == helium["Z"], stopped), 3] +x_H = output[np.logical_and(Z == hydrogen["Z"], stopped), 3] + +num_bins = 50 +bins = np.linspace(0.0, 400.0, num_bins) +plt.hist(x_He, bins=bins, histtype='step', label='He') +plt.hist(x_H, bins=bins, histtype='step', label='H') +plt.plot([dx[0], dx[0]], [0.0, number_ions], linestyle='--', color='gray') +plt.ylim([0.0, 375]) +plt.xlabel('x [A]') +plt.ylabel('f(x) [counts]') +plt.title('H/He implantation in TiB2 on Al') +plt.legend() +plt.show() \ No newline at end of file From 377c2c5e1beb00a141fa13b500cd05920eb3cce0 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 23 Sep 2025 13:51:21 -0700 Subject: [PATCH 27/31] Fix to multiple interaction potentials; somehow was using 1.0 geometry format.: --- examples/multiple_interaction_potentials.toml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index ba1aea1..aa8dc2d 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -43,13 +43,14 @@ pos = [ [ -1.7453292519934434e-8, 0.0, 0.0,], [ -1.7453292519934434e-8, 0.0, 0.0 dir = [ [ 1.0, 0.0, 0.0,], [ 1.0, 0.0, 0.0,],] [geometry_input] -energy_barrier_thickness = 2.2E-4 length_unit = "MICRON" -triangles = [ [ 0.0, 0.01, 0.0, 0.5, 0.5, -0.5,], [ 0.0, 0.01, 0.01, -0.5, 0.5, -0.5,], [ 0.01, 0.01, 0.04, -0.5, 0.5, -0.5,], [ 0.01, 0.04, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.04, 0.5, 0.5, -0.5,], [ 0.04, 0.5, 0.5, -0.5, 0.5, -0.5,],] +points = [[0.0, -0.5], [0.01, -0.5], [0.04, -0.5], [0.5, -0.5], [0.5, 0.5], [0.04, 0.5], [0.01, 0.5], [0.0, 0.5]] +triangles = [[0, 1, 6], [0, 6, 7], [1, 2, 5], [1, 5, 6], [2, 3, 4], [2, 4, 5]] densities = [ [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 3.0e+10, 6.0e+10, 0.0, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 6.026e+10, 0.0,], [ 0.0, 0.0, 0.0, 4.996e+10,], [ 0.0, 0.0, 0.0, 4.996e+10,],] -material_boundary_points = [ [ 0.0, 0.5,], [ 0.5, 0.5,], [ 0.5, -0.5,], [ 0.0, -0.5,],] +boundary = [0, 3, 4, 7] simulation_boundary_points = [ [ 0.6, -0.6,], [ -0.1, -0.6,], [ -0.1, 0.6,], [ 0.6, 0.6,], [ 0.6, -0.6,],] electronic_stopping_correction_factors = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0] +energy_barrier_thickness = 2.2E-4 [material_parameters] energy_unit = "EV" From 6c9f46aa52602369d9ab13ef6d7a818357502aaa Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 23 Sep 2025 15:20:20 -0700 Subject: [PATCH 28/31] Add example using tomllib --- examples/make_input_file_and_run.py | 258 ++++++++++++++++++++++++++++ 1 file changed, 258 insertions(+) create mode 100644 examples/make_input_file_and_run.py diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py new file mode 100644 index 0000000..68e5013 --- /dev/null +++ b/examples/make_input_file_and_run.py @@ -0,0 +1,258 @@ +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 formulas import * +import time +from tomlkit import parse, dumps + +''' +This script is a first draft of a comprehensive RustBCA input +file creation script using tomlkit. + +It includes two geometry modes, 1D and 0D. + +It simulates the following situations: + +if mode == '1D': + H+ (1 keV) + | + V +__________ +| | +| B | dx = 100 A +_________| +| | +| TiB2 | dx = 100 A +_________| +| | +| Ti | dx = 1000 A + + +if mode == '0D': + H+ (1 keV) + | + V +__________ +| | +| TiB2 | +| | +| | + +And calculates implantation profiles, reflection coefficients, +and sputtering yields. It also uses the ergonomic python functions +to compare the result of using the default values for H on B with +the custom values of this input file. + +It creates an input file as a nested dictionary which is written to +a TOML file using tomlkit. + +It runs the input file with cargo run --release and reads the output files. +''' +run_sim = True +mode = '1D' + +# species definitions +ion = hydrogen +target1 = boron +target2 = titanium + +# geometry definitions +n_i = 0.0328 # 1 / A^3 from n_i = rho_TiB2 / (mB * 2 + mTi) +layer_thicknesses = [100.0, 100.0, 1000.0] # A +layer_1_densities = [boron["n"]/10**30, 0.0] # 1/A^3 +layer_2_densities = [n_i * 2, n_i] # 1/A^3 +layer_3_densities = [0.0, titanium["n"]/10**30] # 1/A^3 + +# displacement energies +target1["Ed"] = 25.0 # eV +target2["Ed"] = 19.0 # eV + +incident_energy = 1000.0 # eV +number_ions = 100000 +angle = 45.0 # degrees; measured from surface normal + +options = { + 'name': 'input_file', + 'track_trajectories': False, # whether to track trajectories for plotting; memory intensive + 'track_recoils': True, # whether to track recoils; must enable for sputtering + 'track_recoil_trajectories': False, # whether to track recoil trajectories for plotting + 'track_displacements': False, # whether to track collisions with T > Ed for each species + 'track_energy_losses': False, # whether to track detailed collision energies; memory intensive + 'write_buffer_size': 8192, # how big the buffer is for file writing + 'weak_collision_order': 0, # weak collisions at radii (k + 1)*r; enable only when required + 'suppress_deep_recoils': False, # suppress recoils too deep to ever sputter + 'high_energy_free_flight_paths': False, # SRIM-style high energy free flight distances; use with caution + 'num_threads': 4, # number of threads to run in parallel + 'num_chunks': 10, # code will write to file every nth chunk; for very large simulations, increase num_chunks + 'electronic_stopping_mode': 'LOW_ENERGY_NONLOCAL', + 'mean_free_path_model': 'LIQUID', # liquid is amorphous (constant mean free path); gas is exponentially-distributed mean free paths + 'interaction_potential': [['KR_C']], + 'scattering_integral': [ + [ + { + 'GAUSS_MEHLER': {'n_points': 6} + } + ] + ], + + 'root_finder': [ + [ + { + 'NEWTON': { + 'max_iterations': 100, + 'tolerance': 1e-6 + } + } + ] + ], +} + +# material parameters are per-species +material_parameters = { + 'energy_unit': 'EV', + 'mass_unit': 'AMU', + # bulk binding energy; typically zero as a model choice + 'Eb': [ + target1["Eb"], + target2["Eb"] + ], + # surface binding energy + 'Es': [ + target1["Es"], + target2["Es"] + ], + # cutoff energy - particles with E < Ec stop + 'Ec': [ + target1["Ec"], + target2["Ec"] + ], + # displacement energy - only used to track displacements + 'Ed': [ + target1["Ed"], + target2["Ed"] + ], + # atomic number + 'Z': [ + target1["Z"], + target2["Z"] + ], + # atomic mass + 'm': [ + target1["m"], + target2["m"] + ], + # used to pick interaction potential from matrix in [options] + 'interaction_index': [0, 0], + 'surface_binding_model': { + "PLANAR": {'calculation': "INDIVIDUAL"} + }, + 'bulk_binding_model': 'INDIVIDUAL' +} + +particle_parameters = { + 'length_unit': 'ANGSTROM', + 'energy_unit': 'EV', + 'mass_unit': 'AMU', + # number of computational ions of this species to run at this energy + 'N': [number_ions], + # atomic mass + 'm': [ion["m"]], + # atomic number + 'Z': [ion["Z"]], + # incidenet energy + 'E': [incident_energy], + # cutoff energy - if E < Ec, particle stops + 'Ec': [ion["Ec"]], + # surface binding energy + 'Es': [ion["Es"]], + # initial position - if Es significant and E low, start (n)^(-1/3) above surface + # otherwise 0, 0, 0 is fine; most geometry modes have surface at x=0 with target x>0 + 'pos': [[0.0, 0.0, 0.0]], + # initial direction unit vector; most geometry modes have x-axis into the surface + 'dir': [ + [ + np.cos(angle*np.pi/180.0), + np.sin(angle*np.pi/180.0), + 0.0 + ] + ], +} + +geometry_0D = { + 'length_unit': 'ANGSTROM', + # used to correct nonlocal stopping for known compound discrpancies + 'electronic_stopping_correction_factor': 1.0, + # number densities of each species + 'densities': [2 * n_i, n_i] +} + +geometry_1D = { + 'length_unit': 'ANGSTROM', + # used to correct nonlocal stopping for known compound discrpancies + 'electronic_stopping_correction_factors': [1.0, 1.0, 1.0], + # thickness of each layer in order from top (x=0) to bottom + 'layer_thicknesses': layer_thicknesses, + # number densitiy of each layer in order from top to bottom + 'densities': [ + layer_1_densities, + layer_2_densities, + layer_3_densities, + ] +} + +if mode == '1D': + input_data = { + 'options': options, + 'material_parameters': material_parameters, + 'particle_parameters': particle_parameters, + 'geometry_input': geometry_1D + } +elif mode == '0D': + input_data = { + 'options': options, + 'material_parameters': material_parameters, + 'particle_parameters': particle_parameters, + 'geometry_input': geometry_0D +} + +# Attempt to cleanup line endings +input_string = dumps(input_data).replace('\r', '') +with open('examples/input_file.toml', 'w') as input_file: + input_file.write(input_string) + +if run_sim: + os.system(f'cargo run --release {mode} examples/input_file.toml') + +# Read output files - ensure arrays are at least 2D for indexing +sputtered = np.atleast_2d(np.genfromtxt('input_filesputtered.output', delimiter=',')) +reflected = np.atleast_2d(np.genfromtxt('input_filereflected.output', delimiter=',')) +implanted = np.atleast_2d(np.genfromtxt('input_filedeposited.output', delimiter=',')) + +print('H on B-TiB2-Ti') +print(f'R_N, R_E: {np.size(reflected[:, 0])/number_ions}, {np.sum(reflected[:, 2]/incident_energy/number_ions)}') +print(f'Y: {np.size(sputtered[:, 0])/number_ions} [at./ion]') +print() +print('H on B (default settings)') +print(f'Y: {sputtering_yield(hydrogen, boron, incident_energy, angle, 10000)} [at./ion]') +print(f'R_N, R_E: {reflection_coefficient(hydrogen, boron, incident_energy, angle, 10000)}') + +x = implanted[:, 2] +num_bins=100 +bins = np.linspace(0.0, 500.0, num_bins) +plt.hist(x, bins=bins, histtype='step') +if mode == '1D': + plt.plot([layer_thicknesses[0], layer_thicknesses[0]], [0.0, number_ions], color='gray') + plt.plot([layer_thicknesses[0] + layer_thicknesses[1], layer_thicknesses[0] + layer_thicknesses[1]], [0.0, number_ions], linestyle='--', color='gray') + +plt.ylim([0.0, 3000]) +plt.xlim([0.0, 500.0]) +plt.xlabel('x [A]') +plt.ylabel('f(x) [counts]') +plt.title('Implantation Distribution H on B-TiB2-Ti') +plt.show() \ No newline at end of file From 61909286eb7f0e4540ba5877fac2237056ca258a Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 23 Sep 2025 15:23:00 -0700 Subject: [PATCH 29/31] Add example using tomllib --- examples/make_input_file_and_run.py | 38 +++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py index 68e5013..f32ad06 100644 --- a/examples/make_input_file_and_run.py +++ b/examples/make_input_file_and_run.py @@ -6,8 +6,6 @@ #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 formulas import * import time from tomlkit import parse, dumps @@ -57,6 +55,42 @@ run_sim = True mode = '1D' +# For organizational purposes, species are commonly defined in dictionaries. +# Additional examples can be found in scripts/materials.py, but values +# should be checked for correctness before use. +hydrogen = { + 'symbol': 'H', + 'name': 'hydrogen', + 'Z': 1.0, + 'm': 1.008, # AMU + 'Ec': 0.95, # eV + 'Es': 1.5, # eV +} + +titanium = { + 'symbol': 'Ti', + 'name': 'titanium', + 'Z': 22.0, + 'm': 47.867, # AMU + 'Es': 4.84, # eV + 'Ec': 3.5, # eV + 'Eb': 0., # eV + 'Ed': 19.0, # eV + 'n': 5.67e28, # 1/m^3 +} + +boron = { + 'symbol': 'B', + 'name': 'boron', + 'Z': 5.0, + 'm': 10.811, # AMU + 'n': 1.309E29, # 1/m^3 + 'Es': 5.77, # eV + 'Eb': 0., # eV + 'Ec': 5., # eV + 'Ed': 25.0 # eV +} + # species definitions ion = hydrogen target1 = boron From 51673c04d2d0f1249f6b468eddfc229f1d93d73a Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 23 Sep 2025 15:24:35 -0700 Subject: [PATCH 30/31] Adjust comments on tomllib example. --- examples/make_input_file_and_run.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/examples/make_input_file_and_run.py b/examples/make_input_file_and_run.py index f32ad06..9a63f8b 100644 --- a/examples/make_input_file_and_run.py +++ b/examples/make_input_file_and_run.py @@ -52,12 +52,19 @@ It runs the input file with cargo run --release and reads the output files. ''' + run_sim = True mode = '1D' +incident_energy = 1000.0 # eV +number_ions = 100000 # at least 10k are typically needed for decent results +angle = 45.0 # degrees; measured from surface normal -# For organizational purposes, species are commonly defined in dictionaries. -# Additional examples can be found in scripts/materials.py, but values -# should be checked for correctness before use. +''' +For organizational purposes, species are commonly defined in dictionaries. +Additional examples can be found in scripts/materials.py, but values +should be checked for correctness before use. Values are explained +in the relevant sections below. +''' hydrogen = { 'symbol': 'H', 'name': 'hydrogen', @@ -103,14 +110,6 @@ layer_2_densities = [n_i * 2, n_i] # 1/A^3 layer_3_densities = [0.0, titanium["n"]/10**30] # 1/A^3 -# displacement energies -target1["Ed"] = 25.0 # eV -target2["Ed"] = 19.0 # eV - -incident_energy = 1000.0 # eV -number_ions = 100000 -angle = 45.0 # degrees; measured from surface normal - options = { 'name': 'input_file', 'track_trajectories': False, # whether to track trajectories for plotting; memory intensive From e01474c8075f29662bd7b8238e499ce8d6f9ad6c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 23 Sep 2025 15:25:56 -0700 Subject: [PATCH 31/31] Add tomlkit example to workflow.: --- .github/workflows/rustbca_compile_check.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 99de991..2830b7f 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -30,7 +30,7 @@ jobs: sudo apt-get install python3-pip python3-dev - name: Install Python libraries run: | - python3 -m pip install numpy shapely scipy matplotlib toml + python3 -m pip install numpy shapely scipy matplotlib toml tomlkit - name: Install HDF5 Libraries run: | sudo apt install libhdf5-dev @@ -42,6 +42,7 @@ jobs: python3 -c "from libRustBCA import *;" python3 examples/test_rustbca.py python3 examples/test_cube.py + python3 examples/make_input_file_and_run.py - name: Test Fortran and C bindings run : | cargo build --release --lib --features parry3d