Do not add constraints involving extra particles (#3506)

* Do not add constraints involving extra particles

* Added test case

wrappers/python/openmm/app/internal/amber_file_parser.py

@@ -768,15 +768,18 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
 
     # Add constraints.
     isWater = [prmtop.getResidueLabel(i) in ('WAT', 'HOH', 'TP4', 'TP5', 'T4E') for i in range(prmtop.getNumAtoms())]
+    isEP = [a.element is None for a in topology.atoms()]
     if shake in ('h-bonds', 'all-bonds', 'h-angles'):
         for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
-            system.addConstraint(iAtom, jAtom, rMin)
+            if not (isEP[iAtom] or isEP[jAtom]):
+                system.addConstraint(iAtom, jAtom, rMin)
     if shake in ('all-bonds', 'h-angles'):
         for (iAtom, jAtom, k, rMin) in prmtop.getBondsNoH():
-            system.addConstraint(iAtom, jAtom, rMin)
+            if not (isEP[iAtom] or isEP[jAtom]):
+                system.addConstraint(iAtom, jAtom, rMin)
     if rigidWater and shake is None:
         for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
-            if isWater[iAtom] and isWater[jAtom]:
+            if isWater[iAtom] and isWater[jAtom] and not (isEP[iAtom] or isEP[jAtom]):
                 system.addConstraint(iAtom, jAtom, rMin)
 
     # Add harmonic bonds.

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -448,5 +448,14 @@ class TestAmberPrmtopFile(unittest.TestCase):
                 OpenMM_CMAP_E = simulation.context.getState(getEnergy=True, groups=1<<i).getPotentialEnergy().value_in_unit(kilojoules_per_mole)/conversion
                 self.assertAlmostEqual(OpenMM_CMAP_E, sander_CMAP_E, places=4)
 
+    def testEPConstraints(self):
+        """Test different types of constraints when using extra particles"""
+        prmtop = AmberPrmtopFile('systems/peptide_with_tip4p.prmtop')
+        for constraints in (HBonds, AllBonds):
+            system = prmtop.createSystem(constraints=constraints)
+            integrator = VerletIntegrator(0.001*picoseconds)
+            # If a constraint was added to a massless particle, this will throw an exception.
+            context = Context(system, integrator, Platform.getPlatformByName('Reference'))
+
 if __name__ == '__main__':
     unittest.main()

wrappers/python/tests/systems/peptide_with_tip4p.inpcrd

+default_name
+    85
+   6.1027258   9.4114116   5.1039335   6.5815895  10.1277022   5.6290157
+   6.6678073   9.1230379   4.3177487   5.2086444   9.7657205   4.7621357
+   5.8010133   8.2419298   5.9556322   5.2338195   8.5904194   6.8193516
+   4.9060379   7.2614871   5.1978254   5.2290030   7.1815685   4.1573041
+   4.9489259   6.2745706   5.6619402   3.4965842   7.8237994   5.2678187
+   3.3776621   8.9780360   4.7980651   2.6472232   7.1669020   5.8997283
+   7.0376422   7.5370608   6.4781338   7.9157036   7.1622066   5.7047041
+   7.0652775   7.3158858   7.7952036   6.3169315   7.6633401   8.3789911
+   8.0764491   6.4982978   8.4537730   9.0634999   6.8670277   8.1688606
+   7.9197212   6.6666247   9.9695324   8.6690479   6.0648171  10.4845540
+   8.0657682   7.7128438  10.2415669   6.9244697   6.3515483  10.2875817
+   7.9931502   5.0209894   8.0274701   9.0297108   4.4106384   7.7836349
+   6.9039267   4.4598555   7.9317597   4.8638450  11.1831355   9.8736302
+   4.1346516  10.5621898   9.8644378   4.8794334  11.5468847   8.9890396
+   4.7726236  11.1496591   9.7590789   2.6329636   2.9920417   7.3605727
+   2.8338815   2.3735456   6.6584157   3.4695441   3.4182530   7.5500773
+   2.7657541   2.9679551   7.2948118  11.6032547   4.2363825   3.6652327
+  11.0489985   4.5413440   2.9475194  11.9080957   3.3741621   3.3823900
+  11.5709231   4.1649390   3.5372464   2.6699140   9.6787405   9.4896712
+   2.7570939   9.4003152   8.5789176   2.7420827   8.8695924   9.9959596
+   2.6901212   9.5395278   9.4380524   4.2027790   3.2218855   3.8639297
+   3.7166501   3.1773864   3.0401503   5.0820834   2.9173323   3.6390697
+   4.2524311   3.1769782   3.7295794   0.5083146   4.1661638   6.2786993
+   1.2617567   3.6970865   6.6372025   0.8431777   4.5870678   5.4874468
+   0.6474560   4.1596318   6.2235450   5.4377311   1.1104254   9.6028083
+   5.7316024   1.7256555   8.9303498   4.4874126   1.2190194   9.6268500
+   5.3534380   1.2030977   9.5200769  12.5431809   8.9055182   5.2286458
+  11.6760008   9.0337082   4.8468081  12.3913378   8.3331541   5.9810074
+  12.4132771   8.8479552   5.2760220  10.2580286  10.2519211   7.1265204
+  10.2459042  10.2237519   6.1698049  10.7147126  11.0667847   7.3358240
+  10.3146089  10.3527584   7.0310610   3.0018903   7.5970825  11.3670397
+   2.8021272   6.7560910  10.9540893   2.7438148   7.4819565  12.2820359
+   2.9430006   7.4746080  11.4306793   6.9630906  12.5058599   7.3124894
+   6.0324019  12.6912046   7.4397686   7.3972580  13.3452184   7.4595676
+   6.8995820  12.6364993   7.3478448   2.2680782   5.1300381  10.2420328
+   1.8500766   5.0630853  11.1004605   1.5463903   5.2949704   9.6353352
+   2.1225859   5.1422855  10.2745597  11.0871036   7.0226771   2.3181943
+  10.9127438   7.6505629   1.6181586  10.4907008   7.2778323   3.0231798
+  10.9883767   7.1353536   2.3196085   8.9439794   4.4033560   2.5175984
+   8.7517218   4.6458555   1.6110875   8.3729800   3.6546287   2.6915467
+   8.8458299   4.3392611   2.4235532  11.9946982   8.0588113   7.7954434
+  12.4305976   8.4597111   8.5478050  11.2851347   8.6646517   7.5790687
+  11.9594798   8.1878177   7.8647399
+  13.0633205  13.0633205  13.0633205 109.4712190 109.4712190 109.4712190