Skip to content

GitLab

Commit 52042c41 authored by peastman's avatar peastman
Browse files

Merge pull request #922 from swails/gmxdefines 

Fix Gromacs #define replacement and add support for Gromacs 5
parents 78e8a786 3e84be61
Expand all Show whitespace changes
Inline Side-by-side
Showing with 1638 additions and 5 deletions
wrappers/python/simtk/openmm/app/gromacstopfile.py
View file @ 52042c41
...@@ -8,7 +8,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org. ...@@ -8,7 +8,7 @@ Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2015 Stanford University and the Authors. Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
Contributors: Contributors: Jason Swails
Permission is hereby granted, free of charge, to any person obtaining a Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"), copy of this software and associated documentation files (the "Software"),
...@@ -40,7 +40,9 @@ import simtk.unit as unit ...@@ -40,7 +40,9 @@ import simtk.unit as unit
import simtk.openmm as mm import simtk.openmm as mm
import math import math
import os import os
import re
import distutils.spawn import distutils.spawn
from collections import OrderedDict
HBonds = ff.HBonds HBonds = ff.HBonds
AllBonds = ff.AllBonds AllBonds = ff.AllBonds
...@@ -48,6 +50,57 @@ HAngles = ff.HAngles ...@@ -48,6 +50,57 @@ HAngles = ff.HAngles
OBC2 = prmtop.OBC2 OBC2 = prmtop.OBC2
novarcharre = re.compile(r'\W')
def _find_all_instances_in_string(string, substr):
""" Find indices of all instances of substr in string """
indices = []
idx = string.find(substr, 0)
while idx > -1:
indices.append(idx)
idx = string.find(substr, idx+1)
return indices
def _replace_defines(line, defines):
""" Replaces defined tokens in a given line """
if not defines: return line
for define in reversed(defines):
value = defines[define]
indices = _find_all_instances_in_string(line, define)
if not indices: continue
# Check to see if it's inside of quotes
inside = ''
idx = 0
n_to_skip = 0
new_line = []
for i, char in enumerate(line):
if n_to_skip:
n_to_skip -= 1
continue
if char in ('\'"'):
if not inside:
inside = char
else:
if inside == char:
inside = ''
if idx < len(indices) and i == indices[idx]:
if inside:
new_line.append(char)
idx += 1
continue
if i == 0 or novarcharre.match(line[i-1]):
endidx = indices[idx] + len(define)
if endidx >= len(line) or novarcharre.match(line[endidx]):
new_line.extend(list(value))
n_to_skip = len(define) - 1
idx += 1
continue
idx += 1
new_line.append(char)
line = ''.join(new_line)
return line
class GromacsTopFile(object): class GromacsTopFile(object):
"""GromacsTopFile parses a Gromacs top file and constructs a Topology and (optionally) an OpenMM System from it.""" """GromacsTopFile parses a Gromacs top file and constructs a Topology and (optionally) an OpenMM System from it."""
...@@ -113,6 +166,7 @@ class GromacsTopFile(object): ...@@ -113,6 +166,7 @@ class GromacsTopFile(object):
name = fields[1] name = fields[1]
valueStart = stripped.find(name, len(command))+len(name)+1 valueStart = stripped.find(name, len(command))+len(name)+1
value = line[valueStart:].strip() value = line[valueStart:].strip()
value = value or '1' # Default define is 1
self._defines[name] = value self._defines[name] = value
elif command == '#ifdef': elif command == '#ifdef':
# See whether this block should be ignored. # See whether this block should be ignored.
...@@ -121,6 +175,12 @@ class GromacsTopFile(object): ...@@ -121,6 +175,12 @@ class GromacsTopFile(object):
name = fields[1] name = fields[1]
self._ifStack.append(name in self._defines) self._ifStack.append(name in self._defines)
self._elseStack.append(False) self._elseStack.append(False)
elif command == '#undef':
# Un-define a variable
if len(fields) < 2:
raise ValueError('Illegal line in .top file: '+line)
if fields[1] in self._defines:
self._defines.pop(fields[1])
elif command == '#ifndef': elif command == '#ifndef':
# See whether this block should be ignored. # See whether this block should be ignored.
if len(fields) < 2: if len(fields) < 2:
...@@ -145,6 +205,11 @@ class GromacsTopFile(object): ...@@ -145,6 +205,11 @@ class GromacsTopFile(object):
self._elseStack[-1] = True self._elseStack[-1] = True
elif not ignore: elif not ignore:
# Gromacs occasionally uses #define's to introduce specific
# parameters for individual terms (for instance, this is how
# ff99SB-ILDN is implemented). So make sure we do the appropriate
# pre-processor replacements necessary
line = _replace_defines(line, self._defines)
# A line of data for the current category # A line of data for the current category
if self._currentCategory is None: if self._currentCategory is None:
raise ValueError('Unexpected line in .top file: '+line) raise ValueError('Unexpected line in .top file: '+line)
...@@ -379,9 +444,11 @@ class GromacsTopFile(object): ...@@ -379,9 +444,11 @@ class GromacsTopFile(object):
# unless the preprocessor #define FLEXIBLE is given, don't define # unless the preprocessor #define FLEXIBLE is given, don't define
# bonds between the water hydrogen and oxygens, but only give the # bonds between the water hydrogen and oxygens, but only give the
# constraint distances and exclusions. # constraint distances and exclusions.
self._defines = {'FLEXIBLE': True} self._defines = OrderedDict()
self._defines['FLEXIBLE'] = True
if defines is not None: if defines is not None:
self._defines.update(defines) for define, value in defines.iteritems():
self._defines[define] = value
# Parse the file. # Parse the file.
...@@ -796,8 +863,8 @@ class GromacsTopFile(object): ...@@ -796,8 +863,8 @@ class GromacsTopFile(object):
def _defaultGromacsIncludeDir(): def _defaultGromacsIncludeDir():
"""Find the location where gromacs #include files are referenced from, by """Find the location where gromacs #include files are referenced from, by
searching for (1) gromacs environment variables, (2) for the gromacs binary searching for (1) gromacs environment variables, (2) for the gromacs binary
'pdb2gmx' in the PATH, or (3) just using the default gromacs install 'pdb2gmx' or 'gmx' in the PATH, or (3) just using the default gromacs
location, /usr/local/gromacs/share/gromacs/top """ install location, /usr/local/gromacs/share/gromacs/top """
if 'GMXDATA' in os.environ: if 'GMXDATA' in os.environ:
return os.path.join(os.environ['GMXDATA'], 'top') return os.path.join(os.environ['GMXDATA'], 'top')
if 'GMXBIN' in os.environ: if 'GMXBIN' in os.environ:
...@@ -806,5 +873,9 @@ def _defaultGromacsIncludeDir(): ...@@ -806,5 +873,9 @@ def _defaultGromacsIncludeDir():
pdb2gmx_path = distutils.spawn.find_executable('pdb2gmx') pdb2gmx_path = distutils.spawn.find_executable('pdb2gmx')
if pdb2gmx_path is not None: if pdb2gmx_path is not None:
return os.path.abspath(os.path.join(os.path.dirname(pdb2gmx_path), '..', 'share', 'gromacs', 'top')) return os.path.abspath(os.path.join(os.path.dirname(pdb2gmx_path), '..', 'share', 'gromacs', 'top'))
else:
gmx_path = distutils.spawn.find_executable('gmx')
if gmx_path is not None:
return os.path.abspath(os.path.join(os.path.dirname(gmx_path), '..', 'share', 'gromacs', 'top'))
return '/usr/local/gromacs/share/gromacs/top' return '/usr/local/gromacs/share/gromacs/top'
wrappers/python/tests/TestGromacsTopFile.py
View file @ 52042c41
...@@ -31,6 +31,20 @@ class TestGromacsTopFile(unittest.TestCase): ...@@ -31,6 +31,20 @@ class TestGromacsTopFile(unittest.TestCase):
f.getNonbondedMethod()==methodMap[method] f.getNonbondedMethod()==methodMap[method]
for f in forces)) for f in forces))
def test_ff99SBILDN(self):
""" Test Gromacs topology #define replacement as used in ff99SB-ILDN """
top = GromacsTopFile('systems/aidilnaaaaa.top')
gro = GromacsGroFile('systems/aidilnaaaaa.gro')
system = top.createSystem()
for force in system.getForces():
if isinstance(force, PeriodicTorsionForce):
force.setForceGroup(1)
context = Context(system, VerletIntegrator(1*femtosecond),
Platform.getPlatformByName('Reference'))
context.setPositions(gro.positions)
ene = context.getState(getEnergy=True, groups=2**1).getPotentialEnergy()
self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), 341.6905133582857)
def test_Cutoff(self): def test_Cutoff(self):
"""Test to make sure the nonbondedCutoff parameter is passed correctly.""" """Test to make sure the nonbondedCutoff parameter is passed correctly."""
......
wrappers/python/tests/systems/aidilnaaaaa.gro 0 → 100644
View file @ 52042c41
S C A M O R G
146
1ALA N 1 0.333 0.155 -0.000
1ALA H1 2 0.360 0.104 -0.082
1ALA H2 3 0.360 0.104 0.082
1ALA H3 4 0.233 0.167 0.000
1ALA CA 5 0.397 0.285 -0.000
1ALA HA 6 0.364 0.328 -0.084
1ALA CB 7 0.358 0.365 0.123
1ALA HB1 8 0.403 0.454 0.121
1ALA HB2 9 0.259 0.379 0.124
1ALA HB3 10 0.385 0.316 0.206
1ALA C 11 0.549 0.271 -0.000
1ALA O 12 0.601 0.159 -0.000
2ILE N 13 0.619 0.384 -0.000
2ILE H 14 0.571 0.471 0.000
2ILE CA 15 0.764 0.384 -0.000
2ILE HA 16 0.791 0.337 0.084
2ILE CB 17 0.815 0.312 -0.124
2ILE HB 18 0.781 0.218 -0.124
2ILE CG2 19 0.764 0.384 -0.249
2ILE HG21 20 0.797 0.337 -0.331
2ILE HG22 21 0.664 0.384 -0.249
2ILE HG23 22 0.797 0.478 -0.249
2ILE CG1 23 0.967 0.312 -0.124
2ILE HG11 24 1.001 0.265 -0.043
2ILE HG12 25 1.001 0.406 -0.125
2ILE CD 26 1.018 0.240 -0.249
2ILE HD1 27 1.118 0.240 -0.249
2ILE HD2 28 0.985 0.146 -0.249
2ILE HD3 29 0.985 0.287 -0.331
2ILE C 30 0.819 0.526 -0.000
2ILE O 31 0.743 0.622 0.000
3ASP N 32 0.952 0.539 -0.000
3ASP H 33 1.008 0.456 0.000
3ASP CA 34 1.016 0.668 0.000
3ASP HA 35 0.983 0.711 -0.084
3ASP CB 36 0.977 0.749 0.123
3ASP HB1 37 1.004 0.700 0.206
3ASP HB2 38 0.878 0.763 0.124
3ASP CG 39 1.047 0.885 0.119
3ASP OD1 40 0.986 0.979 0.062
3ASP OD2 41 1.160 0.893 0.174
3ASP C 42 1.168 0.654 -0.000
3ASP O 43 1.220 0.543 -0.000
4ILE N 44 1.238 0.768 -0.000
4ILE H 45 1.190 0.855 0.000
4ILE CA 46 1.383 0.768 -0.000
4ILE HA 47 1.410 0.720 0.084
4ILE CB 48 1.434 0.696 -0.124
4ILE HB 49 1.401 0.602 -0.124
4ILE CG2 50 1.383 0.768 -0.249
4ILE HG21 51 1.416 0.721 -0.331
4ILE HG22 52 1.283 0.768 -0.249
4ILE HG23 53 1.416 0.862 -0.249
4ILE CG1 54 1.586 0.696 -0.124
4ILE HG11 55 1.620 0.649 -0.043
4ILE HG12 56 1.620 0.790 -0.125
4ILE CD 57 1.637 0.624 -0.249
4ILE HD1 58 1.737 0.624 -0.249
4ILE HD2 59 1.604 0.530 -0.249
4ILE HD3 60 1.604 0.671 -0.331
4ILE C 61 1.438 0.910 -0.000
4ILE O 62 1.362 1.006 0.000
5LEU N 63 1.571 0.923 -0.000
5LEU H 64 1.628 0.840 0.000
5LEU CA 65 1.635 1.052 0.000
5LEU HA 66 1.602 1.095 -0.084
5LEU CB 67 1.596 1.133 0.123
5LEU HB1 68 1.623 1.084 0.206
5LEU HB2 69 1.497 1.147 0.124
5LEU CG 70 1.666 1.269 0.119
5LEU HG 71 1.765 1.255 0.119
5LEU CD1 72 1.623 1.344 -0.006
5LEU HD11 73 1.669 1.433 -0.009
5LEU HD12 74 1.649 1.291 -0.087
5LEU HD13 75 1.524 1.357 -0.006
5LEU CD2 76 1.626 1.349 0.243
5LEU HD21 77 1.672 1.438 0.240
5LEU HD22 78 1.527 1.363 0.244
5LEU HD23 79 1.654 1.300 0.325
5LEU C 80 1.787 1.038 -0.000
5LEU O 81 1.839 0.927 -0.000
6ASN N 82 1.857 1.152 -0.000
6ASN H 83 1.809 1.239 0.000
6ASN CA 84 2.002 1.152 -0.000
6ASN HA 85 2.026 1.103 0.084
6ASN CB 86 2.057 1.080 -0.123
6ASN HB1 87 2.025 1.127 -0.205
6ASN HB2 88 2.025 0.986 -0.124
6ASN CG 89 2.209 1.080 -0.123
6ASN OD1 90 2.272 1.133 -0.032
6ASN ND2 91 2.269 1.021 -0.227
6ASN HD21 92 2.214 0.979 -0.299
6ASN HD22 93 2.369 1.018 -0.232
6ASN C 94 2.057 1.294 -0.000
6ASN O 95 1.981 1.390 0.000
7ALA N 96 2.190 1.306 -0.000
7ALA H 97 2.247 1.224 0.000
7ALA CA 98 2.254 1.436 0.000
7ALA HA 99 2.222 1.479 -0.084
7ALA CB 100 2.215 1.517 0.123
7ALA HB1 101 2.261 1.606 0.121
7ALA HB2 102 2.116 1.530 0.124
7ALA HB3 103 2.243 1.468 0.206
7ALA C 104 2.406 1.422 -0.000
7ALA O 105 2.458 1.311 -0.000
8ALA N 106 2.476 1.535 -0.000
8ALA H 107 2.428 1.623 0.000
8ALA CA 108 2.621 1.535 -0.000
8ALA HA 109 2.645 1.487 0.084
8ALA CB 110 2.676 1.464 -0.123
8ALA HB1 111 2.776 1.466 -0.121
8ALA HB2 112 2.644 1.370 -0.124
8ALA HB3 113 2.644 1.511 -0.206
8ALA C 114 2.676 1.677 0.000
8ALA O 115 2.600 1.774 0.000
9ALA N 116 2.809 1.690 -0.000
9ALA H 117 2.866 1.608 -0.000
9ALA CA 118 2.873 1.820 0.000
9ALA HA 119 2.841 1.863 -0.084
9ALA CB 120 2.834 1.901 0.123
9ALA HB1 121 2.880 1.990 0.121
9ALA HB2 122 2.735 1.914 0.124
9ALA HB3 123 2.862 1.851 0.206
9ALA C 124 3.025 1.806 -0.000
9ALA O 125 3.077 1.695 -0.000
10ALA N 126 3.095 1.919 -0.000
10ALA H 127 3.047 2.007 0.000
10ALA CA 128 3.240 1.919 -0.000
10ALA HA 129 3.264 1.871 0.084
10ALA CB 130 3.295 1.848 -0.123
10ALA HB1 131 3.395 1.850 -0.121
10ALA HB2 132 3.263 1.753 -0.124
10ALA HB3 133 3.263 1.895 -0.206
10ALA C 134 3.295 2.061 0.000
10ALA O 135 3.219 2.158 0.000
11ALA N 136 3.428 2.074 -0.000
11ALA H 137 3.485 1.992 -0.000
11ALA CA 138 3.493 2.204 0.000
11ALA HA 139 3.460 2.247 -0.084
11ALA CB 140 3.453 2.285 0.123
11ALA HB1 141 3.499 2.374 0.121
11ALA HB2 142 3.354 2.298 0.124
11ALA HB3 143 3.481 2.235 0.206
11ALA C 144 3.644 2.190 -0.000
11ALA OC1 145 3.716 2.289 -0.000
11ALA OC2 146 3.696 2.079 -0.000
3.48257 2.26962 0.65582
wrappers/python/tests/systems/aidilnaaaaa.top 0 → 100644
View file @ 52042c41
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment