New rcpr

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}

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')

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,40 +620,41 @@ 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,
     interval_limit: f64, derivative_free: bool) -> Result <f64, anyhow::Error> {
 
     //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)));

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.

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;

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")]