Skip to content

GitLab

Commit b0c7239b authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #1560 from rafwiewiora/master 

New code for impropers assignments
parents 915b84f6 3850a251
Show whitespace changes
Inline Side-by-side
Showing with 224 additions and 104 deletions
wrappers/python/simtk/openmm/app/forcefield.py
View file @ b0c7239b
......@@ -538,6 +538,7 @@ class ForceField(object):
def __init__(self):
self.atomType = {}
self.atomParameters = {}
self.atomTemplateIndexes = {}
self.atoms = []
self.excludeAtomWith = []
self.virtualSites = {}
......@@ -565,6 +566,7 @@ class ForceField(object):
for atom, match in zip(residue.atoms(), matches):
self.atomType[atom] = template.atoms[match].type
self.atomParameters[atom] = template.atoms[match].parameters
self.atomTemplateIndexes[atom] = match
for site in template.virtualSites:
if match == site.index:
self.virtualSites[atom] = (site, [matchAtoms[i].index for i in site.atoms], matchAtoms[site.excludeWith].index)
......@@ -1734,6 +1736,108 @@ def _createResidueTemplate(residue):
template.addExternalBondByName(atom2.name)
return template
def _matchImproper(data, torsion, generator):
type1 = data.atomType[data.atoms[torsion[0]]]
type2 = data.atomType[data.atoms[torsion[1]]]
type3 = data.atomType[data.atoms[torsion[2]]]
type4 = data.atomType[data.atoms[torsion[3]]]
wildcard = generator.ff._atomClasses['']
match = None
for tordef in generator.improper:
types1 = tordef.types1
types2 = tordef.types2
types3 = tordef.types3
types4 = tordef.types4
hasWildcard = (wildcard in (types1, types2, types3, types4))
if match is not None and hasWildcard:
# Prefer specific definitions over ones with wildcards
continue
if type1 in types1:
for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
if tordef.ordering == 'default':
# Workaround to be more consistent with AMBER. It uses wildcards to define most of its
# impropers, which leaves the ordering ambiguous. It then follows some bizarre rules
# to pick the order.
a1 = torsion[t2[1]]
a2 = torsion[t3[1]]
e1 = data.atoms[a1].element
e2 = data.atoms[a2].element
if e1 == e2 and a1 > a2:
(a1, a2) = (a2, a1)
elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
(a1, a2) = (a2, a1)
match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
break
elif tordef.ordering == 'charmm':
if hasWildcard:
# Workaround to be more consistent with AMBER. It uses wildcards to define most of its
# impropers, which leaves the ordering ambiguous. It then follows some bizarre rules
# to pick the order.
a1 = torsion[t2[1]]
a2 = torsion[t3[1]]
e1 = data.atoms[a1].element
e2 = data.atoms[a2].element
if e1 == e2 and a1 > a2:
(a1, a2) = (a2, a1)
elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
(a1, a2) = (a2, a1)
match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
else:
# There are no wildcards, so the order is unambiguous.
match = (torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef)
break
elif tordef.ordering == 'amber':
# topology atom indexes
a2 = torsion[t2[1]]
a3 = torsion[t3[1]]
a4 = torsion[t4[1]]
# residue indexes
r2 = data.atoms[a2].residue.index
r3 = data.atoms[a3].residue.index
r4 = data.atoms[a4].residue.index
# template atom indexes
ta2 = data.atomTemplateIndexes[data.atoms[a2]]
ta3 = data.atomTemplateIndexes[data.atoms[a3]]
ta4 = data.atomTemplateIndexes[data.atoms[a4]]
# elements
e2 = data.atoms[a2].element
e3 = data.atoms[a3].element
e4 = data.atoms[a4].element
if not hasWildcard:
if t2[0] == t4[0] and (r2 > r4 or (r2 == r4 and ta2 > ta4)):
(a2, a4) = (a4, a2)
r2 = data.atoms[a2].residue.index
r4 = data.atoms[a4].residue.index
ta2 = data.atomTemplateIndexes[data.atoms[a2]]
ta4 = data.atomTemplateIndexes[data.atoms[a4]]
if t3[0] == t4[0] and (r3 > r4 or (r3 == r4 and ta3 > ta4)):
(a3, a4) = (a4, a3)
r3 = data.atoms[a3].residue.index
r4 = data.atoms[a4].residue.index
ta3 = data.atomTemplateIndexes[data.atoms[a3]]
ta4 = data.atomTemplateIndexes[data.atoms[a4]]
if t2[0] == t3[0] and (r2 > r3 or (r2 == r3 and ta2 > ta3)):
(a2, a3) = (a3, a2)
else:
if e2 == e4 and (r2 > r4 or (r2 == r4 and ta2 > ta4)):
(a2, a4) = (a4, a2)
r2 = data.atoms[a2].residue.index
r4 = data.atoms[a4].residue.index
ta2 = data.atomTemplateIndexes[data.atoms[a2]]
ta4 = data.atomTemplateIndexes[data.atoms[a4]]
if e3 == e4 and (r3 > r4 or (r3 == r4 and ta3 > ta4)):
(a3, a4) = (a4, a3)
r3 = data.atoms[a3].residue.index
r4 = data.atoms[a4].residue.index
ta3 = data.atomTemplateIndexes[data.atoms[a3]]
ta4 = data.atomTemplateIndexes[data.atoms[a4]]
if r2 > r3 or (r2 == r3 and ta2 > ta3):
(a2, a3) = (a3, a2)
match = (a2, a3, torsion[0], a4, tordef)
break
return match
# The following classes are generators that know how to create Force subclasses and add them to a System that is being
# created. Each generator class must define two methods: 1) a static method that takes an etree Element and a ForceField,
# and returns the corresponding generator object; 2) a createForce() method that constructs the Force object and adds it
......@@ -1893,6 +1997,7 @@ class PeriodicTorsion(object):
self.periodicity = []
self.phase = []
self.k = []
self.ordering = 'default'
## @private
class PeriodicTorsionGenerator(object):
......@@ -1908,9 +2013,13 @@ class PeriodicTorsionGenerator(object):
if torsion is not None:
self.proper.append(torsion)
def registerImproperTorsion(self, parameters):
def registerImproperTorsion(self, parameters, ordering='default'):
torsion = self.ff._parseTorsion(parameters)
if torsion is not None:
if ordering in ['default', 'charmm', 'amber']:
torsion.ordering = ordering
else:
raise ValueError('Illegal ordering type %s for improper torsion %s' % (ordering, torsion))
self.improper.append(torsion)
@staticmethod
......@@ -1924,6 +2033,9 @@ class PeriodicTorsionGenerator(object):
for torsion in element.findall('Proper'):
generator.registerProperTorsion(torsion.attrib)
for torsion in element.findall('Improper'):
if 'ordering' in element.attrib:
generator.registerImproperTorsion(torsion.attrib, element.attrib['ordering'])
else:
generator.registerImproperTorsion(torsion.attrib)
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
......@@ -1957,36 +2069,7 @@ class PeriodicTorsionGenerator(object):
if match.k[i] != 0:
force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.periodicity[i], match.phase[i], match.k[i])
for torsion in data.impropers:
type1 = data.atomType[data.atoms[torsion[0]]]
type2 = data.atomType[data.atoms[torsion[1]]]
type3 = data.atomType[data.atoms[torsion[2]]]
type4 = data.atomType[data.atoms[torsion[3]]]
match = None
for tordef in self.improper:
types1 = tordef.types1
types2 = tordef.types2
types3 = tordef.types3
types4 = tordef.types4
hasWildcard = (wildcard in (types1, types2, types3, types4))
if match is not None and hasWildcard:
# Prefer specific definitions over ones with wildcards
continue
if type1 in types1:
for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
# Workaround to be more consistent with AMBER. It uses wildcards to define most of its
# impropers, which leaves the ordering ambiguous. It then follows some bizarre rules
# to pick the order.
a1 = torsion[t2[1]]
a2 = torsion[t3[1]]
e1 = data.atoms[a1].element
e2 = data.atoms[a2].element
if e1 == e2 and a1 > a2:
(a1, a2) = (a2, a1)
elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
(a1, a2) = (a2, a1)
match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
break
match = _matchImproper(data, torsion, self)
if match is not None:
(a1, a2, a3, a4, tordef) = match
for i in range(len(tordef.phase)):
......@@ -2000,12 +2083,16 @@ parsers["PeriodicTorsionForce"] = PeriodicTorsionGenerator.parseElement
class RBTorsion(object):
"""An RBTorsion records the information for a Ryckaert-Bellemans torsion definition."""
def __init__(self, types, c):
def __init__(self, types, c, ordering='charmm'):
self.types1 = types[0]
self.types2 = types[1]
self.types3 = types[2]
self.types4 = types[3]
self.c = c
if ordering in ['default', 'charmm', 'amber']:
self.ordering = ordering
else:
raise ValueError('Illegal ordering type %s for RBTorsion (%s,%s,%s,%s)' % (ordering, types[0], types[1], types[2], types[3]))
## @private
class RBTorsionGenerator(object):
......@@ -2031,6 +2118,9 @@ class RBTorsionGenerator(object):
for torsion in element.findall('Improper'):
types = ff._findAtomTypes(torsion.attrib, 4)
if None not in types:
if 'ordering' in element.attrib:
generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)], element.attrib['ordering']))
else:
generator.improper.append(RBTorsion(types, [float(torsion.attrib['c'+str(i)]) for i in range(6)]))
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
......@@ -2062,40 +2152,7 @@ class RBTorsionGenerator(object):
if match is not None:
force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.c[0], match.c[1], match.c[2], match.c[3], match.c[4], match.c[5])
for torsion in data.impropers:
type1 = data.atomType[data.atoms[torsion[0]]]
type2 = data.atomType[data.atoms[torsion[1]]]
type3 = data.atomType[data.atoms[torsion[2]]]
type4 = data.atomType[data.atoms[torsion[3]]]
match = None
for tordef in self.improper:
types1 = tordef.types1
types2 = tordef.types2
types3 = tordef.types3
types4 = tordef.types4
hasWildcard = (wildcard in (types1, types2, types3, types4))
if match is not None and hasWildcard:
# Prefer specific definitions over ones with wildcards
continue
if type1 in types1:
for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
if hasWildcard:
# Workaround to be more consistent with AMBER. It uses wildcards to define most of its
# impropers, which leaves the ordering ambiguous. It then follows some bizarre rules
# to pick the order.
a1 = torsion[t2[1]]
a2 = torsion[t3[1]]
e1 = data.atoms[a1].element
e2 = data.atoms[a2].element
if e1 == e2 and a1 > a2:
(a1, a2) = (a2, a1)
elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
(a1, a2) = (a2, a1)
match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
else:
# There are no wildcards, so the order is unambiguous.
match = (torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef)
break
match = _matchImproper(data, torsion, self)
if match is not None:
(a1, a2, a3, a4, tordef) = match
force.addTorsion(a1, a2, a3, a4, tordef.c[0], tordef.c[1], tordef.c[2], tordef.c[3], tordef.c[4], tordef.c[5])
......@@ -2556,12 +2613,16 @@ parsers["CustomAngleForce"] = CustomAngleGenerator.parseElement
class CustomTorsion(object):
"""A CustomTorsion records the information for a custom torsion definition."""
def __init__(self, types, paramValues):
def __init__(self, types, paramValues, ordering='charmm'):
self.types1 = types[0]
self.types2 = types[1]
self.types3 = types[2]
self.types4 = types[3]
self.paramValues = paramValues
if ordering in ['default', 'charmm', 'amber']:
self.ordering = ordering
else:
raise ValueError('Illegal ordering type %s for CustomTorsion (%s,%s,%s,%s)' % (ordering, types[0], types[1], types[2], types[3]))
## @private
class CustomTorsionGenerator(object):
......@@ -2590,6 +2651,9 @@ class CustomTorsionGenerator(object):
for torsion in element.findall('Improper'):
types = ff._findAtomTypes(torsion.attrib, 4)
if None not in types:
if 'ordering' in element.attrib:
generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams], element.attrib['ordering']))
else:
generator.improper.append(CustomTorsion(types, [float(torsion.attrib[param]) for param in generator.perTorsionParams]))
def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
......@@ -2620,40 +2684,7 @@ class CustomTorsionGenerator(object):
if match is not None:
force.addTorsion(torsion[0], torsion[1], torsion[2], torsion[3], match.paramValues)
for torsion in data.impropers:
type1 = data.atomType[data.atoms[torsion[0]]]
type2 = data.atomType[data.atoms[torsion[1]]]
type3 = data.atomType[data.atoms[torsion[2]]]
type4 = data.atomType[data.atoms[torsion[3]]]
match = None
for tordef in self.improper:
types1 = tordef.types1
types2 = tordef.types2
types3 = tordef.types3
types4 = tordef.types4
hasWildcard = (wildcard in (types1, types2, types3, types4))
if match is not None and hasWildcard:
# Prefer specific definitions over ones with wildcards
continue
if type1 in types1:
for (t2, t3, t4) in itertools.permutations(((type2, 1), (type3, 2), (type4, 3))):
if t2[0] in types2 and t3[0] in types3 and t4[0] in types4:
if hasWildcard:
# Workaround to be more consistent with AMBER. It uses wildcards to define most of its
# impropers, which leaves the ordering ambiguous. It then follows some bizarre rules
# to pick the order.
a1 = torsion[t2[1]]
a2 = torsion[t3[1]]
e1 = data.atoms[a1].element
e2 = data.atoms[a2].element
if e1 == e2 and a1 > a2:
(a1, a2) = (a2, a1)
elif e1 != elem.carbon and (e2 == elem.carbon or e1.mass < e2.mass):
(a1, a2) = (a2, a1)
match = (a1, a2, torsion[0], torsion[t4[1]], tordef)
else:
# There are no wildcards, so the order is unambiguous.
match = (torsion[0], torsion[t2[1]], torsion[t3[1]], torsion[t4[1]], tordef)
break
match = _matchImproper(data, torsion, self)
if match is not None:
(a1, a2, a3, a4, tordef) = match
force.addTorsion(a1, a2, a3, a4, tordef.paramValues)
......
wrappers/python/tests/TestForceField.py
View file @ b0c7239b
......@@ -769,6 +769,35 @@ class TestForceField(unittest.TestCase):
self.assertTrue('spce-O' in forcefield._atomTypes)
self.assertTrue('HOH' in forcefield._templates)
def test_ImpropersOrdering(self):
"""Test correctness of the ordering of atom indexes in improper torsions
and the torsion.ordering parameter.
"""
xml = """
<ForceField>
<PeriodicTorsionForce ordering="amber">
<Improper class1="C" class2="" class3="O2" class4="O2" periodicity1="2" phase1="3.14159265359" k1="43.932"/>
</PeriodicTorsionForce>
</ForceField>
"""
pdb = PDBFile('systems/impropers_ordering_tetrapeptide.pdb')
# ff1 uses default ordering of impropers, ff2 uses "amber" for the one
# problematic improper
ff1 = ForceField('amber99sbildn.xml')
ff2 = ForceField(StringIO(xml), 'amber99sbildn.xml')
system1 = ff1.createSystem(pdb.topology)
system2 = ff2.createSystem(pdb.topology)
imp1 = system1.getForce(2).getTorsionParameters(158)
imp2 = system2.getForce(0).getTorsionParameters(158)
system1_indexes = [imp1[0], imp1[1], imp1[2], imp1[3]]
system2_indexes = [imp2[0], imp2[1], imp2[2], imp2[3]]
self.assertEqual(system1_indexes, [51, 56, 54, 55])
self.assertEqual(system2_indexes, [51, 55, 54, 56])
class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......
wrappers/python/tests/systems/impropers_ordering_tetrapeptide.pdb 0 → 100644
View file @ b0c7239b
ATOM 1 N LEU A 1 25.160 14.160 19.440 1.00 0.00 N
ATOM 2 H1 LEU A 1 24.224 13.904 20.146 1.00 0.00 H
ATOM 3 H2 LEU A 1 25.993 13.474 19.964 1.00 0.00 H
ATOM 4 H3 LEU A 1 25.448 15.264 19.811 1.00 0.00 H
ATOM 5 CA LEU A 1 25.090 13.920 17.980 1.00 0.00 C
ATOM 6 HA LEU A 1 24.679 12.800 17.995 1.00 0.00 H
ATOM 7 CB LEU A 1 24.420 14.970 17.100 1.00 0.00 C
ATOM 8 HB3 LEU A 1 24.592 14.703 15.949 1.00 0.00 H
ATOM 9 HB2 LEU A 1 24.907 16.046 17.290 1.00 0.00 H
ATOM 10 CG LEU A 1 22.930 15.010 17.400 1.00 0.00 C
ATOM 11 HG LEU A 1 22.677 15.678 18.357 1.00 0.00 H
ATOM 12 CD1 LEU A 1 22.410 15.830 16.210 1.00 0.00 C
ATOM 13 HD11 LEU A 1 22.229 15.260 15.170 1.00 0.00 H
ATOM 14 HD12 LEU A 1 21.323 16.283 16.456 1.00 0.00 H
ATOM 15 HD13 LEU A 1 22.952 16.853 15.894 1.00 0.00 H
ATOM 16 CD2 LEU A 1 22.100 13.730 17.590 1.00 0.00 C
ATOM 17 HD21 LEU A 1 22.228 12.870 16.765 1.00 0.00 H
ATOM 18 HD22 LEU A 1 22.102 13.167 18.648 1.00 0.00 H
ATOM 19 HD23 LEU A 1 20.924 13.966 17.534 1.00 0.00 H
ATOM 20 C LEU A 1 26.520 13.970 17.430 1.00 0.00 C
ATOM 21 O LEU A 1 27.360 14.710 17.880 1.00 0.00 O
ATOM 22 N SER A 2 26.720 13.080 16.460 1.00 0.00 N
ATOM 23 H SER A 2 26.006 12.134 16.365 1.00 0.00 H
ATOM 24 CA SER A 2 27.950 12.790 15.610 1.00 0.00 C
ATOM 25 HA SER A 2 28.762 12.484 16.429 1.00 0.00 H
ATOM 26 CB SER A 2 27.740 11.640 14.610 1.00 0.00 C
ATOM 27 HB3 SER A 2 28.674 11.459 13.887 1.00 0.00 H
ATOM 28 HB2 SER A 2 26.789 11.677 13.883 1.00 0.00 H
ATOM 29 OG SER A 2 27.520 10.520 15.410 1.00 0.00 O
ATOM 30 HG SER A 2 27.939 9.540 14.882 1.00 0.00 H
ATOM 31 C SER A 2 28.360 14.010 14.820 1.00 0.00 C
ATOM 32 O SER A 2 27.440 14.680 14.350 1.00 0.00 O
ATOM 33 N ASP A 3 29.650 14.360 14.620 1.00 0.00 N
ATOM 34 H ASP A 3 30.534 13.595 14.830 1.00 0.00 H
ATOM 35 CA ASP A 3 29.940 15.810 14.230 1.00 0.00 C
ATOM 36 HA ASP A 3 29.441 16.464 15.093 1.00 0.00 H
ATOM 37 CB ASP A 3 31.420 16.030 14.240 1.00 0.00 C
ATOM 38 HB3 ASP A 3 32.161 15.358 13.586 1.00 0.00 H
ATOM 39 HB2 ASP A 3 31.910 15.981 15.330 1.00 0.00 H
ATOM 40 CG ASP A 3 31.690 17.460 13.840 1.00 0.00 C
ATOM 41 OD1 ASP A 3 31.390 18.400 14.660 1.00 0.00 O
ATOM 42 OD2 ASP A 3 32.230 17.650 12.700 1.00 0.00 O
ATOM 43 C ASP A 3 29.350 16.300 12.880 1.00 0.00 C
ATOM 44 O ASP A 3 28.860 17.390 12.790 1.00 0.00 O
ATOM 45 N GLU A 4 29.370 15.470 11.800 1.00 0.00 N
ATOM 46 H GLU A 4 29.943 14.435 11.896 1.00 0.00 H
ATOM 47 CA GLU A 4 28.630 15.770 10.590 1.00 0.00 C
ATOM 48 HA GLU A 4 28.813 16.889 10.221 1.00 0.00 H
ATOM 49 CB GLU A 4 29.000 14.740 9.500 1.00 0.00 C
ATOM 50 HB3 GLU A 4 28.209 14.877 8.613 1.00 0.00 H
ATOM 51 HB2 GLU A 4 28.927 13.570 9.734 1.00 0.00 H
ATOM 52 CG GLU A 4 30.400 15.140 9.010 1.00 0.00 C
ATOM 53 HG3 GLU A 4 31.338 14.917 9.713 1.00 0.00 H
ATOM 54 HG2 GLU A 4 30.559 16.245 8.583 1.00 0.00 H
ATOM 55 CD GLU A 4 30.820 14.370 7.750 1.00 0.00 C
ATOM 56 OE1 GLU A 4 31.770 13.490 7.830 1.00 0.00 O
ATOM 57 OE2 GLU A 4 30.220 14.580 6.660 1.00 0.00 O
ATOM 58 C GLU A 4 27.080 15.870 10.880 1.00 0.00 C
ATOM 59 O GLU A 4 26.440 16.810 10.390 1.00 0.00 O
ATOM 60 OXT GLU A 4 26.692 14.850 11.569 1.00 0.00 O
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