diff --git a/.github/ISSUE_TEMPLATE/benchmark_request.md b/.github/ISSUE_TEMPLATE/benchmark_request.md new file mode 100644 index 0000000..be36f65 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/benchmark_request.md @@ -0,0 +1,14 @@ +--- +name: Benchmark request +about: Suggest a benchmark for validation +title: "[benchmark]" +labels: benchmark +assignees: '' + +--- + +**Description** +A clear and concise description of the requested benchmark. + +**Source** +Source of data for benchmark; a paper, code, or database with sputtering yields, reflection coefficients, etc. diff --git a/examples/test_rustbca.py b/examples/test_rustbca.py index 8a9580d..cee6673 100644 --- a/examples/test_rustbca.py +++ b/examples/test_rustbca.py @@ -81,6 +81,49 @@ def main(): print(f'RustBCA R: {len(reflected[:, 0])/number_ions} Thomas R: {thomas}') print(f'Time per ion: {delta_time/number_ions*1e3} us/{ion["symbol"]}') + #Next up is the layered target version. I'll add a 50 Angstrom layer of W-H to the top of the target. + + #1 keV is above the He on W sputtering threshold of ~150 eV + energies_eV = 1000.0*np.ones(number_ions) + + #Working with angles of exactly 0 is problematic due to gimbal lock + angle = 0.0001 + + #In RustBCA's 0D geometry, +x -> into the surface + ux = np.cos(angle*np.pi/180.)*np.ones(number_ions) + uy = np.sin(angle*np.pi/180.)*np.ones(number_ions) + uz = np.zeros(number_ions) + + print(f'Running RustBCA for {number_ions} {ion["symbol"]} ions on {target["symbol"]} with hydrogenated layer at {energies_eV[0]/1000.} keV...') + print(f'This may take several minutes.') + #Not the different argument order; when a breaking change is due, this will + #be back-ported to the other bindings as well for consistency. + output, incident, stopped = compound_bca_list_1D_py( + ux, uy, uz, energies_eV, [ion['Z']]*number_ions, + [ion['m']]*number_ions, [ion['Ec']]*number_ions, [ion['Es']]*number_ions, [target['Z'], 1.0], [target['m'], 1.008], + [target['Ec'], 1.0], [target['Es'], 1.5], [target['Eb'], 0.0], [[target['n']/10**30, target['n']/10**30], [target['n']/10**30, 0.0]], [50.0, 1e6] + ) + + output = np.array(output) + + Z = output[:, 0] + m = output[:, 1] + E = output[:, 2] + x = output[:, 3] + y = output[:, 4] + z = output[:, 5] + ux = output[:, 6] + uy = output[:, 7] + uz = output[:, 8] + + plt.figure(3) + plt.title(f'Implanted {ion["symbol"]} Depth Distribution with 50A {target["symbol"]}-H layer') + plt.xlabel('x [A]') + plt.ylabel('f(x) [A.U.]') + heights, _, _ = plt.hist(x[np.logical_and(incident, stopped)], bins=100, density=True, histtype='step') + plt.plot([50.0, 50.0], [0.0, np.max(heights)*1.1]) + plt.gca().set_ylim([0.0, np.max(heights)*1.1]) + plt.show() if __name__ == '__main__': diff --git a/src/RustBCA.egg-info/PKG-INFO b/src/RustBCA.egg-info/PKG-INFO new file mode 100644 index 0000000..f9a0c9c --- /dev/null +++ b/src/RustBCA.egg-info/PKG-INFO @@ -0,0 +1,4 @@ +Metadata-Version: 2.1 +Name: RustBCA +Version: 1.2.0 +License-File: LICENSE diff --git a/src/RustBCA.egg-info/SOURCES.txt b/src/RustBCA.egg-info/SOURCES.txt new file mode 100644 index 0000000..e08c132 --- /dev/null +++ b/src/RustBCA.egg-info/SOURCES.txt @@ -0,0 +1,28 @@ +Cargo.toml +LICENSE +MANIFEST.in +README.md +pyproject.toml +setup.py +src/bca.rs +src/consts.rs +src/enums.rs +src/geometry.rs +src/gui.rs +src/input.rs +src/interactions.rs +src/lib.rs +src/main.rs +src/material.rs +src/output.rs +src/parry.rs +src/particle.rs +src/physics.rs +src/sphere.rs +src/structs.rs +src/tests.rs +src/RustBCA.egg-info/PKG-INFO +src/RustBCA.egg-info/SOURCES.txt +src/RustBCA.egg-info/dependency_links.txt +src/RustBCA.egg-info/not-zip-safe +src/RustBCA.egg-info/top_level.txt \ No newline at end of file diff --git a/src/RustBCA.egg-info/dependency_links.txt b/src/RustBCA.egg-info/dependency_links.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/RustBCA.egg-info/dependency_links.txt @@ -0,0 +1 @@ + diff --git a/src/RustBCA.egg-info/not-zip-safe b/src/RustBCA.egg-info/not-zip-safe new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/RustBCA.egg-info/not-zip-safe @@ -0,0 +1 @@ + diff --git a/src/RustBCA.egg-info/top_level.txt b/src/RustBCA.egg-info/top_level.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/src/RustBCA.egg-info/top_level.txt @@ -0,0 +1 @@ + diff --git a/src/bca.rs b/src/bca.rs index 544f4c3..32150af 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -75,9 +75,6 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate //Remove particle from top of vector as particle_1 let mut particle_1 = particles.pop().unwrap(); - - //println!("Particle start Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); - //BCA loop while !particle_1.stopped & !particle_1.left { @@ -114,8 +111,6 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate particle_1.pos.x, particle_2.pos.x, &binary_collision_geometry)) .unwrap(); - //println!("{}", binary_collision_result); - //Only use 0th order collision for local electronic stopping if k == 0 { normalized_distance_of_closest_approach = binary_collision_result.normalized_distance_of_closest_approach; @@ -127,17 +122,14 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate particle_2.E = binary_collision_result.recoil_energy - material.average_bulk_binding_energy(particle_2.pos.x, particle_2.pos.y, particle_2.pos.z); particle_2.energy_origin = particle_2.E; - //Accumulate asymptotic deflections for primary particle + //Accumulate energy losses and asymptotic deflections for primary particle total_energy_loss += binary_collision_result.recoil_energy; - - //total_deflection_angle += psi; total_asymptotic_deflection += binary_collision_result.asymptotic_deflection; - //Rotate particle 1, 2 by lab frame scattering angles - particle::rotate_particle(&mut particle_1, binary_collision_result.psi, + particle_1.rotate(binary_collision_result.psi, binary_collision_geometry.phi_azimuthal); - particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil, + particle_2.rotate(-binary_collision_result.psi_recoil, binary_collision_geometry.phi_azimuthal); particle_2.dir_old.x = particle_2.dir.x; @@ -184,20 +176,17 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate //Advance particle in space and track total distance traveled #[cfg(not(feature = "accelerated_ions"))] - let distance_traveled = particle::particle_advance(&mut particle_1, + let distance_traveled = particle_1.advance( binary_collision_geometries[0].mfp, total_asymptotic_deflection); #[cfg(feature = "accelerated_ions")] - let distance_traveled = particle::particle_advance(&mut particle_1, + let distance_traveled = particle_1.advance( binary_collision_geometries[0].mfp + distance_to_target - material.geometry.get_energy_barrier_thickness(), total_asymptotic_deflection); //Subtract total energy from all simultaneous collisions and electronic stopping bca::update_particle_energy(&mut particle_1, &material, distance_traveled, total_energy_loss, normalized_distance_of_closest_approach, strong_collision_Z, strong_collision_index, &options); - //println!("Particle finished collision loop Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); - - //println!("{} {} {}", energy_0/EV, energy_1/EV, (energy_1 - energy_0)/EV); //Check boundary conditions on leaving and stopping material::boundary_condition_planar(&mut particle_1, &material); @@ -205,7 +194,6 @@ pub fn single_ion_bca(particle: particle::Particle, material: &mate //Set particle index to topmost particle particle_index = particles.len(); } - //println!("Particle stopped or left Z: {} E: {} ({}, {}, {})", particle_1.Z, particle_1.E/Q, particle_1.pos.x/ANGSTROM, particle_1.pos.y/ANGSTROM, particle_1.pos.z/ANGSTROM); particle_output.push(particle_1); } particle_output diff --git a/src/consts.rs b/src/consts.rs index a0dfcd2..052627d 100644 --- a/src/consts.rs +++ b/src/consts.rs @@ -15,6 +15,8 @@ pub const MICRON: f64 = 1E-6; pub const NM: f64 = 1E-9; /// One centimeter in meters. pub const CM: f64 = 1E-2; +/// One milimter in meters. +pub const MM: f64 = 1E-3; /// Vacuum permitivity in Farads/meter. pub const EPS0: f64 = 8.8541878128E-12; /// Bohr radius in meters. diff --git a/src/geometry.rs b/src/geometry.rs index 462bad5..5d8cdd7 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -52,6 +52,7 @@ impl Geometry for Mesh0D { let length_unit: f64 = match input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., @@ -97,7 +98,6 @@ impl Geometry for Mesh0D { } fn inside_simulation_boundary(&self, x: f64, y: f64, z: f64) -> bool { - //println!("x: {} energy_barrier_thickness: {}", x/ANGSTROM, self.energy_barrier_thickness/ANGSTROM); x > -10.*self.energy_barrier_thickness } @@ -143,6 +143,7 @@ impl Geometry for Mesh1D { let length_unit: f64 = match geometry_input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., @@ -324,6 +325,7 @@ impl Geometry for Mesh2D { let length_unit: f64 = match geometry_input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., diff --git a/src/input.rs b/src/input.rs index 714c638..91693aa 100644 --- a/src/input.rs +++ b/src/input.rs @@ -381,6 +381,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static { let length_unit: f64 = match particle_parameters.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., diff --git a/src/interactions.rs b/src/interactions.rs index 3c81cb7..db965a9 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -8,7 +8,6 @@ pub fn crossing_point_doca(interaction_potential: InteractionPotential) -> f64 { InteractionPotential::MORSE{D, alpha, r0} => (alpha*r0 - (2.0_f64).ln())/alpha, InteractionPotential::WW => 50.*ANGSTROM, _ => 10.*ANGSTROM, - //_ => panic!("Input error: potential never crosses zero for r > 0. Consider using the Newton rootfinder.") } } diff --git a/src/lib.rs b/src/lib.rs index 32e0a46..095b5f3 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -80,6 +80,8 @@ pub use crate::parry::{ParryBall, ParryBallInput, InputParryBall, ParryTriMesh, pub fn pybca(py: Python, m: &PyModule) -> PyResult<()> { m.add_function(wrap_pyfunction!(simple_bca_py, m)?)?; m.add_function(wrap_pyfunction!(simple_bca_list_py, m)?)?; + m.add_function(wrap_pyfunction!(compound_bca_list_py, m)?)?; + m.add_function(wrap_pyfunction!(compound_bca_list_1D_py, m)?)?; Ok(()) } @@ -979,6 +981,404 @@ pub extern "C" fn simple_bca_c(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64 } } +#[cfg(feature = "python")] +///compound_tagged_bca_list_py(ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2) +/// 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 +/// uy (list(f64)): initial ion directions y. +/// uz (list(f64)): initial ion directions z. +/// Z1 (list(f64)): initial ion atomic numbers. +/// m1 (list(f64)): initial ion masses in amu. +/// Ec1 (list(f64)): ion cutoff energies in eV. If ion energy < Ec1, it stops in the material. +/// Es1 (list(f64)): ion surface binding energies. Assumed planar. +/// Z2 (list(f64)): target material species atomic numbers. +/// m2 (list(f64)): target material species masses in amu. +/// Ec2 (list(f64)): target material species cutoff energies in eV. If recoil energy < Ec2, it stops in the material. +/// Es2 (list(f64)): target species surface binding energies. Assumed planar. +/// n2 (list(f64)): target material species atomic number densities in inverse cubic Angstroms. +/// Eb2 (list(f64)): target material species bulk binding energies in eV. +/// Returns: +/// output (NX9 list of f64): each row in the list represents an output particle (implanted, +/// sputtered, or reflected). Each row consists of: +/// [Z, m (amu), E (eV), x, y, z, (angstrom), ux, uy, uz] +/// incident (list(bool)): whether each row of output was an incident ion or originated in the target +#[pyfunction] +pub fn compound_bca_list_py(energies: Vec, ux: Vec, uy: Vec, uz: Vec, Z1: Vec, m1: Vec, Ec1: Vec, Es1: Vec, Z2: Vec, m2: Vec, Ec2: Vec, Es2: Vec, n2: Vec, Eb2: Vec) -> (Vec<[f64; 9]>, Vec) { + let mut total_output = vec![]; + let mut incident = vec![]; + let num_species_target = Z2.len(); + let num_incident_ions = energies.len(); + + assert_eq!(ux.len(), num_incident_ions, "Input error: list of x-directions is not the same length as list of incident energies."); + assert_eq!(uy.len(), num_incident_ions, "Input error: list of y-directions is not the same length as list of incident energies."); + assert_eq!(uz.len(), num_incident_ions, "Input error: list of z-directions is not the same length as list of incident energies."); + assert_eq!(Z1.len(), num_incident_ions, "Input error: list of incident atomic numbers is not the same length as list of incident energies."); + assert_eq!(m1.len(), num_incident_ions, "Input error: list of incident atomic masses is not the same length as list of incident energies."); + assert_eq!(Es1.len(), num_incident_ions, "Input error: list of incident surface binding energies is not the same length as list of incident energies."); + assert_eq!(Ec1.len(), num_incident_ions, "Input error: list of incident cutoff energies is not the same length as list of incident energies."); + + assert_eq!(m2.len(), num_species_target, "Input error: list of target atomic masses is not the same length as atomic numbers."); + assert_eq!(Ec2.len(), num_species_target, "Input error: list of target cutoff energies is not the same length as atomic numbers."); + assert_eq!(Es2.len(), num_species_target, "Input error: list of target surface binding energies is not the same length as atomic numbers."); + assert_eq!(Eb2.len(), num_species_target, "Input error: list of target bulk binding energies is not the same length as atomic numbers."); + assert_eq!(n2.len(), num_species_target, "Input error: list of target number densities is not the same length as atomic numbers."); + + #[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: 3, + suppress_deep_recoils: true, + high_energy_free_flight_paths: true, + electronic_stopping_mode: ElectronicStoppingMode::LOW_ENERGY_NONLOCAL, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::KR_C]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::DEFAULTNEWTON]], + 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, + }; + + #[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: 3, + suppress_deep_recoils: true, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::LOW_ENERGY_NONLOCAL, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::KR_C]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::DEFAULTNEWTON]], + track_displacements: false, + track_energy_losses: false, + }; + + let x = -2.*(n2.iter().sum::()*10E30).powf(-1./3.); + let y = 0.0; + let z = 0.0; + + let material_parameters = material::MaterialParameters { + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: Eb2, + Es: Es2, + Ec: Ec2, + Z: Z2, + m: m2, + interaction_index: vec![0; num_species_target], + surface_binding_model: SurfaceBindingModel::INDIVIDUAL, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let geometry_input = geometry::Mesh0DInput { + length_unit: "ANGSTROM".to_string(), + densities: n2, + electronic_stopping_correction_factor: 1.0 + }; + + let m = material::Material::::new(&material_parameters, &geometry_input); + + let mut index: usize = 0; + for (((((((E1_, ux_), uy_), uz_), Z1_), Ec1_), Es1_), m1_) in energies.iter().zip(ux).zip(uy).zip(uz).zip(Z1).zip(Ec1).zip(Es1).zip(m1) { + + let mut energy_out; + let p = particle::Particle { + m: m1_*AMU, + Z: Z1_, + E: E1_*EV, + Ec: Ec1_*EV, + Es: Es1_*EV, + pos: Vector::new(x, y, z), + dir: Vector::new(ux_, uy_, uz_), + pos_origin: Vector::new(x, y, z), + pos_old: Vector::new(x, y, z), + dir_old: Vector::new(ux_, uy_, uz_), + energy_origin: E1_*EV, + asymptotic_deflection: 0.0, + stopped: false, + left: false, + incident: true, + first_step: true, + trajectory: vec![], + energies: vec![], + track_trajectories: false, + number_collision_events: 0, + backreflected: false, + interaction_index : 0, + weight: 0.0, + tag: 0, + tracked_vector: Vector::new(0., 0., 0.), + }; + + let output = bca::single_ion_bca(p, &m, &options); + + for particle in output { + if (particle.left) | (particle.incident) { + + incident.push(particle.incident); + + if particle.stopped { + energy_out = 0. + } else { + energy_out = particle.E/EV + } + total_output.push( + [ + particle.Z, + particle.m/AMU, + energy_out, + particle.pos.x, + particle.pos.y, + particle.pos.z, + particle.dir.x, + particle.dir.y, + particle.dir.z, + ] + ); + } + } + index += 1; + } + (total_output, incident) +} + +#[cfg(feature = "python")] +///compound_bca_list_1D_py(ux, uy, uz, energy, 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: +/// energies (list(f64)): initial ion energies in eV. +/// ux (list(f64)): initial ion directions x. ux != 0.0 to avoid gimbal lock +/// uy (list(f64)): initial ion directions y. +/// uz (list(f64)): initial ion directions z. +/// Z1 (list(f64)): initial ion atomic numbers. +/// m1 (list(f64)): initial ion masses in amu. +/// Ec1 (list(f64)): ion cutoff energies in eV. If ion energy < Ec1, it stops in the material. +/// Es1 (list(f64)): ion surface binding energies. Assumed planar. +/// Z2 (list(f64)): target material species atomic numbers. +/// m2 (list(f64)): target material species masses in amu. +/// Ec2 (list(f64)): target material species cutoff energies in eV. If recoil energy < Ec2, it stops in the material. +/// Es2 (list(f64)): target species surface binding energies. Assumed planar. +/// Eb2 (list(f64)): target material species bulk binding energies in eV. +/// n2 (list(list(f64))): target material species atomic number densities in inverse cubic Angstroms. +/// dx (list(f64)): target material layer thicknesses starting at surface. +/// Returns: +/// output (NX9 list of f64): each row in the list represents an output particle (implanted, +/// sputtered, or reflected). Each row consists of: +/// [Z, m (amu), E (eV), x, y, z, (angstrom), ux, uy, uz] +/// incident (list(bool)): whether each row of output was an incident ion or originated in the target +/// stopped (list(bool)): whether each row of output is associated with a particle that stopped in the target +#[pyfunction] +pub fn compound_bca_list_1D_py(ux: Vec, uy: Vec, uz: Vec, energies: Vec, Z1: Vec, m1: Vec, Ec1: Vec, Es1: Vec, Z2: Vec, m2: Vec, Ec2: Vec, Es2: Vec, Eb2: Vec, n2: Vec>, dx: Vec) -> (Vec<[f64; 9]>, Vec, Vec) { + let mut total_output = vec![]; + let mut incident = vec![]; + let mut stopped = vec![]; + + let num_layers_target = n2.len(); + let num_species = Z2.len(); + let num_incident_ions = energies.len(); + + assert_eq!(ux.len(), num_incident_ions, "Input error: list of x-directions is not the same length as list of incident energies."); + assert_eq!(uy.len(), num_incident_ions, "Input error: list of y-directions is not the same length as list of incident energies."); + assert_eq!(uz.len(), num_incident_ions, "Input error: list of z-directions is not the same length as list of incident energies."); + assert_eq!(Z1.len(), num_incident_ions, "Input error: list of incident atomic numbers is not the same length as list of incident energies."); + assert_eq!(m1.len(), num_incident_ions, "Input error: list of incident atomic masses is not the same length as list of incident energies."); + assert_eq!(Es1.len(), num_incident_ions, "Input error: list of incident surface binding energies is not the same length as list of incident energies."); + assert_eq!(Ec1.len(), num_incident_ions, "Input error: list of incident cutoff energies is not the same length as list of incident energies."); + + assert_eq!(m2.len(), num_species, "Input error: list of target atomic masses is not the same length as atomic numbers."); + assert_eq!(Ec2.len(), num_species, "Input error: list of target cutoff energies is not the same length as atomic numbers."); + assert_eq!(Es2.len(), num_species, "Input error: list of target surface binding energies is not the same length as atomic numbers."); + assert_eq!(Eb2.len(), num_species, "Input error: list of target bulk binding energies is not the same length as atomic numbers."); + assert_eq!(n2[0].len(), num_species, "Input error: first layer list of target number densities is not the same length as atomic numbers."); + + assert_eq!(n2[0].len(), num_species, "Input error: first layer species list of target number densities is not the same length as atomic numbers."); + assert_eq!(dx.len(), num_layers_target, "Input error: number of layer thicknesses not the same as number of layers in atomic densities list."); + + #[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: 3, + suppress_deep_recoils: true, + high_energy_free_flight_paths: true, + electronic_stopping_mode: ElectronicStoppingMode::LOW_ENERGY_NONLOCAL, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::KR_C]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::DEFAULTNEWTON]], + 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, + }; + + #[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: 3, + suppress_deep_recoils: true, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::LOW_ENERGY_NONLOCAL, + mean_free_path_model: MeanFreePathModel::LIQUID, + interaction_potential: vec![vec![InteractionPotential::KR_C]], + scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], + num_threads: 1, + num_chunks: 1, + use_hdf5: false, + root_finder: vec![vec![Rootfinder::DEFAULTNEWTON]], + track_displacements: false, + track_energy_losses: false, + }; + + let y = 0.0; + let z = 0.0; + + let material_parameters = material::MaterialParameters { + energy_unit: "EV".to_string(), + mass_unit: "AMU".to_string(), + Eb: Eb2, + Es: Es2, + Ec: Ec2, + Z: Z2, + m: m2, + interaction_index: vec![0; num_species], + surface_binding_model: SurfaceBindingModel::INDIVIDUAL, + bulk_binding_model: BulkBindingModel::INDIVIDUAL, + }; + + let geometry_input = geometry::Mesh1DInput { + length_unit: "ANGSTROM".to_string(), + densities: n2, + layer_thicknesses: dx, + electronic_stopping_correction_factors: vec![1.0; num_layers_target] + }; + + let m = material::Material::::new(&material_parameters, &geometry_input); + + let x = -m.geometry.top_energy_barrier_thickness/2.; + + let mut index: usize = 0; + for (((((((E1_, ux_), uy_), uz_), Z1_), Ec1_), Es1_), m1_) in energies.iter().zip(ux).zip(uy).zip(uz).zip(Z1).zip(Ec1).zip(Es1).zip(m1) { + + let mut dir = Vector::new(ux_, uy_, uz_); + dir.normalize(); + + let mut energy_out; + let p = particle::Particle { + m: m1_*AMU, + Z: Z1_, + E: E1_*EV, + Ec: Ec1_*EV, + Es: Es1_*EV, + pos: Vector::new(x, y, z), + dir: dir, + pos_origin: Vector::new(x, y, z), + pos_old: Vector::new(x, y, z), + dir_old: dir, + energy_origin: E1_*EV, + asymptotic_deflection: 0.0, + stopped: false, + left: false, + incident: true, + first_step: true, + trajectory: vec![], + energies: vec![], + track_trajectories: false, + number_collision_events: 0, + backreflected: false, + interaction_index : 0, + weight: 0.0, + tag: 0, + tracked_vector: Vector::new(0., 0., 0.), + }; + + let output = bca::single_ion_bca(p, &m, &options); + + for particle in output { + if (particle.left) | (particle.incident) { + + incident.push(particle.incident); + stopped.push(particle.stopped); + + if particle.stopped { + energy_out = 0. + } else { + energy_out = particle.E/EV + } + total_output.push( + [ + particle.Z, + particle.m/AMU, + energy_out, + particle.pos.x/ANGSTROM, + particle.pos.y/ANGSTROM, + particle.pos.z/ANGSTROM, + particle.dir.x, + particle.dir.y, + particle.dir.z, + ] + ); + } + } + index += 1; + } + (total_output, incident, stopped) +} + #[cfg(feature = "python")] /// simple_bca_py( x, y, z, ux, uy, uz, energy, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2) /// -- diff --git a/src/material.rs b/src/material.rs index 8be2149..9b3eb86 100644 --- a/src/material.rs +++ b/src/material.rs @@ -323,7 +323,6 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, //Actual surface binding energies let Es = material.actual_surface_binding_energy(particle_1, x_old, y_old, z_old); - //println!("Actual Es: {}", Es); let Ec = particle_1.Ec; let inside_now = material.inside_energy_barrier(x, y, z); @@ -380,7 +379,6 @@ pub fn surface_binding_energy(particle_1: &mut particle::Particle, particle_1.dir.x = -2.*(costheta)*dx/mag + cosx; particle_1.dir.y = -2.*(costheta)*dy/mag + cosy; particle_1.dir.z = -2.*(costheta)*dz/mag + cosz; - particle_1.backreflected = true; } } diff --git a/src/parry.rs b/src/parry.rs index 75d9a24..8686683 100644 --- a/src/parry.rs +++ b/src/parry.rs @@ -63,6 +63,7 @@ impl Geometry for ParryBall { let length_unit: f64 = match input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., @@ -186,6 +187,7 @@ impl Geometry for ParryTriMesh { let length_unit: f64 = match input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., diff --git a/src/particle.rs b/src/particle.rs index 6aff9f6..e82412b 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -83,7 +83,7 @@ pub struct Particle { pub left: bool, pub incident: bool, pub first_step: bool, - pub trajectory: Vec, + pub trajectory: Vec, pub energies: Vec, pub track_trajectories: bool, pub number_collision_events: usize, @@ -122,7 +122,7 @@ impl Particle { left: false, incident: true, first_step: true, - trajectory: vec![Vector4::new(input.E, input.x, input.y, input.z)], + trajectory: vec![TrajectoryElement::new(input.E, input.x, input.y, input.z)], energies: vec![], track_trajectories: options.track_trajectories, number_collision_events: 0, @@ -170,7 +170,7 @@ impl Particle { /// If `track_trajectories`, add the current (E, x, y, z) to the trajectory. pub fn add_trajectory(&mut self) { if self.track_trajectories { - self.trajectory.push(Vector4 {E: self.E, x: self.pos.x, y: self.pos.y, z: self.pos.z}); + self.trajectory.push(TrajectoryElement {E: self.E, x: self.pos.x, y: self.pos.y, z: self.pos.z}); } } @@ -190,58 +190,57 @@ impl Particle { self.m*speed*self.dir.z, ) } -} -/// Rotate a particle by a deflection psi at an azimuthal angle phi. -pub fn rotate_particle(particle_1: &mut particle::Particle, psi: f64, phi: f64) { - let cosx: f64 = particle_1.dir.x; - let cosy: f64 = particle_1.dir.y; - let cosz: f64 = particle_1.dir.z; - let cphi: f64 = phi.cos(); - let sphi: f64 = phi.sin(); - let sa = (1. - cosx*cosx).sqrt(); + /// Rotate a particle by a deflection psi at an azimuthal angle phi. + pub fn rotate(&mut self, psi: f64, phi: f64) { + 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(); - //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); + //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 dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new}; + let dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new}; - particle_1.dir.assign(&dir_new); - particle_1.dir.normalize(); -} + self.dir.assign(&dir_new); + self.dir.normalize(); + } -/// Push particle in space according to previous direction and return the distance traveled. -pub fn particle_advance(particle_1: &mut particle::Particle, mfp: f64, asymptotic_deflection: f64) -> f64 { + /// Push particle in space according to previous direction and return the distance traveled. + pub fn advance(&mut self, mfp: f64, asymptotic_deflection: f64) -> f64 { - if particle_1.E > particle_1.Ec { - particle_1.add_trajectory(); - } + if self.E > self.Ec { + self.add_trajectory(); + } - //Update previous position - particle_1.pos_old.x = particle_1.pos.x; - particle_1.pos_old.y = particle_1.pos.y; - particle_1.pos_old.z = particle_1.pos.z; + //Update previous position + self.pos_old.x = self.pos.x; + self.pos_old.y = self.pos.y; + self.pos_old.z = self.pos.z; - //In order to keep average denisty constant, must add back previous asymptotic deflection - let distance_traveled = mfp + particle_1.asymptotic_deflection - asymptotic_deflection; - //let distance_traveled = mfp - asymptotic_deflection; + //In order to keep average denisty constant, must add back previous asymptotic deflection + let distance_traveled = mfp + self.asymptotic_deflection - asymptotic_deflection; - //dir has been updated, so use previous direction to advance in space - particle_1.pos.x += particle_1.dir_old.x*distance_traveled; - particle_1.pos.y += particle_1.dir_old.y*distance_traveled; - particle_1.pos.z += particle_1.dir_old.z*distance_traveled; - particle_1.asymptotic_deflection = asymptotic_deflection; + //dir has been updated, so use previous direction to advance in space + self.pos.x += self.dir_old.x*distance_traveled; + self.pos.y += self.dir_old.y*distance_traveled; + self.pos.z += self.dir_old.z*distance_traveled; + self.asymptotic_deflection = asymptotic_deflection; - //Update previous direction - particle_1.dir_old.x = particle_1.dir.x; - particle_1.dir_old.y = particle_1.dir.y; - particle_1.dir_old.z = particle_1.dir.z; + //Update previous direction + self.dir_old.x = self.dir.x; + self.dir_old.y = self.dir.y; + self.dir_old.z = self.dir.z; - return distance_traveled; + return distance_traveled; + } } pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) { @@ -249,10 +248,6 @@ pub fn surface_refraction(particle: &mut Particle, normal: Vector, Es: f64) { let costheta = particle.dir.dot(&normal); - let a = (E/(E + Es)).sqrt(); - let b = -(E).sqrt()*costheta; - let c = (E*costheta.powi(2) + Es).sqrt(); - let u1x = (E/(E + Es)).sqrt()*particle.dir.x + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.x; let u1y = (E/(E + Es)).sqrt()*particle.dir.y + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.y; let u1z = (E/(E + Es)).sqrt()*particle.dir.z + ((-(E).sqrt()*costheta + (E*costheta.powi(2) + Es).sqrt())/(E + Es).sqrt())*normal.z; diff --git a/src/sphere.rs b/src/sphere.rs index c4994f3..92553b6 100644 --- a/src/sphere.rs +++ b/src/sphere.rs @@ -58,6 +58,7 @@ impl Geometry for Sphere { let length_unit: f64 = match input.length_unit.as_str() { "MICRON" => MICRON, "CM" => CM, + "MM" => MM, "ANGSTROM" => ANGSTROM, "NM" => NM, "M" => 1., diff --git a/src/structs.rs b/src/structs.rs index 2b9eb90..47f20a0 100644 --- a/src/structs.rs +++ b/src/structs.rs @@ -44,18 +44,18 @@ impl Vector { } } -/// Vector4 is a trajectory-tracking object that includes x, y, z, and the current energy. +/// TrajectoryElement is a trajectory-tracking object that includes x, y, z, and the current energy. #[derive(Clone)] -pub struct Vector4 { +pub struct TrajectoryElement { pub E: f64, pub x: f64, pub y: f64, pub z: f64, } -impl Vector4 { - pub fn new(E: f64, x: f64, y: f64, z: f64) -> Vector4 { - Vector4 { +impl TrajectoryElement { + pub fn new(E: f64, x: f64, y: f64, z: f64) -> TrajectoryElement { + TrajectoryElement { E, x, y, diff --git a/src/tests.rs b/src/tests.rs index ea34548..a1f4a8c 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -660,9 +660,6 @@ fn test_surface_refraction() { let cosy_new = cosy_new/dir_mag; let cosz_new = cosz_new/dir_mag; - //let delta_theta = particle::refraction_angle(cosx, E, E + Es); - //particle::rotate_particle(&mut particle_1, delta_theta, 0.); - let normal = Vector::new(1.0, 0.0, 0.0); particle::surface_refraction(&mut particle_1, normal, Es); @@ -690,9 +687,6 @@ fn test_surface_refraction() { println!("{} {} {}", particle_1.dir.x, particle_1.dir.y, particle_1.dir.z); } - //let delta_theta = particle::refraction_angle(particle_1.dir.x, particle_1.E, particle_1.E - Es); - //particle::rotate_particle(&mut particle_1, delta_theta, 0.); - let normal = Vector::new(1.0, 0.0, 0.0); particle::surface_refraction(&mut particle_1, normal, -Es); @@ -856,10 +850,10 @@ fn test_momentum_conservation() { 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::rotate_particle(&mut particle_1, binary_collision_result.psi, + particle_1.rotate(binary_collision_result.psi, binary_collision_geometries[0].phi_azimuthal); - particle::rotate_particle(&mut particle_2, -binary_collision_result.psi_recoil, + particle_2.rotate(-binary_collision_result.psi_recoil, binary_collision_geometries[0].phi_azimuthal); //Subtract total energy from all simultaneous collisions and electronic stopping @@ -894,7 +888,7 @@ fn test_momentum_conservation() { } #[test] -fn test_rotate_particle() { +fn test_rotate() { let mass = 1.; let Z = 1.; let E = 1.; @@ -912,18 +906,18 @@ fn test_rotate_particle() { let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); //Check that rotation in 2D works - particle::rotate_particle(&mut particle, psi, phi); + particle.rotate(psi, phi); assert!(approx_eq!(f64, particle.dir.x, 0., epsilon = 1E-12), "particle.dir.x: {} Should be ~0.", particle.dir.x); assert!(approx_eq!(f64, particle.dir.y, 1., epsilon = 1E-12), "particle.dir.y: {} Should be ~1.", particle.dir.y); //Check that rotating back by negative psi returns to the previous values - particle::rotate_particle(&mut particle, -psi, phi); + particle.rotate(-psi, phi); assert!(approx_eq!(f64, particle.dir.x, cosx, epsilon = 1E-12), "particle.dir.x: {} Should be ~{}", particle.dir.x, cosx); assert!(approx_eq!(f64, particle.dir.y, cosy, epsilon = 1E-12), "particle.dir.y: {} Should be ~{}", particle.dir.y, cosy); //Check that azimuthal rotation by 180 degrees works correctly let phi = PI; - particle::rotate_particle(&mut particle, psi, phi); + particle.rotate(psi, phi); assert!(approx_eq!(f64, particle.dir.x, 1., epsilon = 1E-12), "particle.dir.x: {} Should be ~1.", particle.dir.x); assert!(approx_eq!(f64, particle.dir.y, 0., epsilon = 1E-12), "particle.dir.y: {} Should be ~0.", particle.dir.y); @@ -950,7 +944,7 @@ fn test_particle_advance() { let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); - let distance_traveled = particle::particle_advance(&mut particle, mfp, asymptotic_deflection); + let distance_traveled = particle.advance(mfp, asymptotic_deflection); assert_eq!(particle.pos.x, (1. - 0.5)*cosx); assert_eq!(particle.pos.y, (1. - 0.5)*cosy);