New rcpr

examples/f_kporha_et_al_benchmark.py

@@ -27,7 +27,7 @@ def input_file(ion, target, incident_energy, angle, number_ions=1000, pot="WW",
     else:
         pot_name = "MORSE"
 
-    cpr = {'CPR': {'n0': 2, 'nmax': 32, 'epsilon': 5e-4, 'complex_threshold': 1E-9, 'truncation_threshold': 1E-9, 'far_from_zero': 1e3, 'interval_limit': 1E-4, 'derivative_free': True}}
+    cpr = {'CPR': {'n0': 2, 'nmax': 32, 'epsilon': 5e-4, 'complex_threshold': 1E-9, 'far_from_zero': 1e3, 'interval_limit': 1E-4, 'derivative_free': True}}
     options = {
         'name': f'input_file_{ion["symbol"]}_{target["symbol"]}_{np.round(angle, 1)}_{np.round(incident_energy/1000, 4)}_{pot_name}_{Es}',
         'track_trajectories': False, # whether to track trajectories for plotting; memory intensive
@@ -190,7 +190,7 @@ def tungsten_sputtering(ion, energy, angle, num_ions=1000, run_sim=False, pot="W
 colors = []
 
 num_ions = 10000
-run_sim = False
+run_sim = True
 
 tungsten["Es"] = 8.79
 tungsten["Eb"] = 0.0

examples/test_morse.py

@@ -29,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={n0}, 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, far_from_zero=1E22, interval_limit={interval_limit}, derivative_free=true}}}}]]
     num_threads = 4
     num_chunks = 10
 
@@ -91,7 +91,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=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}}}}]]
+    root_finder = [[{{"CPR"={{n0=3, nmax={nmax}, epsilon={epsilon}, complex_threshold=1E-9, far_from_zero=1E22, interval_limit={interval_limit}, derivative_free=true}}}}]]
     num_threads = 6
     num_chunks = 1
 

examples/xenon_sputtering_benchmark.py

@@ -14,7 +14,7 @@ def input_file(ion, target, incident_energy, angle, number_ions=1000):
 
     mfp = (target["n"]/10**30)**(-1./3.)
 
-    cpr = {'CPR': {'n0': 2, 'nmax': 32, 'epsilon': 1e-3, 'complex_threshold': 1E-9, 'truncation_threshold': 1E-12, 'far_from_zero': 1e3, 'interval_limit': 1E-3, 'derivative_free': True}}
+    cpr = {'CPR': {'n0': 2, 'nmax': 32, 'epsilon': 1e-3, 'complex_threshold': 1E-9, 'far_from_zero': 1e3, 'interval_limit': 1E-3, 'derivative_free': True}}
     options = {
         'name': f'input_file_{ion["symbol"]}_{target["symbol"]}_{np.round(angle, 1)}_{np.round(incident_energy/1000, 4)}',
         'track_trajectories': False, # whether to track trajectories for plotting; memory intensive
@@ -196,7 +196,7 @@ def input_file(ion, target, incident_energy, angle, number_ions=1000):
 
 num_ions = 10000
 num_energies = 25
-run_sim = True
+run_sim = False
 energies = np.logspace(np.log10(25), np.log10(1600), num_energies)
 angle = 0.0
 
@@ -231,8 +231,57 @@ def input_file(ion, target, incident_energy, angle, number_ions=1000):
 
 plt.plot(energies, Y_Xe_Mo, label='RustBCA Default')
 plt.gca().set_xscale('log')
+plt.gca().set_yscale('log')
 plt.legend()
 plt.xlabel('E [eV]')
 plt.ylabel('Y [at/ion]')
 plt.title('Xe on Mo Sputtering Yields')
+
+ # Xe on Ti
+
+data = np.array([
+[0.06339144215530546, 0.1428571428571308],
+[29.920760697305866, 0.39030612244896945],
+[44.88114104595877, 0.5459183673469292],
+[59.968304278922346, 0.6530612244897864],
+[69.79397781299527, 0.5586734693877453],
+[0.1267828843106198, 0.33928571428570375],
+[30.047543581616488, 0.7755102040816233],
+[44.8811410459588, 1.0714285714285632],
+[59.96830427892236, 1.2602040816326459],
+[69.79397781299525, 1.1938775510204005],
+[0.25356576862123426, 0.4846938775510101],
+[29.920760697305866, 1.030612244897951],
+[45.00792393026942, 1.4285714285714222],
+[59.96830427892236, 1.6173469387755044],
+[69.98415213946117, 1.6632653061224434],
+[0.06339144215530546, 0.8061224489795826],
+[29.984152139461166, 1.5255102040816264],
+[45.00792393026942, 2.056122448979588],
+[59.96830427892236, 2.4693877551020385],
+[69.85736925515056, 2.6377551020408143],
+[0.12678288431062157, 1.0076530612244818],
+[29.98415213946116, 1.834183673469383],
+[44.94453248811412, 2.461734693877548],
+[59.96830427892236, 3.229591836734694],
+[69.92076069730587, 3.357142857142858],
+])
+
+run_sim = True
+num_angles = 15
+num_ions = 1000
+angles = np.linspace(0, 89, num_angles)
+energies = [200, 400, 600, 1000, 1400]
+titanium['Eb'] = 3.0
+
+plt.figure()
+for energy in energies:
+    Y_Xe_Ti = np.zeros(num_angles)
+    for index, angle in enumerate(angles):
+            input_data = input_file(xenon, titanium, energy, angle, num_ions)
+            if run_sim: rustbca_py(input_data, geometry_mode="0D")
+            sputtered = np.genfromtxt(f'{input_data["options"]["name"]}sputtered.output', delimiter=',')
+            Y_Xe_Ti[index] = np.shape(sputtered)[0]/num_ions
+    plt.plot(angles, Y_Xe_Ti, label=f'{energy} eV')
+plt.scatter(data[:, 0], data[:, 1])
 plt.show()

src/bca.rs

@@ -412,8 +412,8 @@ fn distance_of_closest_approach(particle_1: &particle::Particle, particle_2: &pa
         Rootfinder::POLYNOMIAL{complex_threshold} => polynomial_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, complex_threshold)
             .with_context(|| format!("Numerical error: Polynomial rootfinder failed for {} at {} eV with p = {} A.", interaction_potential, E0/EV, p/ANGSTROM))
             .unwrap(),
-        Rootfinder::CPR{n0, nmax, epsilon, complex_threshold, truncation_threshold, far_from_zero, interval_limit, derivative_free} =>
-            cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, n0, nmax, epsilon, complex_threshold, truncation_threshold, far_from_zero, interval_limit, derivative_free)
+        Rootfinder::CPR{n0, nmax, epsilon, complex_threshold, far_from_zero, interval_limit, derivative_free} =>
+            cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, n0, nmax, epsilon, complex_threshold, far_from_zero, interval_limit, derivative_free)
             .with_context(|| format!("Numerical error: CPR rootfinder failed for {} at {} eV with p = {} A.", interaction_potential, E0/EV, p/ANGSTROM))
             .unwrap(),
         Rootfinder::NEWTON{max_iterations, tolerance} => newton_rootfinder(Za, Zb, Ma, Mb, E0, p, interaction_potential, max_iterations, tolerance)
@@ -620,15 +620,14 @@ fn transform(x: f64) -> f64 {
 /// `epsilon`: absolute tolerance of Chebyshev interpolant.
 /// `complex_threshold`: slightly-complex roots with an imaginary part below this value are considered real.
 /// `far_from_zero`: if the distance of closest approach function, evaluated over an interval [a, b] on the Lobatto grid, is always greater than this value, it is assumed that there are no roots in the interval [a, b].
-/// `truncation_threshold`: trailing terms of the Chebyshev interpolant with coefficients smaller than this value are ignored.
 /// `interval_limit`: if subdivision produces an interval smaller than this value, the root-finder will panic.
 /// `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 (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,
+    complex_threshold: f64, far_from_zero: f64,
     interval_limit: f64, derivative_free: bool) -> Result <f64, anyhow::Error> {
 
     //Lindhard screening length and reduced energy
@@ -650,7 +649,6 @@ pub fn cpr_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact_parame
         n0,
         nmax,
         complex_threshold,
-        truncation_threshold,
         far_from_zero,
         interval_limit
     );

src/enums.rs

@@ -241,7 +241,7 @@ impl PartialEq for ScatteringIntegral {
 pub enum Rootfinder {
     /// Newton root-finder with user-defined `max_iterations` and `tolerance`.
     NEWTON{max_iterations: usize, tolerance: f64},
-    CPR{n0: usize, nmax: usize, epsilon: f64, complex_threshold: f64, truncation_threshold: f64,
+    CPR{n0: usize, nmax: usize, epsilon: f64, complex_threshold: f64,
         far_from_zero: f64, interval_limit: f64, derivative_free: bool},
     POLYNOMIAL{complex_threshold: f64},
     DEFAULTNEWTON,
@@ -251,7 +251,7 @@ impl fmt::Display for Rootfinder {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         match *self {
             Rootfinder::NEWTON{max_iterations, tolerance} => write!(f, "Newton-Raphson Rootfinder with maximum {} iterations and toleance = {}", max_iterations, tolerance),
-            Rootfinder::CPR{n0, nmax, epsilon, complex_threshold, truncation_threshold, far_from_zero, interval_limit, derivative_free} =>
+            Rootfinder::CPR{n0, nmax, epsilon, complex_threshold, far_from_zero, interval_limit, derivative_free} =>
                 write!(f, "Chebyshev-Proxy Rootfinder with {}-polishing", match derivative_free { true => "Secant", false => "Newton"}),
             Rootfinder::POLYNOMIAL{complex_threshold} => write!(f, "Frobenius Companion Matrix Polynomial Real Rootfinder with a complex tolerance of {}", complex_threshold),
             Rootfinder::DEFAULTNEWTON => write!(f, "Newton-Raphson Rootfinder with maximum {} iterations and toleance = {}", 100, 1E-3),

src/lib.rs

@@ -145,6 +145,9 @@ mod libRustBCA {
 
     #[pymodule_export]
     use super::rustbca_local_py;
+
+    #[pymodule_export]
+    use super::scattering_integral_morse;
 }
 
 #[cfg(feature = "python")]
@@ -2129,6 +2132,7 @@ pub fn compound_reflection_coefficient<'py>(ion: &Bound<'py, PyDict>, targets: V
 /// E. Kadric et al., IEEE Transactions on Computers 65 11
 /// doi: 10.1109/TC.2016.2532874
 #[cfg(feature = "python")]
+#[allow(dead_code)]
 fn moller_knuth_two_sum(a: f64, b: f64) -> (f64, f64) {
     let s = a + b;
     let b_prime = s - a;
@@ -2169,6 +2173,36 @@ fn scattering_integrals(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, n_g
     Ok((theta_gm, theta_gl, theta_mw, theta_magic))
 }
 
+#[cfg(feature = "python")]
+#[pyfunction]
+#[pyo3(signature = (Za, Zb, Ma, Mb, E0, p, d, alpha, r0, n0=2, nmax=64, epsilon=1e-6, complex_threshold=0.0, far_from_zero=1e22, interval_limit=1e-3, n_gl_points=100, interaction_potential="KR_C"))]
+fn scattering_integral_morse(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, p: f64, d: f64, alpha: f64, r0: f64, n0: usize, nmax: usize, epsilon: f64, complex_threshold: f64, far_from_zero: f64, interval_limit: f64, n_gl_points: usize, interaction_potential: &str) -> PyResult<(f64, f64, f64)> {
+    let E0 = E0*EV;
+    let p = p*ANGSTROM;
+    let d = d*EV;
+    let alpha = alpha/ANGSTROM;
+    let r0 = r0*ANGSTROM;
+
+    let potential = InteractionPotential::MORSE{D: d, alpha: alpha, r0: r0};
+    let screening_length = interactions::screening_length(Za, Zb, potential);
+
+    let x0 = bca::cpr_rootfinder(
+        Za, Zb, Ma, Mb, E0, p,
+        potential, n0, nmax, epsilon,
+        complex_threshold,
+        far_from_zero, interval_limit, true
+    ).map_err(|e| PyValueError::new_err(format!(
+        "Failed to find root: {}; p={}, E={}", e, p/ANGSTROM, E0/EV
+    )))?;
+
+    //Compute center of mass deflection angle with each algorithm
+    let theta_gm = bca::gauss_mehler(Za, Zb, Ma, Mb, E0, p, x0, screening_length, potential, n_gl_points);
+    let theta_gl = bca::gauss_legendre(Za, Zb, Ma, Mb, E0, p, x0, screening_length, potential);
+
+
+    Ok((x0, theta_gm, theta_gl))
+}
+
 #[cfg(feature = "python")]
 #[pyfunction]
 #[pyo3(signature=(input, geometry_mode="1D"))]

src/tests.rs

@@ -788,15 +788,13 @@ fn test_momentum_conservation() {
     #[cfg(not(feature = "cpr_rootfinder"))]
     let  rootfinders = vec![Rootfinder::NEWTON{max_iterations: 100, tolerance: 1E-3}; 4];
 
-    //[[{"CPR"={n0=2, nmax=100, epsilon=1E-9, complex_threshold=1E-3, truncation_threshold=1E-9, far_from_zero=1E9, interval_limit=1E-12, derivative_free=true}}]]
     #[cfg(feature = "cpr_rootfinder")]
     rootfinders.push(
         Rootfinder::CPR{
             n0: 2,
             nmax: 100,
             epsilon: 1e-9,
             complex_threshold: 1e-3,
-            truncation_threshold: 1e-9,
             far_from_zero: 1e9,
             interval_limit:1e-12,
             derivative_free: true
@@ -1171,7 +1169,7 @@ fn test_quadrature() {
 
     //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, interaction_potential, 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, true) {
         println!("CPR: {} Newton: {}", x0_cpr, x0_newton);
         assert!(approx_eq!(f64, x0_newton, x0_cpr, epsilon=1E-3));
     };