Merge pull request #233 from lcpp-org/4_8_potential

4 8 potential

Author: drobnyjt

scripts/materials.py

@@ -12,6 +12,29 @@
     's': 2.5
 }
 
+carbon = {
+    'symbol': 'C',
+    'name': 'carbon',
+    'Z': 6,
+    'm': 12.011,
+    'n': 1.1331E29,
+    'Es': 7.37,
+    'Ec': 1.0,
+    'Eb': 3.0,
+}
+
+chromium = {
+ 'symbol': 'Cr',
+ 'name': 'chromium',
+ 'Z': 24,
+ 'm': 51.9961,
+ 'Es': 4.10,
+ 'Ed': 28.0,
+ 'Ec': 3.0,
+ 'Eb': 3.0,
+ 'n': 8.327E28
+}
+
 hydrogen = {
     'symbol': 'H',
     'name': 'hydrogen',

src/bca.rs

@@ -576,7 +576,21 @@ pub fn polynomial_rootfinder(Za: f64, Zb: f64, Ma: f64, Mb: f64, E0: f64, impact
 
     let coefficients = interactions::polynomial_coefficients(relative_energy, impact_parameter, interaction_potential);
     let roots = real_polynomial_roots(coefficients.clone(), polynom_complex_threshold).unwrap();
+
+    /*
+    println!("p={}", impact_parameter/ANGSTROM);
+
+    for coefficient in &coefficients {
+        println!("{}", {coefficient});
+    }
+
+    for root in &roots {
+        println!("{} A", {root});
+    }
+    */
+
     let max_root = roots.iter().cloned().fold(f64::NAN, f64::max);
+
     let inverse_transformed_root = interactions::inverse_transform(max_root, interaction_potential);
 
     if roots.is_empty() || inverse_transformed_root.is_nan() {

src/enums.rs

@@ -173,7 +173,9 @@ pub enum InteractionPotential {
     /// Unscreened Coulombic interatomic potential between ions with charges Za and Zb.
     COULOMB{Za: f64, Zb: f64},
     /// Morse potential that smoothly interpolates to Kr-C at short ranges.
-    KRC_MORSE{D: f64, alpha: f64, r0: f64, k: f64, x0: f64}
+    KRC_MORSE{D: f64, alpha: f64, r0: f64, k: f64, x0: f64},
+    /// Born repulsive potential.
+    FOUR_EIGHT{alpha: f64, beta: f64}
 }
 
 impl fmt::Display for InteractionPotential {
@@ -190,6 +192,7 @@ impl fmt::Display for InteractionPotential {
             InteractionPotential::WW => write!(f, "W-W cubic spline interaction potential."),
             InteractionPotential::COULOMB{Za, Zb} => write!(f, "Coulombic interaction with Za = {} and Zb = {}", Za, Zb),
             InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => write!(f, "Morse potential with D = {} eV, alpha = {} 1/A, and r0 = {} A with short-range Kr-C with interpolation width k = {} 1/A and transition point x0 = {} A", D/EV, alpha*ANGSTROM, r0/ANGSTROM, k*ANGSTROM, x0/ANGSTROM),
+            InteractionPotential::FOUR_EIGHT{alpha, beta} => write!(f, "Inverse polynomial potential with beta/r^8 (beta = {} A) Born repulsive term and alpha/r^4 (alpha = {} A) attractive term. Used to describe hydrogen interactions with transition metals.", alpha/ANGSTROM, beta/ANGSTROM),
         }
     }
 }

src/input.rs

@@ -450,6 +450,7 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
 
             assert!(
                 match (interaction_potential, root_finder) {
+                    (InteractionPotential::FOUR_EIGHT{..}, Rootfinder::POLYNOMIAL{..}) => true,
                     (InteractionPotential::LENNARD_JONES_12_6{..}, Rootfinder::CPR{..}) => true,
                     (InteractionPotential::LENNARD_JONES_12_6{..}, Rootfinder::POLYNOMIAL{..}) => true,
                     (InteractionPotential::LENNARD_JONES_12_6{..}, _) => false,

src/interactions.rs

@@ -40,6 +40,10 @@ pub fn interaction_potential(r: f64, a: f64, Za: f64, Zb: f64, interaction_poten
         InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => {
             krc_morse(r, a, Za, Zb, D, alpha, r0, k, x0)
         }
+        InteractionPotential::FOUR_EIGHT{alpha, beta} => {
+            four_eight(r, alpha, beta)
+        }
+
     }
 }
 
@@ -48,6 +52,7 @@ pub fn energy_threshold_single_root(interaction_potential: InteractionPotential)
     match interaction_potential{
         InteractionPotential::LENNARD_JONES_12_6{..} | InteractionPotential::LENNARD_JONES_65_6{..} => f64::INFINITY,
         InteractionPotential::MORSE{..} | InteractionPotential::KRC_MORSE{..} => f64::INFINITY,
+        InteractionPotential::FOUR_EIGHT{..} => f64::INFINITY,
         InteractionPotential::WW => f64::INFINITY,
         InteractionPotential::COULOMB{..} => f64::INFINITY,
         InteractionPotential::MOLIERE | InteractionPotential::KR_C | InteractionPotential::LENZ_JENSEN | InteractionPotential::ZBL | InteractionPotential::TRIDYN => 0.,
@@ -103,6 +108,9 @@ pub fn distance_of_closest_approach_function(r: f64, a: f64, Za: f64, Zb: f64, r
         InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => {
             doca_krc_morse(r, impact_parameter, relative_energy, a, Za, Zb, D, alpha, r0, k, x0)
         },
+        InteractionPotential::FOUR_EIGHT{alpha, beta} => {
+            doca_four_eight(r, impact_parameter, relative_energy, alpha, beta)
+        }
         InteractionPotential::COULOMB{..} => panic!("Coulombic potential cannot be used with rootfinder.")
     }
 }
@@ -131,6 +139,9 @@ pub fn distance_of_closest_approach_function_singularity_free(r: f64, a: f64, Za
         InteractionPotential::KRC_MORSE{D, alpha, r0, k, x0} => {
             doca_krc_morse(r, impact_parameter, relative_energy, a, Za, Zb, D, alpha, r0, k, x0)
         },
+        InteractionPotential::FOUR_EIGHT{alpha, beta} => {
+            doca_four_eight(r, impact_parameter, relative_energy, alpha, beta)
+        }
         InteractionPotential::WW => {
             doca_tungsten_tungsten_cubic_spline(r, impact_parameter, relative_energy)
         },
@@ -152,6 +163,10 @@ pub fn scaling_function(r: f64, a: f64, interaction_potential: InteractionPotent
             let n = 6.;
             1./(1. + (r/sigma).powf(n))
         },
+        InteractionPotential::FOUR_EIGHT{alpha, beta} => {
+            let n = 8.;
+            1./(1. + r.powi(8)/beta)
+        }
         InteractionPotential::MORSE{D, alpha, r0} => {
             1./(1. + (r*alpha).powi(2))
         }
@@ -256,35 +271,65 @@ pub fn screening_length(Za: f64, Zb: f64, interaction_potential: InteractionPote
         InteractionPotential::MORSE{D, alpha, r0} => alpha,
         InteractionPotential::COULOMB{Za: Z1, Zb: Z2} => 0.88534*A0/(Z1.powf(0.23) + Z2.powf(0.23)),
         InteractionPotential::KRC_MORSE{..} => 0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.),
+        InteractionPotential::FOUR_EIGHT{..} =>0.8853*A0*(Za.sqrt() + Zb.sqrt()).powf(-2./3.),
     }
 }
 
 /// Coefficients of inverse-polynomial interaction potentials.
+/// Note: these go in order of decreasing orders of r, e.g., [r^n, r^n-1, r^n-2, ... r, 1].
 pub fn polynomial_coefficients(relative_energy: f64, impact_parameter: f64, interaction_potential: InteractionPotential) -> Vec<f64> {
     match interaction_potential {
         InteractionPotential::LENNARD_JONES_12_6{sigma, epsilon} => {
-            vec![1., 0., -impact_parameter.powi(2), 0., 0., 0., 4.*epsilon*sigma.powf(6.)/relative_energy, 0., 0., 0., 0., 0., -4.*epsilon*sigma.powf(12.)/relative_energy]
+            let impact_parameteter_angstroms = impact_parameter/ANGSTROM;
+            let epsilon_ev = epsilon/EV;
+            let sigma_angstroms = sigma/ANGSTROM;
+            let relative_energy_ev = relative_energy/EV;
+            vec![1., 0., -impact_parameteter_angstroms.powi(2), 0., 0., 0., 4.*epsilon_ev*sigma_angstroms.powf(6.)/relative_energy_ev, 0., 0., 0., 0., 0., -4.*epsilon_ev*sigma_angstroms.powf(12.)/relative_energy_ev]
         },
         InteractionPotential::LENNARD_JONES_65_6{sigma, epsilon} => {
-            vec![1., 0., 0., 0., -impact_parameter.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon*sigma.powf(6.)/relative_energy, -4.*epsilon*sigma.powf(6.5)/relative_energy]
-        },
+            let impact_parameteter_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_parameteter_angstroms.powi(2), 0., 0., 0., 0., 0., 0., 0., 4.*epsilon_ev*sigma_angstroms.powf(6.)/relative_energy, -4.*epsilon_ev*sigma_angstroms.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.
+            //The units on alpha [m^4 J] and beta [m^8 J] make them unreasonably small to use numerically.
+            let alpha_angstroms4_joule = alpha/ANGSTROM.powi(4);
+            let beta_angstroms8_joule = beta/ANGSTROM.powi(8);
+            let impact_parameteter_angstroms = impact_parameter/ANGSTROM;
+            vec![1., 0., -(impact_parameteter_angstroms).powi(2), 0., alpha_angstroms4_joule/relative_energy, 0., 0., 0., -beta_angstroms8_joule/relative_energy,]
+        }
         _ => panic!("Input error: non-polynomial interaction potential used with polynomial root-finder.")
     }
 }
 
-/// Inverse-polynomial interaction potentials transformation to remove singularities for the CPR root-finder.
+/// Inverse-polynomial interaction potentials transformation to reduce numercial issues for the polynomial rootfinder.
+/// This is for if, for example, you have non-integer powers and need to multiply by a factor to get to integer powers.
+/// E.g., the 6.5-6 potential needs to be squared (x*x) giving you a 13th order polynomial in polynomial_coefficients().
 pub fn inverse_transform(x: f64, interaction_potential: InteractionPotential) -> f64 {
     match interaction_potential {
         InteractionPotential::LENNARD_JONES_12_6{..} => {
-            x
+            x*ANGSTROM
         },
         InteractionPotential::LENNARD_JONES_65_6{..} => {
             x*x
         },
+        | InteractionPotential::FOUR_EIGHT{..} => {
+            x*ANGSTROM
+        }
         _ => panic!("Input error: non-polynomial interaction potential used with polynomial root-finder transformation.")
     }
 }
 
+/// Four-Eight potential (Born repulsion)
+pub fn four_eight(r: f64, alpha: f64, beta: f64) -> f64 {
+    let a = alpha.powf(1./4.);
+    let b = beta.powf(1./8.);
+    -(a/r).powi(4) + (b/r).powi(8)
+}
+
 /// 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.))
@@ -305,6 +350,14 @@ pub fn krc_morse(r: f64, a: f64, Za: f64, Zb: f64, D: f64, alpha: f64, r0: f64,
     sigmoid(r, k, x0)*morse(r, D, alpha, r0) + sigmoid(r, -k, x0)*screened_coulomb(r, a, Za, Zb, InteractionPotential::KR_C)
 }
 
+/// Distance of closest approach function for four-eight potential (Born repulsion)
+pub fn doca_four_eight(r: f64, impact_parameter: f64, relative_energy: f64, alpha: f64, beta: f64) -> f64 {
+    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)
+}
+
 /// Distance of closest approach function for Morse potential.
 pub fn doca_morse(r: f64, impact_parameter: f64, relative_energy: f64, D: f64, alpha: f64, r0: f64) -> f64 {
     (r*alpha).powi(2) - (r*alpha).powi(2)*D/relative_energy*((-2.*alpha*(r - r0)).exp() - 2.*(-alpha*(r - r0)).exp()) - (impact_parameter*alpha).powi(2)
@@ -463,5 +516,6 @@ pub fn first_screening_radius(interaction_potential: InteractionPotential) -> f6
         InteractionPotential::MORSE{..} | InteractionPotential::KRC_MORSE{..} => 1.,
         InteractionPotential::WW => 1.,
         InteractionPotential::COULOMB{..} => 0.3,
+        InteractionPotential::FOUR_EIGHT{..} => 1.,
     }
 }

src/tests.rs

@@ -3,6 +3,21 @@ use super::*;
 #[cfg(test)]
 use float_cmp::*;
 
+
+#[test]
+#[cfg(feature = "cpr_rootfinder")]
+fn test_polynom() {
+    use rcpr::chebyshev::*;
+    let interaction_potential = InteractionPotential::FOUR_EIGHT{alpha: 1., beta: 1.};
+    let coefficients = interactions::polynomial_coefficients(1., 1., interaction_potential);
+    println!("{} {} {}", coefficients[0], coefficients[4], coefficients[8]);
+    let roots = real_polynomial_roots(coefficients.clone(), 1e-9).unwrap();
+
+    let max_root = roots.iter().cloned().fold(f64::NAN, f64::max);
+    println!("{}", max_root);
+    let inverse_transformed_root = interactions::inverse_transform(max_root, interaction_potential);
+}
+
 #[test]
 #[cfg(feature = "parry3d")]
 fn test_parry_cuboid() {
@@ -1072,7 +1087,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, InteractionPotential::KR_C, 2, 10000, 1E-6, 1E-6, 1E-9, 1E9, 1E-13, false) {
+    if let Ok(x0_cpr) = bca::cpr_rootfinder(Za, Zb, Ma, Mb, E0, p, InteractionPotential::KR_C, 2, 10000, 1E-6, 1E-6, 1E-9, 1E9, 1E-13, true) {
         println!("CPR: {} Newton: {}", x0_cpr, x0_newton);
         assert!(approx_eq!(f64, x0_newton, x0_cpr, epsilon=1E-3));
     };
@@ -1090,4 +1105,4 @@ fn test_quadrature() {
 
     println!("Gauss-Mehler: {} Gauss-Legendre: {} Mendenhall-Weller: {} MAGIC: {}",
         theta_gm, theta_gl, theta_mw, theta_magic);
-}
+}