From edde2773124aa2f02eb49180d7fdc813b5aaf1bb Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 18 Aug 2026 10:36:09 -0700 Subject: [PATCH 1/9] Refactored rcpr version --- src/bca.rs | 25 ++++++++++++++++++++----- src/tests.rs | 2 ++ 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index eb9a8ab..fdf1941 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -2,7 +2,12 @@ use super::*; use rand::RngExt; #[cfg(feature = "cpr_rootfinder")] -use rcpr::chebyshev::*; +use rcpr::rootfinders::{ + real_polynomial_roots, + find_roots_with_newton_polishing, + find_roots_with_secant_polishing, + Config +}; /// Geometrical quantities of binary collision. pub struct BinaryCollisionGeometry { @@ -625,16 +630,26 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame let upper_bound = impact_parameter + interactions::crossing_point_doca(interaction_potential); let lower_bound = impact_parameter / 1000.0; + let delta = 1e-5; + let config = Config::new( + epsilon, + delta, + n0, + nmax, + complex_threshold, + truncation_threshold, + far_from_zero, + interval_limit + ); + 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), + config), 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) + config) } }.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))?; 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")] From 19f5d050f4d9a9c985fd1be43a3752b9ee27b44b Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 18 Aug 2026 10:37:02 -0700 Subject: [PATCH 2/9] wip cargo.toml --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index ee9637a..39bd31f 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 = { path="../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} From e62919dd84268946bd674c4e4f3e1804846f3d17 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 18 Aug 2026 10:52:43 -0700 Subject: [PATCH 3/9] I suspected root polishing was a waste of time for rustbca - and I think I'm right. Deprecated. --- src/bca.rs | 13 ++----------- 1 file changed, 2 insertions(+), 11 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index fdf1941..32058e6 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -3,6 +3,7 @@ use rand::RngExt; #[cfg(feature = "cpr_rootfinder")] use rcpr::rootfinders::{ + find_roots, real_polynomial_roots, find_roots_with_newton_polishing, find_roots_with_secant_polishing, @@ -642,17 +643,7 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame interval_limit ); - let roots = match derivative_free { - true => find_roots_with_secant_polishing(&g, &f, lower_bound, upper_bound, - config), - - 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, - config) - } - }.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 roots = find_roots(&g, vec![(lower_bound, upper_bound)], config)?; let max_root = roots.iter().cloned().fold(f64::NAN, f64::max)/a; From b6fcf1ade56eebce86f820b8672bf1b66684277d Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Tue, 18 Aug 2026 12:10:05 -0700 Subject: [PATCH 4/9] Morse screening length was in 1/A - should be 1/alpha, not alpha --- src/bca.rs | 1 - src/interactions.rs | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/bca.rs b/src/bca.rs index 32058e6..20b3b9a 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -624,7 +624,6 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame 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)}; diff --git a/src/interactions.rs b/src/interactions.rs index 5ba58ee..3d5488f 100644 --- a/src/interactions.rs +++ b/src/interactions.rs @@ -286,7 +286,7 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote //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::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), From 9f2cf35f5aa22b9fe80736b55d83a9d96b03125e Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 19 Aug 2026 15:45:14 -0700 Subject: [PATCH 5/9] Working version of Boyd's semi-infinite map - this removes the magic number on the upper bound and should significantly improve numerical stability of the cpr rootfinder. --- Cargo.toml | 2 +- src/bca.rs | 26 ++++++++++++++++++-------- src/interactions.rs | 14 ++++++++++---- 3 files changed, 29 insertions(+), 13 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 39bd31f..6110a90 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 = { path="../rcpr/", optional = true} +rcpr = { git = "https://github.com/drobnyjt/rcpr" , tag="0.3.0", 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/src/bca.rs b/src/bca.rs index 20b3b9a..1cbc9e8 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -5,8 +5,6 @@ use rand::RngExt; use rcpr::rootfinders::{ find_roots, real_polynomial_roots, - find_roots_with_newton_polishing, - find_roots_with_secant_polishing, Config }; @@ -593,6 +591,15 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact } } +const L: f64 = 8.0; +fn transform(x: f64) -> f64 { + L/(x*PI/2.).tan().powi(2) +} + +fn inverse_transform(x: f64) -> f64 { + 2./PI*((L/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. /// @@ -612,7 +619,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, @@ -624,11 +631,13 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame let relative_energy = E0*Mb/(Ma + Mb); let p = impact_parameter; - 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 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 = impact_parameter + interactions::crossing_point_doca(interaction_potential); - let lower_bound = impact_parameter / 1000.0; + let upper_bound = 1.0; + let lower_bound = 1e-5; let delta = 1e-5; let config = Config::new( @@ -644,7 +653,8 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame 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); 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 3d5488f..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,9 +282,9 @@ 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} => 1./alpha, InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => zbl_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. From f640ef50d5cc5cd96af5a8a6fcf8509aff3400be Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 19 Aug 2026 16:25:23 -0700 Subject: [PATCH 6/9] Update morse params. --- examples/test_morse.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/examples/test_morse.py b/examples/test_morse.py index 2716eb5..416c242 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -13,6 +13,9 @@ hydrogen['Ec'] = 0.1 hydrogen['Es'] = 1.5 +epsilon = 1e-3 +interval_limit = 1e-6 +nmax = 64 #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 +28,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=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 = 4 num_chunks = 10 @@ -87,7 +90,7 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): 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}}}}]] 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 @@ -195,8 +198,8 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): 0.46, 0.23, 0.05, 0.00, 0.00 ] -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, rtol=0.5) +np.testing.assert_allclose(R_N_2, R_N_2_test, rtol=0.5) 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') From 4e0ebea3c13f2da31fd2679858306398e60bf851 Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 19 Aug 2026 21:34:58 -0700 Subject: [PATCH 7/9] Updated test_morse with new values for new rcpr --- examples/test_morse.py | 37 +++++++++++++++++++------------------ src/bca.rs | 14 +++++++------- src/lib.rs | 1 + 3 files changed, 27 insertions(+), 25 deletions(-) diff --git a/examples/test_morse.py b/examples/test_morse.py index 416c242..0ada221 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -13,9 +13,10 @@ hydrogen['Ec'] = 0.1 hydrogen['Es'] = 1.5 -epsilon = 1e-3 -interval_limit = 1e-6 -nmax = 64 +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): @@ -28,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={nmax}, epsilon={epsilon}, complex_threshold=1E-9, truncation_threshold=1E-9, far_from_zero=1E22, interval_limit={interval_limit}, 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 @@ -86,9 +87,9 @@ 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=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 @@ -173,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 +run_sim = False +num_samples = 1000 R_N = np.zeros(num_energies) R_E = np.zeros(num_energies) R_N_2 = np.zeros(num_energies) @@ -187,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, rtol=0.5) -np.testing.assert_allclose(R_N_2, R_N_2_test, rtol=0.5) +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 1cbc9e8..3e057b1 100644 --- a/src/bca.rs +++ b/src/bca.rs @@ -591,13 +591,14 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact } } -const L: f64 = 8.0; +#[cfg(feature = "cpr_rootfinder")] fn transform(x: f64) -> f64 { - L/(x*PI/2.).tan().powi(2) + 2.0/(x*PI/2.).tan().powi(2) } +#[cfg(feature = "cpr_rootfinder")] fn inverse_transform(x: f64) -> f64 { - 2./PI*((L/x).sqrt()).atan() + 2./PI*((2.0/x).sqrt()).atan() } #[cfg(feature = "cpr_rootfinder")] @@ -627,9 +628,7 @@ 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 g = |r: f64| -> f64 { interactions::distance_of_closest_approach_function_singularity_free(transform(r)*a, a, Za, Zb, relative_energy, impact_parameter, interaction_potential)* @@ -637,7 +636,7 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame }; let upper_bound = 1.0; - let lower_bound = 1e-5; + let lower_bound = 1e-4; let delta = 1e-5; let config = Config::new( @@ -654,7 +653,8 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame let roots = find_roots(&g, vec![(lower_bound, upper_bound)], config)?; // 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)).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/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; From 29e63cfb547069dbcf6ea49de3005925fbe36c6f Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 19 Aug 2026 22:00:19 -0700 Subject: [PATCH 8/9] Left run_sim off in test_morse. fixed. --- examples/test_morse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/test_morse.py b/examples/test_morse.py index 0ada221..3f17c54 100644 --- a/examples/test_morse.py +++ b/examples/test_morse.py @@ -176,7 +176,7 @@ def run_krc_morse_potential(energy, index, num_samples=10000, run_sim=True): #Running and plotting the H-Ni simulations with the Morse potential and updated Es num_energies = 20 energies = np.logspace(-1, 4, num_energies) -run_sim = False +run_sim = True num_samples = 1000 R_N = np.zeros(num_energies) R_E = np.zeros(num_energies) From a6ec5337e86d8ab239dd5bbf0f377e75d812a44c Mon Sep 17 00:00:00 2001 From: Jon Drobny Date: Wed, 19 Aug 2026 23:00:16 -0700 Subject: [PATCH 9/9] Need to update cargo to use most recent rcpr with Schur fix. --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 6110a90..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" , tag="0.3.0", 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}