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/test_cube.py b/examples/test_cube.py index 13e6f7b..6d41713 100644 --- a/examples/test_cube.py +++ b/examples/test_cube.py @@ -115,10 +115,10 @@ def main(): run_sim = True track_trajectories = False a = 1000 - num_samples = 100000 + num_samples = 10000 num_angles = 10 - energy = 100 - angles = np.linspace(0.01, 89.9, num_angles) + energy = 200 + angles = np.linspace(0.0, 89.9, num_angles) sputtering_yields_x = np.zeros(num_angles) sputtering_yields_z = np.zeros(num_angles) @@ -138,43 +138,48 @@ def main(): sim_index = 0 impl_index = num_angles*3//4 + bins = np.linspace(0, 500, 100) + for angle_index, angle in enumerate(angles): R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, ux=np.cos(angle*np.pi/180.), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2)) sim_index += 1 plt.figure(1) - if angle_index == impl_index: plt.hist(implanted[:, 2], histtype='step', bins=100, density=False, label='x') - + if angle_index == impl_index: plt.hist(implanted[:, 2], histtype='step', bins=bins, density=False, label='x') sputtering_yields_x[angle_index] = Y R_x[angle_index] = R R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a, ux=-np.cos(angle*np.pi/180.), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2)) sim_index += 1 - if angle_index == impl_index: plt.hist(a - implanted[:, 2], histtype='step', bins=100, density=False, label='-x') + if angle_index == impl_index: plt.hist(a - implanted[:, 2], histtype='step', bins=bins, density=False, label='-x') sputtering_yields_x_minus[angle_index] = Y R_x_minus[angle_index] = R R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=a/2., z0=a, ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=-np.cos(angle*np.pi/180.)) sim_index += 1 - if angle_index == impl_index: plt.hist(a - implanted[:, 4], histtype='step', bins=100, density=False, label='-z') + if angle_index == impl_index: plt.hist(a - implanted[:, 4], histtype='step', bins=bins, density=False, label='-z') + sputtering_yields_z_minus[angle_index] = Y R_z_minus[angle_index] = R R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=a/2., z0=0.0, ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.cos(angle*np.pi/180.)) sim_index += 1 - if angle_index == impl_index: plt.hist(implanted[:, 4], histtype='step', bins=100, density=False, label='z') + if angle_index == impl_index: plt.hist(implanted[:, 4], histtype='step', bins=bins, density=False, label='z') + sputtering_yields_z[angle_index] = Y R_z[angle_index] = R R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=0.0, z0=a/2., ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=np.cos(angle*np.pi/180.)) sim_index += 1 - if angle_index == impl_index: plt.hist(implanted[:, 3], histtype='step', bins=100, density=False, label='y') + if angle_index == impl_index: plt.hist(implanted[:, 3], histtype='step', bins=bins, density=False, label='y') + sputtering_yields_y[angle_index] = Y R_y[angle_index] = R R, Y, reflected, sputtered, implanted = run(energy, sim_index, num_samples, track_trajectories=track_trajectories, run_sim=run_sim, x0=a/2., y0=a, z0=a/2., ux=np.sin(angle*np.pi/180.)/np.sqrt(2), uz=np.sin(angle*np.pi/180.)/np.sqrt(2), uy=-np.cos(angle*np.pi/180.)) sim_index += 1 - if angle_index == impl_index: plt.hist(a - implanted[:, 3], histtype='step', bins=100, density=False, label='-y') + if angle_index == impl_index: plt.hist(a - implanted[:, 3], histtype='step', bins=bins, density=False, label='-y') + sputtering_yields_y_minus[angle_index] = Y R_y_minus[angle_index] = R @@ -206,7 +211,6 @@ def main(): plt.show() if track_trajectories and do_plots: - do_trajectory_plot('cube_19', boundary=[[0.0, 0.0], [0.0, a], [a, a], [a, 0.0], [0.0, 0.0]]) do_trajectory_plot('cube_20', boundary=[[0.0, 0.0], [0.0, a], [a, a], [a, 0.0], [0.0, 0.0]]) diff --git a/examples/tungsten_twist_trimesh.toml b/examples/tungsten_twist_trimesh.toml 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] diff --git a/src/bca.rs b/src/bca.rs index 7d3f013..353c24d 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -127,9 +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); + // Since psi and psi_r are absolute-valued in the BC calculation, note + // that it is negated here to correctly rotate the recoil in the correct + // direction (that is, away from the direction the incident atom is deflected) particle_2.rotate(-binary_collision_result.psi_recoil, binary_collision_geometry.phi_azimuthal); @@ -360,10 +364,28 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma let sinx: f64 = (1. - cosx*cosx).sqrt(); let cosphi: f64 = phi_azimuthal.cos(); - //Find recoil location - let x_recoil: f64 = x + mfp*cosx - impact_parameter*cosphi*sinx; - let y_recoil: f64 = y + mfp*cosy - impact_parameter*(sinphi*cosz - cosphi*cosy*cosx)/sinx; - let z_recoil: f64 = z + mfp*cosz + impact_parameter*(sinphi*cosy - cosphi*cosx*cosz)/sinx; + // These formulas find the recoil one mfp away at an impact parameter p at angle phi + // To resolve the singularity, a different set of rotations is used when cosx == -1 + // Because of this, the recoil location is not consistent between the two formulas at a given phi + // Since phi is sampled uniformly from (0, 2pi), this does not matter + // However, if a crystalline structure is ever added, this needs to be considered + let x_recoil = if cosx > -1. { + x + mfp*cosx - impact_parameter*(cosz*sinphi + cosy*cosphi) + } else { + x + mfp*cosx - impact_parameter*((1. + cosz - cosx*cosx)*cosphi - cosx*cosy*sinphi)/(1. + cosz) + }; + + let y_recoil = if cosx > -1. { + y + mfp*cosy + impact_parameter*((1. + cosx - cosy*cosy)*cosphi - cosy*cosz*sinphi)/(1. + cosx) + } else { + y + mfp*cosy + impact_parameter*((1. + cosz - cosy*cosy)*sinphi - cosx*cosy*cosphi)/(1. + cosz) + }; + + let z_recoil = if cosx > -1. { + z + mfp*cosz + impact_parameter*((1. + cosx - cosz*cosz)*sinphi - cosy*cosz*cosphi)/(1. + cosx) + } else { + z + mfp*cosz + impact_parameter*(cosx*cosphi + cosy*sinphi) + }; //Choose recoil Z, M let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, Ed_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil); @@ -504,8 +526,9 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: & InteractionPotential::COULOMB{..} => 0., _ => x0*a*(theta/2.).sin() }; - let psi = (theta.sin().atan2(Ma/Mb + theta.cos())).abs(); - let psi_recoil = (theta.sin().atan2(1. - theta.cos())).abs(); + + let psi = theta.sin().atan2(Ma/Mb + theta.cos());//.abs(); + let psi_recoil = theta.sin().atan2(1. - theta.cos());//.abs(); let recoil_energy = 4.*(Ma*Mb)/(Ma + Mb).powi(2)*E0*(theta/2.).sin().powi(2); Ok(BinaryCollisionResult::new(theta, psi, psi_recoil, recoil_energy, asymptotic_deflection, x0)) diff --git a/src/input.rs b/src/input.rs index af04e1f..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!(ux != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); ux} + Distributions::POINT(ux) => ux, }, uy: match cosy { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, @@ -654,7 +654,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static { ux: match cosx { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit}, - Distributions::POINT(x) => {assert!(x != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); x*length_unit}, + Distributions::POINT(x) => x*length_unit }, uy: match cosy { Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit}, diff --git a/src/lib.rs b/src/lib.rs index 4e52076..41ea254 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. @@ -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,12 +1511,9 @@ 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 rotation_matrix = if (c + 1.0).abs() > DELTA { - 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 @@ -1528,15 +1522,8 @@ 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 + 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 - ); - - *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(); @@ -1613,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 @@ -1623,18 +1608,15 @@ 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 Rotation3::from_axis_angle(&Vector3::y_axis(), PI).into() }; + // Note: transpose of R == R^-1 let u = rotation_matrix.transpose()*direction; *ux = u.x; @@ -1715,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.")); @@ -1734,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 { @@ -1806,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.")); @@ -1827,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); @@ -1908,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.")); @@ -1930,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); diff --git a/src/particle.rs b/src/particle.rs index 9de3daa..79cc09d 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 { @@ -239,16 +236,33 @@ 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 sa = (1. - cosx*cosx).sqrt(); + let cosphi: f64 = (phi + PI).cos(); + let sinphi: f64 = (phi + PI).sin(); - //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); + + // To resolve the singularity, a different set of rotations is used when cosx == -1 + // Because of this, the recoil location is not consistent between the two formulas at a given phi + // Since phi is sampled uniformly from (0, 2pi), this does not matter + // However, if a crystalline structure is ever added, this needs to be considered + let 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 fd263f0..b926984 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -761,62 +761,101 @@ 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(); - let cosz = 0.; - - 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}] { + + #[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 + 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.; + + 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, 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] { + + //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); @@ -834,7 +873,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, @@ -856,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, @@ -930,9 +969,10 @@ 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)); + //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.)); assert!(!particle_1.E.is_nan()); assert!(!particle_2.E.is_nan());