Add new 'smirnoff' mode to ForceField improper torsion handling behaviors to...

Add new 'smirnoff' mode to ForceField improper torsion handling behaviors to allow SMIRNOFF trefoil torsions to be added correctly

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

@@ -632,29 +632,6 @@ class ForceField(object):
             atom = self.getAtomIndexByName(atom_name)
             self.addExternalBond(atom)
 
-        def areParametersIdentical(self, template2, matchingAtoms, matchingAtoms2):
-            """Get whether this template and another one both assign identical atom types and parameters to all atoms.
-            
-            Parameters
-            ----------
-            template2: _TemplateData
-                the template to compare this one to
-            matchingAtoms: list
-                the indices of atoms in this template that atoms of the residue are matched to
-            matchingAtoms2: list
-                the indices of atoms in template2 that atoms of the residue are matched to
-            """
-            atoms1 = [self.atoms[m] for m in matchingAtoms]
-            atoms2 = [template2.atoms[m] for m in matchingAtoms2]
-            if any(a1.type != a2.type or a1.parameters != a2.parameters for a1,a2 in zip(atoms1, atoms2)):
-                return False
-            # Properly comparing virtual sites really needs a much more complicated analysis.  This simple check
-            # could easily fail for templates containing vsites, even if they're actually identical.  Since we
-            # currently have no force fields that include both patches and vsites, I'm not going to worry about it now.
-            if self.virtualSites != template2.virtualSites:
-                return False
-            return True
-
     class _TemplateAtomData(object):
         """Inner class used to encapsulate data about an atom in a residue template definition."""
         def __init__(self, name, type, element, parameters={}):
@@ -889,7 +866,7 @@ class ForceField(object):
 
         Returns
         -------
-        template : _TemplateData
+        template : _ForceFieldTemplate
             The matching forcefield residue template, or None if no matches are found.
         matches : list
             a list specifying which atom of the template each atom of the residue
@@ -911,13 +888,7 @@ class ForceField(object):
                 template = allMatches[0][0]
                 matches = allMatches[0][1]
             elif len(allMatches) > 1:
-                # We found multiple matches.  This is OK if and only if they assign identical types and parameters to all atoms.
-                t1, m1 = allMatches[0]
-                for t2, m2 in allMatches[1:]:
-                    if not t1.areParametersIdentical(t2, m1, m2):
-                        raise Exception('Multiple non-identical matching templates found for residue %d (%s): %s.' % (res.index+1, res.name, ', '.join(match[0].name for match in allMatches)))
-                template = allMatches[0][0]
-                matches = allMatches[0][1]
+                raise Exception('Multiple matching templates found for residue %d (%s): %s.' % (res.index+1, res.name, ', '.join(match[0].name for match in allMatches)))
         return [template, matches]
 
     def _buildBondedToAtomList(self, topology):
@@ -1841,6 +1812,16 @@ def _matchImproper(data, torsion, generator):
                                 (a2, a3) = (a3, a2)
                         match = (a2, a3, torsion[0], a4, tordef)
                         break
+                    elif tordef.ordering == 'smirnoff':
+                        # topology atom indexes
+                        a1 = torsion[0]
+                        a2 = torsion[t2[1]]
+                        a3 = torsion[t3[1]]
+                        a4 = torsion[t4[1]]
+                        # enforce exact match
+                        match = (a1, a2, a3, a4, tordef)
+                        break
+
     return match
 
 
@@ -2032,7 +2013,7 @@ class PeriodicTorsionGenerator(object):
     def registerImproperTorsion(self, parameters, ordering='default'):
         torsion = self.ff._parseTorsion(parameters)
         if torsion is not None:
-            if ordering in ['default', 'charmm', 'amber']:
+            if ordering in ['default', 'charmm', 'amber', 'smirnoff']:
                 torsion.ordering = ordering
             else:
                 raise ValueError('Illegal ordering type %s for improper torsion %s' % (ordering, torsion))
@@ -2116,7 +2097,13 @@ class PeriodicTorsionGenerator(object):
                 (a1, a2, a3, a4, tordef) = match
                 for i in range(len(tordef.phase)):
                     if tordef.k[i] != 0:
-                        force.addTorsion(a1, a2, a3, a4, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                        if tordef.ordering == 'smirnoff':
+                            # Add all torsions in trefoil
+                            force.addTorsion(a1, a2, a3, a4, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                            force.addTorsion(a1, a3, a4, a2, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                            force.addTorsion(a1, a4, a2, a3, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
+                        else:
+                            force.addTorsion(a1, a2, a3, a4, tordef.periodicity[i], tordef.phase[i], tordef.k[i])
 parsers["PeriodicTorsionForce"] = PeriodicTorsionGenerator.parseElement
 
 
@@ -4460,7 +4447,7 @@ class AmoebaVdwGenerator(object):
     def createForce(self, sys, data, nonbondedMethod, nonbondedCutoff, args):
 
         sigmaMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'CUBIC-MEAN':1}
-        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'W-H':1, 'HHG':1}
+        epsilonMap = {'ARITHMETIC':1, 'GEOMETRIC':1, 'HARMONIC':1, 'HHG':1}
 
         force = mm.AmoebaVdwForce()
         sys.addForce(force)