Merge pull request #42 from sansona/master

Added sigma epsilon version of Lennard-Jones potential and the square well potential although only for an energy calculation

.zenodo.json

@@ -25,6 +25,11 @@
 			"affiliation": "Department of Chemistry, University of Bath, Claverton Down, Bath, BA2 7AY, UK.",
 			"orcid": "0000-0003-3804-0975"
 		}
+		{
+			"name": "Jiaming Chen",
+			"affiliation": "Department of Chemistry & Biochemistry, University of California, Los Angeles."
+
+		}
 	],
 	"description": "An open source teaching tool for classical atomistic simulation",
 	"access_right": "open",

pylj/forcefields.py

@@ -19,21 +19,57 @@ def lennard_jones(dr, constants, force=False):
         The distances between the all pairs of particles.
     constants: float, array_like
         An array of length two consisting of the A and B parameters for the
-        12-6 Lennard-Jones function.
+        12-6 Lennard-Jones function
     force: bool (optional)
         If true, the negative first derivative will be found.
 
     Returns
     -------
-    float:
+    float: array_like
         The potential energy or force between the particles.
     """
     if force:
         return 12 * constants[0] * np.power(dr, -13) - (
-            6 * constants[1] * np.power(dr, -7)
-        )
+            6 * constants[1] * np.power(dr, -7))
+
+    else:
+        return constants[0] * np.power(dr, -12) - (
+            constants[1] * np.power(dr, -6))
+
+
+@jit
+def lennard_jones_sigma_epsilon(dr, constants, force=False):
+    r"""Calculate the energy or force for a pair of particles using the
+    Lennard-Jones (sigma/epsilon variant) forcefield.
+
+    .. math::
+        E = \frac{4e*a^{12}}{dr^{12}} - \frac{4e*a^{6}}{dr^6}
+
+    .. math::
+        f = \frac{48e*a^{12}}{dr^{13}} - \frac{24e*a^{6}}{dr^7}
+
+    Parameters
+    ----------
+    dr: float, array_like
+        The distances between the all pairs of particles.
+    constants: float, array_like
+        An array of length two consisting of the sigma (a) and epsilon (e)
+        parameters for the 12-6 Lennard-Jones function
+    force: bool (optional)
+        If true, the negative first derivative will be found.
+
+    Returns
+    -------
+    float: array_like
+        The potential energy or force between the particles.
+    """
+    if force:
+        return 48 * constants[1] * np.power(constants[0], 12) * np.power(
+            dr, -13) - (24 * constants[1] * np.power(
+                constants[0], 6) * np.power(dr, -7))
     else:
-        return constants[0] * np.power(dr, -12) - (constants[1] * np.power(dr, -6))
+        return 4 * constants[1] * np.power(dr, -12) - (
+            4 * constants[1] * np.power(constants[0], 6) * np.power(dr, -6))
 
 
 @jit
@@ -47,26 +83,86 @@ def buckingham(dr, constants, force=False):
     .. math::
         f = ABe^{(-Bdr)} - \frac{6C}{dr^7}
 
-    Paramters
-    ---------
+    Parameters
+    ----------
     dr: float, array_like
         The distances between all the pairs of particles.
     constants: float, array_like
-        An array of lenght three consisting of the A, B and C parameters for
+        An array of length three consisting of the A, B and C parameters for
         the Buckingham function.
     force: bool (optional)
         If true, the negative first derivative will be found.
 
     Returns
     -------
-    float:
-        the potential energy or force between the particles.
+    float: array_like
+        The potential energy or force between the particles.
     """
     if force:
-        return constants[0] * constants[1] * np.exp(-constants[1] * dr) - 6 * constants[
-            2
-        ] / np.power(dr, 7)
+        return constants[0] * constants[1] * np.exp(
+            -constants[1] * dr) - 6 * constants[2] / np.power(dr, 7)
+    else:
+        return constants[0] * np.exp(
+            -constants[1] * dr) - constants[2] / np.power(dr, 6)
+
+
+def square_well(dr, constants, max_val=np.inf, force=False):
+    r'''Calculate the energy or force for a pair of particles using a
+    square well model.
+
+    .. math::
+        E = {
+        if dr < sigma:
+            E = max_val
+        elif sigma <= dr < lambda * sigma:
+            E = -epsilon
+        elif r >= lambda * sigma:
+            E = 0
+        }
+    .. math::
+        f = {
+        if sigma <= dr < lambda * sigma:
+            f = inf
+        else:
+            f = 0
+        }
+    Parameters
+    ----------
+    dr: float, array_like
+        The distances between all the pairs of particles.
+    constants: float, array_like
+        An array of length three consisting of the epsilon, sigma, and lambda
+        parameters for the square well model.
+    max_val: int (optional)
+        Upper bound for values in square well - replaces usual infinite values
+    force: bool (optional)
+        If true, the negative first derivative will be found.
+
+    Returns
+    -------
+    float: array_like
+        The potential energy between the particles.
+    '''
+    if not isinstance(dr, np.ndarray):
+        if isinstance(dr, list):
+            dr = np.array(dr, dtype='float')
+        elif isinstance(dr, float):
+            dr = np.array([dr], dtype='float')
+
+    if force:
+        raise ValueError("Force is infinite at sigma <= dr < lambda * sigma")
+
     else:
-        return constants[0] * np.exp(-constants[1] * dr) - constants[2] / np.power(
-            dr, 6
-        )
+        E = np.zeros_like(dr)
+        E[np.where(dr < constants[0])] = max_val
+        E[np.where(dr >= constants[2] * constants[1])] = 0
+
+        # apply mask for sigma <= dr < lambda * sigma
+        a = constants[1] <= dr
+        b = dr < constants[2] * constants[1]
+        E[np.where(a & b)] = -constants[0]
+
+        if len(E) == 1:
+            return float(E[0])
+        else:
+            return np.array(E, dtype='float')

pylj/tests/test_forcefields.py

@@ -1,4 +1,4 @@
-from numpy.testing import assert_almost_equal
+from numpy.testing import assert_almost_equal, assert_equal
 from pylj import forcefields
 import unittest
 
@@ -12,10 +12,46 @@ def test_lennard_jones_force(self):
         a = forcefields.lennard_jones(2.0, [1.0, 1.0], force=True)
         assert_almost_equal(a, -0.045410156)
 
+    def test_lennard_jones_sigma_epsilon_energy(self):
+        a = forcefields.lennard_jones_sigma_epsilon(2.0, [1.0, 0.25])
+        assert_almost_equal(a, -0.015380859)
+
+    def test_lennard_jones_sigma_epsilon_force(self):
+        a = forcefields.lennard_jones_sigma_epsilon(
+            2.0, [1.0, 0.25], force=True)
+        assert_almost_equal(a, -0.045410156)
+
     def test_buckingham_energy(self):
         a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0])
         assert_almost_equal(a, 0.1197103832)
 
     def test_buckingham_force(self):
         a = forcefields.buckingham(2.0, [1.0, 1.0, 1.0], force=True)
         assert_almost_equal(a, 0.08846028324)
+
+    def test_square_well_energy(self):
+        a = forcefields.square_well(2.0, [1.0, 1.5, 2.0])
+        assert_equal(a, -1.0)
+        b = forcefields.square_well(0.5, [1.0, 2.0, 1.25])
+        assert_equal(b, float('inf'))
+        c = forcefields.square_well(3.0, [0.5, 1.5, 1.25])
+        assert_equal(c, 0)
+        d = forcefields.square_well([2.0, 0.5], [1.0, 1.5, 2.0])
+        assert_equal(d, [-1.0, float('inf')])
+        e = forcefields.square_well([3.0, 3.0, 0.25], [1.0, 1.5, 1.25])
+        assert_equal(e, [0, 0, float('inf')])
+        f = forcefields.square_well(
+            [3.0, 3.0, 0.25], [1.0, 1.5, 1.25], max_val=5000)
+        assert_equal(f, [0, 0, 5000])
+
+    def test_square_well_force(self):
+        with self.assertRaises(ValueError):
+            forcefields.square_well(
+                2.0, [1.0, 1.5, 2.0], force=True)
+        with self.assertRaises(ValueError):
+            forcefields.square_well(
+                [2.0], [1.0, 1.5, 2.0], force=True)
+
+
+if __name__ == '__main__':
+    unittest.main(exit=False)