Skip to content

Commit a6cd8c2

Browse files
peastmanpeastman
authored
Do not add constraints involving extra particles (openmm#3506)
* Do not add constraints involving extra particles * Added test case
1 parent 4dc111b commit a6cd8c2

File tree

4 files changed

+492
-3
lines changed

4 files changed

+492
-3
lines changed

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

Lines changed: 6 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -768,15 +768,18 @@ def readAmberSystem(topology, prmtop_filename=None, prmtop_loader=None, shake=No
768768

769769
# Add constraints.
770770
isWater = [prmtop.getResidueLabel(i) in ('WAT', 'HOH', 'TP4', 'TP5', 'T4E') for i in range(prmtop.getNumAtoms())]
771+
isEP = [a.element is None for a in topology.atoms()]
771772
if shake in ('h-bonds', 'all-bonds', 'h-angles'):
772773
for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
773-
system.addConstraint(iAtom, jAtom, rMin)
774+
if not (isEP[iAtom] or isEP[jAtom]):
775+
system.addConstraint(iAtom, jAtom, rMin)
774776
if shake in ('all-bonds', 'h-angles'):
775777
for (iAtom, jAtom, k, rMin) in prmtop.getBondsNoH():
776-
system.addConstraint(iAtom, jAtom, rMin)
778+
if not (isEP[iAtom] or isEP[jAtom]):
779+
system.addConstraint(iAtom, jAtom, rMin)
777780
if rigidWater and shake is None:
778781
for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
779-
if isWater[iAtom] and isWater[jAtom]:
782+
if isWater[iAtom] and isWater[jAtom] and not (isEP[iAtom] or isEP[jAtom]):
780783
system.addConstraint(iAtom, jAtom, rMin)
781784

782785
# Add harmonic bonds.

‎wrappers/python/tests/TestAmberPrmtopFile.py‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -448,5 +448,14 @@ def testAmberCMAP(self):
448448
OpenMM_CMAP_E = simulation.context.getState(getEnergy=True, groups=1<<i).getPotentialEnergy().value_in_unit(kilojoules_per_mole)/conversion
449449
self.assertAlmostEqual(OpenMM_CMAP_E, sander_CMAP_E, places=4)
450450

451+
def testEPConstraints(self):
452+
"""Test different types of constraints when using extra particles"""
453+
prmtop = AmberPrmtopFile('systems/peptide_with_tip4p.prmtop')
454+
for constraints in (HBonds, AllBonds):
455+
system = prmtop.createSystem(constraints=constraints)
456+
integrator = VerletIntegrator(0.001*picoseconds)
457+
# If a constraint was added to a massless particle, this will throw an exception.
458+
context = Context(system, integrator, Platform.getPlatformByName('Reference'))
459+
451460
if __name__ == '__main__':
452461
unittest.main()

‎wrappers/python/tests/systems/peptide_with_tip4p.inpcrd‎

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

0 commit comments

Comments
 (0)