sci_form/
potentials.rs

1//! Analytical potentials for reaction path benchmarking and testing.
2
3/// 1D double well potential: V(x) = (x^2 - 1)^2
4/// Symmetric, minima at +/- 1, barrier at x=0
5pub fn double_well_1d(x: f64) -> f64 {
6    (x * x - 1.0).powi(2)
7}
8
9/// 1D asymmetric potential: V(x) = (x^2 - 1)^2 - 0.3x
10/// Exothermic with product at x=+1 being more stable
11pub fn asymmetric_1d(x: f64) -> f64 {
12    (x * x - 1.0).powi(2) - 0.3 * x
13}
14
15/// 1D SN2 model potential: V(x) = -2*exp(-2(x+1.5)^2) - 3*exp(-2(x-1.5)^2) + 1.5*exp(-8x^2)
16pub fn sn2_model_1d(x: f64) -> f64 {
17    -2.0 * (-(2.0 * (x + 1.5).powi(2))).exp() - 3.0 * (-(2.0 * (x - 1.5).powi(2))).exp()
18        + 1.5 * (-(8.0 * x.powi(2))).exp()
19}
20
21/// 2D Müller-Brown potential
22pub fn muller_brown_2d(x: f64, y: f64) -> f64 {
23    let a = [-200.0_f64, -100., -170., 15.];
24    let b = [-1.0, -1., -6.5, 0.7];
25    let c = [0.0, 0., 11.0, 0.6];
26    let d = [-10.0, -10., -6.5, 0.7];
27    let x0 = [1.0, 0., -0.5, -1.0];
28    let y0 = [0.0, 0.5, 1.5, 1.0];
29    (0..4)
30        .map(|k| {
31            let dx = x - x0[k];
32            let dy = y - y0[k];
33            a[k] * (b[k] * dx * dx + c[k] * dx * dy + d[k] * dy * dy).exp()
34        })
35        .sum::<f64>()
36        * 0.001
37}
38
39/// 1D Marcus/EVB adiabatic potential for proton transfer
40/// s in [-1, 1], symmetric TS exactly at s=0
41pub fn proton_transfer_evb(s: f64, k: f64, coupling: f64) -> f64 {
42    let eps_d = k * (s + 1.0).powi(2);
43    let eps_a = k * (s - 1.0).powi(2);
44    let sum2 = (eps_d + eps_a) / 2.0;
45    let diff2 = (eps_d - eps_a) / 2.0;
46    sum2 - (diff2 * diff2 + coupling * coupling).sqrt()
47}
48
49/// 2D Morse-based potential for proton transfer A-H...B
50pub fn proton_transfer_2d_morse(r_ah: f64, r_hb: f64) -> f64 {
51    let de = 2.0;
52    let alpha = 2.0;
53    let r0 = 1.0;
54    let morse_ah = de * (1.0 - (-alpha * (r_ah - r0)).exp()).powi(2);
55    let morse_hb = de * (1.0 - (-alpha * (r_hb - r0)).exp()).powi(2);
56    let constraint = 50.0 * (r_ah + r_hb - 2.5).powi(2);
57    morse_ah + morse_hb + constraint
58}
sci_form/potentials.rs

sci_form/
potentials.rs
