Add DampedShiftedForce potential

Author: seatonullberg

crates/velvet-core/Cargo.toml

@@ -9,6 +9,7 @@ edition = "2018"
 
 [dependencies]
 indicatif = "0.15"
+libm = "0.2"
 nalgebra = "0.26"
 rand = "0.7"
 rand_distr = "0.3"

crates/velvet-core/src/potentials/coulomb.rs

@@ -1,8 +1,15 @@
 //! Potentials which describe Coulombic electrostatic interactions.
 
+#[cfg(feature = "f64")]
+use libm::erfc as erfc;
+
+#[cfg(not(feature = "f64"))]
+use libm::erfcf as erfc;
+
 use crate::internal::consts::COULOMB;
 use crate::internal::Float;
-use crate::potentials::types::StandardCoulombic;
+use crate::internal::consts::FRAC_2_SQRT_PI;
+use crate::potentials::types::{DampedShiftedForce, StandardCoulombic};
 use crate::potentials::Potential;
 use crate::selection::{setup_pairs_with_charge, update_pairs_by_cutoff_radius, Selection};
 use crate::system::System;
@@ -15,6 +22,36 @@ pub trait CoulombPotential: Potential {
     fn force(&self, qi: Float, qj: Float, r: Float) -> Float;
 }
 
+impl CoulombPotential for DampedShiftedForce {
+    fn energy(&self, qi: Float, qj: Float, r: Float) -> Float {
+        let factor = FRAC_2_SQRT_PI * self.alpha;
+        let alpha2 = self.alpha.powi(2);
+        let cutoff2 = self.cutoff.powi(2);
+
+        let term_a = erfc(self.alpha * r) / r;
+        let term_b = erfc(self.alpha * self.cutoff) / self.cutoff;
+        let term_c = erfc(self.alpha * self.cutoff) / cutoff2;
+        let term_d = factor * (Float::exp(-alpha2 * cutoff2) / self.cutoff);
+        let term_e = r - self.cutoff;
+
+        qi * qj * (term_a - term_b + (term_c + term_d) * term_e)
+    }
+
+    fn force(&self, qi: Float, qj: Float, r: Float) -> Float {
+        let factor = FRAC_2_SQRT_PI * self.alpha;
+        let r2 = r.powi(2);
+        let alpha2 = self.alpha.powi(2);
+        let cutoff2 = self.cutoff.powi(2);
+
+        let term_a = erfc(self.alpha * r) / r2;
+        let term_b = factor * Float::exp(-alpha2 * r2) / r;
+        let term_c = erfc(self.alpha * self.cutoff) / cutoff2;
+        let term_d = factor * Float::exp(-alpha2 * cutoff2) / self.cutoff;
+
+        qi * qj * ((term_a + term_b) - (term_c + term_d))
+    }
+}
+
 impl CoulombPotential for StandardCoulombic {
     fn energy(&self, qi: Float, qj: Float, r: Float) -> Float {
         (COULOMB * qi * qj) / (self.dielectric * r)

crates/velvet-core/src/potentials/types.rs

@@ -23,6 +23,26 @@ impl Buckingham {
 
 impl Potential for Buckingham {}
 
+/// [Damped Shifted Force](https://lammps.sandia.gov/doc/pair_coul.html#description) potential.
+#[derive(Clone, Copy, Debug)]
+pub struct DampedShiftedForce {
+    /// Damping parameter.
+    pub alpha: Float,
+    /// Cutoff radius
+    pub cutoff: Float,
+}
+
+impl DampedShiftedForce {
+    /// Returns a new [`DampedShiftedForce`] potential.
+    pub fn new(alpha: Float, cutoff: Float) -> DampedShiftedForce {
+        DampedShiftedForce {alpha, cutoff}
+    }
+}
+
+impl Potential for DampedShiftedForce {}
+
+
+
 /// [Harmonic](https://lammps.sandia.gov/doc/bond_harmonic.html#description) oscillator potential.
 #[derive(Clone, Copy, Debug)]
 pub struct Harmonic {

examples/magnesium-oxide.rs

@@ -16,25 +16,21 @@ fn main() {
     let buck_mg_o = Buckingham::new(18946.9176, 0.32, 0.0);
     let buck_o_o = Buckingham::new(524960.604, 0.149, 642.94068);
 
-    // Initialize standard Coulombic potentials between all charged particles.
-    let coul = StandardCoulombic::new(1.0);
+    // Initialize a DSF potential to evaluate electrostatic interactions.
+    let dsf = DampedShiftedForce::new(0.1, 10.0);
+    // let coul = StandardCoulombic::new(1.0);
 
     // Store all of the system's potentials in a Potentials struct.
-    let cutoff_buck = 5.0;
-    let cutoff_coul = 10.0;
-    let thickness = 3.0;
     let potentials = PotentialsBuilder::new()
-        .update_frequency(3)
-        .pair(buck_mg_o, (magnesium, oxygen), cutoff_buck, thickness)
-        .pair(buck_o_o, (oxygen, oxygen), cutoff_buck, thickness)
-        .coulomb(coul, cutoff_coul, thickness)
+        .update_frequency(1)
+        .coulomb(dsf, 10.0, 3.0)
         .build();
 
     // Initialize a velocity Verlet style integrator.
-    let velocity_verlet = VelocityVerlet::new(1e-6);
+    let velocity_verlet = VelocityVerlet::new(0.001);
 
     // Initialize a Nose-Hoover style thermostat.
-    let nose_hoover = NoseHoover::new(300.0, 1.25, 1e-4);
+    let nose_hoover = NoseHoover::new(300.0, 1.25, 0.1);
 
     // Run MD with a Nose-Hoover thermostat to simulate the NVT ensemble.
     let md = MolecularDynamics::new(velocity_verlet, nose_hoover);
@@ -56,5 +52,5 @@ fn main() {
 
     // Run the simulation.
     let mut sim = Simulation::new(system, potentials, md, config);
-    sim.run(250_000);
+    sim.run(50_000);
 }

resources/lammps/magnesium-oxide/NVE.dat

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:97738dfdc2e5afcae0e5fd552e7d931dba5e3c63753bfd20f2e1fd0b6331ead0
-size 4017
+oid sha256:c65504f513f0af286c3fdda2927dfada301266ef7d6014c5d2bc8d5e6c3e7ca6
+size 396

resources/lammps/magnesium-oxide/lammps.in

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:fd7994cdc5ee50c39f9955d72bc451decbc11958ea2bf963a4cb525a2cdfd186
-size 890
+oid sha256:9a6e41063b0b29ea9bf13cff6b56989de865aef742887ccbf81747f87ad3fc53
+size 939

resources/lammps/magnesium-oxide/log.lammps

@@ -1,3 +0,0 @@
-version https://git-lfs.github.com/spec/v1
-oid sha256:f37d8c76b6889ab5d32418e01f823346cab3d866c7bc4578d951e3434e635d96
-size 11123