Fix glycam ffxml and conversion script (#180)

* fix glycam

* fix imports

* fix OLP external bonds

* remove diffs in atom ordering

* remove atom ordering diffs introduced by last commit

* fix remaining atom ordering diffs

amber/convert_amber.py

@@ -3,8 +3,7 @@
 from __future__ import print_function, division
 from io import StringIO
 import parmed
-from parmed.utils.six import iteritems
-from parmed.utils.six.moves import StringIO, zip
+import openmm
 import openmm.app as app
 import openmm.unit as u
 import os
@@ -1493,19 +1492,37 @@ def modify_glycan_ffxml(input_ffxml_path):
 
     # Add initialization script for setting up GlycamTemplateMatcher
     initialization_script = etree.SubElement(root, 'InitializationScript')
-    initialization_script.text = """
+    initialization_script.text = """ 
+from openmm.app.internal import compiled
+
 class GlycamTemplateMatcher(object):
+
   def __init__(self, glycam_residues):
     self.glycam_residues = glycam_residues
-  def __call__(self, ff, residue):
+
+  def __call__(self, ff, residue, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles):
     if residue.name in self.glycam_residues:
-      return ff._templates[residue.name]
+      template = ff._templates[residue.name]
+      if compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles) is not None:
+        return template
+
+      # The residue doesn't actually match the template with the same name.  Try the terminal variants.
+
+      if 'N'+residue.name in self.glycam_residues:
+        template = ff._templates['N'+residue.name]
+        if compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles) is not None:
+          return template
+      if 'C'+residue.name in self.glycam_residues:
+        template = ff._templates['C'+residue.name]
+        if compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles) is not None:
+          return template
     return None
 
 glycam_residues = set()
 for residue in tree.getroot().find('Residues').findall('Residue'):
   glycam_residues.add(residue.get('name'))
 self.registerTemplateMatcher(GlycamTemplateMatcher(glycam_residues))
+
 """
 
     tree.write(input_ffxml_path)

amber/ffxml/GLYCAM_06j-1.xml

@@ -1,6 +1,6 @@
 <ForceField>
   <Info>
-    <DateGenerated>2021-07-13</DateGenerated>
+    <DateGenerated>2021-10-29</DateGenerated>
     <Source Source="glycam/leaprc.GLYCAM_06j-1" md5hash="a7359f8d717f92dba02d57d6e0c8f3f7" sourcePackage="AmberTools" sourcePackageVersion="20.15">glycam/leaprc.GLYCAM_06j-1</Source>
     <Reference>R. Kadirvelraj; O. C. Grant; I. J. Goldstein; H. C. Winter; H. Tateno; E. Fadda; R. J. Woods. Structure and binding analysis of Polyporus squamosus lectin in complex with the Neu5Ac&#945;2-6Gal&#946;1-4GlcNAc human- type influenza receptor. Glycobiology, 2011, 21, 973&#8211;984. M. L. DeMarco; R. J. Woods. From agonist to antagonist: Structure and dynamics of innate immune glycoprotein MD-2 upon recognition of variably acylated bacterial endotoxins. Mol. Immunol., 2011, 49, 124&#8211;133. B. L. Foley; M. B. Tessier; R. J. Woods. Carbohydrate force fields. WIREs Comput. Mol. Sci., 2012, 2, 652&#8211;697. E. Ficko-Blean; C. P. Stuart; M. D. Suits; M. Cid; M. Tessier; R. J. Woods; A. B. Boraston. Carbohy- drate Recognition by an Architecturally Complex &#945;-N-Acetylglucosaminidase from Clostridium perfrin- gens. PLoS ONE, 2012, 7, e33524.</Reference>
   </Info>
@@ -40829,9 +40829,7 @@
       <Bond atomName1="CH3" atomName2="H2" />
       <Bond atomName1="CH3" atomName2="H3" />
       <Bond atomName1="CH3" atomName2="O" />
-      <ExternalBond atomName="H1" />
       <ExternalBond atomName="O" />
-      <ExternalBond atomName="CH3" />
     </Residue>
     <Residue name="PEA" override="1">
       <Atom name="C1" type="glycam-Cg" charge="0.48" />
@@ -45177,7 +45175,6 @@
       <Atom name="HO1" type="glycam-Ho" charge="0.445" />
       <Atom name="O1" type="glycam-Oh" charge="-0.639" />
       <Bond atomName1="HO1" atomName2="O1" />
-      <ExternalBond atomName="HO1" />
       <ExternalBond atomName="O1" />
     </Residue>
     <Residue name="RTA" override="1">
@@ -46421,9 +46418,6 @@
       <Bond atomName1="S1" atomName2="O2" />
       <Bond atomName1="S1" atomName2="O3" />
       <ExternalBond atomName="S1" />
-      <ExternalBond atomName="O1" />
-      <ExternalBond atomName="O2" />
-      <ExternalBond atomName="O3" />
     </Residue>
     <Residue name="STA" override="1">
       <Atom name="C1" type="glycam-Cg" charge="0.456" />
@@ -47229,7 +47223,6 @@
       <Bond atomName1="C4" atomName2="H7" />
       <Bond atomName1="C4" atomName2="H8" />
       <Bond atomName1="C4" atomName2="H9" />
-      <ExternalBond atomName="C2" />
       <ExternalBond atomName="O1" />
     </Residue>
     <Residue name="TDA" override="1">
@@ -66389,9 +66382,6 @@
       <Bond atomName1="C" atomName2="O" />
       <ExternalBond atomName="N" />
       <ExternalBond atomName="C" />
-      <ExternalBond atomName="CG" />
-      <ExternalBond atomName="OD1" />
-      <ExternalBond atomName="HOD" />
     </Residue>
     <Residue name="OLP" override="1">
       <Atom name="N" type="protein-N" charge="-0.2548" />
@@ -66424,7 +66414,7 @@
       <Bond atomName1="C" atomName2="O" />
       <ExternalBond atomName="N" />
       <ExternalBond atomName="C" />
-      <ExternalBond atomName="CG" />
+      <ExternalBond atomName="OD1" />
     </Residue>
     <Residue name="CHYP" override="1">
       <Atom name="N" type="protein-N" charge="-0.28" />
@@ -66460,10 +66450,6 @@
       <Bond atomName1="C" atomName2="O" />
       <Bond atomName1="C" atomName2="OXT" />
       <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
-      <ExternalBond atomName="CG" />
-      <ExternalBond atomName="OD1" />
-      <ExternalBond atomName="HD1" />
     </Residue>
     <Residue name="CNLN" override="1">
       <Atom name="N" type="protein-N" charge="-0.3821" />
@@ -66494,7 +66480,6 @@
       <Bond atomName1="C" atomName2="O" />
       <Bond atomName1="C" atomName2="OXT" />
       <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
       <ExternalBond atomName="ND2" />
     </Residue>
     <Residue name="COLP" override="1">
@@ -66529,8 +66514,7 @@
       <Bond atomName1="C" atomName2="O" />
       <Bond atomName1="C" atomName2="OXT" />
       <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
-      <ExternalBond atomName="CG" />
+      <ExternalBond atomName="OD1" />
     </Residue>
     <Residue name="COLS" override="1">
       <Atom name="N" type="protein-N" charge="-0.3821" />
@@ -66555,7 +66539,6 @@
       <Bond atomName1="C" atomName2="O" />
       <Bond atomName1="C" atomName2="OXT" />
       <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
       <ExternalBond atomName="OG" />
     </Residue>
     <Residue name="COLT" override="1">
@@ -66587,7 +66570,6 @@
       <Bond atomName1="C" atomName2="O" />
       <Bond atomName1="C" atomName2="OXT" />
       <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
       <ExternalBond atomName="OG1" />
     </Residue>
     <Residue name="NHYP" override="1">
@@ -66625,11 +66607,7 @@
       <Bond atomName1="CA" atomName2="HA" />
       <Bond atomName1="CA" atomName2="C" />
       <Bond atomName1="C" atomName2="O" />
-      <ExternalBond atomName="N" />
       <ExternalBond atomName="C" />
-      <ExternalBond atomName="CG" />
-      <ExternalBond atomName="OD1" />
-      <ExternalBond atomName="HD1" />
     </Residue>
     <Residue name="NNLN" override="1">
       <Atom name="N" type="protein-N3" charge="0.1801" />
@@ -66661,9 +66639,8 @@
       <Bond atomName1="CG" atomName2="ND2" />
       <Bond atomName1="ND2" atomName2="HD21" />
       <Bond atomName1="C" atomName2="O" />
-      <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
       <ExternalBond atomName="ND2" />
+      <ExternalBond atomName="C" />
     </Residue>
     <Residue name="NOLP" override="1">
       <Atom name="N" type="protein-N3" charge="-0.202" />
@@ -66698,9 +66675,8 @@
       <Bond atomName1="CA" atomName2="HA" />
       <Bond atomName1="CA" atomName2="C" />
       <Bond atomName1="C" atomName2="O" />
-      <ExternalBond atomName="N" />
+      <ExternalBond atomName="OD1" />
       <ExternalBond atomName="C" />
-      <ExternalBond atomName="CG" />
     </Residue>
     <Residue name="NOLS" override="1">
       <Atom name="N" type="protein-N3" charge="0.1849" />
@@ -66726,9 +66702,8 @@
       <Bond atomName1="CB" atomName2="HB3" />
       <Bond atomName1="CB" atomName2="OG" />
       <Bond atomName1="C" atomName2="O" />
-      <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
       <ExternalBond atomName="OG" />
+      <ExternalBond atomName="C" />
     </Residue>
     <Residue name="NOLT" override="1">
       <Atom name="N" type="protein-N3" charge="0.1812" />
@@ -66760,9 +66735,8 @@
       <Bond atomName1="CG2" atomName2="HG22" />
       <Bond atomName1="CG2" atomName2="HG23" />
       <Bond atomName1="C" atomName2="O" />
-      <ExternalBond atomName="N" />
-      <ExternalBond atomName="C" />
       <ExternalBond atomName="OG1" />
+      <ExternalBond atomName="C" />
     </Residue>
   </Residues>
   <HarmonicBondForce>
@@ -67918,17 +67892,34 @@ for force in sys.getForces():
           force.setExceptionParameters(i, p1, p2, atom_charges[p1]*atom_charges[p2], (atom_sigmas[p1]+atom_sigmas[p2])/2, unit.sqrt(atom_epsilons[p1]*atom_epsilons[p2]))
 </Script>
 <InitializationScript> 
+from openmm.app.internal import compiled
+
 class GlycamTemplateMatcher(object):
+
   def __init__(self, glycam_residues):
     self.glycam_residues = glycam_residues
-  def __call__(self, ff, residue):
+
+  def __call__(self, ff, residue, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles):
     if residue.name in self.glycam_residues:
-      return ff._templates[residue.name]
+      template = ff._templates[residue.name]
+      if compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles) is not None:
+        return template
+
+      # The residue doesn't actually match the template with the same name.  Try the terminal variants.
+
+      if 'N'+residue.name in self.glycam_residues:
+        template = ff._templates['N'+residue.name]
+        if compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles) is not None:
+          return template
+      if 'C'+residue.name in self.glycam_residues:
+        template = ff._templates['C'+residue.name]
+        if compiled.matchResidueToTemplate(residue, template, bondedToAtom, ignoreExternalBonds, ignoreExtraParticles) is not None:
+          return template
     return None
 
 glycam_residues = set()
 for residue in tree.getroot().find('Residues').findall('Residue'):
   glycam_residues.add(residue.get('name'))
 self.registerTemplateMatcher(GlycamTemplateMatcher(glycam_residues))
 
-    </InitializationScript></ForceField>
+</InitializationScript></ForceField>