Merge pull request #239 from open-forcefield-group/constraints

Initial support for bond constraints

The-SMIRFF-force-field-format.md

@@ -172,6 +172,38 @@ Here is an example not intended for actual use:
 ```
 The charge model specified must be a method understood by the OpenEye toolkits, and the charge correction `increment` (in units of proton charge) will be applied on top of this by subtracting `increment` from the atom tagged as 1 and adding it to the atom tagged as 2.
 
+### CONSTRAINTS
+
+Bond length constraints can be specified through a `<Constraints/>` block, which can constrain bonds to their equilibrium lengths or specify an interatomic constraint distance.
+Two atoms must be tagged in the `smirks` attribute of each `<Constraint/>` record.
+
+To constrain two atoms to their equilibrium bond length, it is critical that a `<HarmonicBondForce/>` record be specified for those atoms:
+```XML
+<Constraints>
+  <!-- constrain all bonds to hydrogen to their equilibrium bond length -->
+  <Constraint smirks="[#1:1]-[*:2]" />
+</Constraints>
+```
+However, this constraint distance can be overridden, or two atoms that are not directly bonded constrained, by specifying the `distance` attribute (and optional `distance_unit` attribute for the `<Constraints/>` tag):
+```XML
+<Constraints distance_unit="angstroms">
+  <!-- constrain water O-H bond to equilibrium bond length (overrides earlier constraint) -->
+  <Constraint smirks="[#1:1]-[#8X2H2:2]-[#1]" distance="0.9572"/>
+  <!-- constrain water H...H, calculating equilibrium length from H-O-H equilibrium angle and H-O equilibrium bond lengths -->
+  <Constraint smirks="[#1:1]-[#8X2H2]-[#1:2]" distance="1.8532"/>
+</Constraints>
+```
+Typical molecular simulation practice is to constrain all bonds to hydrogen to their equilibrium bond lengths and enforce rigid TIP3P geometry on water molecules:
+```XML
+<Constraints distance_unit="angstroms">
+  <!-- constrain all bonds to hydrogen to their equilibrium bond length -->
+  <Constraint smirks="[#1:1]-[*:2]" />
+  <!-- TIP3P rigid water -->
+  <Constraint smirks="[#1:1]-[#8X2H2:2]-[#1]" distance="0.9572"/>
+  <Constraint smirks="[#1:1]-[#8X2H2]-[#1:2]" distance="1.8532"/>
+</Constraints>
+```
+
 ## Advanced features
 
 Standard usage is expected to rely primarily on the features documented above and potentially new features. However, some advanced features are also currently supported.

smarty/data/systems/README.md

@@ -0,0 +1,3 @@
+Solvated systems here were built using SolvationToolkit (github.com/mobleylab/solvationtoolkit) version 0.41 and 0.42 (equivalent except for DOI assignment for the latter), except as noted.
+
+Because of a bug in handling of insertion of CONECT entries in SolvationToolkit/openmoltools for systems containing water, water boxes were generated by taking GROMACS coordinate/topology files for solvated cyclohexane ('mobley_2689721') and ethanol ('mobley_2310185') from FreeSolv 0.5 and converting to PDB format via OpenMM 7.1 (reading via GromacsGroFile and GromacsTopFile and writing via PDBFile). 

smarty/data/systems/monomers/water.mol2

@@ -0,0 +1,14 @@
+@<TRIPOS>MOLECULE
+ZZU
+3 2 1 0 1
+SMALL
+USER_CHARGES
+@<TRIPOS>ATOM
+       1 O1          -0.2950    -0.2180     0.1540 oh            1 ZZU       -0.785000
+       2 H1          -0.0170     0.6750     0.4080 ho            1 ZZU        0.392500
+       3 H2           0.3120    -0.4570    -0.5630 ho            1 ZZU        0.392500
+@<TRIPOS>BOND
+       1        1        2 1
+       2        1        3 1
+@<TRIPOS>SUBSTRUCTURE
+       1 ZZU             1 RESIDUE    0 **** ROOT      0

smarty/data/systems/monomers/water.sdf

@@ -0,0 +1,28 @@
+water
+  -OEChem-03221709533D
+
+  3  2  0     0  0  0  0  0  0999 V2000
+   -0.2950   -0.2180    0.1540 O   0  0  0  0  0  0  0  0  0  0  0  0
+   -0.0170    0.6750    0.4080 H   0  0  0  0  0  0  0  0  0  0  0  0
+    0.3120   -0.4570   -0.5630 H   0  0  0  0  0  0  0  0  0  0  0  0
+  1  2  1  0  0  0  0
+  1  3  1  0  0  0  0
+M  END
+> <partial_charges>
+-0.7850000262260437
+0.39250001311302185
+0.39250001311302185
+
+
+> <partial_bond_orders>
+1.0
+1.0
+
+
+> <atom_types>
+oh
+ho
+ho
+
+
+$$$$

smarty/forcefield.py

@@ -232,6 +232,53 @@ def __init__(self, topology, reference_molecules):
         # Get/initialize bond orders
         self._updateBondOrders()
 
+        # Track constraints
+        self._constrainedAtomPairs = dict()
+
+    def constrainAtomPair(self, iatom, jatom, distance=True):
+        """
+        Mark a pair of atoms as constrained.
+
+        Parameters
+        ----------
+        iatom, jatom : int
+            Indices of atoms to mark as constrained.
+        distance : simtk.unit.Quantity, optional, default=True
+            Constraint distance if constraint has been applied,
+            or True if no constraint has yet been applied
+        """
+        # Check that constraint hasn't already been specified.
+        if (iatom,jatom) in self._constrainedAtomPairs:
+            existing_distance = self._constrainedAtomPairs[(iatom,jatom)]
+            if unit.is_quantity(existing_distance) and (distance is True):
+                raise Exception('Atoms (%d,%d) already constrained with distance %s but attempting to override with unspecified distance' % (iatom, jatom, existing_distance))
+            if (existing_distance is True) and (distance is True):
+                raise Exception('Atoms (%d,%d) already constrained with unspecified distance but attempting to override with unspecified distance' % (iatom, jatom))
+
+        self._constrainedAtomPairs[(iatom,jatom)] = distance
+        self._constrainedAtomPairs[(jatom,iatom)] = distance
+
+    def atomPairIsConstrained(self, iatom, jatom):
+        """
+        Check if a pair of atoms are marked as constrained.
+
+        Parameters
+        ----------
+        iatom, jatom : int
+            Indices of atoms to mark as constrained.
+
+        Returns
+        -------
+        distance : simtk.unit.Quantity or bool
+            True if constrained but constraints have not yet been applied
+            Distance if constraint has already been added to System
+
+        """
+        if (iatom,jatom) in self._constrainedAtomPairs:
+            return self._constrainedAtomPairs[(iatom,jatom)]
+        else:
+            return False
+
     def angles(self):
         """
         Get an iterator over all i-j-k angles.
@@ -618,7 +665,12 @@ def getGenerators(self):
 
     def registerGenerator(self, generator):
         """Register a new generator."""
-        self._forces.append(generator)
+        # Special case: ConstraintGenerator has to come before HarmonicBondGenerator and HarmonicAngleGenerator.
+        # TODO: Figure out a more general way to allow generators to specify enforced orderings.
+        if isinstance(generator, ConstraintGenerator):
+            self._forces.insert(0, generator)
+        else:
+            self._forces.append(generator)
 
     def getParameter(self, smirks = None, paramID=None, force_type='Implied'):
         """Get info associated with a particular parameter as specified by SMIRKS or parameter ID, and optionally force term.
@@ -1243,6 +1295,117 @@ def __keytransform__(self,key):
 # Force generators
 #=============================================================================================
 
+#=============================================================================================
+# Force generators
+#=============================================================================================
+
+## @private
+class ConstraintGenerator(object):
+    """A ConstraintGenerator adds constraints.
+
+    The ConstraintGenerator must be applied before HarmonicBondGenerator and HarmonicAngleGenerator if constrained bonds are to not also have harmonic bond terms added.
+
+    ConstraintGenerator will mark bonds as being constrained for HarmonicBondGenerator to assign constraints to equilibrium bond lengths,
+    while constraints with distances specified will be assigned by ConstraintGenerator.
+
+    """
+
+    class ConstraintType(object):
+        """A SMIRFF constraint type."""
+        def __init__(self, node, parent):
+            self.smirks = _validateSMIRKS(node.attrib['smirks'], node=node)
+            self.pid = _extractQuantity(node, parent, 'id')
+            if 'distance' in node.attrib:
+                # Constraint distance is specified, will be handled by ConstraintGenerator
+                self.distance = _extractQuantity(node, parent, 'distance')
+            else:
+                # Constraint to equilibrium bond length, handled by HarmonicBondGenerator
+                self.distance = True
+
+    def __init__(self, forcefield):
+        self.ff = forcefield
+        self._constraint_types = list()
+
+    def registerConstraint(self, node, parent):
+        """Register a SMIRFF constraint type definition."""
+        constraint = ConstraintGenerator.ConstraintType(node, parent)
+        self._constraint_types.append(constraint)
+
+    @staticmethod
+    def parseElement(element, ff):
+        # Find existing force generator or create new one.
+        existing = [f for f in ff._forces if isinstance(f, ConstraintGenerator)]
+        if len(existing) == 0:
+            generator = ConstraintGenerator(ff)
+            ff.registerGenerator(generator)
+        else:
+            generator = existing[0]
+
+        # Register all SMIRFF constraint definitions.
+        for constraint in element.findall('Constraint'):
+            generator.registerConstraint(constraint, element)
+
+    def createForce(self, system, topology, verbose=False, **kwargs):
+        # Iterate over all defined constraint types, allowing later matches to override earlier ones.
+        constraints = ValenceDict()
+        for constraint in self._constraint_types:
+            for atom_indices in topology.unrollSMIRKSMatches(constraint.smirks, aromaticity_model = self.ff._aromaticity_model):
+                constraints[atom_indices] = constraint
+
+        if verbose:
+            print('')
+            print('ConstraintGenerator:')
+            print('')
+            for constraint in self._constraint_types:
+                print('%64s : %8d matches' % (constraint.smirks, len(topology.unrollSMIRKSMatches(constraint.smirks, aromaticity_model = self.ff._aromaticity_model))))
+            print('')
+
+        for (atom_indices, constraint) in constraints.items():
+            # Update constrained atom pairs in topology
+            topology.constrainAtomPair(atom_indices[0], atom_indices[1], constraint.distance)
+            # If a distance is specified, add the constraint here.
+            # Otherwise, the equilibrium bond length will be used to constrain the atoms in HarmonicBondGenerator
+            if constraint.distance is not True:
+                system.addConstraint(atom_indices[0], atom_indices[1], constraint.distance)
+
+        if verbose: print('%d constraints added' % (len(constraints)))
+
+    def labelForce(self, oemol, verbose=False, **kwargs):
+        """Take a provided OEMol and parse Constraint terms for this molecule.
+
+        Parameters
+        ----------
+            oemol : OEChem OEMol object for molecule to be examined
+
+        Returns
+        ---------
+            force_terms: list
+                Returns a list of tuples, [ ([atom id 1, ... atom id N], parameter id, smirks) , (....), ... ] for all forces of this type which would be applied.
+        """
+
+        # Iterate over all defined constraint SMIRKS, allowing later matches to override earlier ones.
+        constraints = ValenceDict()
+        for constraint in self._constraint_types:
+            for atom_indices in topology.unrollSMIRKSMatches(constraint.smirks, aromaticity_model = self.ff._aromaticity_model):
+                constraints[atom_indices] = constraint
+
+        if verbose:
+            print('')
+            print('ConstraintGenerator:')
+            print('')
+            for constraint in self._constraint_types:
+                print('%64s : %8d matches' % (constraint.smirks, len(topology.unrollSMIRKSMatches(constraint.smirks, aromaticity_model = self.ff._aromaticity_model))))
+            print('')
+
+        # Add all bonds to the output list
+        force_terms = []
+        for (atom_indices, constraint) in constraints.items():
+            force_terms.append( ([atom_indices[0], atom_indices[1]], constraint.pid, constraint.smirks) )
+
+        return force_terms
+
+parsers["Constraints"] = ConstraintGenerator.parseElement
+
 ## @private
 class HarmonicBondGenerator(object):
     """A HarmonicBondGenerator constructs a HarmonicBondForce."""
@@ -1327,14 +1490,16 @@ def createForce(self, system, topology, verbose=False, **kwargs):
             print('')
 
         # Add all bonds to the system.
+        skipped_constrained_bonds = 0 # keep track of how many bonds were constrained (and hence skipped)
         for (atom_indices, bond) in bonds.items():
             # Ensure atoms are actually bonded correct pattern in Topology
             assert topology._isBonded(atom_indices[0], atom_indices[1]), 'Atom indices %d and %d are not bonded in topology' % (atom_indices[0], atom_indices[1])
 
+            # Compute equilibrium bond length and spring constant.
             if bond.fractional_bondorder==None:
-                force.addBond(atom_indices[0], atom_indices[1], bond.length, bond.k)
-            # If this bond uses partial bond orders
+                [k, length] = [bond.k, bond.length]
             else:
+                # This bond uses partial bond orders
                 # Make sure forcefield asks for fractional bond orders
                 if not self.ff._use_fractional_bondorder:
                     raise ValueError("Error: your forcefield file does not request to use fractional bond orders in its header, but a harmonic bond attempts to use them.")
@@ -1343,12 +1508,25 @@ def createForce(self, system, topology, verbose=False, **kwargs):
                 if bond.fractional_bondorder=='interpolate-linear':
                     k = bond.k[0] + (bond.k[1]-bond.k[0])*(order-1.)
                     length = bond.length[0] + (bond.length[1]-bond.length[0])*(order-1.)
-                    force.addBond(atom_indices[0], atom_indices[1], length, k)
-                    if verbose: print("%64s" % "Added %s bond, order %.2f; length=%.2g; k=%.2g" % (bond.smirks, order, length, k))
                 else:
                     raise Exception("Partial bondorder treatment %s is not implemented." % bond.fractional_bondorder)
 
-        if verbose: print('%d bonds added' % (len(bonds)))
+            # Handle constraints.
+            if topology.atomPairIsConstrained(*atom_indices):
+                # Atom pair is constrained; we don't need to add a bond term.
+                skipped_constrained_bonds += 1
+                # Check if we need to add the constraint here to the equilibrium bond length.
+                if topology.atomPairIsConstrained(*atom_indices) is True:
+                    # Mark that we have now assigned a specific constraint distance to this constraint.
+                    topology.constrainAtomPair(atom_indices[0], atom_indices[1], length)
+                    # Add the constraint to the System.
+                    system.addConstraint(atom_indices[0], atom_indices[1], length)
+                continue
+
+            # Add bond
+            force.addBond(atom_indices[0], atom_indices[1], length, k)
+
+        if verbose: print('%d bonds added (%d skipped due to constraints)' % (len(bonds) - skipped_constrained_bonds, skipped_constrained_bonds))
 
         # Check that no topological bonds are missing force parameters
         _check_for_missing_valence_terms('HarmonicBondForce', topology, bonds.keys(), topology.bonds())
@@ -1455,14 +1633,20 @@ def createForce(self, system, topology, verbose=False, **kwargs):
             print('')
 
         # Add all angles to the system.
+        skipped_constrained_angles = 0 # keep track of how many angles were constrained (and hence skipped)
         for (atom_indices, angle) in angles.items():
             # Ensure atoms are actually bonded correct pattern in Topology
             assert topology._isBonded(atom_indices[0], atom_indices[1]), 'Atom indices %d and %d are not bonded in topology' % (atom_indices[0], atom_indices[1])
             assert topology._isBonded(atom_indices[1], atom_indices[2]), 'Atom indices %d and %d are not bonded in topology' % (atom_indices[1], atom_indices[2])
 
+            if topology.atomPairIsConstrained(atom_indices[0], atom_indices[1]) and topology.atomPairIsConstrained(atom_indices[1], atom_indices[2]) and topology.atomPairIsConstrained(atom_indices[0], atom_indices[2]):
+                # Angle is constrained; we don't need to add an angle term.
+                skipped_constrained_angles += 1
+                continue
+
             force.addAngle(atom_indices[0], atom_indices[1], atom_indices[2], angle.angle, angle.k)
 
-        if verbose: print('%d angles added' % (len(angles)))
+        if verbose: print('%d angles added (%d skipped due to constraints)' % (len(angles) - skipped_constrained_angles, skipped_constrained_angles))
 
         # Check that no topological angles are missing force parameters
         _check_for_missing_valence_terms('HarmonicAngleForce', topology, angles.keys(), topology.angles())