Fixed bug that caused ForceField to generate an incorrect System if used multiple times

wrappers/python/simtk/openmm/app/forcefield.py

@@ -1435,7 +1435,6 @@ class AmoebaBondGenerator:
         self.types2 = []
         self.length = []
         self.k = []
-        self.hasBeenCalled = 0
     
     #=============================================================================================
 
@@ -1464,11 +1463,6 @@ class AmoebaBondGenerator:
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
-        if (self.hasBeenCalled):
-             return
-
-        self.hasBeenCalled = 1
-
         #countConstraint(data)
 
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
@@ -1554,8 +1548,6 @@ class AmoebaAngleGenerator:
 
         self.angle = []
         self.k = []
-
-        self.hasBeenCalled = 0
     
     #=============================================================================================
 
@@ -1610,10 +1602,6 @@ class AmoebaAngleGenerator:
     
     def createForcePostOpBendAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
 
-        if (self.hasBeenCalled):
-             return
-        self.hasBeenCalled += 1
-
         # get force
 
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
@@ -1688,8 +1676,6 @@ class AmoebaAngleGenerator:
     
     def createForcePostOpBendInPlaneAngle(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
 
-        self.hasBeenCalled += 1
-
         # get force
 
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
@@ -1773,7 +1759,6 @@ class AmoebaOutOfPlaneBendGenerator:
         self.types4 = []
 
         self.ks = []
-        self.hasBeenCalled = 0
 
     #=============================================================================================
     # Local version of findAtomTypes needed since class indices are 0 (i.e., not recognized)
@@ -1864,10 +1849,6 @@ class AmoebaOutOfPlaneBendGenerator:
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
-        if (self.hasBeenCalled):
-             return
-        self.hasBeenCalled = 1
-
         # get force
 
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
@@ -2017,8 +1998,6 @@ class AmoebaOutOfPlaneBendGenerator:
             if (force.__class__.__name__ == 'AmoebaAngleGenerator'): 
                 force.createForcePostOpBendAngle(sys, data, nonbondedMethod, nonbondedCutoff, nonInPlaneAngles, args)
                 force.createForcePostOpBendInPlaneAngle(sys, data, nonbondedMethod, nonbondedCutoff, inPlaneAngles, args)
-            if (force.__class__.__name__ == 'AmoebaBondGenerator'): 
-                force.createForce(sys, data, nonbondedMethod, nonbondedCutoff, args)
 
         for force in self.forceField._forces:
             if (force.__class__.__name__ == 'AmoebaStretchBendGenerator'): 
@@ -2049,7 +2028,6 @@ class AmoebaTorsionGenerator:
         self.t1 = []
         self.t2 = []
         self.t3 = []
-        self.hasBeenCalled = 0
     
     #=============================================================================================
 
@@ -2103,10 +2081,6 @@ class AmoebaTorsionGenerator:
 
     def createForce(self, sys, data, nontorsionedMethod, nontorsionedCutoff, args):
 
-        if (self.hasBeenCalled):
-             return
-        self.hasBeenCalled = 1
-
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
         existing = [f for f in existing if type(f) == mm.PeriodicTorsionForce]
         if len(existing) == 0:
@@ -2169,7 +2143,6 @@ class AmoebaPiTorsionGenerator:
         self.types1 = []
         self.types2 = []
         self.k = []
-        self.hasBeenCalled = 0
     
     #=============================================================================================
 
@@ -2198,10 +2171,6 @@ class AmoebaPiTorsionGenerator:
 
     def createForce(self, sys, data, nonpiTorsionedMethod, nonpiTorsionedCutoff, args):
 
-        if (self.hasBeenCalled):
-             return
-        self.hasBeenCalled = 1
-
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
         existing = [f for f in existing if type(f) == mm.AmoebaPiTorsionForce]
 
@@ -2292,8 +2261,6 @@ class AmoebaTorsionTorsionGenerator:
         self.gridIndex = []
 
         self.grids = []
-
-        self.hasBeenCalled = 0
     
     #=============================================================================================
 
@@ -2443,10 +2410,6 @@ class AmoebaTorsionTorsionGenerator:
  
     def createForce(self, sys, data, nonpiTorsionedMethod, nonpiTorsionedCutoff, args):
 
-        if (self.hasBeenCalled):
-             return
-        self.hasBeenCalled = 1
-
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
         existing = [f for f in existing if type(f) == mm.AmoebaTorsionTorsionForce]
 
@@ -2536,7 +2499,6 @@ class AmoebaStretchBendGenerator:
 
         self.k1 = []
         self.k2 = []
-        self.hasBeenCalled = 0
     
     #=============================================================================================
 
@@ -2589,10 +2551,6 @@ class AmoebaStretchBendGenerator:
 
     def createForcePostAmoebaBondForce(self, sys, data, nonbondedMethod, nonbondedCutoff, angleList, args):
 
-        if (self.hasBeenCalled):
-             return
-        self.hasBeenCalled = 1
-
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
         existing = [f for f in existing if type(f) == mm.AmoebaStretchBendForce]
         if len(existing) == 0:
@@ -2900,7 +2858,6 @@ class AmoebaMultipoleGenerator:
         self.polar15Scale = polar15Scale 
 
         self.typeMap = {}
-        self.hasBeenCalled = 0
 
     #=============================================================================================
     # Set axis type
@@ -3208,10 +3165,6 @@ class AmoebaMultipoleGenerator:
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
-        if (self.hasBeenCalled ):
-            return
-        self.hasBeenCalled = 1
-
         methodMap = {NoCutoff:mm.AmoebaMultipoleForce.NoCutoff,
                      PME:mm.AmoebaMultipoleForce.PME}
 
@@ -3638,7 +3591,6 @@ class AmoebaGeneralizedKirkwoodGenerator:
         self.includeCavityTerm = includeCavityTerm
         self.probeRadius = probeRadius
         self.surfaceAreaFactor = surfaceAreaFactor
-        self.hasBeenCalled = 0
 
         self.radiusTypeMap = {}
         self.radiusTypeMap['Bondi'] = {}
@@ -3820,9 +3772,6 @@ class AmoebaGeneralizedKirkwoodGenerator:
     
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
-        if (self.hasBeenCalled ):
-            return
-
         if( nonbondedMethod != NoCutoff ):
             raise ValueError( "Only the nonbondedMethod=NoCutoff option is available for implicit solvent simulations." )
 
@@ -3866,8 +3815,6 @@ class AmoebaGeneralizedKirkwoodGenerator:
         else:
             force = existing[0]
 
-        self.hasBeenCalled = 1
-
         # add particles to force
         # throw error if particle type not available 
 
@@ -3913,8 +3860,6 @@ class AmoebaUreyBradleyGenerator:
         self.length = []
         self.k = []
 
-        self.hasBeenCalled = 0
-    
     #=============================================================================================
 
     @staticmethod
@@ -3945,11 +3890,6 @@ class AmoebaUreyBradleyGenerator:
 
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
-        if (self.hasBeenCalled):
-             return
- 
-        self.hasBeenCalled = 1
-
         existing = [sys.getForce(i) for i in range(sys.getNumForces())]
         existing = [f for f in existing if type(f) == mm.HarmonicBondForce]