Merge pull request #491 from jchodera/fix-amoeba

Fix error in setting up AMOEBA systems when useDispersionCorrection is specified (retry)

plugins/amoeba/openmmapi/include/openmm/AmoebaMultipoleForce.h

@@ -126,7 +126,7 @@ public:
      *
      * @return the cutoff distance, measured in nm
      */
-    double getCutoffDistance(void) const;
+    double getCutoffDistance() const;
 
     /**
      * Set the cutoff distance (in nm) being used for nonbonded interactions.  If the NonbondedMethod in use

plugins/amoeba/openmmapi/include/openmm/AmoebaVdwForce.h

@@ -191,7 +191,7 @@ public:
     /**
      * Get the cutoff distance.
      */
-    double getCutoff(void) const;
+    double getCutoff() const;
 
     /**
      * Get the method used for handling long range nonbonded interactions.

plugins/amoeba/openmmapi/src/AmoebaMultipoleForce.cpp

@@ -61,7 +61,7 @@ void AmoebaMultipoleForce::setPolarizationType( AmoebaMultipoleForce::Polarizati
     polarizationType = type;
 }
 
-double AmoebaMultipoleForce::getCutoffDistance( void ) const {
+double AmoebaMultipoleForce::getCutoffDistance() const {
     return cutoffDistance;
 }
 

plugins/amoeba/openmmapi/src/AmoebaVdwForce.cpp

@@ -106,7 +106,7 @@ void AmoebaVdwForce::setCutoff( double inputCutoff ){
     cutoff = inputCutoff;
 }
 
-double AmoebaVdwForce::getCutoff( void ) const {
+double AmoebaVdwForce::getCutoff() const {
     return cutoff;
 }
 

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

@@ -1107,6 +1107,8 @@ class NonbondedGenerator:
         force.setCutoffDistance(nonbondedCutoff)
         if 'ewaldErrorTolerance' in args:
             force.setEwaldErrorTolerance(args['ewaldErrorTolerance'])
+        if 'useDispersionCorrection' in args:
+            force.setUseDispersionCorrection(bool(args['useDispersionCorrection']))
         sys.addForce(force)
 
     def postprocessSystem(self, sys, data, args):
@@ -3037,7 +3039,7 @@ class AmoebaVdwGenerator:
             # dispersion correction
 
             if ('useDispersionCorrection' in args):
-                force.setUseDispersionCorrection(int(args['useDispersionCorrection']))
+                force.setUseDispersionCorrection(bool(args['useDispersionCorrection']))
 
             if (nonbondedMethod == PME):
                 force.setNonbondedMethod(mm.AmoebaVdwForce.CutoffPeriodic)

wrappers/python/tests/TestForceField.py

@@ -19,6 +19,7 @@ class TestForceField(unittest.TestCase):
         self.topology1 = self.pdb1.topology
         self.topology1.setUnitCellDimensions(Vec3(2, 2, 2))
 
+        # alalnine dipeptide with implicit water
         self.pdb2 = PDBFile('systems/alanine-dipeptide-implicit.pdb')
         self.forcefield2 = ForceField('amber99sb.xml', 'amber99_obc.xml')
 
@@ -38,6 +39,18 @@ class TestForceField(unittest.TestCase):
                                 f.getNonbondedMethod()==methodMap[method] 
                                 for f in forces))
 
+    def test_DispersionCorrection(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        for useDispersionCorrection in [True, False]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedCutoff=2*nanometer, 
+                                                   useDispersionCorrection=useDispersionCorrection)
+
+            for force in system.getForces():
+                if isinstance(force, NonbondedForce):
+                    self.assertEqual(useDispersionCorrection, force.getUseDispersionCorrection())
+
     def test_Cutoff(self):
         """Test to make sure the nonbondedCutoff parameter is passed correctly."""
 
@@ -120,6 +133,61 @@ class TestForceField(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
 
+class AmoebaTestForceField(unittest.TestCase):
+    """Test the ForceField.createSystem() method with the AMOEBA forcefield."""
+ 
+    def setUp(self):
+        """Set up the tests by loading the input pdb files and force field 
+        xml files.
+
+        """
+
+        self.pdb1 = PDBFile('systems/amoeba-ion-in-water.pdb')
+        self.forcefield1 = ForceField('amoeba2009.xml')
+        self.topology1 = self.pdb1.topology
+
+
+    def test_NonbondedMethod(self):
+        """Test all five options for the nonbondedMethod parameter."""
+
+        methodMap = {NoCutoff:AmoebaMultipoleForce.NoCutoff,
+                     PME:AmoebaMultipoleForce.PME}
+
+        for method in methodMap:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                  nonbondedMethod=method)
+            forces = system.getForces()
+            self.assertTrue(any(isinstance(f, AmoebaMultipoleForce) and
+                                f.getNonbondedMethod()==methodMap[method]
+                                for f in forces))
+    def test_Cutoff(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        cutoff_distance = 0.7*nanometer
+        for method in [NoCutoff, PME]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedMethod=method,
+                                                   nonbondedCutoff=cutoff_distance,
+                                                   constraints=None)
+
+            for force in system.getForces():
+                if isinstance(force, AmoebaVdwForce):
+                    self.assertEqual(force.getCutoff(), cutoff_distance)
+                if isinstance(force, AmoebaMultipoleForce):
+                    self.assertEqual(force.getCutoffDistance(), cutoff_distance)
+
+    def test_DispersionCorrection(self):
+        """Test to make sure the nonbondedCutoff parameter is passed correctly."""
+
+        for useDispersionCorrection in [True, False]:
+            system = self.forcefield1.createSystem(self.pdb1.topology,
+                                                   nonbondedMethod=PME,
+                                                   useDispersionCorrection=useDispersionCorrection)
+
+            for force in system.getForces():
+                if isinstance(force, AmoebaVdwForce):
+                    self.assertEqual(useDispersionCorrection, force.getUseDispersionCorrection())
+
 
 if __name__ == '__main__':
     unittest.main()