Unverified Commit 7432ad8c authored by Peter Eastman's avatar Peter Eastman Committed by GitHub
Browse files

Ensure bond lengths are set before creating forces that depend on them (#3493)

* Ensure bond lengths are set before creating forces that depend on them

* Minor formatting
parent a7e6f0a9
...@@ -2057,6 +2057,11 @@ class HarmonicAngleGenerator(object): ...@@ -2057,6 +2057,11 @@ class HarmonicAngleGenerator(object):
generator.registerAngle(angle.attrib) generator.registerAngle(angle.attrib)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
pass
def postprocessSystem(self, sys, data, args):
# We need to wait until after all bonds have been added so their lengths will be set correctly.
existing = [f for f in sys.getForces() if type(f) == mm.HarmonicAngleForce] existing = [f for f in sys.getForces() if type(f) == mm.HarmonicAngleForce]
if len(existing) == 0: if len(existing) == 0:
force = mm.HarmonicAngleForce() force = mm.HarmonicAngleForce()
...@@ -3733,8 +3738,11 @@ class AmoebaOutOfPlaneBendGenerator(object): ...@@ -3733,8 +3738,11 @@ class AmoebaOutOfPlaneBendGenerator(object):
#============================================================================================= #=============================================================================================
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args): def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
self._nonbondedMethod = nonbondedMethod
self._nonbondedCutoff = nonbondedCutoff
# get force def postprocessSystem(self, sys, data, args):
# We need to wait until after all bonds have been added so their lengths will be set correctly.
energy = """k*(theta^2 + %s*theta^3 + %s*theta^4 + %s*theta^5 + %s*theta^6); energy = """k*(theta^2 + %s*theta^3 + %s*theta^4 + %s*theta^5 + %s*theta^6);
theta = %.15g*pointangle(x2, y2, z2, x4, y4, z4, projx, projy, projz); theta = %.15g*pointangle(x2, y2, z2, x4, y4, z4, projx, projy, projz);
...@@ -3857,12 +3865,12 @@ class AmoebaOutOfPlaneBendGenerator(object): ...@@ -3857,12 +3865,12 @@ class AmoebaOutOfPlaneBendGenerator(object):
for force in self.forceField._forces: for force in self.forceField._forces:
if (force.__class__.__name__ == 'AmoebaAngleGenerator'): if (force.__class__.__name__ == 'AmoebaAngleGenerator'):
force.createForcePostOpBendAngle(sys, data, nonbondedMethod, nonbondedCutoff, angles, args) force.createForcePostOpBendAngle(sys, data, self._nonbondedMethod, self._nonbondedCutoff, angles, args)
force.createForcePostOpBendInPlaneAngle(sys, data, nonbondedMethod, nonbondedCutoff, angles, args) force.createForcePostOpBendInPlaneAngle(sys, data, self._nonbondedMethod, self._nonbondedCutoff, angles, args)
for force in self.forceField._forces: for force in self.forceField._forces:
if (force.__class__.__name__ == 'AmoebaStretchBendGenerator'): if (force.__class__.__name__ == 'AmoebaStretchBendGenerator'):
force.createForcePostAmoebaBondForce(sys, data, nonbondedMethod, nonbondedCutoff, angles, args) force.createForcePostAmoebaBondForce(sys, data, self._nonbondedMethod, self._nonbondedCutoff, angles, args)
parsers["AmoebaOutOfPlaneBendForce"] = AmoebaOutOfPlaneBendGenerator.parseElement parsers["AmoebaOutOfPlaneBendForce"] = AmoebaOutOfPlaneBendGenerator.parseElement
...@@ -4672,7 +4680,6 @@ class AmoebaStretchBendGenerator(object): ...@@ -4672,7 +4680,6 @@ class AmoebaStretchBendGenerator(object):
k1, k2 = self.k1[i], self.k2[i] k1, k2 = self.k1[i], self.k2[i]
else: else:
k1, k2 = self.k2[i], self.k1[i] k1, k2 = self.k2[i], self.k1[i]
force.addBond((angle[0], angle[1], angle[2]), (bondAB, bondCB, angleDict['idealAngle']/radian, k1, k2)) force.addBond((angle[0], angle[1], angle[2]), (bondAB, bondCB, angleDict['idealAngle']/radian, k1, k2))
break break
......
...@@ -1020,7 +1020,7 @@ class TestForceField(unittest.TestCase): ...@@ -1020,7 +1020,7 @@ class TestForceField(unittest.TestCase):
system1 = ff1.createSystem(pdb.topology) system1 = ff1.createSystem(pdb.topology)
system2 = ff2.createSystem(pdb.topology) system2 = ff2.createSystem(pdb.topology)
imp1 = system1.getForce(2).getTorsionParameters(158) imp1 = system1.getForce(1).getTorsionParameters(158)
imp2 = system2.getForce(0).getTorsionParameters(158) imp2 = system2.getForce(0).getTorsionParameters(158)
system1_indexes = [imp1[0], imp1[1], imp1[2], imp1[3]] system1_indexes = [imp1[0], imp1[1], imp1[2], imp1[3]]
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment