diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index a50af72..f174909 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -32,9 +32,6 @@ jobs: - name: Install Python libraries run: | python3 -m pip install numpy shapely scipy matplotlib toml tomlkit - - name: Install HDF5 Libraries - run: | - sudo apt install libhdf5-dev - name: test Python Bindings and run testing scripts run: | python3 -m pip install setuptools_rust testresources setuptools wheel @@ -59,7 +56,7 @@ jobs: ./a.out - name: Test RustBCA run: | - cargo test --features cpr_rootfinder,hdf5_input,distributions,parry3d + cargo test --features cpr_rootfinder,distributions,parry3d - name: Run Examples run: | cargo run --release 0D examples/boron_nitride_0D.toml diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..4c20965 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +target/ +debug/ +Cargo.lock +*.output +**/*.rs.bk +*.pdb +src/RustBCA.egg-info +src/RustBCA.egg-info/* diff --git a/Cargo.toml b/Cargo.toml index 6933a13..c9379be 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,50 +1,48 @@ -[package] -name = "RustBCA" -version = "2.9.0" -default-run = "RustBCA" -authors = ["Jon Drobny ", "Jon Drobny "] -edition = "2021" - -[[bin]] -name = "RustBCA" -path = "src/main.rs" - -[lib] -name = "libRustBCA" -path = "src/lib.rs" -crate-type = ["cdylib", "lib"] - -[dependencies] -rand = {version = "0.10.2", features=["chacha"]} -rand_distr = "0.6.0" -toml = "1.1.4" -anyhow = "1.0.71" -itertools = "0.15.0" -rayon = "1.10.0" -geo = {version = "0.25", optional = false} -indicatif = {version = "0.18.6", features=["rayon"]} -serde = { version = "1.0.163", features = ["derive"] } -hdf5 = {version = "0.7.1", optional = true} -rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} -ndarray = {version = "0.17.2", features = ["serde"], optional = true} -parry3d-f64 = {optional = true, version="0.2.0"} -pyo3 = {version = "0.29.0", optional=true} -pythonize = {version = "0.29.0", optional=true} - -[dev-dependencies] -float-cmp = "0.10.0" - -[profile.release] -lto = "fat" -codegen-units = 1 -opt-level = 3 -debug = false - -[features] -hdf5_input = ["hdf5"] -cpr_rootfinder = ["rcpr"] -distributions = ["ndarray"] -no_list_output = [] -parry3d = ["parry3d-f64"] -python = ["pyo3", "pythonize"] +[package] +name = "RustBCA" +version = "3.0.0" +default-run = "RustBCA" +authors = ["Jon Drobny ", "Jon Drobny "] +edition = "2024" + +[[bin]] +name = "RustBCA" +path = "src/main.rs" + +[lib] +name = "libRustBCA" +path = "src/lib.rs" +crate-type = ["cdylib", "lib"] + +[dependencies] +rand = {version = "0.10.2", features=["chacha"]} +rand_distr = "0.6.0" +toml = "1.1.4" +anyhow = "1.0.71" +itertools = "0.15.0" +rayon = "1.10.0" +geo = {version = "0.25", optional = false} +indicatif = {version = "0.18.6", features=["rayon"]} +serde = { version = "1.0.163", features = ["derive"] } +rcpr = { git = "https://github.com/drobnyjt/rcpr", optional = true} +ndarray = {version = "0.17.2", features = ["serde"], optional = true} +parry3d-f64 = {optional = true, version="0.2.0"} +pyo3 = {version = "0.29.0", optional=true} +pythonize = {version = "0.29.0", optional=true} + +[dev-dependencies] +float-cmp = "0.10.0" + +[profile.release] +lto = "fat" +codegen-units = 1 +opt-level = 3 +debug = false + +[features] +cpr_rootfinder = ["rcpr"] +distributions = ["ndarray"] +no_list_output = [] +parry3d = ["parry3d-f64"] +python = ["pyo3", "pythonize"] extended_max_z = [] \ No newline at end of file diff --git a/examples/layered_geometry.toml b/examples/layered_geometry.toml index aaaa434..7590a59 100644 --- a/examples/layered_geometry.toml +++ b/examples/layered_geometry.toml @@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" mean_free_path_model = "LIQUID" interaction_potential = [["KR_C"]] scattering_integral = [["MENDENHALL_WELLER"]] -use_hdf5 = false root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]] num_threads = 4 num_chunks = 10 diff --git a/examples/layered_geometry_1D.toml b/examples/layered_geometry_1D.toml index 0c6938b..0bd2db7 100644 --- a/examples/layered_geometry_1D.toml +++ b/examples/layered_geometry_1D.toml @@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" mean_free_path_model = "LIQUID" interaction_potential = [["KR_C"]] scattering_integral = [["MENDENHALL_WELLER"]] -use_hdf5 = false root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]] num_threads = 4 num_chunks = 10 diff --git a/examples/lithium_vapor_shield.toml b/examples/lithium_vapor_shield.toml index ea887b9..981f14e 100644 --- a/examples/lithium_vapor_shield.toml +++ b/examples/lithium_vapor_shield.toml @@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_LOCAL" mean_free_path_model = {"THRESHOLD"={density=1e28}} interaction_potential = [["KR_C"]] scattering_integral = [["MENDENHALL_WELLER"]] -use_hdf5 = false root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]] num_threads = 4 num_chunks = 2 diff --git a/examples/multiple_interaction_potentials.toml b/examples/multiple_interaction_potentials.toml index aa8dc2d..bc821be 100644 --- a/examples/multiple_interaction_potentials.toml +++ b/examples/multiple_interaction_potentials.toml @@ -22,7 +22,6 @@ root_finder = [ [{"NEWTON"={max_iterations=100, tolerance=1E-6}}, {"POLYNOMIAL"={complex_threshold=1E-12}}], [{"POLYNOMIAL"={complex_threshold=1E-12}}, {"NEWTON"={max_iterations=100, tolerance=1E-6}}] ] -use_hdf5 = false num_threads = 4 num_chunks = 10 track_displacements = false diff --git a/examples/titanium_dioxide_0D.toml b/examples/titanium_dioxide_0D.toml index 178f855..7ebf8f7 100644 --- a/examples/titanium_dioxide_0D.toml +++ b/examples/titanium_dioxide_0D.toml @@ -12,7 +12,6 @@ electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" mean_free_path_model = {"THRESHOLD"={density=1E28}} interaction_potential = [["KR_C"]] scattering_integral = [["MENDENHALL_WELLER"]] -use_hdf5 = false root_finder = [[{"NEWTON"={max_iterations=100, tolerance=1E-6}}]] num_threads = 4 num_chunks = 100 diff --git a/examples/tungsten_tiles_trimesh.toml b/examples/tungsten_tiles_trimesh.toml index 793590c..e60dcf2 100644 --- a/examples/tungsten_tiles_trimesh.toml +++ b/examples/tungsten_tiles_trimesh.toml @@ -10,7 +10,6 @@ suppress_deep_recoils = false high_energy_free_flight_paths = false num_threads = 4 num_chunks = 5 -use_hdf5 = false electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" mean_free_path_model = "LIQUID" interaction_potential = [["KR_C"]] diff --git a/examples/tungsten_twist_trimesh.toml b/examples/tungsten_twist_trimesh.toml index e41cf9b..abb23e7 100644 --- a/examples/tungsten_twist_trimesh.toml +++ b/examples/tungsten_twist_trimesh.toml @@ -10,7 +10,6 @@ suppress_deep_recoils = false high_energy_free_flight_paths = false num_threads = 4 num_chunks = 5 -use_hdf5 = false electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" mean_free_path_model = "LIQUID" interaction_potential = [["KR_C"]] diff --git a/setup.py b/setup.py index 21ff36d..7d2962f 100644 --- a/setup.py +++ b/setup.py @@ -3,15 +3,12 @@ setup( name="RustBCA", - version="2.8.4", + version="3.0.0", rust_extensions=[ RustExtension( "libRustBCA", binding=Binding.PyO3, - features=["python", "parry3d"], - #args=["+nightly", "--edition 2021", "-Z unstable-options"], - #optional=True, - #rust_version="1.57.0" + features=["python", "parry3d", "pythonize"], ) ], # rust extensions are not zip safe, just like C-extensions. diff --git a/src/bca.rs b/src/bca.rs index 993c101..2915d1a 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -425,8 +425,6 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa pub fn subtract_electronic_stopping_energy(particle_1: &mut particle::Particle, material: &material::Material, distance_traveled: f64, x0: f64, strong_collision_Z: f64, strong_collision_index: usize, options: &Options) -> f64 { - //assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following collision."); - let x = particle_1.pos.x; let y = particle_1.pos.y; let z = particle_1.pos.z; @@ -452,9 +450,7 @@ pub fn subtract_electronic_stopping_energy(particle_1: &mut particl particle_1.E += -delta_energy_electronic; //Make sure particle energy doesn't become negative again - /* - assert!(!particle_1.E.is_nan(), "Numerical error: particle energy is NaN following electronic stopping."); - */ + if particle_1.E < 0. { particle_1.E = 0.; } @@ -479,17 +475,12 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: & let x0 = distance_of_closest_approach(particle_1, particle_2, binary_collision_geometry, options); let theta: f64 = match scattering_integral { - ScatteringIntegral::MENDENHALL_WELLER => mendenhall_weller(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential), - ScatteringIntegral::GAUSS_MEHLER{n_points} => gauss_mehler(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential, n_points), - ScatteringIntegral::GAUSS_LEGENDRE => gauss_legendre(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential), - ScatteringIntegral::MAGIC => magic(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, interaction_potential), + ScatteringIntegral::MENDENHALL_WELLER => mendenhall_weller(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential), + ScatteringIntegral::GAUSS_MEHLER{n_points} => gauss_mehler(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential, n_points), + ScatteringIntegral::GAUSS_LEGENDRE => gauss_legendre(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential), + ScatteringIntegral::MAGIC => magic(Za, Zb, Ma, Mb, E0, binary_collision_geometry.impact_parameter, x0, a, interaction_potential), }; - /* - if theta.is_nan() { - return Err(anyhow!("Numerical error: CoM deflection angle is NaN for {}. Check input parameters.", binary_collision_geometry)); - } - */ let (sin_theta, cos_theta) = theta.sin_cos(); let sin_2_theta = (theta/2.).sin(); @@ -695,8 +686,7 @@ pub fn newton_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par } /// Gauss-Mehler quadrature. -pub fn gauss_mehler(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential, n_points: usize) -> f64 { - let a: f64 = interactions::screening_length(Za, Zb, interaction_potential); +pub fn gauss_mehler(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential, n_points: usize) -> f64 { let r0 = x0*a; let V = |r| {interactions::interaction_potential(r, a, Za, Zb, interaction_potential)}; let relative_energy = E0*Mb/(Ma + Mb); @@ -704,8 +694,7 @@ pub fn gauss_mehler(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_paramete } /// Gauss-Legendre quadrature. -pub fn gauss_legendre(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential) -> f64 { - let a: f64 = interactions::screening_length(Za, Zb, interaction_potential); +pub fn gauss_legendre(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential) -> f64 { let r0 = x0*a; let V = |r| {interactions::interaction_potential(r, a, Za, Zb, interaction_potential)}; let relative_energy = E0*Mb/(Ma + Mb); @@ -713,7 +702,7 @@ pub fn gauss_legendre(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame } /// Ziegler's MAGIC algorithm for approximating the scattering integral. -pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential) -> f64 { +pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential) -> f64 { //MAGIC algorithm //Since this is legacy code I don't think I will clean this up let C_ = match interaction_potential { @@ -723,7 +712,6 @@ pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, InteractionPotential::TRIDYN => [1.0144, 0.235809, 0.126, 69350., 83550.], //Undocumented Tridyn constants _ => panic!("Input error: unimplemented interaction potential {} for MAGIC algorithm. Use a screened Coulomb potential.", interaction_potential) }; - let a: f64 = interactions::screening_length(Za, Zb, interaction_potential); let beta: f64 = impact_parameter/a; let V0 = Za*Zb*Q*Q/4.0/PI/EPS0/a; let relative_energy = E0*Mb/(Ma + Mb); @@ -740,9 +728,8 @@ pub fn magic(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, } /// Mendenhall-Weller quadrature. -pub fn mendenhall_weller(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, interaction_potential: InteractionPotential) -> f64 { +pub fn mendenhall_weller(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, x0: f64, a: f64, interaction_potential: InteractionPotential) -> f64 { //Lindhard screening length and reduced energy - let a: f64 = interactions::screening_length(Za, Zb, interaction_potential); let reduced_energy: f64 = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E0; let beta: f64 = impact_parameter/a; diff --git a/src/input.rs b/src/input.rs index e6c607f..0a95d6e 100644 --- a/src/input.rs +++ b/src/input.rs @@ -246,8 +246,6 @@ pub struct Options { #[serde(default = "one_u64")] pub num_chunks: u64, #[serde(default = "default_false")] - pub use_hdf5: bool, - #[serde(default = "default_false")] pub track_displacements: bool, #[serde(default = "default_false")] pub track_energy_losses: bool, @@ -274,7 +272,6 @@ impl Options { root_finder: default_rootfinder(), num_threads: 1, num_chunks: 1, - use_hdf5: false, track_displacements: false, track_energy_losses: false, seed: default_seed(), @@ -315,8 +312,6 @@ pub struct Options { #[serde(default = "one_u64")] pub num_chunks: u64, #[serde(default = "default_false")] - pub use_hdf5: bool, - #[serde(default = "default_false")] pub track_displacements: bool, #[serde(default = "default_false")] pub track_energy_losses: bool, @@ -358,7 +353,6 @@ impl Options { root_finder: default_rootfinder(), num_threads: 1, num_chunks: 1, - use_hdf5: false, track_displacements: false, track_energy_losses: false, energy_min: 0.0, @@ -552,160 +546,73 @@ pub fn process_input_file(input: ::InputFileFormat) .unwrap_or_else(|_| panic!("Input errror: could nor parse mass unit {}. Use a valid float or one of AMU, KG", &particle_parameters.mass_unit.as_str())), }; - //HDF5 - #[cfg(feature = "hdf5_input")] let particle_input_array: Vec = { - if options.use_hdf5 { - let particle_input_filename = particle_parameters.particle_input_filename.as_str(); - let _e = hdf5::silence_errors(); - let particle_input_file = hdf5::File::open(particle_input_filename) - .context("Input error: cannot open HDF5 file.") - .unwrap(); - let particle_input = particle_input_file.dataset("particles") - .context("Input error: cannot read from HDF5 file.") - .unwrap(); - particle_input.read_raw::().unwrap() - - } else { - let mut particle_input: Vec = Vec::new(); - - for particle_index in 0..N { - let N_ = particle_parameters.N[particle_index]; - let m = particle_parameters.m[particle_index]; - let Z = particle_parameters.Z[particle_index]; - let E = particle_parameters.E[particle_index]; - 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]; - - 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: match E { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*energy_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*energy_unit}, - Distributions::POINT(E) => E*energy_unit, - }, - Ec: Ec*energy_unit, - Es: Es*energy_unit, - x: match x { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut 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 rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(y) => y*length_unit, - }, - z: match z { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(z) => z*length_unit, - }, - ux: match cosx { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(ux) => ux, - }, - uy: match cosy { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(uy) => uy, - }, - uz: match cosz { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(uz) => uz, - }, - interaction_index: interaction_index, - tag: 0, - weight: 1.0, - } - ); - } - } - particle_input - } - }; - #[cfg(not(feature = "hdf5_input"))] - let particle_input_array: Vec = { - if options.use_hdf5 { - panic!("HDF5 particle input not enabled. Enable with: cargo build --features hdf5_input") - } else { - let mut particle_input: Vec = Vec::new(); - - for particle_index in 0..N { - let N_ = particle_parameters.N[particle_index]; - let m = particle_parameters.m[particle_index]; - let Z = particle_parameters.Z[particle_index]; - let E = particle_parameters.E[particle_index]; - 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]; - - for sub_particle_index in 0..N_ { - - //Add new particle to particle vector - particle_input.push( - particle::ParticleInput{ - m: m*mass_unit, - Z, - E: match E { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*energy_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*energy_unit}, - Distributions::POINT(x) => x*energy_unit, - }, - Ec: Ec*energy_unit, - Es: Es*energy_unit, - x: match x { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut 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 rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(y) => y*length_unit, - }, - z: match z { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(z) => z*length_unit, - }, - ux: match cosx { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(ux) => ux - }, - uy: match cosy { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(uy) => uy, - }, - uz: match cosz { - Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, - Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, - Distributions::POINT(uz) => uz, - }, - interaction_index, - tag: 0, - weight: 1.0, - } - ); - } + let mut particle_input: Vec = Vec::new(); + + for particle_index in 0..N { + let N_ = particle_parameters.N[particle_index]; + let m = particle_parameters.m[particle_index]; + let Z = particle_parameters.Z[particle_index]; + let E = particle_parameters.E[particle_index]; + 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]; + + for sub_particle_index in 0..N_ { + + //Add new particle to particle vector + particle_input.push( + particle::ParticleInput{ + m: m*mass_unit, + Z, + E: match E { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*energy_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*energy_unit}, + Distributions::POINT(x) => x*energy_unit, + }, + Ec: Ec*energy_unit, + Es: Es*energy_unit, + x: match x { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut 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 rng)*length_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, + Distributions::POINT(y) => y*length_unit, + }, + z: match z { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, + Distributions::POINT(z) => z*length_unit, + }, + ux: match cosx { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, + Distributions::POINT(ux) => ux + }, + uy: match cosy { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, + Distributions::POINT(uy) => uy, + }, + uz: match cosz { + Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rng)*length_unit}, + Distributions::UNIFORM{min, max} => {let uniform = Uniform::new(min, max).unwrap(); uniform.sample(&mut rng)*length_unit}, + Distributions::POINT(uz) => uz, + }, + interaction_index, + tag: 0, + weight: 1.0, + } + ); } - particle_input } + particle_input }; (particle_input_array, material, options, OutputUnits {length_unit, energy_unit, mass_unit}) } diff --git a/src/interactions.rs b/src/interactions.rs index 2c22b78..2def3ac 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -541,7 +541,6 @@ fn kr_c(xi: f64) -> f64 { 0.190945*(-0.278544*xi).exp() + 0.473674*(-0.637174*xi).exp() + 0.335381*(-1.919249*xi).exp() } - fn zbl(xi: f64) -> f64 { 0.02817*(-0.20162*xi).exp() + 0.28022*(-0.40290*xi).exp() + 0.50986*(-0.94229*xi).exp() + 0.18175*(-3.1998*xi).exp() } @@ -558,7 +557,6 @@ fn diff_kr_c(xi: f64) -> f64 { -0.278544*0.190945*(-0.278544*xi).exp() - 0.637174*0.473674*(-0.637174*xi).exp() - 0.335381*1.919249*(-1.919249*xi).exp() } - fn diff_zbl(xi: f64) -> f64 { -0.20162*0.02817*(-0.20162*xi).exp() -0.40290*0.28022*(-0.40290*xi).exp() -0.94229*0.50986*(-0.94229*xi).exp() -3.1998*0.18175*(-3.1998*xi).exp() } diff --git a/src/lib.rs b/src/lib.rs index 67ebc4a..1e282db 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -20,10 +20,6 @@ use anyhow::{Result, Context, anyhow}; //Serializing/Deserializing crate use serde::*; -//Array input via hdf5 -#[cfg(feature = "hdf5_input")] -use hdf5::*; - //I/O use std::fs::OpenOptions; use std::io::prelude::*; @@ -229,7 +225,7 @@ pub struct OutputTaggedBCA { pub incident: *mut bool, } -#[no_mangle] +#[unsafe(no_mangle)] pub extern "C" fn drop_output_tagged_bca(output: OutputTaggedBCA) { let length = output.len; @@ -249,7 +245,7 @@ pub extern "C" fn drop_output_tagged_bca(output: OutputTaggedBCA) { } } -#[no_mangle] +#[unsafe(no_mangle)] pub extern "C" fn drop_output_bca(output: OutputBCA) { let length = output.len; @@ -262,7 +258,7 @@ pub extern "C" fn drop_output_bca(output: OutputBCA) { } } -#[no_mangle] +#[unsafe(no_mangle)] pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTaggedBCA { let mut total_output = vec![]; @@ -402,7 +398,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag } } -#[no_mangle] +#[unsafe(no_mangle)] pub unsafe 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); @@ -484,7 +480,7 @@ pub unsafe extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux } } -#[no_mangle] +#[unsafe(no_mangle)] pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { let x = -2.*(input.n2*1E30).powf(-1./3.); @@ -602,7 +598,7 @@ pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { } } -#[no_mangle] +#[unsafe(no_mangle)] pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { let mut total_output = vec![]; @@ -727,7 +723,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { } } -#[no_mangle] +#[unsafe(no_mangle)] pub unsafe 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, @@ -856,7 +852,7 @@ pub unsafe extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int particles } -#[no_mangle] +#[unsafe(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); @@ -1623,7 +1619,7 @@ pub fn simple_compound_bca(x: f64, y: f64, z: f64, ux: f64, uy: f64, uz: f64, E1 } #[cfg(feature = "parry3d")] -#[no_mangle] +#[unsafe(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) { let direction = Vector3::new(*ux, *uy, *uz); @@ -1711,7 +1707,7 @@ pub fn rotate_given_surface_normal_vec_py<'py>(nx: Vec, ny: Vec, nz: V } #[cfg(feature = "parry3d")] -#[no_mangle] +#[unsafe(no_mangle)] pub extern "C" fn rotate_back(nx: f64, ny: f64, nz: f64, ux: &mut f64, uy: &mut f64, uz: &mut f64) { let direction = Vector3::new(*ux, *uy, *uz); @@ -2188,15 +2184,17 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g _ => return Err(PyValueError::new_err(format!("Unimplemented interaction potential {}; try 'KR_C'", interaction_potential))) }; + let screening_length = interactions::screening_length(Za, Zb, potential); + let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, potential, 1000, 1E-12).map_err( |error| PyRuntimeError::new_err(format!("Rootfinder failed to find distance of closest approach; check input values.")) )?; //Compute center of mass deflection angle with each algorithm - let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, potential, n_gl_points); - let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, potential); - let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, potential); - let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, potential); + let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, screening_length, potential, n_gl_points); + let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, screening_length, potential); + let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, screening_length, potential); + let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, screening_length, potential); Ok((theta_gm, theta_gl, theta_mw, theta_magic)) } diff --git a/src/main.rs b/src/main.rs index 0920812..9422bd3 100644 --- a/src/main.rs +++ b/src/main.rs @@ -11,10 +11,6 @@ use anyhow::{Result, Context, anyhow}; //Serializing/Deserializing crate use serde::*; -//Array input via hdf5 -#[cfg(feature = "hdf5_input")] -use hdf5::*; - //Parallelization //use rayon::prelude::*; //use rayon::ThreadPoolBuilder; diff --git a/src/particle.rs b/src/particle.rs index 3132ba3..f64f3b2 100644 --- a/src/particle.rs +++ b/src/particle.rs @@ -4,10 +4,8 @@ fn default_vec_zero() -> Vec { vec![0] } -#[cfg(feature = "hdf5_input")] #[derive(Deserialize, Clone)] pub struct ParticleParameters { - pub particle_input_filename: String, pub length_unit: String, pub energy_unit: String, pub mass_unit: String, @@ -23,27 +21,7 @@ pub struct ParticleParameters { pub interaction_index: Vec, } -#[cfg(not(feature = "hdf5_input"))] -#[derive(Deserialize, Clone)] -pub struct ParticleParameters { - pub length_unit: String, - pub energy_unit: String, - pub mass_unit: String, - pub N: Vec, - pub m: Vec, - pub Z: Vec, - pub E: Vec, - pub Ec: Vec, - pub Es: Vec, - pub pos: Vec<(Distributions, Distributions, Distributions)>, - pub dir: Vec<(Distributions, Distributions, Distributions)>, - #[serde(default = "default_vec_zero")] - pub interaction_index: Vec, -} - -/// HDF5 version of particle input. #[derive(Clone, PartialEq, Debug, Copy)] -#[cfg_attr(feature = "hdf5_input", derive(hdf5::H5Type))] #[repr(C)] pub struct ParticleInput { pub m: f64, @@ -239,13 +217,9 @@ impl Particle { let (sinpsi, cospsi) = psi.sin_cos(); let (e1, e2) = math::duff_orthonormal_basis(self.dir); - let cosx_new = cospsi*cosx - sinpsi*(cosphi*e1.x + sinphi*e2.x); - let cosy_new = cospsi*cosy - sinpsi*(cosphi*e1.y + sinphi*e2.y); - let cosz_new = cospsi*cosz - sinpsi*(cosphi*e1.z + sinphi*e2.z); - - let dir_new = Vector {x: cosx_new, y: cosy_new, z: cosz_new}; - - self.dir.assign(&dir_new); + self.dir.x = cospsi*cosx - sinpsi*(cosphi*e1.x + sinphi*e2.x); + self.dir.y = cospsi*cosy - sinpsi*(cosphi*e1.y + sinphi*e2.y); + self.dir.z = cospsi*cosz - sinpsi*(cosphi*e1.z + sinphi*e2.z); self.dir.normalize(); } diff --git a/src/tests.rs b/src/tests.rs index 8d68d80..d346226 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -226,7 +226,6 @@ fn test_distributions() { scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], num_threads: 1, num_chunks: 1, - use_hdf5: false, root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}]], track_displacements: false, track_energy_losses: true, @@ -882,7 +881,6 @@ fn test_momentum_conservation() { scattering_integral: vec![vec![scattering_integral]], num_threads: 1, num_chunks: 1, - use_hdf5: false, root_finder: vec![vec![root_finder]], track_displacements: false, track_energy_losses: false, @@ -905,7 +903,6 @@ fn test_momentum_conservation() { scattering_integral: vec![vec![scattering_integral]], num_threads: 1, num_chunks: 1, - use_hdf5: false, root_finder: vec![vec![root_finder]], track_displacements: false, track_energy_losses: false, @@ -1118,7 +1115,6 @@ fn test_quadrature() { scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], num_threads: 1, num_chunks: 1, - use_hdf5: false, root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-14}]], track_displacements: false, track_energy_losses: false, @@ -1141,7 +1137,6 @@ fn test_quadrature() { scattering_integral: vec![vec![ScatteringIntegral::MENDENHALL_WELLER]], num_threads: 1, num_chunks: 1, - use_hdf5: false, root_finder: vec![vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-14}]], track_displacements: false, track_energy_losses: false, @@ -1163,20 +1158,24 @@ fn test_quadrature() { z_num: 11, }; - let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 100, 1E-12).unwrap(); + let interaction_potential = InteractionPotential::KR_C; + + let x0_newton = bca::newton_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, 100, 1E-12).unwrap(); //If cpr_rootfinder is enabled, compare Newton to CPR - they should be nearly identical #[cfg(feature = "cpr_rootfinder")] - if let Ok(x0_cpr) = bca::cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 2, 10000, 1E-6, 1E-6, 1E-9, 1E9, 1E-13, true) { + if let Ok(x0_cpr) = bca::cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, 2, 10000, 1E-6, 1E-6, 1E-9, 1E9, 1E-13, true) { println!("CPR: {} Newton: {}", x0_cpr, x0_newton); assert!(approx_eq!(f64, x0_newton, x0_cpr, epsilon=1E-3)); }; + let a = interactions::screening_length(Za, Zb, interaction_potential); + //Compute center of mass deflection angle with each algorithm - let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C, 10); - let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); - let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); - let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, InteractionPotential::KR_C); + let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0_newton, a, interaction_potential, 10); + let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0_newton, a, interaction_potential); + let theta_mw = bca::mendenhall_weller(Za, Zb, Ma, Mb, E0, p, x0_newton, a, interaction_potential); + let theta_magic = bca::magic(Za, Zb, Ma, Mb, E0, p, x0_newton, a, interaction_potential); //Gauss-Mehler and Gauss-Legendre should be very close to each other assert!(approx_eq!(f64, theta_gm, theta_gl, epsilon=0.001));