diff --git a/Cargo.toml b/Cargo.toml index ee9637a..1183c7b 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -24,7 +24,7 @@ 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} +rcpr = { git = "https://github.com/drobnyjt/rcpr", branch="main", 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} diff --git a/examples/test_morse.py b/examples/test_morse.py index 2716eb5..3f17c54 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -13,6 +13,10 @@ hydrogen['Ec'] = 0.1 hydrogen['Es'] = 1.5 +epsilon = 1e-4 +interval_limit = 1e-3 +nmax = 32 +n0=3 #This function simply contains an entire input file as a multi-line f-string to modify some inputs. def run_morse_potential(energy, index, num_samples=10000, run_sim=True): @@ -25,7 +29,7 @@ def run_morse_potential(energy, index, num_samples=10000, run_sim=True): mean_free_path_model = "LIQUID" interaction_potential = [[{{"MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10}}}}]] scattering_integral = [["GAUSS_LEGENDRE"]] - root_finder = [[{{"CPR"={{n0=3, nmax=100, epsilon=1E-9, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-13, derivative_free=true}}}}]] + root_finder = [[{{"CPR"={{n0={n0}, nmax={nmax}, epsilon={epsilon}, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E22, interval_limit={interval_limit}, derivative_free=true}}}}]] num_threads = 4 num_chunks = 10 @@ -83,11 +87,11 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): name = "krc_morse_{index}" track_recoils = false weak_collision_order = 0 - electronic_stopping_mode = "LOW_ENERGY_NONLOCAL" + electronic_stopping_mode = "INTERPOLATED" mean_free_path_model = "LIQUID" - interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=7E10, x0=0.75E-10}}}}]] + interaction_potential = [[{{"KRC_MORSE"={{D=5.4971E-20, r0=2.782E-10, alpha=1.4198E10, k=8E10, x0=0.75E-10}}}}]] scattering_integral = [["GAUSS_LEGENDRE"]] - root_finder = [[{{"CPR"={{n0=2, nmax=200, epsilon=1E-9, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-13, derivative_free=true}}}}]] + root_finder = [[{{"CPR"={{n0=3, nmax={nmax}, epsilon={epsilon}, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E22, interval_limit={interval_limit}, derivative_free=true}}}}]] num_threads = 6 num_chunks = 1 @@ -170,10 +174,10 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): plt.semilogx(energies, r_benchmark, marker='^', linestyle='', label='Exp.') #Running and plotting the H-Ni simulations with the Morse potential and updated Es -num_energies = 15 +num_energies = 20 energies = np.logspace(-1, 4, num_energies) run_sim = True -num_samples = 100 +num_samples = 1000 R_N = np.zeros(num_energies) R_E = np.zeros(num_energies) R_N_2 = np.zeros(num_energies) @@ -184,26 +188,26 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): R_N_2[index], R_E_2[index] = run_morse_potential(energy, index, num_samples=num_samples, run_sim=run_sim) R_N_test = [ - 0.00, 0.01, 0.28, 0.60, 0.90, - 0.95, 0.88, 0.81, 0.70, 0.49, - 0.30, 0.24, 0.17, 0.11, 0.10 + 0.0, 0.028, 0.141, 0.327, 0.65, 0.832, 0.913, 0.926, + 0.889, 0.844, 0.783, 0.652, 0.47, 0.38, 0.344, 0.292, + 0.253, 0.197, 0.139, 0.084 ] R_N_2_test = [ - 0.00, 0.01, 0.28, 0.60, 0.90, - 0.95, 0.88, 0.81, 0.74, 0.60, - 0.46, 0.23, 0.05, 0.00, 0.00 + 0.0, 0.028, 0.141, 0.327, 0.65, 0.832, 0.913, 0.926, + 0.889, 0.844, 0.786, 0.722, 0.622, 0.531, 0.411, 0.254, + 0.104, 0.022, 0.003, 0.0 ] -np.testing.assert_allclose(R_N, R_N_test) -np.testing.assert_allclose(R_N_2, R_N_2_test) +np.testing.assert_allclose(R_N, R_N_test, atol=0.1) +np.testing.assert_allclose(R_N_2, R_N_2_test, atol=0.1) plt.semilogx(energies, R_N, label='R_N Morse-Kr-C H-Ni, Es=1.5eV', color='purple') plt.semilogx(energies, R_N_2, label='R_N Morse H-Ni, Es=1.5eV', color='green') #Plotting RustBCA data points, using the ergonomic helper function reflection_coefficient(). energies = np.logspace(-1, 4, 50) -r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 10000) for energy in energies]) +r_rustbca = np.array([reflection_coefficient(hydrogen, nickel, energy, 0.0, 1000) for energy in energies]) r_n = r_rustbca[:, 0] r_e = r_rustbca[:, 1] plt.semilogx(energies, r_n, label='R_N, Default Settings', color='black') diff --git a/src/bca.rs b/src/bca.rs index eb9a8ab..3e057b1 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -2,7 +2,11 @@ use super::*; use rand::RngExt; #[cfg(feature = "cpr_rootfinder")] -use rcpr::chebyshev::*; +use rcpr::rootfinders::{ + find_roots, + real_polynomial_roots, + Config +}; /// Geometrical quantities of binary collision. pub struct BinaryCollisionGeometry { @@ -587,6 +591,16 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact } } +#[cfg(feature = "cpr_rootfinder")] +fn transform(x: f64) -> f64 { + 2.0/(x*PI/2.).tan().powi(2) +} + +#[cfg(feature = "cpr_rootfinder")] +fn inverse_transform(x: f64) -> f64 { + 2./PI*((2.0/x).sqrt()).atan() +} + #[cfg(feature = "cpr_rootfinder")] /// Computes the distance of closest approach of two particles with atomic numbers `Za`, `Zb` and masses `Ma`, `Mb` for an arbitrary interaction potential (e.g., Morse) for a given impact parameter and incident energy `E0` using the Chebyshev-Proxy Root-Finder method. /// @@ -606,7 +620,7 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact /// `derivative_free`: if false, use Newton's method to polish roots from the CPR. If true, use the secant method. /// /// # Returns -/// Returns the distance of closest approach or an error if the root-finder failed. +/// Returns the distance of closest approach (reduced by a) or an error if the root-finder failed. pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parameter: f64, interaction_potential: InteractionPotential, n0: usize, nmax: usize, epsilon: f64, complex_threshold: f64, truncation_threshold: f64, far_from_zero: f64, @@ -614,32 +628,33 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame //Lindhard screening length and reduced energy let a = interactions::screening_length(Za, Zb, interaction_potential); - let reduced_energy = LINDHARD_REDUCED_ENERGY_PREFACTOR*a*Mb/(Ma+Mb)/Za/Zb*E0; let relative_energy = E0*Mb/(Ma + Mb); - let p = impact_parameter; - let f = |r: f64| -> f64 {interactions::distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)}; - let g = |r: f64| -> f64 {interactions::distance_of_closest_approach_function_singularity_free(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)* - interactions::scaling_function(r, impact_parameter, interaction_potential)}; - - let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential); - let lower_bound = impact_parameter / 1000.0; - - let roots = match derivative_free { - true => find_roots_with_secant_polishing(&g, &f, lower_bound, upper_bound, - n0, epsilon, nmax, complex_threshold, - truncation_threshold, interval_limit, far_from_zero), - - false => { - let df = |r: f64| -> f64 {interactions::diff_distance_of_closest_approach_function(r, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)}; - find_roots_with_newton_polishing(&g, &f, &df, lower_bound, upper_bound, - n0, epsilon, nmax, complex_threshold, - truncation_threshold, interval_limit, far_from_zero) - } - }.with_context(|| format!("Numerical error: CPR Rootfinder failed to converge when calculating distance of closest approach for Er = {} eV p = {} A using {}.", - relative_energy/EV, impact_parameter/ANGSTROM, interaction_potential))?; + let g = |r: f64| -> f64 { + interactions::distance_of_closest_approach_function_singularity_free(transform(r)*a, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)* + interactions::scaling_function(transform(r)*a, a, interaction_potential) + }; + + let upper_bound = 1.0; + let lower_bound = 1e-4; + + let delta = 1e-5; + let config = Config::new( + epsilon, + delta, + n0, + nmax, + complex_threshold, + truncation_threshold, + far_from_zero, + interval_limit + ); + + let roots = find_roots(&g, vec![(lower_bound, upper_bound)], config)?; - let max_root = roots.iter().cloned().fold(f64::NAN, f64::max)/a; + // Since above the arg to doca is transform(r)*a, this is already scaled as output + //let max_root = roots.iter().map(|&x| transform(x)).fold(f64::NAN, f64::max); + let max_root = roots.iter().map(|&x| transform(x)).max_by(f64::total_cmp).expect("Numerical error: failed to find maximum root."); if roots.is_empty() || max_root.is_nan() { return Err(anyhow!("Numerical error: CPR rootfinder failed to find root. x0: {}, F(a): {}, F(b): {};", max_root, g(0.), g(upper_bound))); diff --git a/src/interactions.rs b/src/interactions.rs index 5ba58ee..5ffb125 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -179,7 +179,7 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent 1./(1. + (r/ANGSTROM).powi(2)) }, InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => { - 1./(1. + (r*alpha).powi(2)) + 1. } InteractionPotential::COULOMB{..} => panic!("Coulombic potential cannot be used with rootfinder.") } @@ -282,11 +282,11 @@ pub fn dphi(xi: f64, interaction_potential: InteractionPotential) -> f64 { pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPotential) -> f64 { match interaction_potential { //ZBL screening length, Eckstein (4.1.8) - InteractionPotential::ZBL => zbl_screening_length_lookup(Za as u64, Zb as u64), + InteractionPotential::ZBL | InteractionPotential::WW => zbl_screening_length_lookup(Za as u64, Zb as u64), //Lindhard/Firsov screening length, Eckstein (4.1.5) - InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN | InteractionPotential::WW => lindhard_screening_length_lookup(Za as u64, Zb as u64), + InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::TRIDYN => lindhard_screening_length_lookup(Za as u64, Zb as u64), InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64), - InteractionPotential::MORSE{D, alpha, r0} => alpha, + InteractionPotential::MORSE{D, alpha, r0} => 1./alpha, InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => zbl_screening_length_lookup(Za as u64, Zb as u64), InteractionPotential::KRC_MORSE{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64), InteractionPotential::FOUR_EIGHT{..} => lindhard_screening_length_lookup(Za as u64, Zb as u64), @@ -434,7 +434,13 @@ pub fn doca_morse(r: f64, impact_parameter: f64, relative_energy: f64, D: f64, a /// Distance of closest approach function for Morse potential. pub fn doca_krc_morse(r: f64, impact_parameter: f64, relative_energy: f64, a: f64, Za: f64, Zb: f64, D: f64, alpha: f64, r0: f64, k: f64, x0: f64) -> f64 { - (r*alpha).powi(2) - (r*alpha).powi(2)/relative_energy*krc_morse(r, a, Za, Zb, D, alpha, r0, k, x0) - (impact_parameter*alpha).powi(2) + let K = coulomb_constant(Za, Zb); + let ralpha = r*alpha; + let term_1 = (ralpha).powi(2) - (impact_parameter*alpha).powi(2); + let term_2 = -(ralpha)*alpha*(K/relative_energy)*phi(r/a, InteractionPotential::KR_C)*smootherstep(r, -k, x0); + let term_3 = -(ralpha).powi(2)*(morse(r, D, alpha, r0)/relative_energy)*smootherstep(r, k, x0); + let scale = 1./(1. + ralpha).powi(2); + term_1*scale + term_2*scale + term_3*scale } /// Distance of closest approach function for LJ 6.5-6 potential. diff --git a/src/lib.rs b/src/lib.rs index 4ebb80f..6b0f24e 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -54,6 +54,7 @@ use pyo3::exceptions::{PyValueError, PyRuntimeError}; //Load internal modules pub mod material; pub mod particle; +#[cfg(test)] pub mod tests; pub mod interactions; pub mod bca; diff --git a/src/tests.rs b/src/tests.rs index afba5de..3c2b4b4 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -4,6 +4,8 @@ use super::*; use float_cmp::*; #[cfg(test)] use rand::RngExt; +#[cfg(feature = "cpr_rootfinder")] +use rcpr::rootfinders::*; #[test] #[cfg(feature = "cpr_rootfinder")]