Residue templates can specify atoms by name instead of index

docs-source/usersguide/application.rst

@@ -1889,12 +1889,12 @@ The residue template definitions look like this:
       <Atom name="HH33" type="710"/>
       <Atom name="C" type="712"/>
       <Atom name="O" type="713"/>
-      <Bond from="0" to="1"/>
-      <Bond from="1" to="2"/>
-      <Bond from="1" to="3"/>
-      <Bond from="1" to="4"/>
-      <Bond from="4" to="5"/>
-      <ExternalBond from="4"/>
+      <Bond atomName1="HH31" atomName2="CH3"/>
+      <Bond atomName1="CH3" atomName2="HH32"/>
+      <Bond atomName1="CH3" atomName2="HH33"/>
+      <Bond atomName1="CH3" atomName2="C"/>
+      <Bond atomName1="C" atomName2="O"/>
+      <ExternalBond atomName="C"/>
      </Residue>
      <Residue name="ALA">
       ...
@@ -1908,10 +1908,12 @@ contains the following tags:
 * An :code:`<Atom>` tag for each atom in the residue.  This specifies the
   name of the atom and its atom type.
 * A :code:`<Bond>` tag for each pair of atoms that are bonded to each
-  other.  The :code:`to` and :code:`from` attributes are the indices of
-  the two bonded atoms (starting from 0) in the order they were listed.  For
-  example, :code:`<Bond from="1" to="3"/>` describes a bond between atom CH3
-  and atom HH33.
+  other.  The :code:`atomName1` and :code:`atomName2` attributes are the names
+  of the two bonded atoms.  (Some older force fields use the alternate tags
+  :code:`to` and :code:`from` to specify the atoms by index instead of name.
+  This is still supported for backward compatibility, but specifying atoms by
+  name is recommended, since it makes the residue definition much easier to
+  understand.)
 * An :code:`<ExternalBond>` tag for each atom that will be bonded to an
   atom of a different residue.
 
@@ -1926,20 +1928,20 @@ as in the following example:
      <Atom name="H1" type="tip4pew-H"/>
      <Atom name="H2" type="tip4pew-H"/>
      <Atom name="M" type="tip4pew-M"/>
-     <VirtualSite type="average3" index="3" atom1="0" atom2="1" atom3="2"
+     <VirtualSite type="average3" siteName="M" atomName1="O" atomName2="H1" atomName3="H2"
          weight1="0.786646558" weight2="0.106676721" weight3="0.106676721"/>
-     <Bond from="0" to="1"/>
-     <Bond from="0" to="2"/>
+     <Bond atomName1="O" atomName2="H1"/>
+     <Bond atomName1="O" atomName2="H2"/>
     </Residue>
 
 Each :code:`<VirtualSite>` tag indicates an atom in the residue that should
 be represented with a virtual site.  The :code:`type` attribute may equal
 :code:`"average2"`\ , :code:`"average3"`\ , :code:`"outOfPlane"`\ , or
 :code:`"localCoords"`\ , which correspond to the TwoParticleAverageSite, ThreeParticleAverageSite,
-OutOfPlaneSite, and LocalCoordinatesSite classes respectively.  The :code:`index` attribute gives the
-index (starting from 0) of the atom to represent with a virtual site.  The atoms
-it is calculated based on are specified by :code:`atom1`\ , :code:`atom2`\ ,
-and (for virtual site classes that involve three atoms) :code:`atom3`\ .  The
+OutOfPlaneSite, and LocalCoordinatesSite classes respectively.  The :code:`siteName`
+attribute gives the name of the atom to represent with a virtual site.  The atoms
+it is calculated based on are specified by :code:`atomName1`\ , :code:`atomName2`\ ,
+and (for virtual site classes that involve three atoms) :code:`atomName3`\ .  The
 remaining attributes are specific to the virtual site class, and specify the
 parameters for calculating the site position.  For a TwoParticleAverageSite,
 they are :code:`weight1` and :code:`weight2`\ .  For a

wrappers/python/simtk/openmm/app/data/spce.xml

@@ -8,8 +8,8 @@
    <Atom name="O" type="spce-O"/>
    <Atom name="H1" type="spce-H"/>
    <Atom name="H2" type="spce-H"/>
-   <Bond from="0" to="1"/>
-   <Bond from="0" to="2"/>
+   <Bond atomName1="O" atomName2="H1"/>
+   <Bond atomName1="O" atomName2="H2"/>
   </Residue>
  </Residues>
  <HarmonicBondForce>

wrappers/python/simtk/openmm/app/data/tip3p.xml

@@ -8,8 +8,8 @@
    <Atom name="O" type="tip3p-O"/>
    <Atom name="H1" type="tip3p-H"/>
    <Atom name="H2" type="tip3p-H"/>
-   <Bond from="0" to="1"/>
-   <Bond from="0" to="2"/>
+   <Bond atomName1="O" atomName2="H1"/>
+   <Bond atomName1="O" atomName2="H2"/>
   </Residue>
  </Residues>
  <HarmonicBondForce>

wrappers/python/simtk/openmm/app/data/tip4pew.xml

@@ -10,9 +10,9 @@
    <Atom name="H1" type="tip4pew-H"/>
    <Atom name="H2" type="tip4pew-H"/>
    <Atom name="M" type="tip4pew-M"/>
-   <VirtualSite type="average3" index="3" atom1="0" atom2="1" atom3="2" weight1="0.786646558" weight2="0.106676721" weight3="0.106676721"/>
-   <Bond from="0" to="1"/>
-   <Bond from="0" to="2"/>
+   <VirtualSite type="average3" siteName="M" atomName1="O" atomName2="H1" atomName3="H2" weight1="0.786646558" weight2="0.106676721" weight3="0.106676721"/>
+   <Bond atomName1="O" atomName2="H1"/>
+   <Bond atomName1="O" atomName2="H2"/>
   </Residue>
  </Residues>
  <HarmonicBondForce>

wrappers/python/simtk/openmm/app/data/tip5p.xml

@@ -11,10 +11,10 @@
    <Atom name="H2" type="tip5p-H"/>
    <Atom name="M1" type="tip5p-M"/>
    <Atom name="M2" type="tip5p-M"/>
-   <VirtualSite type="outOfPlane" index="3" atom1="0" atom2="1" atom3="2" weight12="-0.34490826" weight13="-0.34490826" weightCross="-6.4437903"/>
-   <VirtualSite type="outOfPlane" index="4" atom1="0" atom2="1" atom3="2" weight12="-0.34490826" weight13="-0.34490826" weightCross="6.4437903"/>
-   <Bond from="0" to="1"/>
-   <Bond from="0" to="2"/>
+   <VirtualSite type="outOfPlane" siteName="M1" atomName1="O" atomName2="H1" atomName3="H2" weight12="-0.34490826" weight13="-0.34490826" weightCross="-6.4437903"/>
+   <VirtualSite type="outOfPlane" siteName="M2" atomName1="O" atomName2="H1" atomName3="H2" weight12="-0.34490826" weight13="-0.34490826" weightCross="6.4437903"/>
+   <Bond atomName1="O" atomName2="H1"/>
+   <Bond atomName1="O" atomName2="H2"/>
   </Residue>
  </Residues>
  <HarmonicBondForce>

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

@@ -149,18 +149,29 @@ class ForceField(object):
             for residue in root.find('Residues').findall('Residue'):
                 resName = residue.attrib['name']
                 template = ForceField._TemplateData(resName)
+                atomIndices = {}
                 for atom in residue.findall('Atom'):
                     params = {}
                     for key in atom.attrib:
                         if key not in ('name', 'type'):
                             params[key] = _convertParameterToNumber(atom.attrib[key])
-                    template.atoms.append(ForceField._TemplateAtomData(atom.attrib['name'], atom.attrib['type'], self._atomTypes[atom.attrib['type']][2], params))
+                    atomName = atom.attrib['name']
+                    if atomName in atomIndices:
+                        raise ValueError('Residue '+resName+' contains multiple atoms named '+atomName)
+                    atomIndices[atomName] = len(template.atoms)
+                    template.atoms.append(ForceField._TemplateAtomData(atomName, atom.attrib['type'], self._atomTypes[atom.attrib['type']][2], params))
                 for site in residue.findall('VirtualSite'):
-                    template.virtualSites.append(ForceField._VirtualSiteData(site))
+                    template.virtualSites.append(ForceField._VirtualSiteData(site, atomIndices))
                 for bond in residue.findall('Bond'):
-                    template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
+                    if 'atomName1' in bond.attrib:
+                        template.addBond(atomIndices[bond.attrib['atomName1']], atomIndices[bond.attrib['atomName2']])
+                    else:
+                        template.addBond(int(bond.attrib['from']), int(bond.attrib['to']))
                 for bond in residue.findall('ExternalBond'):
-                    b = int(bond.attrib['from'])
+                    if 'atomName' in bond.attrib:
+                        b = atomIndices[bond.attrib['atomName']]
+                    else:
+                        b = int(bond.attrib['from'])
                     template.externalBonds.append(b)
                     template.atoms[b].externalBonds += 1
                 self.registerResidueTemplate(template)
@@ -316,27 +327,32 @@ class ForceField(object):
 
     class _VirtualSiteData:
         """Inner class used to encapsulate data about a virtual site."""
-        def __init__(self, node):
+        def __init__(self, node, atomIndices):
             attrib = node.attrib
-            self.index = int(attrib['index'])
             self.type = attrib['type']
             if self.type == 'average2':
-                self.atoms = [int(attrib['atom1']), int(attrib['atom2'])]
+                numAtoms = 2
                 self.weights = [float(attrib['weight1']), float(attrib['weight2'])]
             elif self.type == 'average3':
-                self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
+                numAtoms = 3
                 self.weights = [float(attrib['weight1']), float(attrib['weight2']), float(attrib['weight3'])]
             elif self.type == 'outOfPlane':
-                self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
+                numAtoms = 3
                 self.weights = [float(attrib['weight12']), float(attrib['weight13']), float(attrib['weightCross'])]
             elif self.type == 'localCoords':
-                self.atoms = [int(attrib['atom1']), int(attrib['atom2']), int(attrib['atom3'])]
+                numAtoms = 3
                 self.originWeights = [float(attrib['wo1']), float(attrib['wo2']), float(attrib['wo3'])]
                 self.xWeights = [float(attrib['wx1']), float(attrib['wx2']), float(attrib['wx3'])]
                 self.yWeights = [float(attrib['wy1']), float(attrib['wy2']), float(attrib['wy3'])]
                 self.localPos = [float(attrib['p1']), float(attrib['p2']), float(attrib['p3'])]
             else:
                 raise ValueError('Unknown virtual site type: %s' % self.type)
+            if 'siteName' in attrib:
+                self.index = atomIndices[attrib['siteName']]
+                self.atoms = [atomIndices[attrib['atomName%d'%(i+1)]] for i in range(numAtoms)]
+            else:
+                self.index = int(attrib['index'])
+                self.atoms = [int(attrib['atom%d'%(i+1)]) for i in range(numAtoms)]
             if 'excludeWith' in attrib:
                 self.excludeWith = int(attrib['excludeWith'])
             else: