diff --git a/.github/workflows/rustbca_compile_check.yml b/.github/workflows/rustbca_compile_check.yml index 2830b7f..c42b665 100644 --- a/.github/workflows/rustbca_compile_check.yml +++ b/.github/workflows/rustbca_compile_check.yml @@ -1,6 +1,7 @@ name: RustBCA Compile check on: + workflow_dispatch: push: branches: [ main, dev ] pull_request: diff --git a/src/bca.rs b/src/bca.rs index 353c24d..60e30a6 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -538,7 +538,7 @@ pub fn calculate_binary_collision(particle_1: &particle::Particle, particle_2: & /// Mendenhall-Weller scattering integrand. fn scattering_integral_mw(x: f64, beta: f64, reduced_energy: f64, interaction_potential: InteractionPotential) -> f64 { //Function for scattering integral - see Mendenhall and Weller, 1991 & 2005 - return (1. - interactions::phi(x, interaction_potential)/x/reduced_energy - beta*beta/x/x).powf(-0.5); + return 1./(1. - interactions::phi(x, interaction_potential)/x/reduced_energy - beta*beta/x/x).sqrt(); } /// Gauss-Legendre scattering integrand. @@ -776,7 +776,7 @@ pub fn mendenhall_weller(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_par let beta: f64 = impact_parameter/a; //Scattering integral quadrature from Mendenhall and Weller 2005 - let lambda0 = (0.5 + beta*beta/x0/x0/2. - interactions::dphi(x0, interaction_potential)/2./reduced_energy).powf(-1./2.); + let lambda0 = 1./(0.5 + beta*beta/x0/x0/2. - interactions::dphi(x0, interaction_potential)/2./reduced_energy).sqrt(); let alpha = 1./12.*(1. + lambda0 + 5.*(0.4206*scattering_integral_mw(x0/0.9072, beta, reduced_energy, interaction_potential) + 0.9072*scattering_integral_mw(x0/0.4206, beta, reduced_energy, interaction_potential))); PI*(1. - beta*alpha/x0) } diff --git a/src/geometry.rs b/src/geometry.rs index a915629..38cfb2d 100644 --- a/src/geometry.rs +++ b/src/geometry.rs @@ -70,7 +70,7 @@ impl Geometry for Mesh0D { let total_density: f64 = densities.iter().sum(); - let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; + let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.; let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); @@ -180,8 +180,8 @@ impl Geometry for Mesh1D { layer_top = layer_bottom; } - let top_energy_barrier_thickness = layers[0].densities.iter().sum::().powf(-1./3.)/SQRTPI*2.; - let bottom_energy_barrier_thickness = layers[layers.len() - 1].densities.iter().sum::().powf(-1./3.)/SQRTPI*2.; + let top_energy_barrier_thickness = 1./layers[0].densities.iter().sum::().cbrt()/SQRTPI*2.; + let bottom_energy_barrier_thickness = 1./layers[layers.len() - 1].densities.iter().sum::().cbrt()/SQRTPI*2.; Mesh1D { layers, @@ -315,7 +315,7 @@ impl Geometry for HomogeneousMesh2D { let total_density: f64 = densities.iter().sum(); - let energy_barrier_thickness = total_density.powf(-1./3.)/SQRTPI*2.; + let energy_barrier_thickness = 1./total_density.cbrt()/SQRTPI*2.; let concentrations: Vec = densities.iter().map(|&density| density/total_density).collect::>(); @@ -361,7 +361,7 @@ impl Geometry for HomogeneousMesh2D { true } else { if let Closest::SinglePoint(p) = self.boundary.closest_point(&point!(x: x, y: y)) { - let distance = ((x - p.x()).powf(2.) + (y - p.y()).powf(2.)).sqrt(); + let distance = ((x - p.x()).powi(2) + (y - p.y()).powi(2)).sqrt(); distance < self.energy_barrier_thickness } else if let Closest::Intersection(p) = self.boundary.closest_point(&point!(x: x, y: y)) { true diff --git a/src/interactions.rs b/src/interactions.rs index 921aa7e..e15ab94 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -166,15 +166,12 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent 1./(1. + (r/a).powi(2)) }, InteractionPotential::LENNARD_JONES_12_6{sigma, ..} => { - let n = 11.; - 1./(1. + (r/sigma).powf(n)) + 1./(1. + (r/sigma).powi(11)) }, InteractionPotential::LENNARD_JONES_65_6{sigma, ..} => { - let n = 6.; - 1./(1. + (r/sigma).powf(n)) + 1./(1. + (r/sigma).powi(6)) }, InteractionPotential::FOUR_EIGHT{alpha, beta} => { - let n = 8.; 1./(1. + r.powi(8)/beta) } InteractionPotential::MORSE{D, alpha, r0} => { @@ -294,15 +291,14 @@ pub fn polynomial_coefficients(relative_energy: f64, impact_parameter: f64, inte let epsilon_ev = epsilon/EV; let sigma_angstroms = sigma/ANGSTROM; let relative_energy_ev = relative_energy/EV; - //vec![1., 0., -impact_parameter.powi(2), 0., 0., 0., 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, 0., 0., 0., 0., 0., -4.*epsilon_ev*sigma.powf(12.)/relative_energy_ev] - vec![1.0, -impact_parameter.powi(2), 0.0, 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, 0.0, 0.0, -4.*epsilon_ev*sigma.powf(12.)/relative_energy_ev] + vec![1.0, -impact_parameter.powi(2), 0.0, 4.*epsilon_ev*sigma.powi(6)/relative_energy_ev, 0.0, 0.0, -4.*epsilon_ev*sigma.powi(12)/relative_energy_ev] }, InteractionPotential::LENNARD_JONES_65_6{sigma, epsilon} => { let impact_parameter_angstroms = impact_parameter/ANGSTROM; let epsilon_ev = epsilon/EV; let sigma_angstroms = sigma/ANGSTROM; let relative_energy_ev = relative_energy/EV; - vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon_ev*sigma.powf(6.)/relative_energy_ev, -4.*epsilon_ev*sigma.powf(6.5)/relative_energy_ev] + vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon_ev*sigma.powi(6)/relative_energy_ev, -4.*epsilon_ev*sigma.powf(6.5)/relative_energy_ev] }, InteractionPotential::FOUR_EIGHT{alpha, beta} => { //Note: I've transformed to angstroms here to help the rootfinder with numerical issues. @@ -343,12 +339,12 @@ pub fn four_eight(r: f64, alpha: f64, beta: f64) -> f64 { /// Lennard-Jones 12-6 pub fn lennard_jones(r: f64, sigma: f64, epsilon: f64) -> f64 { - 4.*epsilon*((sigma/r).powf(12.) - (sigma/r).powf(6.)) + 4.*epsilon*((sigma/r).powi(12) - (sigma/r).powi(6)) } /// Lennard-Jones 6.5-6 pub fn lennard_jones_65_6(r: f64, sigma: f64, epsilon: f64) -> f64 { - 4.*epsilon*((sigma/r).powf(6.5) - (sigma/r).powf(6.)) + 4.*epsilon*((sigma/r).powf(6.5) - (sigma/r).powi(6)) } /// Morse potential @@ -366,7 +362,7 @@ pub fn doca_four_eight(r: f64, impact_parameter: f64, relative_energy: f64, alph let a = alpha.powf(1./4.); let b = beta.powf(1./8.); let b4 = beta.sqrt(); - (r/b).powf(8.) - (-(a*r/b/b) + 1.)/relative_energy - (impact_parameter*r.powf(3.)/b4) + (r/b).powi(8) - (-(a*r/b/b) + 1.)/relative_energy - (impact_parameter*r.powi(3)/b4) } /// Distance of closest approach function for Morse potential. @@ -391,17 +387,17 @@ pub fn doca_lennard_jones_65_6(r: f64, p: f64, relative_energy: f64, sigma: f64, /// Distance of closest approach function for LJ 12-6 potential. pub fn doca_lennard_jones(r: f64, p: f64, relative_energy: f64, sigma: f64, epsilon: f64) -> f64 { - (r/sigma).powf(12.) - 4.*epsilon/relative_energy*(1. - (r/sigma).powf(6.)) - p.powi(2)*r.powf(10.)/sigma.powf(12.) + (r/sigma).powi(12) - 4.*epsilon/relative_energy*(1. - (r/sigma).powi(6)) - p.powi(2)*r.powi(10)/sigma.powi(12) } /// First derivative w.r.t. `r` of the distance of closest approach function for LJ 12-6 potential. pub fn diff_doca_lennard_jones(r: f64, p: f64, relative_energy: f64, sigma: f64, epsilon: f64) -> f64 { - 12.*(r/sigma).powf(11.)/sigma + 4.*epsilon/relative_energy*6.*(r/sigma).powf(5.)/sigma - 10.*p.powi(2)*r.powf(9.)/sigma.powf(12.) + 12.*(r/sigma).powi(11)/sigma + 4.*epsilon/relative_energy*6.*(r/sigma).powi(5)/sigma - 10.*p.powi(2)*r.powi(9)/sigma.powi(12) } /// First derivative w.r.t. `r` of the distance of closest approach function for LJ 6.5-6 potential. pub fn diff_doca_lennard_jones_65_6(r: f64, p: f64, relative_energy: f64, sigma: f64, epsilon: f64) -> f64 { - 6.5*(r/sigma).powf(5.5)/sigma + 4.*epsilon/relative_energy*0.5*(sigma*r).powf(-0.5) - (p/sigma).powi(2)*4.5*(r/sigma).powf(3.5)/sigma + 6.5*(r/sigma).powf(5.5)/sigma + 4.*epsilon/relative_energy*0.5/(sigma*r).sqrt() - (p/sigma).powi(2)*4.5*(r/sigma).powf(3.5)/sigma } /// W-W cublic spline potential from Bjorkas et al. @@ -427,7 +423,7 @@ pub fn tungsten_tungsten_cubic_spline(r: f64) -> f64 { -0.050264585985867E4 ]; - (a[0] + a[1]*x + a[2]*x.powi(2) + a[3]*x.powi(3) + a[4]*x.powf(4.) + a[5]*x.powf(5.))*EV + (a[0] + a[1]*x + a[2]*x.powi(2) + a[3]*x.powi(3) + a[4]*x.powi(4) + a[5]*x.powi(5))*EV } else { diff --git a/src/lib.rs b/src/lib.rs index 41ea254..75fc018 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -235,7 +235,7 @@ pub extern "C" fn compound_tagged_bca_list_c(input: InputTaggedBCA) -> OutputTag let tags = unsafe { slice::from_raw_parts(input.tags, input.len).to_vec() }; let weights = unsafe { slice::from_raw_parts(input.weights, input.len).to_vec() }; - let x = -2.*(n2.iter().sum::()*10E30).powf(-1./3.); + let x = -2.*(n2.iter().sum::()*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -368,7 +368,7 @@ pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut 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 x = -2.*(n2.iter().sum::()*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -438,7 +438,7 @@ pub extern "C" fn reflect_single_ion_c(num_species_target: &mut c_int, ux: &mut #[no_mangle] pub extern "C" fn simple_bca_list_c(input: InputSimpleBCA) -> OutputBCA { - let x = -2.*(input.n2*10E30).powf(-1./3.); + let x = -2.*(input.n2*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -560,7 +560,7 @@ pub extern "C" fn compound_bca_list_c(input: InputCompoundBCA) -> OutputBCA { 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 x = -2.*(n2.iter().sum::()*10E30).powf(-1./3.); + let x = -2.*(n2.iter().sum::()*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -699,7 +699,7 @@ pub extern "C" fn compound_bca_list_fortran(num_incident_ions: &mut c_int, track //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 x = -2.*(n2.iter().sum::()*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -844,7 +844,7 @@ pub fn compound_bca_list_py(energies: Vec, ux: Vec, uy: Vec, uz: let options = Options::default_options(true); - let x = -2.*(n2.iter().sum::()*10E30).powf(-1./3.); + let x = -2.*(n2.iter().sum::()*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -968,7 +968,7 @@ pub fn compound_bca_list_tracked_py(energies: Vec, ux: Vec, uy: Vec()*10E30).powf(-1./3.); + let x = -2.*(n2.iter().sum::()*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -1333,7 +1333,7 @@ pub fn simple_bca_list_py(energies: Vec, usx: Vec, usy: Vec, usz: assert_eq!(energies.len(), usy.len()); assert_eq!(energies.len(), usz.len()); - let x = -2.*(n2*10E30).powf(-1./3.); + let x = -2.*(n2*1E30).powf(-1./3.); let y = 0.0; let z = 0.0; @@ -1525,7 +1525,7 @@ pub extern "C" fn rotate_given_surface_normal(nx: f64, ny: f64, nz: f64, ux: &mu *ux = incident.x; *uy = incident.y; *uz = incident.z; - let mag = (ux.powf(2.) + uy.powf(2.) + uz.powf(2.)).sqrt(); + let mag = (ux*ux + uy*uy + uz*uz).sqrt(); *ux /= mag; *uy /= mag; diff --git a/src/material.rs b/src/material.rs index e6931fa..ae8bd18 100644 --- a/src/material.rs +++ b/src/material.rs @@ -121,7 +121,7 @@ impl Material { /// Determines the local mean free path from the formula sum(n(x, y))^(-1/3) pub fn mfp(&self, x: f64, y: f64, z: f64) -> f64 { - return self.total_number_density(x, y, z).powf(-1./3.); + return 1./self.total_number_density(x, y, z).cbrt(); } /// Total number density of triangle that contains or is nearest to (x, y). @@ -274,7 +274,7 @@ impl Material { for (n, Zb) in self.number_densities(x, y, z).iter().zip(&self.Z) { - let beta = (1. - (1. + E/Ma/C.powi(2)).powf(-2.)).sqrt(); + let beta = (1. - 1./(1. + E/Ma/C.powi(2)).powi(2)).sqrt(); let v = beta*C; // This term is an empirical fit to the mean ionization potential @@ -296,11 +296,12 @@ impl Material { let S_high = prefactor*(eb + 1. + B/eb).ln(); //Lindhard-Scharff electronic stopping - let S_low = LINDHARD_SCHARFF_PREFACTOR*(Za.powf(7./6.)*Zb)/(Za.powf(2./3.) + Zb.powf(2./3.)).powf(3./2.)*(E/Q/Ma*AMU).sqrt(); + //let S_low = LINDHARD_SCHARFF_PREFACTOR*(Za.powf(7./6.)*Zb)/(Za.powf(2./3.) + Zb.powf(2./3.)).powf(3./2.)*(E/Q/Ma*AMU).sqrt(); + let S_low = LINDHARD_SCHARFF_PREFACTOR*(Za*Za.cbrt().sqrt()*Zb)/(Za.cbrt().powi(2) + Zb.cbrt().powi(2)).powi(3).sqrt()*(E/Q/Ma*AMU).sqrt(); let stopping_power = match electronic_stopping_mode { //Biersack-Varelas Interpolation - ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low*ck), + ElectronicStoppingMode::INTERPOLATED => 1./(1./S_high + 1./S_low)*ck, //Oen-Robinson ElectronicStoppingMode::LOW_ENERGY_LOCAL => S_low*ck, //Lindhard-Scharff