Allow epsilon and/or sigma to be zero in FENE bonds (#63)

* Added requested changes to azplugins/BondEvaluatorFENE.h * Removed error check on epsilon and sigma * Removed test which caused error when eps=0 and sigma=0 * Added new tests for sigma=0 && epsilon=0, sigma=0 && epsilon=nonzero , sigma=nonzero && epsilon=zero * made changes to calculated test values U and F for the new three cases * Added changes to azplugins/test-py/test_bond_fene.py * changed comments and value in azplugins/test-py/test_bond_fene.py * Added the suggested changes in azplugins/BondEvaluatorFENE.h * made requested changes in azplugins/bond.py * removed test for k=0 and added test to check force and energy when k=0,3 cases * changed a mistake in documentation * added a check for force and energy when k=0,delta=0,sigma=0 * Apply suggestions from code review Co-authored-by: Michael Howard <[email protected]>
mphowardlab · Mar 5, 2022 · 44411ca · 44411ca
1 parent 239cd81
commit 44411ca
Show file tree

Hide file tree

Showing 3 changed files with 121 additions and 28 deletions.
diff --git a/azplugins/BondEvaluatorFENE.h b/azplugins/BondEvaluatorFENE.h
@@ -108,7 +108,7 @@ class BondEvaluatorFENE
         DEVICE bool evalForceAndEnergy(Scalar& force_divr, Scalar& bond_eng)
             {
             // check for invalid parameters
-            if (lj1 == Scalar(0.0) || r_0 == Scalar(0.0) || K == Scalar(0.0) ) return false;
+            if (r_0 == Scalar(0.0)) return false;
 
             // Check if bond length restriction is violated
             if (rsq >= r_0*r_0) return false;
@@ -118,13 +118,18 @@ class BondEvaluatorFENE
             Scalar sigma6inv = lj2/lj1;
 
             // wca cutoff: r < 2^(1/6)*sigma
-            // wca cutoff: 1 / r^6 > 1/((2^(1/6))^6 sigma^6)
-            if (r6inv > sigma6inv/Scalar(2.0))
+            // if epsilon or sigma is zero OR r is beyond cutoff, the force and energy are zero
+            if (lj1 != Scalar(0) && r6inv > sigma6inv/Scalar(2.0))
                 {
                 Scalar epsilon = lj2*lj2/Scalar(4.0)/lj1;
                 force_divr = r2inv * r6inv * (Scalar(12.0)*lj1*r6inv - Scalar(6.0)*lj2);
                 bond_eng = r6inv * (lj1*r6inv - lj2) + epsilon;
                 }
+            else
+                {
+                force_divr = 0;
+                bond_eng = 0;
+                }
             // Check if delta is non -zero,if non-zero sqrt of rsq is calculated
             if (delta != Scalar(0))
                 {

diff --git a/azplugins/bond.py b/azplugins/bond.py
@@ -120,7 +120,7 @@ class fene(hoomd.md.bond._bond):
     - :math:`r_0` - size parameter ``r0`` (in distance units)
     - :math:`\varepsilon` - repulsive force strength ``epsilon`` (in energy units)
     - :math:`\sigma` - repulsive force interaction distance ``sigma`` (in distance units)
-    - :math:`\delta` - extra shift parameter for FENE bonds ``delta`` (in distance units),default value of  is ``delta`` zero
+    - :math:`\delta` - bond-length shift ``delta`` (in distance units), default value is  zero
 
     Examples::
 
@@ -161,15 +161,6 @@ def process_coeff(self, coeff):
         lj1 = 4.0 * epsilon * math.pow(sigma, 12.0)
         lj2 = 4.0 * epsilon * math.pow(sigma, 6.0)
 
-        if epsilon==0:
-            hoomd.context.msg.error("azplugins.bond.fene(): epsilon must be non-zero.\n")
-            raise ValueError('epsilon must be non-zero')
-        if sigma==0:
-            hoomd.context.msg.error("azplugins.bond.fene(): sigma must be non-zero.\n")
-            raise ValueError('sigma must be non-zero')
-        if k==0:
-            hoomd.context.msg.error("azplugins.bond.fene(): k must be non-zero.\n")
-            raise ValueError('k must be non-zero')
         if r0==0:
             hoomd.context.msg.error("azplugins.bond.fene(): r0 must be non-zero.\n")
             raise ValueError('r0 must be non-zero')

diff --git a/azplugins/test-py/test_bond_fene.py b/azplugins/test-py/test_bond_fene.py
@@ -52,21 +52,6 @@ def test_set_missing_r0(self):
         self.assertRaises(RuntimeError, fene.update_coeffs)
 
      # test coefficients = 0
-    def test_set_zero_epsilon(self):
-        fene = azplugins.bond.fene()
-        fene.bond_coeff.set('bond', epsilon=0, sigma=1.0, k=30,r0=1.5)
-        self.assertRaises(ValueError, fene.update_coeffs)
-
-    def test_set_zero_sigma(self):
-        fene = azplugins.bond.fene()
-        fene.bond_coeff.set('bond', epsilon=1.0, sigma=0.0, k=30, r0=1.5)
-        self.assertRaises(ValueError, fene.update_coeffs)
-
-    def test_set_zero_k(self):
-        fene = azplugins.bond.fene()
-        fene.bond_coeff.set('bond', epsilon=1.0, sigma=1.0, k=0, r0=1.5)
-        self.assertRaises(ValueError, fene.update_coeffs)
-
     def test_set_zero_r0(self):
         fene = azplugins.bond.fene()
         fene.bond_coeff.set('bond', epsilon=1.0, sigma=1.0, k=30, r0=0)
@@ -144,6 +129,118 @@ def test_potential_delta_nonzero(self):
         self.assertAlmostEqual(f1[1],0)
         self.assertAlmostEqual(f1[2],0)
 
+    # test the calculation of foce and potential when sigma=0 and epsilon=0
+    def test_potential_sigma_zero_epsilon_zero(self):
+        fene = azplugins.bond.fene()
+        fene.bond_coeff.set('bond', epsilon=0.0, sigma=0.0, k=30,r0=1.5,delta=1.8)
+
+        md.integrate.mode_standard(dt=0)
+        nve = md.integrate.nve(group = hoomd.group.all())
+        hoomd.run(1)
+        #values of F and  U are caluclated using a calculator, by substituting
+        #r0=1.5,delta=1.8,sigma=0.0,epsilon=0.0, with r=1.0
+        F = 33.540372671
+        U = 11.29599912562
+        f0 = fene.forces[0].force
+        f1 = fene.forces[1].force
+        e0 = fene.forces[0].energy
+        e1 = fene.forces[1].energy
+
+        self.assertAlmostEqual(e0,0.5*U,3)
+        self.assertAlmostEqual(e1,0.5*U,3)
+
+        self.assertAlmostEqual(f0[0],-F,3)
+        self.assertAlmostEqual(f0[1],0)
+        self.assertAlmostEqual(f0[2],0)
+
+        self.assertAlmostEqual(f1[0],F,3)
+        self.assertAlmostEqual(f1[1],0)
+        self.assertAlmostEqual(f1[2],0)
+
+    # test the calculation of foce and potential when sigma=nonzero and epsilon=0
+    def test_potential_sigma_nonzero_epsilon_zero(self):
+        fene = azplugins.bond.fene()
+        fene.bond_coeff.set('bond', epsilon=0.0, sigma=1.0, k=30,r0=1.5,delta=1.8)
+
+        md.integrate.mode_standard(dt=0)
+        nve = md.integrate.nve(group = hoomd.group.all())
+        hoomd.run(1)
+        #values of F and  U are caluclated using a calculator, by substituting
+        #r0=1.5,delta=1.8,sigma=1.0,epsilon=0.0, with r=1.0
+        F = 33.540372671
+        U = 11.29599912562
+        f0 = fene.forces[0].force
+        f1 = fene.forces[1].force
+        e0 = fene.forces[0].energy
+        e1 = fene.forces[1].energy
+
+        self.assertAlmostEqual(e0,0.5*U,3)
+        self.assertAlmostEqual(e1,0.5*U,3)
+
+        self.assertAlmostEqual(f0[0],-F,3)
+        self.assertAlmostEqual(f0[1],0)
+        self.assertAlmostEqual(f0[2],0)
+
+        self.assertAlmostEqual(f1[0],F,3)
+        self.assertAlmostEqual(f1[1],0)
+        self.assertAlmostEqual(f1[2],0)
+
+    # test the calculation of foce and potential when sigma=0 and epsilon=nonzero
+    def test_potential_sigma_zero_epsilon_nonzero(self):
+        fene = azplugins.bond.fene()
+        fene.bond_coeff.set('bond', epsilon=1.0, sigma=0.0, k=30,r0=1.5,delta=1.8)
+
+        md.integrate.mode_standard(dt=0)
+        nve = md.integrate.nve(group = hoomd.group.all())
+        hoomd.run(1)
+        #values of F and  U are caluclated using a calculator, by substituting
+        #r0=1.5,delta=1.8,sigma=0.0,epsilon=1.0, with r=1.0
+        F = 33.540372671
+        U = 11.29599912562 #no contribution from WCA as 2^1/6 sigma < r
+        f0 = fene.forces[0].force
+        f1 = fene.forces[1].force
+        e0 = fene.forces[0].energy
+        e1 = fene.forces[1].energy
+
+        self.assertAlmostEqual(e0,0.5*U,3)
+        self.assertAlmostEqual(e1,0.5*U,3)
+
+        self.assertAlmostEqual(f0[0],-F,3)
+        self.assertAlmostEqual(f0[1],0)
+        self.assertAlmostEqual(f0[2],0)
+
+        self.assertAlmostEqual(f1[0],F,3)
+        self.assertAlmostEqual(f1[1],0)
+        self.assertAlmostEqual(f1[2],0)
+
+    # test the calculation of foce and potential when k=0 and epsilon=nonzero sigma=nonzero
+    def test_potential_k_zero(self):
+        fene = azplugins.bond.fene()
+        fene.bond_coeff.set('bond', epsilon=1.0, sigma=1.0, k=0,r0=1.5,delta=1.8)
+
+        md.integrate.mode_standard(dt=0)
+        nve = md.integrate.nve(group = hoomd.group.all())
+        hoomd.run(1)
+        #values of F and  U are caluclated using a calculator, by substituting
+        #k=0,r0=1.5,delta=1.8,sigma=1.0,epsilon=1.0, with r=1.0
+        F = 24
+        U = 1               #no contribution from FENE bonds as k=0
+        f0 = fene.forces[0].force
+        f1 = fene.forces[1].force
+        e0 = fene.forces[0].energy
+        e1 = fene.forces[1].energy
+
+        self.assertAlmostEqual(e0,0.5*U,3)
+        self.assertAlmostEqual(e1,0.5*U,3)
+
+        self.assertAlmostEqual(f0[0],-F,3)
+        self.assertAlmostEqual(f0[1],0)
+        self.assertAlmostEqual(f0[2],0)
+
+        self.assertAlmostEqual(f1[0],F,3)
+        self.assertAlmostEqual(f1[1],0)
+        self.assertAlmostEqual(f1[2],0)
+
     # test streching beyond maximal bond length of r0
     def test_potential_strech_beyond_r0(self):
         fene = azplugins.bond.fene()