diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 50c2ba7..12713e0 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -27,7 +27,7 @@ jobs: sudo apt-get install rustc - name: Install pip for Python-3 run: | - sudo apt-get install python3-pip + sudo apt-get install python3-pip python3-dev - name: Install Python libraries run: | python3 -m pip install numpy shapely scipy @@ -35,15 +35,25 @@ jobs: run: | git clone https://github.com/uiri/toml.git cd toml - sudo python3 setup.py install + python3 setup.py install --root . - name: Install HDF5 Libraries run: | sudo apt install libhdf5-dev - name: test Python Bindings run: | - sudo python3 -m pip install setuptools_rust testresources - sudo python3 setup.py install - python3 -c "from libRustBCA.pybca import *; print(simple_bca_py)" + python3 -m pip install setuptools_rust testresources + python3 -m pip install . + python3 -c "from libRustBCA.pybca import *;" + - name: Test Fortran and C bindings + run : | + cargo build --release + cp examples/test_rustbca.f90 . + gfortran -c rustbca.f90 target/release/liblibRustBCA.so + gfortran test_rustbca.f90 rustbca.f90 target/release/liblibRustBCA.so + ./a.out + cp examples/RustBCA.c . + g++ RustBCA.c RustBCA.h target/release/liblibRustBCA.so -Iexamples/ -I/usr/include/python3.8 + ./a.out - name: Test RustBCA run: | sudo cargo test --features cpr_rootfinder_netlib,hdf5_input,distributions,parry3d diff --git a/Cargo.toml b/Cargo.toml index 92e6997..a908e10 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,7 @@ [package] name = "RustBCA" -version = "1.1.1" +version = "1.2.0" + authors = ["Jon Drobny "] edition = "2018" @@ -14,6 +15,7 @@ crate-type = ["cdylib"] [dependencies] rand = "0.8.3" +rand_distr = "0.4.2" toml = "0.5.8" anyhow = "1.0.38" itertools = "0.10.0" diff --git a/RustBCA.h b/RustBCA.h index f8cd52d..e89188b 100644 --- a/RustBCA.h +++ b/RustBCA.h @@ -97,8 +97,38 @@ struct InputCompoundBCA { double *Eb2; }; +struct OutputTaggedBCA { + uintptr_t len; + double (*particles)[9]; + double *weights; + int32_t *tags; +}; + +struct InputTaggedBCA { + uintptr_t len; + /// x y z + double (*positions)[3]; + /// vx, vy, vz + double (*velocities)[3]; + double Z1; + double m1; + double Ec1; + double Es1; + uintptr_t num_species_target; + double *Z2; + double *m2; + double *n2; + double *Ec2; + double *Es2; + double *Eb2; + int32_t *tags; + double *weights; +}; + extern "C" { +OutputTaggedBCA compound_tagged_bca_list_c(InputTaggedBCA input); + OutputBCA simple_bca_list_c(InputSimpleBCA input); OutputBCA compound_bca_list_c(InputCompoundBCA input); diff --git a/examples/RustBCA.c b/examples/RustBCA.c new file mode 100644 index 0000000..a218531 --- /dev/null +++ b/examples/RustBCA.c @@ -0,0 +1,54 @@ +#include "RustBCA.h" +#include +#include + +int main(int argc, char * argv[]) { + OutputTaggedBCA output; + double velocities[2][3] = {{500000.0, 0.1, 0.0}, {500000.0, 0.1, 0.0}}; + double positions[2][3] = {{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}}; + int tags[2] = {0, 1}; + double weights[2] = {1.0, 1.0}; + double Z[3] = {74.0, 74.0}; + double m[2] = {184.0, 184.0}; + double n[2] = {0.06306, 0.06306}; + double Ec[2] = {1.0, 1.0}; + double Es[2] = {8.79, 8.79}; + double Eb[2] = {0.0, 0.0}; + + InputTaggedBCA input = { + 2, + positions, + velocities, + 1.0, + 1.0, + 1.0, + 1.0, + 2, + Z, + m, + n, + Ec, + Es, + Eb, + tags, + weights + }; + + //output = simple_bca_c(0., 0., 0., 0.5, 0.5, 0.00, 2000.0, 2.0, 4.0, 1.0, 0.0, 74.0, 184.0, 1.0, 8.79, 0.06306, 0.0); + //output = compound_bca_list_c(input); + output = compound_tagged_bca_list_c(input); + + std::cout << "Particle 1 Z: "; + std::cout << output.particles[0][0]; + std::cout << std::endl; + std::cout << "Particle 1 E [eV]: "; + std::cout << output.particles[0][2]; + std::cout << std::endl; + std::cout << "Particle 2 Z: "; + std::cout << output.particles[1][0]; + std::cout << std::endl; + std::cout << "Particle 2 E [eV]: "; + std::cout << output.particles[1][2]; + std::cout << std::endl; + return 0; +} diff --git a/examples/test_rustbca.f90 b/examples/test_rustbca.f90 new file mode 100644 index 0000000..c0b8cae --- /dev/null +++ b/examples/test_rustbca.f90 @@ -0,0 +1,75 @@ + +program test_rustbca + + use rustbca + use, intrinsic :: iso_c_binding + + integer :: N_ions + real(c_double), allocatable, dimension(:) :: ux, uy, uz, E, Z1, m1, Ec1, Es1 + integer(c_int) :: num_species_target, num_emitted_particles + real(c_double), target :: Z2(2), m2(2), Ec2(2), Es2(2), Eb2(2), n2(2) + real(c_double) :: ux1, uy1, uz1, E1 + type(c_ptr) :: bca_output_c + real(c_double), pointer, dimension(:,:) :: bca_output_f + real :: start, stop + logical(c_bool) :: track_recoils + + !Initial ion conditions + N_ions = 100000 + allocate(ux(N_ions), uy(N_ions), uz(N_ions), E(N_ions), Z1(N_ions), m1(N_ions), Ec1(N_ions), Es1(N_ions)) + ux(:) = 0.999 + uy(:) = sqrt(1.0 - 0.999*0.999) + uz(:) = 0.0 + E(:) = 1000.0_8 + + !Hydrogen + Z1(:) = 1.0_8 + m1(:) = 1.008_8 + Ec1(:) = 1.0_8 + Es1(:) = 1.5_8 + + !Titanium Hydride + num_species_target = 2 + Z2(1) = 22.0_8 + m2(1) = 47.867_8 + Ec2(1) = 4.84_8 + Es2(1) = 4.84_8 + Eb2(1) = 3.0_8 + n2(1) = 0.04527_8 + + Z2(2) = 1.0_8 + m2(2) = 1.008_8 + Ec2(2) = 1.5_8 + Es2(2) = 1.5_8 + Eb2(2) = 0.0_8 + n2(2) = 0.09054_8 + + track_recoils = .false. + + call cpu_time(start) + bca_output_c = compound_bca_list_fortran(N_ions, track_recoils, ux, uy, uz, E, & + Z1, m1, Ec1, Es1, & + num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, & + num_emitted_particles) + call c_f_pointer(bca_output_c, bca_output_f, [num_emitted_particles, 6]) + call cpu_time(stop) + + write(*,*) "Elapsed time in seconds per ion per eV: ", (stop - start)/N_ions/1000.0 + + !write(*,*) bca_output_f + + call cpu_time(start) + do i = 0, N_ions + !Test reflect_single_ion routine + ux1 = 0.999 + uy1 = sqrt(1.0 - 0.999*0.999) + uz1 = 0.0 + E1 = 1000.0 + call reflect_single_ion_c(num_species_target, ux1, uy1, uz1, E1, Z1(1), m1(1), Ec1(1), Es1(1), Z2, m2, Ec2, Es2, Eb2, n2) + end do + call cpu_time(stop) + write(*,*) "Elapsed time in ions per eV per s: ", (stop - start)/N_ions/1000.0 + + !call exit(1) + +end program test_rustbca diff --git a/rustbca.f90 b/rustbca.f90 new file mode 100644 index 0000000..681a78c --- /dev/null +++ b/rustbca.f90 @@ -0,0 +1,95 @@ +module rustbca + + use, intrinsic :: iso_c_binding + + interface + + subroutine reflect_single_ion_c(num_species_target, ux, uy, uz, E1, & + Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, Eb2, n2) bind(c) + + !Runs a single ion BCA trajectory with no recoils + !Args: + ! num_species_target (integer): number of species in target + ! ux (real(c_double)): x-direction of incident ion, x-direction of reflected ion + ! uy (real(c_double)): y-direction of incident ion, y-direction of reflected ion + ! uz (real(c_double)): z-direction of incident ion, z-direction of reflected ion + ! E1 (real(c_double)): initial energy of incident ion in eV + ! Z1 (real(c_double)): atomic number of incident ion + ! m1 (real(c_double)): atomic mass of incident ion in eV + ! Ec1 (real(c_double)): cutoff energy of incident ion in eV + ! Es1 (real(c_double)): surface binding energy of incident ion in eV + ! Z2 (real(c_double), dimension(:)): list of atomic numbers of target speciesd + ! m2 (real(c_double), dimension(:)): list of atomic masses of target species in amu + ! Ec2 (real(c_double), dimension(:)): list of cutoff energies of target species in eV + ! Es2 (real(c_double), dimension(:)): list of surface binding energies of target species in eV + ! Eb2 (real(c_double), dimension(:)): list of bulk binding energies of target species in eV + ! n2 (real(c_double), dimension(:)): list of number densities of target species in 1/angstrom^3 + + use, intrinsic :: iso_c_binding + real(c_double), intent(inout) :: ux, uy, uz, E1 + real(c_double), intent(in) :: Z1, m1, Ec1, Es1 + integer(c_int), intent(in) :: num_species_target + real(c_double), intent(in), dimension(*) :: Z2, m2, Ec2, Es2, Eb2, n2 + + end subroutine reflect_single_ion_c + + function compound_bca_list_fortran(num_incident_ions, track_recoils, ux, uy, uz, E1, & + Z1, m1, Ec1, Es1, & + num_species_target, Z2, m2, Ec2, Es2, Eb2, n2, & + num_emitted_particles) bind(c) result(output) + + + !Runs a homogeneous, flat, compound target BCA with an arbitrary list of ions. + !Args: + ! num_incident_ion (integer(c_int)): number of incident ions + ! track_recoils (logical(c_bool)): whether to generate recoils (disable to turn off sputtering) + ! ux (real(c_double), dimension(:)): x-direction of incident ion, x-direction of reflected ion + ! uy (real(c_double), dimension(:)): y-direction of incident ion, y-direction of reflected ion + ! uz (real(c_double), dimension(:)): z-direction of incident ion, z-direction of reflected ion + ! E1 (real(c_double), dimension(:)): initial energy of incident ion in eV + ! Z1 (real(c_double), dimension(:)): atomic number of incident ion + ! m1 (real(c_double), dimension(:)): atomic mass of incident ion in eV + ! Ec1 (real(c_double), dimension(:)): cutoff energy of incident ion in eV + ! num_species_target(integer(c_int)): number of species in target + ! Es1 (real(c_double), dimension(:)): surface binding energy of incident ion in eV + ! Z2 (real(c_double), dimension(:)): list of atomic numbers of target speciesd + ! m2 (real(c_double), dimension(:)): list of atomic masses of target species in amu + ! Ec2 (real(c_double), dimension(:)): list of cutoff energies of target species in eV + ! Es2 (real(c_double), dimension(:)): list of surface binding energies of target species in eV + ! Eb2 (real(c_double), dimension(:)): list of bulk binding energies of target species in eV + ! n2 (real(c_double), dimension(:)): list of number densities of target species in 1/angstrom^3 + ! num_emitted_particles (integer(c_int), intent(out)): NOTE THAT THIS IS INTENT(OUT) number of emitted particles in output + !Returns: + ! output (type(c_ptr)): a c pointer to a 2D array of size (num_emitted_particles, 6) that consists of Z, m, E, ux, uy, uz + + use, intrinsic :: iso_c_binding + logical(c_bool), intent(in) :: track_recoils + integer(c_int), intent(in) :: num_incident_ions, num_species_target + integer(c_int), intent(out) :: num_emitted_particles + real(c_double), intent(in), dimension(*) :: ux, uy, uz, E1, Z1, m1, Ec1, Es1 + real(c_double), intent(in), dimension(*) :: Z2, m2, Ec2, Es2, EB2, n2 + type(c_ptr) :: output + end function compound_bca_list_fortran + + end interface + + contains + + subroutine transform_to_local_angle(ux, uy, uz, alpha) + + !Rotates a vector in 2D + !Args: + ! ux (real(c_double)): x-direction + ! uy (real(c_double)): y-direction + ! uz (real(c_double)): z-direction + ! alpha (real(c_double)): local surface angle measured counter-clockwise from x-axis in radians + + real(8), intent(inout) :: ux, uy, uz + real(8), intent(in) :: alpha + + ux = ux*cos(alpha) - uy*sin(alpha) + uy = ux*sin(alpha) + uy*cos(alpha) + + end subroutine transform_to_local_angle + +end module rustbca \ No newline at end of file diff --git a/setup.py b/setup.py index e2ab17b..c9db496 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,18 @@ setup( name="RustBCA", - rust_extensions=[RustExtension("libRustBCA.pybca", binding=Binding.PyO3, features=["python"])], + version="1.1.2", + rust_extensions=[ + RustExtension( + "libRustBCA.pybca", + binding=Binding.PyO3, + features=["python"], + #args=["+nightly", "--edition 2018", "-Z unstable-options"], + #optional=True, + rust_version="1.56.1" + + ) + ], # rust extensions are not zip safe, just like C-extensions. zip_safe=False, ) diff --git a/src/bca.rs b/src/bca.rs index 7ceaa15..4f911bd 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -377,15 +377,19 @@ pub fn choose_collision_partner(particle_1: &particle::Particle, ma //Choose recoil Z, M let (species_index, Z_recoil, M_recoil, Ec_recoil, Es_recoil, interaction_index) = material.choose(x_recoil, y_recoil, z_recoil); - - return (species_index, - particle::Particle::new( - M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil, - x_recoil, y_recoil, z_recoil, - cosx, cosy, cosz, - false, options.track_recoil_trajectories, interaction_index - ) - ) + let mut new_particle = particle::Particle::new( + M_recoil, Z_recoil, 0., Ec_recoil, Es_recoil, + x_recoil, y_recoil, z_recoil, + cosx, cosy, cosz, + false, options.track_recoil_trajectories, interaction_index + ); + + // Carry through tag, weight, and tracked vector from originating particle + new_particle.weight = particle_1.weight; + new_particle.tag = particle_1.tag; + new_particle.tracked_vector = particle_1.tracked_vector; + + return (species_index, new_particle) } /// Calculate the distance of closest approach of two particles given a particular binary collision geometry. @@ -443,6 +447,7 @@ pub fn update_particle_energy(particle_1: &mut particle::Particle, recoil_energy: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) { //If particle energy drops below zero before electronic stopping calcualtion, it produces NaNs + assert!(!recoil_energy.is_nan(), "Numerical error: recoil Energy is NaN. E0 = {} Za = {} Ma = {} x0 = {} Zb = {} delta-x = {}", particle_1.E, particle_1.Z, particle_1.m, x0, strong_collision_Z, distance_traveled); particle_1.E -= recoil_energy; assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following collision."); if particle_1.E < 0. { @@ -694,8 +699,8 @@ pub fn newton_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par return Ok(x0); } } - return Err(anyhow!("Numerical error: exceeded maximum number of Newton-Raphson iterations, {}. E: {}; x0: {}; Error: {}; Tolerance: {}", - max_iterations, E0, x0, err, tolerance)); + return Err(anyhow!("Numerical error: exceeded maximum number of Newton-Raphson iterations, {}. E: {} eV; x0: {}; Error: {}; Tolerance: {}; Za: {}; Zb: {}; Ma: {} amu; Mb: {} amu; a: {}; p: {} A", + max_iterations, E0/Q, x0, err, tolerance, Za, Zb, Ma/AMU, Mb/AMU, a, impact_parameter/ANGSTROM)); } /// Gauss-Mehler quadrature. diff --git a/src/enums.rs b/src/enums.rs index 07a7cbd..272e923 100644 --- a/src/enums.rs +++ b/src/enums.rs @@ -11,6 +11,14 @@ pub enum MaterialType { TRIMESH(material::Material) } +#[derive(Deserialize, PartialEq, Clone, Copy)] +#[serde(untagged)] +pub enum Distributions { + UNIFORM{min: f64, max: f64}, + NORMAL{mean: f64, std: f64}, + POINT(f64), +} + #[derive(Deserialize)] pub enum GeometryType { MESH0D, diff --git a/src/input.rs b/src/input.rs index 42a6779..a07dfc4 100644 --- a/src/input.rs +++ b/src/input.rs @@ -1,4 +1,5 @@ use super::*; +use rand_distr::{Normal, Distribution, Uniform}; pub trait InputFile: GeometryInput { fn new(string: &str) -> Self; @@ -332,26 +333,56 @@ where ::InputFileFormat: Deserialize<'static> + 'static { let Ec = particle_parameters.Ec[particle_index]; let Es = particle_parameters.Es[particle_index]; let interaction_index = particle_parameters.interaction_index[particle_index]; + let (x, y, z) = particle_parameters.pos[particle_index]; let (cosx, cosy, cosz) = particle_parameters.dir[particle_index]; - assert!(cosx < 1., - "Input error: particle x-direction cannot be exactly equal to 1 to avoid numerical gimbal lock."); + for sub_particle_index in 0..N_ { //Add new particle to particle vector particle_input.push( particle::ParticleInput{ m: m*mass_unit, Z: Z, - E: E*energy_unit, + E: match E { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*energy_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*energy_unit}, + Distributions::POINT(x) => x*energy_unit, + }, Ec: Ec*energy_unit, Es: Es*energy_unit, - x: x*length_unit, - y: y*length_unit, - z: z*length_unit, - ux: cosx, - uy: cosy, - uz: cosz, - interaction_index: interaction_index + x: match x { + 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) => x*length_unit, + }, + y: match y { + 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) => x*length_unit, + }, + z: match z { + 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) => x*length_unit, + }, + 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} + }, + uy: match cosy { + 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) => x*length_unit, + }, + uz: match cosz { + 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) => x*length_unit, + }, + interaction_index: interaction_index, + tag: 0, + weight: 1.0, } ); } @@ -377,8 +408,7 @@ where ::InputFileFormat: Deserialize<'static> + 'static { let interaction_index = particle_parameters.interaction_index[particle_index]; let (x, y, z) = particle_parameters.pos[particle_index]; let (cosx, cosy, cosz) = particle_parameters.dir[particle_index]; - assert!(cosx < 1., - "Input error: particle x-direction cannot be exactly equal to 1 to avoid numerical gimbal lock."); + for sub_particle_index in 0..N_ { //Add new particle to particle vector @@ -386,16 +416,46 @@ where ::InputFileFormat: Deserialize<'static> + 'static { particle::ParticleInput{ m: m*mass_unit, Z: Z, - E: E*energy_unit, + E: match E { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*energy_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*energy_unit}, + Distributions::POINT(x) => x*energy_unit, + }, Ec: Ec*energy_unit, Es: Es*energy_unit, - x: x*length_unit, - y: y*length_unit, - z: z*length_unit, - ux: cosx, - uy: cosy, - uz: cosz, - interaction_index + x: match x { + 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) => x*length_unit, + }, + y: match y { + 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) => x*length_unit, + }, + z: match z { + 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) => x*length_unit, + }, + 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}, + }, + uy: match cosy { + 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) => x*length_unit, + }, + uz: match cosz { + 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) => x*length_unit, + }, + interaction_index: interaction_index, + tag: 0, + weight: 1.0, } ); } diff --git a/src/lib.rs b/src/lib.rs index bb3ba11..98ebe6e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -13,8 +13,8 @@ use std::{fmt}; use std::mem::discriminant; //Error handling crate -use anyhow::Result; -use anyhow::*; +use anyhow::{Result, Context, anyhow}; +//use anyhow::*; //Serializing/Deserializing crate use serde::*; @@ -28,6 +28,9 @@ use std::fs::OpenOptions; use std::io::prelude::*; use std::io::BufWriter; +//C integer +use std::os::raw::c_int; + //standard slice use std::slice; @@ -123,6 +126,325 @@ pub struct InputCompoundBCA { pub Eb2: *mut f64, } +#[derive(Debug)] +#[repr(C)] +pub struct InputTaggedBCA { + pub len: usize, + /// x y z + pub positions: *mut [f64; 3], + /// vx, vy, vz + pub velocities: *mut [f64; 3], + pub Z1: f64, + pub m1: f64, + pub Ec1: f64, + pub Es1: f64, + pub num_species_target: usize, + pub Z2: *mut f64, + pub m2: *mut f64, + pub n2: *mut f64, + pub Ec2: *mut f64, + pub Es2: *mut f64, + pub Eb2: *mut f64, + pub tags: *mut i32, + pub weights: *mut f64, +} + +#[derive(Debug)] +#[repr(C)] +pub struct OutputTaggedBCA { + pub len: usize, + pub particles: *mut [f64; 9], + pub weights: *mut f64, + pub tags: *mut i32 +} + +#[no_mangle] +pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTaggedBCA { + + let mut total_output = vec![]; + let mut output_tags = vec![]; + let mut output_weights = vec![]; + + #[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: false, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + 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: false, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + 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 Z2 = unsafe { slice::from_raw_parts(input.Z2, input.num_species_target).to_vec() }; + let m2 = unsafe { slice::from_raw_parts(input.m2, input.num_species_target).to_vec() }; + let n2 = unsafe { slice::from_raw_parts(input.n2, input.num_species_target).to_vec() }; + let Ec2 = unsafe { slice::from_raw_parts(input.Ec2, input.num_species_target).to_vec() }; + let Es2 = unsafe { slice::from_raw_parts(input.Es2, input.num_species_target).to_vec() }; + let Eb2 = unsafe { slice::from_raw_parts(input.Eb2, input.num_species_target).to_vec() }; + let positions = unsafe { slice::from_raw_parts(input.positions, input.num_species_target).to_vec() }; + let tags = unsafe { slice::from_raw_parts(input.tags, input.num_species_target).to_vec() }; + let weights = unsafe { slice::from_raw_parts(input.weights, input.num_species_target).to_vec() }; + + 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; input.num_species_target], + surface_binding_model: SurfaceBindingModel::AVERAGE, + 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 velocities = unsafe { slice::from_raw_parts(input.velocities, input.len) }; + + let mut index: usize = 0; + for velocity in velocities { + + let vx = velocity[0]; + let vy = velocity[1]; + let vz = velocity[2]; + + let v = (vx*vx + vy*vy + vz*vz).sqrt(); + + let E1 = 0.5*input.m1*AMU*v*v; + + let ux = vx/v; + let uy = vy/v; + let uz = vz/v; + + let p = particle::Particle { + m: input.m1*AMU, + Z: input.Z1, + E: E1, + Ec: input.Ec1*EV, + Es: input.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, + 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: weights[index], + tag: tags[index], + tracked_vector: Vector::new(positions[index][0], positions[index][1], positions[index][2]), + }; + + let output = bca::single_ion_bca(p, &m, &options); + + for particle in output { + + if (particle.left) | (particle.incident) { + total_output.push( + [ + particle.Z, + particle.m/AMU, + particle.E/EV, + + particle.tracked_vector.x/ANGSTROM, + particle.tracked_vector.y/ANGSTROM, + particle.tracked_vector.z/ANGSTROM, + + particle.dir.x, + particle.dir.y, + particle.dir.z + ] + ); + output_tags.push(particle.tag); + output_weights.push(particle.weight); + } + } + index += 1; + } + + let len = total_output.len(); + let particles = total_output.as_mut_ptr(); + let tags_ptr = output_tags.as_mut_ptr(); + let weights_ptr = output_weights.as_mut_ptr(); + + std::mem::forget(total_output); + OutputTaggedBCA { + len, + particles, + tags: tags_ptr, + weights: weights_ptr + } +} + + +#[no_mangle] +#[cfg(not(feature = "distributions"))] +pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut f64, uy: &mut f64, uz: &mut f64, E1: &mut f64, Z1: &mut f64, m1: &mut f64, Ec1: &mut f64, Es1: &mut f64, Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64) { + + assert!(E1 > &mut 0.0); + + let options = Options { + name: "test".to_string(), + track_trajectories: false, + track_recoil_trajectories: false, + track_recoils: false, + write_buffer_size: 8000, + weak_collision_order: 3, + suppress_deep_recoils: false, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + 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 Z2 = unsafe { slice::from_raw_parts(Z2, *num_species_target as usize).to_vec() }; + let m2 = unsafe { slice::from_raw_parts(m2, *num_species_target as usize).to_vec() }; + let n2 = unsafe { slice::from_raw_parts(n2, *num_species_target as usize).to_vec() }; + let Ec2 = unsafe { slice::from_raw_parts(Ec2, *num_species_target as usize).to_vec() }; + let Es2 = unsafe { slice::from_raw_parts(Es2, *num_species_target as usize).to_vec() }; + let Eb2 = unsafe { slice::from_raw_parts(Eb2, *num_species_target as usize).to_vec() }; + + 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 as usize], + surface_binding_model: SurfaceBindingModel::AVERAGE, + 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 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: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0), + }; + + let output = bca::single_ion_bca(p, &m, &options); + + *ux = output[0].dir.x; + *uy = output[0].dir.y; + *uz = output[0].dir.z; + if output[0].pos.x >= 0.0 { + *E1 = 0.0 + } else { + *E1 = output[0].E/EV; + } +} + #[no_mangle] pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { @@ -250,7 +572,10 @@ pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { track_trajectories: false, number_collision_events: 0, backreflected: false, - interaction_index : 0 + interaction_index : 0, + weight: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0), }; @@ -420,7 +745,10 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { track_trajectories: false, number_collision_events: 0, backreflected: false, - interaction_index : 0 + interaction_index : 0, + weight: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0), }; @@ -456,6 +784,182 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { } } + +#[no_mangle] +pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track_recoils: &mut bool, + ux: *mut f64, uy: *mut f64, uz: *mut f64, E1: *mut f64, + Z1: *mut f64, m1: *mut f64, Ec1: *mut f64, Es1: *mut f64, + num_species_target: &mut c_int, + Z2: *mut f64, m2: *mut f64, Ec2: *mut f64, Es2: *mut f64, Eb2: *mut f64, n2: *mut f64, + num_emitted_particles: &mut c_int + ) -> *const [f64; 6] { + + //println!("{} {}", num_incident_ions, num_species_target); + + let mut total_output = vec![]; + + #[cfg(feature = "distributions")] + let options = Options { + name: "test".to_string(), + track_trajectories: false, + track_recoils: *track_recoils, + track_recoil_trajectories: false, + write_buffer_size: 8000, + weak_collision_order: 3, + suppress_deep_recoils: false, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + 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: *track_recoils, + track_recoil_trajectories: false, + write_buffer_size: 8000, + weak_collision_order: 3, + suppress_deep_recoils: false, + high_energy_free_flight_paths: false, + electronic_stopping_mode: ElectronicStoppingMode::INTERPOLATED, + 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 ux = unsafe { slice::from_raw_parts(ux, *num_incident_ions as usize).to_vec() }; + let uy = unsafe { slice::from_raw_parts(uy, *num_incident_ions as usize).to_vec() }; + let uz = unsafe { slice::from_raw_parts(uz, *num_incident_ions as usize).to_vec() }; + let Z1 = unsafe { slice::from_raw_parts(Z1, *num_incident_ions as usize).to_vec() }; + let m1 = unsafe { slice::from_raw_parts(m1, *num_incident_ions as usize).to_vec() }; + let E1 = unsafe { slice::from_raw_parts(E1, *num_incident_ions as usize).to_vec() }; + let Ec1 = unsafe { slice::from_raw_parts(Ec1, *num_incident_ions as usize).to_vec() }; + let Es1 = unsafe { slice::from_raw_parts(Es1, *num_incident_ions as usize).to_vec() }; + + //println!("ux: {} uy: {} uz: {} Z1: {} m1: {} E1: {} Ec1: {} Es1: {}", ux[0], uy[0], uz[0], Z1[0], m1[0], E1[0], Ec1[0], Es1[0]); + + let Z2 = unsafe { slice::from_raw_parts(Z2, *num_species_target as usize).to_vec() }; + let m2 = unsafe { slice::from_raw_parts(m2, *num_species_target as usize).to_vec() }; + let n2 = unsafe { slice::from_raw_parts(n2, *num_species_target as usize).to_vec() }; + let Ec2 = unsafe { slice::from_raw_parts(Ec2, *num_species_target as usize).to_vec() }; + let Es2 = unsafe { slice::from_raw_parts(Es2, *num_species_target as usize).to_vec() }; + let Eb2 = unsafe { slice::from_raw_parts(Eb2, *num_species_target as usize).to_vec() }; + + //println!("Z2: {} m2: {} n2: {} Ec2: {} Es2: {} Eb2: {}", Z2[0], m2[0], n2[0], Ec2[0], Es2[0], Eb2[0]); + + 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 as usize], + 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); + + for (((((((E1_, ux_), uy_), uz_), Z1_), Ec1_), Es1_), m1_) in E1.iter().zip(ux).zip(uy).zip(uz).zip(Z1).zip(Ec1).zip(Es1).zip(m1) { + + 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: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0) + }; + + let output = bca::single_ion_bca(p, &m, &options); + + for particle in output { + + if (particle.left) | (particle.incident) { + total_output.push( + [ + particle.Z, + particle.m/AMU, + particle.E/EV, + particle.dir.x, + particle.dir.y, + particle.dir.z + ] + ); + } + } + } + + let len = total_output.len(); + let particles = total_output.as_mut_ptr(); + + std::mem::forget(total_output); + + *num_emitted_particles = len as c_int; + particles +} + #[no_mangle] pub extern "C" fn simple_bca_c(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, Z1: f64, m1: f64, Ec1: f64, Es1: f64, Z2: f64, m2: f64, Ec2: f64, Es2: f64, n2: f64, Eb2: f64) -> OutputBCA { let mut output = simple_bca(x, y, z, ux, uy, uz, E1, Z1, m1, Ec1, Es1, Z2, m2, Ec2, Es2, n2, Eb2); @@ -585,7 +1089,10 @@ pub fn simple_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1: f64, Z1 track_trajectories: false, number_collision_events: 0, backreflected: false, - interaction_index : 0 + interaction_index : 0, + weight: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0), }; let material_parameters = material::MaterialParameters { diff --git a/src/main.rs b/src/main.rs index 5eeb94a..a56fb99 100644 --- a/src/main.rs +++ b/src/main.rs @@ -16,8 +16,7 @@ use std::mem::discriminant; use indicatif::{ProgressBar, ProgressStyle}; //Error handling crate -use anyhow::Result; -use anyhow::*; +use anyhow::{Result, Context, anyhow}; //Serializing/Deserializing crate use serde::*; @@ -43,9 +42,6 @@ use std::f64::consts::FRAC_2_SQRT_PI; use std::f64::consts::PI; use std::f64::consts::SQRT_2; -//rng -//use rand::{Rng, thread_rng}; - //Load internal modules pub mod material; pub mod particle; diff --git a/src/particle.rs b/src/particle.rs index 7c73721..e26ac9e 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -2,7 +2,6 @@ use super::*; /// Rustbca's internal representation of the particle_parameters input. - #[cfg(feature = "hdf5_input")] #[derive(Deserialize, Clone)] pub struct ParticleParameters { @@ -13,12 +12,12 @@ pub struct ParticleParameters { pub N: Vec, pub m: Vec, pub Z: Vec, - pub E: Vec, + pub E: Vec, pub Ec: Vec, pub Es: Vec, - pub pos: Vec<(f64, f64, f64)>, - pub dir: Vec<(f64, f64, f64)>, - pub interaction_index: Vec + pub pos: Vec<(Distributions, Distributions, Distributions)>, + pub dir: Vec<(Distributions, Distributions, Distributions)>, + pub interaction_index: Vec, } #[cfg(not(feature = "hdf5_input"))] @@ -30,12 +29,12 @@ pub struct ParticleParameters { pub N: Vec, pub m: Vec, pub Z: Vec, - pub E: Vec, + pub E: Vec, pub Ec: Vec, pub Es: Vec, - pub pos: Vec<(f64, f64, f64)>, - pub dir: Vec<(f64, f64, f64)>, - pub interaction_index: Vec + pub pos: Vec<(Distributions, Distributions, Distributions)>, + pub dir: Vec<(Distributions, Distributions, Distributions)>, + pub interaction_index: Vec, } /// HDF5 version of particle input. @@ -55,6 +54,8 @@ pub struct ParticleInput { pub uy: f64, pub uz: f64, pub interaction_index: usize, + pub weight: f64, + pub tag: i32, } /// Particle object. Particles in rustbca include incident ions and material atoms. @@ -82,6 +83,9 @@ pub struct Particle { pub number_collision_events: usize, pub backreflected: bool, pub interaction_index: usize, + pub weight: f64, + pub tag: i32, + pub tracked_vector: Vector, } impl Particle { /// Construct a particle object from input. @@ -118,6 +122,9 @@ impl Particle { number_collision_events: 0, backreflected: false, interaction_index: input.interaction_index, + weight: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0), } } @@ -148,6 +155,9 @@ impl Particle { number_collision_events: 0, backreflected: false, interaction_index, + weight: 1.0, + tag: 0, + tracked_vector: Vector::new(0.0, 0.0, 0.0), } } diff --git a/src/structs.rs b/src/structs.rs index 8aa611a..2b9eb90 100644 --- a/src/structs.rs +++ b/src/structs.rs @@ -1,5 +1,5 @@ /// 3D vector. -#[derive(Clone)] +#[derive(Clone, Copy)] pub struct Vector { pub x: f64, pub y: f64, diff --git a/src/tests.rs b/src/tests.rs index f93667b..ea34548 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -249,8 +249,6 @@ fn test_distributions() { let cosz = 0.0; let mut particle = particle::Particle::new(mass, Z, E, Ec, Es, x, y, z, cosx, cosy, cosz, false, false, 0); - - let mut distributions = output::Distributions::new(&options); assert_eq!(distributions.x_range[0], 0.); assert_eq!(distributions.x_range[distributions.x_range.len() - 1], 10.);