feat: Enhance NonbondedInteraction with OpenMM reference tests
- Update COULOMB_CONSTANT for accuracy based on 2019 SI constants - Add nalgebra as a development dependency - Introduce comprehensive tests for NonbondedInteraction against OpenMM NonbondedForce - Create reference values generation script for validation
This commit is contained in:
@@ -18,3 +18,4 @@ rayon = { workspace = true, optional = true }
|
||||
|
||||
[dev-dependencies]
|
||||
approx = { workspace = true }
|
||||
nalgebra = { workspace = true }
|
||||
|
||||
@@ -43,7 +43,8 @@ pub struct NonbondedInteraction<const D: usize> {
|
||||
|
||||
impl<const D: usize> NonbondedInteraction<D> {
|
||||
/// Coulomb constant in kJ·nm/(mol·e²)
|
||||
const COULOMB_CONSTANT: f64 = 138.935456;
|
||||
/// Derived from: 1/(4πε₀) with 2019 SI constants (matches OpenMM)
|
||||
const COULOMB_CONSTANT: f64 = 138.93545764438198;
|
||||
|
||||
/// Create a new nonbonded interaction
|
||||
pub fn new() -> Self {
|
||||
@@ -122,12 +123,6 @@ impl<const D: usize> ForceInteraction<D> for NonbondedInteraction<D> {
|
||||
forces: &mut [SVector<f64, D>],
|
||||
) -> ForceOutput<D> {
|
||||
let n = positions.len();
|
||||
|
||||
// Zero forces
|
||||
for f in forces.iter_mut() {
|
||||
*f = SVector::zeros();
|
||||
}
|
||||
|
||||
let mut total_energy = 0.0;
|
||||
|
||||
// Cutoff squared for distance check
|
||||
|
||||
449
simul-euclidean/tests/nonbonded_openmm.rs
Normal file
449
simul-euclidean/tests/nonbonded_openmm.rs
Normal file
@@ -0,0 +1,449 @@
|
||||
//! Equivalence tests: NonbondedInteraction vs OpenMM NonbondedForce
|
||||
//!
|
||||
//! Reference values generated by: tests/reference/openmm_nonbonded.py
|
||||
//! OpenMM version: 8.4.0, NonbondedMethod: NoCutoff
|
||||
//!
|
||||
//! OpenMM formulae:
|
||||
//! U_LJ = 4ε[(σ/r)^12 - (σ/r)^6]
|
||||
//! U_Coulomb = k_e * q_i * q_j / r (k_e = 138.935456 kJ·nm/(mol·e²))
|
||||
//!
|
||||
//! Combining rules (Lorentz-Berthelot):
|
||||
//! σ_ij = (σ_i + σ_j) / 2
|
||||
//! ε_ij = sqrt(ε_i * ε_j)
|
||||
|
||||
use nalgebra::SVector;
|
||||
use simul_euclidean::interaction::{ForceInteraction, NonbondedInteraction};
|
||||
|
||||
type Vec3 = SVector<f64, 3>;
|
||||
|
||||
const TOL_ENERGY: f64 = 1e-10;
|
||||
const TOL_FORCE: f64 = 1e-6;
|
||||
|
||||
fn assert_energy_eq(actual: f64, expected: f64, case: &str) {
|
||||
let diff = (actual - expected).abs();
|
||||
let scale = expected.abs().max(1.0);
|
||||
assert!(
|
||||
diff / scale < TOL_ENERGY,
|
||||
"[{case}] energy mismatch: got {actual:.15e}, expected {expected:.15e}, \
|
||||
diff {diff:.2e}, rel {:.2e}",
|
||||
diff / scale,
|
||||
);
|
||||
}
|
||||
|
||||
fn assert_force_eq(actual: &SVector<f64, 3>, expected: &[f64; 3], particle: usize, case: &str) {
|
||||
for d in 0..3 {
|
||||
let diff = (actual[d] - expected[d]).abs();
|
||||
let scale = expected[d].abs().max(1.0);
|
||||
assert!(
|
||||
diff / scale < TOL_FORCE,
|
||||
"[{case}] force mismatch particle {particle} dim {d}: \
|
||||
got {:.15e}, expected {:.15e}, diff {diff:.2e}, rel {:.2e}",
|
||||
actual[d],
|
||||
expected[d],
|
||||
diff / scale,
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
fn make_nb(params: &[(f64, f64, f64)]) -> NonbondedInteraction<3> {
|
||||
let mut nb = NonbondedInteraction::<3>::new();
|
||||
for &(sigma, epsilon, charge) in params {
|
||||
nb.add_particle(sigma, epsilon, charge);
|
||||
}
|
||||
nb
|
||||
}
|
||||
|
||||
fn run(nb: &NonbondedInteraction<3>, positions: &[Vec3]) -> (f64, Vec<Vec3>) {
|
||||
let mut forces = vec![Vec3::zeros(); positions.len()];
|
||||
let output = nb.compute_forces(positions, &mut forces);
|
||||
(output.potential_energy, forces)
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 1: LJ only, 2 identical particles along x-axis
|
||||
// σ=0.3 nm, ε=1.0 kJ/mol, q=0, r=0.4 nm (attractive region)
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_two_identical() {
|
||||
let case = "lj_two_identical";
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.0), (0.3, 1.0, 0.0)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.4, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -5.852086544036865e-01, case);
|
||||
assert_force_eq(&forces[0], &[6.877548694610596e+00, 0.0, 0.0], 0, case);
|
||||
assert_force_eq(&forces[1], &[-6.877548694610596e+00, 0.0, 0.0], 1, case);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 2: LJ only, 2 particles at equilibrium (r = 2^(1/6)*σ)
|
||||
// Force ≈ 0, energy = -ε
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_equilibrium() {
|
||||
let case = "lj_equilibrium";
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.0), (0.3, 1.0, 0.0)]);
|
||||
let r_eq = 0.3 * 2.0_f64.powf(1.0 / 6.0);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(r_eq, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -1.000000000000000e+00, case);
|
||||
assert!(
|
||||
forces[0].norm() < 1e-10,
|
||||
"[{case}] force at equilibrium should be ~zero, got {:.2e}",
|
||||
forces[0].norm()
|
||||
);
|
||||
assert!(
|
||||
forces[1].norm() < 1e-10,
|
||||
"[{case}] force at equilibrium should be ~zero, got {:.2e}",
|
||||
forces[1].norm()
|
||||
);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 3: LJ only, 2 particles with different σ/ε
|
||||
// σ_0=0.25, ε_0=0.5, σ_1=0.35, ε_1=2.0
|
||||
// Lorentz-Berthelot: σ_ij=0.3, ε_ij=1.0
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_different_params() {
|
||||
let case = "lj_different_params";
|
||||
let nb = make_nb(&[(0.25, 0.5, 0.0), (0.35, 2.0, 0.0)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.35, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -9.572084907712046e-01, case);
|
||||
assert_force_eq(&forces[0], &[5.625242659312204e+00, 0.0, 0.0], 0, case);
|
||||
assert_force_eq(&forces[1], &[-5.625242659312204e+00, 0.0, 0.0], 1, case);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 4: Coulomb only, opposite charges (attractive)
|
||||
// q1=+0.5e, q2=-0.5e, ε=0 (no LJ), r=0.5 nm
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_coulomb_opposite_charges() {
|
||||
let case = "coulomb_opposite";
|
||||
let nb = make_nb(&[(0.3, 0.0, 0.5), (0.3, 0.0, -0.5)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.5, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -6.946772882219099e+01, case);
|
||||
assert_force_eq(&forces[0], &[1.389354576443820e+02, 0.0, 0.0], 0, case);
|
||||
assert_force_eq(&forces[1], &[-1.389354576443820e+02, 0.0, 0.0], 1, case);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 5: Coulomb only, same charges (repulsive)
|
||||
// q1=+0.8e, q2=+0.8e, ε=0, r=0.3 nm
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_coulomb_same_charges() {
|
||||
let case = "coulomb_same";
|
||||
let nb = make_nb(&[(0.3, 0.0, 0.8), (0.3, 0.0, 0.8)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.3, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, 2.963956429746816e+02, case);
|
||||
assert_force_eq(&forces[0], &[-9.879854765822721e+02, 0.0, 0.0], 0, case);
|
||||
assert_force_eq(&forces[1], &[9.879854765822721e+02, 0.0, 0.0], 1, case);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 6: LJ + Coulomb combined
|
||||
// σ=0.3, ε=1.0, q1=+0.2, q2=-0.3, r=0.4 nm
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_coulomb_combined() {
|
||||
let case = "lj_coulomb_combined";
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.2), (0.3, 1.0, -0.3)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.4, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -2.142552730106098e+01, case);
|
||||
assert_force_eq(&forces[0], &[5.897834531125384e+01, 0.0, 0.0], 0, case);
|
||||
assert_force_eq(&forces[1], &[-5.897834531125384e+01, 0.0, 0.0], 1, case);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 7: 3 particles in triangle (LJ only)
|
||||
// Tests pairwise force accumulation
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_triangle_3() {
|
||||
let case = "lj_triangle_3";
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.0), (0.3, 1.0, 0.0), (0.3, 1.0, 0.0)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.4, 0.0, 0.0),
|
||||
Vec3::new(0.2, 0.35, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -1.713426964257671e+00, case);
|
||||
assert_force_eq(
|
||||
&forces[0],
|
||||
&[1.019072218881330e+01, 5.798053614854739e+00, 0.0],
|
||||
0,
|
||||
case,
|
||||
);
|
||||
assert_force_eq(
|
||||
&forces[1],
|
||||
&[-1.019072218881330e+01, 5.798053614854739e+00, 0.0],
|
||||
1,
|
||||
case,
|
||||
);
|
||||
assert_force_eq(
|
||||
&forces[2],
|
||||
&[0.0, -1.159610722970948e+01, 0.0],
|
||||
2,
|
||||
case,
|
||||
);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 8: 4 particles in 3D, mixed LJ + Coulomb
|
||||
// Different σ, ε, q for each particle (O, H, H, C -like)
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_mixed_4_particles_3d() {
|
||||
let case = "mixed_4_particles_3d";
|
||||
let nb = make_nb(&[
|
||||
(0.30, 0.8, 0.4),
|
||||
(0.10, 0.1, 0.1),
|
||||
(0.10, 0.1, 0.1),
|
||||
(0.34, 0.4, -0.6),
|
||||
]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.1, 0.0, 0.0),
|
||||
Vec3::new(0.0, 0.1, 0.0),
|
||||
Vec3::new(0.15, 0.15, 0.15),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
// Larger tolerance for this case due to very large force magnitudes
|
||||
let tol_e = 1e-8;
|
||||
let tol_f = 1e-4;
|
||||
|
||||
let diff_e = (energy - 9.059241986638443e+03).abs();
|
||||
assert!(
|
||||
diff_e / 9.059241986638443e+03 < tol_e,
|
||||
"[{case}] energy mismatch: got {energy:.15e}, expected 9.059241986638443e+03"
|
||||
);
|
||||
|
||||
let expected_forces: [[f64; 3]; 4] = [
|
||||
[-5.526476354461313e+05, -5.526476354461313e+05, -3.449575370649803e+02],
|
||||
[5.523849713068031e+05, 5.489578471528779e+01, 1.028918868390872e+02],
|
||||
[5.489578471528779e+01, 5.523849713068031e+05, 1.028918868390872e+02],
|
||||
[2.077683546128640e+02, 2.077683546128640e+02, 1.391737633868059e+02],
|
||||
];
|
||||
for (i, exp) in expected_forces.iter().enumerate() {
|
||||
for d in 0..3 {
|
||||
let diff = (forces[i][d] - exp[d]).abs();
|
||||
let scale = exp[d].abs().max(1.0);
|
||||
assert!(
|
||||
diff / scale < tol_f,
|
||||
"[{case}] force mismatch p{i} d{d}: got {:.10e}, expected {:.10e}, rel {:.2e}",
|
||||
forces[i][d],
|
||||
exp[d],
|
||||
diff / scale,
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 9: LJ in strongly repulsive region (r = 0.25 nm < σ = 0.3 nm)
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_repulsive_close() {
|
||||
let case = "lj_repulsive_close";
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.0), (0.3, 1.0, 0.0)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.25, 0.0, 0.0),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, 2.372046579302399e+01, case);
|
||||
assert_force_eq(
|
||||
&forces[0],
|
||||
&[-1.425236822065152e+03, 0.0, 0.0],
|
||||
0,
|
||||
case,
|
||||
);
|
||||
assert_force_eq(
|
||||
&forces[1],
|
||||
&[1.425236822065152e+03, 0.0, 0.0],
|
||||
1,
|
||||
case,
|
||||
);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Case 10: Argon dimer along 3D diagonal
|
||||
// σ=0.3405, ε=0.996 (real argon), non-axis-aligned
|
||||
// ============================================================
|
||||
#[test]
|
||||
fn test_lj_diagonal_3d() {
|
||||
let case = "lj_diagonal_3d";
|
||||
let nb = make_nb(&[(0.3405, 0.996, 0.0), (0.3405, 0.996, 0.0)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.1, 0.2, 0.3),
|
||||
Vec3::new(0.3, 0.4, 0.5),
|
||||
];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
|
||||
assert_energy_eq(energy, -3.524861309580062e-01, case);
|
||||
assert_force_eq(
|
||||
&forces[0],
|
||||
&[-2.888202074083446e+01, -2.888202074083446e+01, -2.888202074083446e+01],
|
||||
0,
|
||||
case,
|
||||
);
|
||||
assert_force_eq(
|
||||
&forces[1],
|
||||
&[2.888202074083446e+01, 2.888202074083446e+01, 2.888202074083446e+01],
|
||||
1,
|
||||
case,
|
||||
);
|
||||
}
|
||||
|
||||
// ============================================================
|
||||
// Property tests (physical invariants, not OpenMM reference)
|
||||
// ============================================================
|
||||
|
||||
#[test]
|
||||
fn test_newton_third_law_lj() {
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.0), (0.3, 1.5, 0.0), (0.25, 0.8, 0.0)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.35, 0.1, 0.0),
|
||||
Vec3::new(0.1, 0.4, 0.2),
|
||||
];
|
||||
let (_, forces) = run(&nb, &positions);
|
||||
|
||||
let total: SVector<f64, 3> = forces.iter().sum();
|
||||
for d in 0..3 {
|
||||
assert!(
|
||||
total[d].abs() < 1e-10,
|
||||
"Newton's 3rd law violated: total force dim {d} = {:.2e}",
|
||||
total[d]
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_newton_third_law_coulomb() {
|
||||
let nb = make_nb(&[(0.3, 0.0, 0.5), (0.3, 0.0, -0.3), (0.3, 0.0, 0.8)]);
|
||||
let positions = vec![
|
||||
Vec3::new(0.0, 0.0, 0.0),
|
||||
Vec3::new(0.3, 0.0, 0.0),
|
||||
Vec3::new(0.15, 0.25, 0.0),
|
||||
];
|
||||
let (_, forces) = run(&nb, &positions);
|
||||
|
||||
let total: SVector<f64, 3> = forces.iter().sum();
|
||||
for d in 0..3 {
|
||||
assert!(
|
||||
total[d].abs() < 1e-8,
|
||||
"Newton's 3rd law violated (Coulomb): total force dim {d} = {:.2e}",
|
||||
total[d]
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_lj_force_is_negative_gradient() {
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.0), (0.3, 1.0, 0.0)]);
|
||||
let h = 1e-7;
|
||||
let base = Vec3::new(0.38, 0.05, 0.02);
|
||||
|
||||
let positions = vec![SVector::<f64, 3>::zeros(), base];
|
||||
let (_, forces) = run(&nb, &positions);
|
||||
let f_analytical = forces[1];
|
||||
|
||||
for d in 0..3 {
|
||||
let mut pos_plus = positions.clone();
|
||||
let mut pos_minus = positions.clone();
|
||||
pos_plus[1][d] += h;
|
||||
pos_minus[1][d] -= h;
|
||||
|
||||
let (e_plus, _) = run(&nb, &pos_plus);
|
||||
let (e_minus, _) = run(&nb, &pos_minus);
|
||||
|
||||
let f_numerical = -(e_plus - e_minus) / (2.0 * h);
|
||||
let diff = (f_analytical[d] - f_numerical).abs();
|
||||
assert!(
|
||||
diff < 1e-3,
|
||||
"F != -grad(E) for particle 1, dim {d}: analytical={:.10e}, numerical={:.10e}, diff={diff:.2e}",
|
||||
f_analytical[d], f_numerical,
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_coulomb_force_is_negative_gradient() {
|
||||
let nb = make_nb(&[(0.3, 0.0, 0.5), (0.3, 0.0, -0.3)]);
|
||||
let h = 1e-7;
|
||||
let base = Vec3::new(0.4, 0.1, 0.05);
|
||||
|
||||
let positions = vec![SVector::<f64, 3>::zeros(), base];
|
||||
let (_, forces) = run(&nb, &positions);
|
||||
let f_analytical = forces[1];
|
||||
|
||||
for d in 0..3 {
|
||||
let mut pos_plus = positions.clone();
|
||||
let mut pos_minus = positions.clone();
|
||||
pos_plus[1][d] += h;
|
||||
pos_minus[1][d] -= h;
|
||||
|
||||
let (e_plus, _) = run(&nb, &pos_plus);
|
||||
let (e_minus, _) = run(&nb, &pos_minus);
|
||||
|
||||
let f_numerical = -(e_plus - e_minus) / (2.0 * h);
|
||||
let diff = (f_analytical[d] - f_numerical).abs();
|
||||
assert!(
|
||||
diff < 1e-3,
|
||||
"Coulomb F != -grad(E) dim {d}: analytical={:.10e}, numerical={:.10e}, diff={diff:.2e}",
|
||||
f_analytical[d], f_numerical,
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_zero_particles_no_crash() {
|
||||
let nb = NonbondedInteraction::<3>::new();
|
||||
let positions: Vec<Vec3> = vec![];
|
||||
let mut forces: Vec<Vec3> = vec![];
|
||||
let output = nb.compute_forces(&positions, &mut forces);
|
||||
assert_eq!(output.potential_energy, 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_single_particle_no_interaction() {
|
||||
let nb = make_nb(&[(0.3, 1.0, 0.5)]);
|
||||
let positions = vec![Vec3::new(0.1, 0.2, 0.3)];
|
||||
let (energy, forces) = run(&nb, &positions);
|
||||
assert_eq!(energy, 0.0);
|
||||
assert_eq!(forces[0], Vec3::zeros());
|
||||
}
|
||||
257
tests/reference/openmm_nonbonded.py
Normal file
257
tests/reference/openmm_nonbonded.py
Normal file
@@ -0,0 +1,257 @@
|
||||
#!/usr/bin/env python3
|
||||
"""
|
||||
Generate reference values for NonbondedInteraction using OpenMM.
|
||||
|
||||
OpenMM NonbondedForce (NoCutoff):
|
||||
U_LJ = 4ε[(σ/r)^12 - (σ/r)^6]
|
||||
U_Coulomb = k_e * q_i * q_j / r (k_e = 138.935456 kJ·nm/(mol·e²))
|
||||
|
||||
Combining rules (Lorentz-Berthelot):
|
||||
σ_ij = (σ_i + σ_j) / 2
|
||||
ε_ij = sqrt(ε_i * ε_j)
|
||||
|
||||
Units: distance=nm, energy=kJ/mol, force=kJ/mol/nm, charge=e
|
||||
"""
|
||||
|
||||
import openmm
|
||||
from openmm import unit
|
||||
import json
|
||||
import numpy as np
|
||||
|
||||
|
||||
def compute_nonbonded(particles, positions_nm):
|
||||
"""
|
||||
Args:
|
||||
particles: list of (mass_Da, sigma_nm, epsilon_kjmol, charge_e)
|
||||
positions_nm: list of (x, y, z) in nm
|
||||
Returns:
|
||||
(energy_kj_mol, forces_kj_mol_nm)
|
||||
"""
|
||||
system = openmm.System()
|
||||
nb_force = openmm.NonbondedForce()
|
||||
nb_force.setNonbondedMethod(openmm.NonbondedForce.NoCutoff)
|
||||
|
||||
for mass, sigma, epsilon, charge in particles:
|
||||
system.addParticle(mass * unit.daltons)
|
||||
nb_force.addParticle(
|
||||
charge * unit.elementary_charge,
|
||||
sigma * unit.nanometers,
|
||||
epsilon * unit.kilojoules_per_mole,
|
||||
)
|
||||
|
||||
system.addForce(nb_force)
|
||||
|
||||
integrator = openmm.VerletIntegrator(0.001 * unit.picoseconds)
|
||||
platform = openmm.Platform.getPlatformByName("Reference")
|
||||
context = openmm.Context(system, integrator, platform)
|
||||
|
||||
pos = [openmm.Vec3(*p) * unit.nanometers for p in positions_nm]
|
||||
context.setPositions(pos)
|
||||
|
||||
state = context.getState(getEnergy=True, getForces=True)
|
||||
energy = state.getPotentialEnergy().value_in_unit(unit.kilojoules_per_mole)
|
||||
forces = state.getForces(asNumpy=True).value_in_unit(
|
||||
unit.kilojoules_per_mole / unit.nanometers
|
||||
)
|
||||
return energy, forces.tolist()
|
||||
|
||||
|
||||
def print_case(name, energy, forces):
|
||||
print(f"\n{'='*60}")
|
||||
print(f" Test Case: {name}")
|
||||
print(f"{'='*60}")
|
||||
print(f" Energy (kJ/mol): {energy:.15e}")
|
||||
print(f" Forces (kJ/mol/nm):")
|
||||
for i, f in enumerate(forces):
|
||||
print(f" particle {i}: ({f[0]:.15e}, {f[1]:.15e}, {f[2]:.15e})")
|
||||
|
||||
|
||||
def print_rust_const(name, energy, forces):
|
||||
tag = name.upper()
|
||||
print(f"\n// --- {name} ---")
|
||||
print(f"const {tag}_ENERGY: f64 = {energy:.15e};")
|
||||
for i, f in enumerate(forces):
|
||||
print(
|
||||
f"const {tag}_FORCE_{i}: [f64; 3] = "
|
||||
f"[{f[0]:.15e}, {f[1]:.15e}, {f[2]:.15e}];"
|
||||
)
|
||||
|
||||
|
||||
def main():
|
||||
results = {}
|
||||
|
||||
# ================================================================
|
||||
# Case 1: LJ only, 2 identical particles along x-axis
|
||||
# σ=0.3 nm, ε=1.0 kJ/mol, q=0
|
||||
# r = 0.4 nm (repulsive region, r < 2^(1/6)*σ ≈ 0.3366 → actually r=0.4 > that, attractive)
|
||||
# ================================================================
|
||||
name = "lj_two_identical"
|
||||
particles = [
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (0.4, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 2: LJ only, 2 particles at equilibrium (r = 2^(1/6)*σ)
|
||||
# Force should be ~zero, energy = -ε
|
||||
# ================================================================
|
||||
name = "lj_equilibrium"
|
||||
r_eq = 0.3 * 2.0 ** (1.0 / 6.0)
|
||||
particles = [
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (r_eq, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 3: LJ only, 2 particles with different σ/ε
|
||||
# Lorentz-Berthelot: σ_ij = (0.25+0.35)/2 = 0.3, ε_ij = sqrt(0.5*2.0) = 1.0
|
||||
# ================================================================
|
||||
name = "lj_different_params"
|
||||
particles = [
|
||||
(12.0, 0.25, 0.5, 0.0),
|
||||
(16.0, 0.35, 2.0, 0.0),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (0.35, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 4: Coulomb only, opposite charges (attractive)
|
||||
# q1=+0.5e, q2=-0.5e, no LJ (ε=0)
|
||||
# ================================================================
|
||||
name = "coulomb_opposite"
|
||||
particles = [
|
||||
(12.0, 0.3, 0.0, 0.5),
|
||||
(12.0, 0.3, 0.0, -0.5),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (0.5, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 5: Coulomb only, same charges (repulsive)
|
||||
# q1=+0.8e, q2=+0.8e, no LJ (ε=0)
|
||||
# ================================================================
|
||||
name = "coulomb_same"
|
||||
particles = [
|
||||
(12.0, 0.3, 0.0, 0.8),
|
||||
(12.0, 0.3, 0.0, 0.8),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (0.3, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 6: LJ + Coulomb combined
|
||||
# σ=0.3, ε=1.0, q1=+0.2, q2=-0.3
|
||||
# ================================================================
|
||||
name = "lj_coulomb_combined"
|
||||
particles = [
|
||||
(12.0, 0.3, 1.0, 0.2),
|
||||
(16.0, 0.3, 1.0, -0.3),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (0.4, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 7: 3 particles in triangle (LJ only)
|
||||
# Tests force accumulation on all particles
|
||||
# ================================================================
|
||||
name = "lj_triangle_3"
|
||||
particles = [
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
]
|
||||
positions = [
|
||||
(0.0, 0.0, 0.0),
|
||||
(0.4, 0.0, 0.0),
|
||||
(0.2, 0.35, 0.0),
|
||||
]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 8: 4 particles in 3D, mixed LJ + Coulomb
|
||||
# Different σ, ε, q for each particle
|
||||
# ================================================================
|
||||
name = "mixed_4_particles_3d"
|
||||
particles = [
|
||||
(16.0, 0.30, 0.8, +0.4), # O-like
|
||||
(1.0, 0.10, 0.1, +0.1), # H-like
|
||||
(1.0, 0.10, 0.1, +0.1), # H-like
|
||||
(12.0, 0.34, 0.4, -0.6), # C-like
|
||||
]
|
||||
positions = [
|
||||
(0.0, 0.0, 0.0),
|
||||
(0.1, 0.0, 0.0),
|
||||
(0.0, 0.1, 0.0),
|
||||
(0.15, 0.15, 0.15),
|
||||
]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 9: LJ in strongly repulsive region (r << σ)
|
||||
# r = 0.25 nm < σ = 0.3 nm
|
||||
# ================================================================
|
||||
name = "lj_repulsive_close"
|
||||
particles = [
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
(12.0, 0.3, 1.0, 0.0),
|
||||
]
|
||||
positions = [(0.0, 0.0, 0.0), (0.25, 0.0, 0.0)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# Case 10: Argon dimer along 3D diagonal
|
||||
# Tests non-axis-aligned geometry
|
||||
# ================================================================
|
||||
name = "lj_diagonal_3d"
|
||||
particles = [
|
||||
(39.948, 0.3405, 0.996, 0.0),
|
||||
(39.948, 0.3405, 0.996, 0.0),
|
||||
]
|
||||
positions = [(0.1, 0.2, 0.3), (0.3, 0.4, 0.5)]
|
||||
energy, forces = compute_nonbonded(particles, positions)
|
||||
print_case(name, energy, forces)
|
||||
print_rust_const(name, energy, forces)
|
||||
results[name] = {"energy": energy, "forces": forces}
|
||||
|
||||
# ================================================================
|
||||
# JSON output
|
||||
# ================================================================
|
||||
print(f"\n\n{'='*60}")
|
||||
print(" JSON Reference Data")
|
||||
print(f"{'='*60}")
|
||||
print(json.dumps(results, indent=2))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in New Issue
Block a user