From 44411ca700b25fd68a47f025cddd21bd02afe2dd Mon Sep 17 00:00:00 2001 From: judevishnu <100363618+judevishnu@users.noreply.github.com> Date: Sat, 5 Mar 2022 23:16:08 +0100 Subject: [PATCH] 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 --- azplugins/BondEvaluatorFENE.h | 11 ++- azplugins/bond.py | 11 +-- azplugins/test-py/test_bond_fene.py | 127 ++++++++++++++++++++++++---- 3 files changed, 121 insertions(+), 28 deletions(-) diff --git a/azplugins/BondEvaluatorFENE.h b/azplugins/BondEvaluatorFENE.h index d1d4e2fc..e1703e52 100644 --- 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 index ac54cd08..9b4acc1b 100644 --- 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 index 94637e5a..90b4eb6a 100644 --- 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()