Merge remote-tracking branch 'origin/master' into qc

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

@@ -111,7 +111,7 @@ class GromacsTopFile(object):
                 if len(fields) < 2:
                     raise ValueError('Illegal line in .top file: '+line)
                 name = fields[1]
-                valueStart = stripped.find(name, len(command))+len(name)
+                valueStart = stripped.find(name, len(command))+len(name)+1
                 value = line[valueStart:].strip()
                 self._defines[name] = value
             elif command == '#ifdef':

wrappers/python/simtk/openmm/app/internal/charmm/topologyobjects.py

@@ -12,7 +12,7 @@ Copyright (c) 2014 the Authors
 
 Author: Jason M. Swails
 Contributors:
-Date: April 18, 2014
+Date: July 3, 2014
 
 Permission is hereby granted, free of charge, to any person obtaining a
 copy of this software and associated documentation files (the "Software"),
@@ -501,7 +501,7 @@ class Angle(object):
         """ See if a bond or an atom is in this angle """
         if isinstance(thing, Bond):
             return self.atom2 in thing and (self.atom1 in thing or
-                                            self.atom2 in thing)
+                                            self.atom3 in thing)
         # Otherwise assume it's an atom
         return self.atom1 is thing or self.atom2 is thing or self.atom3 is thing
 

wrappers/python/tests/TestBytes.py

@@ -8,7 +8,8 @@ class TestBytes(unittest.TestCase):
         system.addParticle(1.0)
         refPositions = [(0,0,0)]
 
-        context = mm.Context(system, mm.VerletIntegrator(0))
+        platform = mm.Platform.getPlatformByName('Reference')
+        context = mm.Context(system, mm.VerletIntegrator(0), platform)
         context.setPositions(refPositions)
         chk = context.createCheckpoint()
         # check that the return value of createCheckpoint is of type bytes (non-unicode)

wrappers/python/tests/TestForceField.py

@@ -4,6 +4,7 @@ from simtk.openmm.app import *
 from simtk.openmm import *
 from simtk.unit import *
 import simtk.openmm.app.element as elem
+import simtk.openmm.app.forcefield as forcefield
 
 class TestForceField(unittest.TestCase):
     """Test the ForceField.createSystem() method."""
@@ -19,6 +20,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 +40,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 +134,113 @@ class TestForceField(unittest.TestCase):
         totalMass2 = sum([system2.getParticleMass(i) for i in range(system2.getNumParticles())]).value_in_unit(amu)
         self.assertAlmostEqual(totalMass1, totalMass2)
     
+    def test_Forces(self):
+        """Compute forces and compare them to ones generated with a previous version of OpenMM to ensure they haven't changed."""
+        
+        pdb = PDBFile('systems/lysozyme-implicit.pdb')
+        system = self.forcefield2.createSystem(pdb.topology)
+        integrator = VerletIntegrator(0.001)
+        context = Context(system, integrator)
+        context.setPositions(pdb.positions)
+        state1 = context.getState(getForces=True)
+        state2 = XmlSerializer.deserialize(open('systems/lysozyme-implicit-forces.xml').read())
+        numDifferences = 0
+        for f1, f2, in zip(state1.getForces().value_in_unit(kilojoules_per_mole/nanometer), state2.getForces().value_in_unit(kilojoules_per_mole/nanometer)):
+            diff = norm(f1-f2)
+            if diff > 0.1 and diff/norm(f1) > 1e-3:
+                numDifferences += 1
+        self.assertTrue(numDifferences < system.getNumParticles()/20) # Tolerate occasional differences from numerical error
+    
+    def test_ProgrammaticForceField(self):
+        """Test building a ForceField programmatically."""
+        
+        # Build the ForceField for TIP3P programmatically.
+        ff = ForceField()
+        ff.registerAtomType({'name':'tip3p-O', 'class':'OW', 'mass':15.99943*daltons, 'element':elem.oxygen})
+        ff.registerAtomType({'name':'tip3p-H', 'class':'HW', 'mass':1.007947*daltons, 'element':elem.hydrogen})
+        residue = ForceField._TemplateData('HOH')
+        residue.atoms.append(ForceField._TemplateAtomData('O', 'tip3p-O', elem.oxygen))
+        residue.atoms.append(ForceField._TemplateAtomData('H1', 'tip3p-H', elem.hydrogen))
+        residue.atoms.append(ForceField._TemplateAtomData('H2', 'tip3p-H', elem.hydrogen))
+        residue.addBond(0, 1)
+        residue.addBond(0, 2)
+        ff.registerResidueTemplate(residue)
+        bonds = forcefield.HarmonicBondGenerator(ff)
+        bonds.registerBond({'class1':'OW', 'class2':'HW', 'length':0.09572*nanometers, 'k':462750.4*kilojoules_per_mole/nanometer})
+        ff.registerGenerator(bonds)
+        angles = forcefield.HarmonicAngleGenerator(ff)
+        angles.registerAngle({'class1':'HW', 'class2':'OW', 'class3':'HW', 'angle':1.82421813418*radians, 'k':836.8*kilojoules_per_mole/radian})
+        ff.registerGenerator(angles)
+        nonbonded = forcefield.NonbondedGenerator(ff, 0.833333, 0.5)
+        nonbonded.registerAtom({'type':'tip3p-O', 'charge':-0.834, 'sigma':0.31507524065751241*nanometers, 'epsilon':0.635968*kilojoules_per_mole})
+        nonbonded.registerAtom({'type':'tip3p-H', 'charge':0.417, 'sigma':1*nanometers, 'epsilon':0*kilojoules_per_mole})
+        ff.registerGenerator(nonbonded)
+        
+        # Build a water box.
+        modeller = Modeller(Topology(), [])
+        modeller.addSolvent(ff, boxSize=Vec3(3, 3, 3)*nanometers)
+        
+        # Create a system using the programmatic force field as well as one from an XML file.
+        system1 = ff.createSystem(modeller.topology)
+        ff2 = ForceField('tip3p.xml')
+        system2 = ff2.createSystem(modeller.topology)
+        self.assertEqual(XmlSerializer.serialize(system1), XmlSerializer.serialize(system2))
+
+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()