Implement tail corrections to the energy and virial

examples/custom-potential.rs

@@ -21,21 +21,31 @@ struct LJ {
     b: f64
 }
 
-/// All we need to do is to implement the Potential trait
+// All we need to do is to implement the Potential trait
 impl Potential for LJ {
-    /// The energy function give the energy at distance `r`
+    // The energy function give the energy at distance `r`
     fn energy(&self, r: f64) -> f64 {
         self.a / r.powi(12) - self.b / r.powi(6)
     }
 
-    /// The force function give the norm of the force at distance `r`
+    // The force function give the norm of the force at distance `r`
     fn force(&self, r: f64) -> f64 {
         12.0 * self.a / r.powi(13) - 6.0 * self.b / r.powi(7)
     }
 }
 
 // We want to use our LJ potential as a pair potential.
-impl PairPotential for LJ {}
+impl PairPotential for LJ {
+    // The long-range correction to the energy at the given cutoff
+    fn tail_energy(&self, cutoff: f64) -> f64 {
+        - (1.0/9.0 * self.a / cutoff.powi(9) - 1.0/3.0 * self.b / cutoff.powi(3))
+    }
+
+    // The long-range correction to the virial at the given cutoff
+    fn tail_virial(&self, cutoff: f64) -> f64 {
+        - (12.0/9.0 * self.a / cutoff.powi(9) - 6.0/3.0 * self.b / cutoff.powi(3))
+    }
+}
 
 fn main() {
     Logger::stdout();

src/core/src/energy/computations.rs

@@ -7,7 +7,7 @@ use energy::{Potential, PairPotential};
 ///
 /// The `Computation` trait represent an alternative way to compute a given
 /// potential. For example using interpolation on a table or on a grid, from a
-/// Fourrier decompostion, *etc.*
+/// Fourier decomposition, *etc.*
 ///
 /// # Examples
 ///
@@ -16,7 +16,7 @@ use energy::{Potential, PairPotential};
 /// use lumol::energy::Computation;
 /// use lumol::energy::Harmonic;
 ///
-/// /// This is just a thin wrapper logging everytime the `energy/force`
+/// /// This is just a thin wrapper logging every time the `energy/force`
 /// /// methods are called.
 /// #[derive(Clone)]
 /// struct LoggingComputation<T: Potential>(T);
@@ -62,16 +62,19 @@ impl<P: Computation + Clone + 'static> Potential for P {
 /// This can be faster than direct computation for smooth potentials, but will
 /// uses more memory and be less precise than direct computation. Values are
 /// tabulated in the `[0, max)` range, and a cutoff is applied after `max`.
-/// Energy is shifted to ensure `E(max) = 0`
 #[derive(Clone)]
 pub struct TableComputation {
     /// Step for tabulated value. `energy_table[i]`/`force_table[i]` contains
     /// energy/force at `r = i * delta`
     delta: f64,
+    /// Cutoff distance
+    cutoff: f64,
     /// Tabulated potential
     energy_table: Vec<f64>,
     /// Tabulated compute_force
     force_table: Vec<f64>,
+    /// Initial potential, kept around for tail corrections
+    potential: Box<PairPotential>,
 }
 
 
@@ -86,26 +89,28 @@ impl TableComputation {
     /// use lumol::energy::TableComputation;
     /// use lumol::energy::Harmonic;
     ///
-    /// let potential = Harmonic{x0: 0.5, k: 4.2};
-    /// let table = TableComputation::new(&potential, 1000, 2.0);
+    /// let potential = Box::new(Harmonic{x0: 0.5, k: 4.2});
+    /// let table = TableComputation::new(potential, 1000, 2.0);
     ///
-    /// assert_eq!(table.energy(1.0), -4.2);
+    /// assert_eq!(table.energy(1.0), 0.525);
     /// assert_eq!(table.energy(3.0), 0.0);
     /// ```
-    pub fn new(potential: &PairPotential, size: usize, max:f64) -> TableComputation {
+    pub fn new(potential: Box<PairPotential>, size: usize, max: f64) -> TableComputation {
         let delta = max / (size as f64);
-        let energy_shift = potential.energy(max);
         let mut energy_table = Vec::with_capacity(size);
         let mut force_table = Vec::with_capacity(size);
         for i in 0..size {
-            let pos = i as f64 * delta;
-            energy_table.push(potential.energy(pos) - energy_shift);
-            force_table.push(potential.force(pos));
+            let r = i as f64 * delta;
+            energy_table.push(potential.energy(r));
+            force_table.push(potential.force(r));
         }
+
         TableComputation {
             delta: delta,
+            cutoff: max,
             energy_table: energy_table,
-            force_table: force_table
+            force_table: force_table,
+            potential: potential,
         }
     }
 }
@@ -136,23 +141,44 @@ impl Computation for TableComputation {
     }
 }
 
-impl PairPotential for TableComputation {}
+impl PairPotential for TableComputation {
+    fn tail_energy(&self, cutoff: f64) -> f64 {
+        if cutoff > self.cutoff {
+            warn_once!("Cutoff in table computation ({}) is smaller than the \
+            pair interaction cutoff ({}) when computing tail correction. This \
+            may lead to wrong values for energy.", cutoff, self.cutoff);
+        }
+        return self.potential.tail_energy(cutoff);
+    }
+
+    fn tail_virial(&self, cutoff: f64) -> f64 {
+        if cutoff > self.cutoff {
+            warn_once!("Cutoff in table computation ({}) is smaller than the \
+            pair interaction cutoff ({}) when computing tail correction. This \
+            may lead to wrong values for pressure.", cutoff, self.cutoff);
+        }
+        return self.potential.tail_virial(cutoff);
+    }
+}
 
 #[cfg(test)]
 mod test {
     use super::*;
-    use energy::Harmonic;
+    use energy::{Harmonic, LennardJones};
+    use energy::PairPotential;
 
     #[test]
     fn table() {
-        let table = TableComputation::new(&Harmonic{k: 50.0, x0: 2.0}, 1000, 4.0);
+        let table = TableComputation::new(
+            Box::new(Harmonic{k: 50.0, x0: 2.0}), 1000, 4.0
+        );
 
-        assert_eq!(table.compute_energy(2.5), -93.75);
+        assert_eq!(table.compute_energy(2.5), 6.25);
         assert_eq!(table.compute_force(2.5), -25.0);
 
         // Check that the table is defined up to the cutoff value
         let delta = 4.0 / 1000.0;
-        assert_eq!(table.compute_energy(4.0 - 2.0 * delta), -0.7984000000000009);
+        assert_eq!(table.compute_energy(4.0 - 2.0 * delta), 99.2016);
         assert_eq!(table.compute_force(4.0 - 2.0 * delta), -99.6);
 
         assert_eq!(table.compute_energy(4.0 - delta), 0.0);
@@ -163,5 +189,11 @@ mod test {
 
         assert_eq!(table.compute_energy(4.1), 0.0);
         assert_eq!(table.compute_force(4.1), 0.0);
+
+
+        let lj = LennardJones{epsilon: 50.0, sigma: 2.0};
+        let table = TableComputation::new(Box::new(lj.clone()), 1000, 4.0);
+        assert_eq!(table.tail_energy(5.0), lj.tail_energy(5.0));
+        assert_eq!(table.tail_virial(5.0), lj.tail_virial(5.0));
     }
 }

src/core/src/energy/functions.rs

@@ -30,7 +30,10 @@ impl Potential for NullPotential {
     #[inline] fn force(&self, _: f64) -> f64 {0.0}
 }
 
-impl PairPotential for NullPotential {}
+impl PairPotential for NullPotential {
+    fn tail_energy(&self, _: f64) -> f64 {0.0}
+    fn tail_virial(&self, _: f64) -> f64 {0.0}
+}
 impl BondPotential for NullPotential {}
 impl AnglePotential for NullPotential {}
 impl DihedralPotential for NullPotential {}
@@ -75,7 +78,23 @@ impl Potential for LennardJones {
     }
 }
 
-impl PairPotential for LennardJones {}
+impl PairPotential for LennardJones {
+    fn tail_energy(&self, cutoff: f64) -> f64 {
+        let s3 = self.sigma * self.sigma * self.sigma;
+        let rc3 = cutoff * cutoff * cutoff;
+        let s9 = s3 * s3 * s3;
+        let rc9 = rc3 * rc3 * rc3;
+        return 4.0 / 3.0 * self.epsilon * s3 * (1.0 / 3.0 * s9 / rc9 - s3 / rc3);
+    }
+
+    fn tail_virial(&self, cutoff: f64) -> f64 {
+        let s3 = self.sigma * self.sigma * self.sigma;
+        let rc3 = cutoff * cutoff * cutoff;
+        let s9 = s3 * s3 * s3;
+        let rc9 = rc3 * rc3 * rc3;
+        return 8.0 * self.epsilon * s3 * (2.0 / 3.0 * s9 / rc9 - s3 / rc3);
+    }
+}
 
 /// Harmonic potential.
 ///
@@ -116,7 +135,14 @@ impl Potential for Harmonic {
     }
 }
 
-impl PairPotential for Harmonic {}
+impl PairPotential for Harmonic {
+    // These two functions should return infinity, as the Harmonic potential
+    // does not goes to zero at infinity. We use 0 instead to ignore the tail
+    // contribution to the energy/virial.
+    fn tail_energy(&self, _: f64) -> f64 {0.0}
+    fn tail_virial(&self, _: f64) -> f64 {0.0}
+}
+
 impl BondPotential for Harmonic {}
 impl AnglePotential for Harmonic {}
 impl DihedralPotential for Harmonic {}
@@ -173,8 +199,8 @@ impl DihedralPotential for CosineHarmonic {}
 
 /// Torsion potential.
 ///
-/// This potential is intended for use with dihedral angles, using a
-/// customisable periodicity and multiple minima.
+/// This potential is intended for use with dihedral angles, using a custom
+/// periodicity and multiple minima.
 ///
 /// The following potential expression is used: `V(x) = k * (1 + cos(n * x -
 /// delta))` where `k` is the force constant, `n` the periodicity of the
@@ -225,7 +251,7 @@ impl DihedralPotential for Torsion {}
 #[cfg(test)]
 mod tests {
     use super::*;
-    use energy::Potential;
+    use energy::{Potential, PairPotential};
     const EPS: f64 = 1e-9;
 
     #[test]
@@ -237,6 +263,9 @@ mod tests {
         assert_eq!(null.force(2.0), 0.0);
         assert_eq!(null.force(2.5), 0.0);
 
+        assert_eq!(null.tail_energy(1.0), 0.0);
+        assert_eq!(null.tail_virial(1.0), 0.0);
+
         let e0 = null.energy(2.0);
         let e1 = null.energy(2.0 + EPS);
         assert_approx_eq!((e0 - e1)/EPS, null.force(2.0), EPS);
@@ -248,6 +277,14 @@ mod tests {
         assert_eq!(lj.energy(2.0), 0.0);
         assert_eq!(lj.energy(2.5), -0.6189584744448002);
 
+        assert_eq!(lj.tail_energy(1.0), 1388.0888888888887);
+        assert_eq!(lj.tail_energy(2.0), -5.688888888888889);
+        assert_eq!(lj.tail_energy(14.42), -0.022767318648783084);
+
+        assert_eq!(lj.tail_virial(1.0), 17066.666666666668);
+        assert_eq!(lj.tail_virial(2.0), -17.06666666666667);
+        assert_eq!(lj.tail_virial(14.42), -0.1366035877536718);
+
         assert_approx_eq!(lj.force(f64::powf(2.0, 1.0/6.0) * 2.0), 0.0);
         assert_approx_eq!(lj.force(2.5), -0.95773475733504);
 
@@ -265,6 +302,9 @@ mod tests {
         assert_eq!(harm.force(2.0), 0.0);
         assert_eq!(harm.force(2.5), -25.0);
 
+        assert_eq!(harm.tail_energy(1.0), 0.0);
+        assert_eq!(harm.tail_virial(1.0), 0.0);
+
         let e0 = harm.energy(2.1);
         let e1 = harm.energy(2.1 + EPS);
         assert_approx_eq!((e0 - e1)/EPS, harm.force(2.1), 1e-6);

src/core/src/energy/mod.rs

@@ -36,7 +36,11 @@
 //!
 //! // It can be used for pair and dihedral potentials, but not for angles or
 //! // bonds.
-//! impl PairPotential for OnePotential {}
+//! impl PairPotential for OnePotential {
+//!     /* some code omitted */
+//!     # fn tail_energy(&self, cutoff: f64) -> f64 {0.0}
+//!     # fn tail_virial(&self, cutoff: f64) -> f64 {0.0}
+//! }
 //! impl DihedralPotential for OnePotential {}
 //! ```
 //!
@@ -104,6 +108,7 @@ pub trait Potential : Sync + Send {
 /// # Example
 ///
 /// ```
+/// # use lumol::types::{Zero, Matrix3};
 /// use lumol::energy::{Potential, PairPotential};
 ///
 /// // A no-op potential
@@ -115,8 +120,17 @@ pub trait Potential : Sync + Send {
 ///     fn force(&self, x: f64) -> f64 {0.0}
 /// }
 ///
-/// // Now we can use the Null potential for pair interactions
-/// impl PairPotential for Null {}
+/// // By implementing this trait, we can use the Null potential for pair
+/// // interactions
+/// impl PairPotential for Null {
+///     fn tail_energy(&self, cutoff: f64) -> f64 {
+///         return 0.0;
+///     }
+///
+///     fn tail_virial(&self, cutoff: f64) -> f64 {
+///         return 0.0;
+///     }
+/// }
 /// ```
 pub trait PairPotential : Potential + BoxClonePair {
     /// Compute the virial contribution corresponding to the distance `r`
@@ -127,6 +141,26 @@ pub trait PairPotential : Potential + BoxClonePair {
         let force = fact * rn;
         force.tensorial(r)
     }
+
+    /// Compute the tail correction to the energy for the given cutoff.
+    ///
+    /// Calling `V(r)` the `Potential::energy(r)` function corresponding to this
+    /// potential, this function should return the integral from `cutoff` to
+    /// infinity of `r^2 V(r)`: `\int_{cutoff}^\infty r^2 V(r) dr`.
+    ///
+    /// If this integral does not converge for the current potential, this
+    /// function should then return 0 to disable tail corrections.
+    fn tail_energy(&self, cutoff: f64) -> f64;
+
+    /// Compute the tail correction to the virial for the given cutoff.
+    ///
+    /// Calling `f(r)` the `Potential::force(r)` function corresponding to this
+    /// potential, this function should return the integral from `cutoff` to
+    /// infinity of `f(r) r^3`: `\int_{cutoff}^\infty r^3 f(r) dr`.
+    ///
+    /// If this integral does not converge for the current potential, this
+    /// function should then return 0.0 to disable tail corrections.
+    fn tail_virial(&self, cutoff: f64) -> f64;
 }
 impl_box_clone!(PairPotential, BoxClonePair, box_clone_pair);