Merge pull request #386 from epretti/charmm-naming-patching

Fixes and improvements for CHARMM force fields

Author: epretti

charmm/README.md

@@ -86,6 +86,22 @@ generated, the TOPPAR files should be reviewed to locate the inconsistency.  You
 may wish to exclude one or another set of parameters at this point to ensure
 that the resulting force field is correct.
 
+### Notice
+
+Due to various outstanding issues in ParmEd, including:
+
+* [ParmEd/ParmEd#1383](https://github.com/ParmEd/ParmEd/issues/1383)
+* [ParmEd/ParmEd#1384](https://github.com/ParmEd/ParmEd/issues/1384)
+* [ParmEd/ParmEd#1395](https://github.com/ParmEd/ParmEd/issues/1395)
+* [ParmEd/ParmEd#1396](https://github.com/ParmEd/ParmEd/issues/1396)
+* [ParmEd/ParmEd#1397](https://github.com/ParmEd/ParmEd/issues/1397)
+
+it is not currently possible to use a prepackaged version of ParmEd or a version
+from its GitHub repository to perform the CHARMM conversion.  Until these issues
+are fixed, a custom patched version of ParmEd containing workarounds for these
+problems specific to this conversion workflow can be found at
+[`patched-for-charmm-conversion` branch of the epretti/ParmEd repository](https://github.com/epretti/ParmEd/tree/patched-for-charmm-conversion).
+
 CHARMM docker image
 ===================
 

charmm/convert_charmm.py

@@ -119,12 +119,12 @@ def convert_yaml(yaml_filename, ffxml_dir):
         if name_collisions:
             raise ValueError(f"Collision between residue/patch names {sorted(name_collisions)}")
 
-        # Ensure that no atom names contain hyphens, as this will confuse the
-        # improper handling script (by the time these names reach the script,
-        # they will have been prepended with residue or patch names joined with
-        # a hyphen).
+        # Ensure that no names contain hyphens, as this will confuse the
+        # improper/anisotropy handling scripts.
         for templates in (params.residues.values(), params.patches.values()):
             for template in templates:
+                if "-" in template.name:
+                    raise ValueError(f"Forbidden character in template {template.name}")
                 for atom in template.atoms:
                     if "-" in atom.name:
                         raise ValueError(f"Forbidden character in template {template.name} atom {atom.name}")
@@ -485,13 +485,13 @@ def preprocess_improper_types(improper_types, prepare_parameters):
 
     script = etree.SubElement(ffxml_tree.getroot(), "Script")
     script.text = script_template.format(
-        residue_impropers=repr(residue_impropers),
-        patch_impropers=repr(patch_impropers),
-        delete_impropers=repr(delete_impropers),
-        residue_order=repr(residue_order),
-        patch_order=repr(patch_order),
-        periodic_improper_types=repr(periodic_improper_types),
-        harmonic_improper_types=repr(harmonic_improper_types),
+        residue_impropers=format_dict(residue_impropers),
+        patch_impropers=format_dict(patch_impropers),
+        delete_impropers=format_dict(delete_impropers),
+        residue_order=format_dict(residue_order),
+        patch_order=format_dict(patch_order),
+        periodic_improper_types=format_dict(periodic_improper_types),
+        harmonic_improper_types=format_dict(harmonic_improper_types),
     )
 
 
@@ -533,13 +533,17 @@ def preprocess_template_anisotropies(template_anisotropies):
 
     script = etree.SubElement(ffxml_tree.getroot(), "Script")
     script.text = script_template.format(
-        residue_anisotropies=repr(residue_anisotropies),
-        patch_anisotropies=repr(patch_anisotropies),
-        delete_anisotropies=repr(delete_anisotropies),
-        residue_order=repr(residue_order),
-        patch_order=repr(patch_order),
+        residue_anisotropies=format_dict(residue_anisotropies),
+        patch_anisotropies=format_dict(patch_anisotropies),
+        delete_anisotropies=format_dict(delete_anisotropies),
+        residue_order=format_dict(residue_order),
+        patch_order=format_dict(patch_order),
     )
 
 
+def format_dict(d):
+    return "\n".join(["{"] + [f"    {k!r}: {v!r}," for k, v in d.items()] + ["}"])
+
+
 if __name__ == "__main__":
     main()

charmm/convert_charmm_anisotropy_script.txt

@@ -7,25 +7,54 @@ delete_anisotropies = {delete_anisotropies}
 residue_order = {residue_order}
 patch_order = {patch_order}
 
+def extract_atom_name(raw_atom_name, templates):
+    # Process Drude particles.
+    atom_name_prefix = ""
+    if raw_atom_name.startswith("Drude-"):
+        raw_atom_name = raw_atom_name[len("Drude-"):]
+        atom_name_prefix = "D"
+
+    # Return None to indicate an atom without template name information.
+    atom_name_parts = raw_atom_name.rsplit("-", maxsplit=1)
+    if len(atom_name_parts) < 2:
+        return None
+    templates.add(atom_name_parts[0])
+    return f"{{atom_name_prefix}}{{atom_name_parts[1]}}"
+
 residue_data = []
 
-for residue in topology.residues():
-    # Extract the residue or patch names and atom names from the atom type names.
-    templates = []
-    atom_names = []
-    for atom in residue.atoms():
-        raw_atom_name = data.atomType[atom]
-        atom_is_drude = raw_atom_name.startswith("Drude-")
-        if atom_is_drude:
-            atom_name_parts = raw_atom_name[len("Drude-"):].rsplit("-", maxsplit=1)
-            atom_names.append(f"D{{atom_name_parts[1]}}")
-        else:
-            atom_name_parts = raw_atom_name.rsplit("-", maxsplit=1)
-            atom_names.append(atom_name_parts[1])
-        templates.append(atom_name_parts[0])
+for residue_index, residue in enumerate(topology.residues()):
+    # Extract the residue or patch names and atom names.
+    template_data = templateForResidue[residue_index]
+    templates = set()
+    if template_data is None:
+        # Multi-residue patch; fall back to parsing atom names because OpenMM
+        # leaves None in templateForResidue in this case.
+        atom_names = []
+        skip_residue = False
+        for atom in residue.atoms():
+            atom_name = extract_atom_name(data.atomType[atom], templates)
+            if atom_name is None:
+                # Skip residues containing atoms without template name data
+                # (they might be from, e.g., another solvent force field file).
+                skip_residue = True
+                break
+            atom_names.append(atom_name)
+        if skip_residue:
+            residue_data.append(([], [], []))
+            continue
+    else:
+        # Extract names from template name.
+        for template_index, template_name in enumerate(template_data.name.split("-")):
+            if template_index:
+                # Patch name.
+                templates.add(template_name.rsplit("_", maxsplit=1)[0])
+            else:
+                # Residue name.
+                templates.add(template_name)
+        atom_names = [template_data.atoms[data.atomTemplateIndexes[atom]].name for atom in residue.atoms()]
 
     # Get all unique residue or patch names that have anisotropies.
-    templates = set(templates)
     residue_templates = sorted(
         (template for template in templates if template in residue_anisotropies), key=residue_order.get
     )

charmm/convert_charmm_improper_script.txt

@@ -12,6 +12,20 @@ patch_order = {patch_order}
 periodic_improper_types = {periodic_improper_types}
 harmonic_improper_types = {harmonic_improper_types}
 
+def extract_atom_name(raw_atom_name, templates):
+    # Process Drude particles.
+    atom_name_prefix = ""
+    if raw_atom_name.startswith("Drude-"):
+        raw_atom_name = raw_atom_name[len("Drude-"):]
+        atom_name_prefix = "D"
+
+    # Return None to indicate an atom without template name information.
+    atom_name_parts = raw_atom_name.rsplit("-", maxsplit=1)
+    if len(atom_name_parts) < 2:
+        return None
+    templates.add(atom_name_parts[0])
+    return f"{{atom_name_prefix}}{{atom_name_parts[1]}}"
+
 def is_valid_improper(atom_indices):
     index_1, index_2, index_3, index_4 = atom_indices
     bonded_to_1 = data.bondedToAtom[index_1]
@@ -34,23 +48,38 @@ def find_improper(improper_types, atom_classes):
 
 residue_data = []
 
-for residue in topology.residues():
-    # Extract the residue or patch names and atom names from the atom type names.
-    templates = []
-    atom_names = []
-    for atom in residue.atoms():
-        raw_atom_name = data.atomType[atom]
-        atom_is_drude = raw_atom_name.startswith("Drude-")
-        if atom_is_drude:
-            atom_name_parts = raw_atom_name[len("Drude-"):].rsplit("-", maxsplit=1)
-            atom_names.append(f"D{{atom_name_parts[1]}}")
-        else:
-            atom_name_parts = raw_atom_name.rsplit("-", maxsplit=1)
-            atom_names.append(atom_name_parts[1])
-        templates.append(atom_name_parts[0])
+for residue_index, residue in enumerate(topology.residues()):
+    # Extract the residue or patch names and atom names.
+    template_data = templateForResidue[residue_index]
+    templates = set()
+    if template_data is None:
+        # Multi-residue patch; fall back to parsing atom names because OpenMM
+        # leaves None in templateForResidue in this case.
+        atom_names = []
+        skip_residue = False
+        for atom in residue.atoms():
+            atom_name = extract_atom_name(data.atomType[atom], templates)
+            if atom_name is None:
+                # Skip residues containing atoms without template name data
+                # (they might be from, e.g., another solvent force field file).
+                skip_residue = True
+                break
+            atom_names.append(atom_name)
+        if skip_residue:
+            residue_data.append(([], [], []))
+            continue
+    else:
+        # Extract names from template name.
+        for template_index, template_name in enumerate(template_data.name.split("-")):
+            if template_index:
+                # Patch name.
+                templates.add(template_name.rsplit("_", maxsplit=1)[0])
+            else:
+                # Residue name.
+                templates.add(template_name)
+        atom_names = [template_data.atoms[data.atomTemplateIndexes[atom]].name for atom in residue.atoms()]
 
     # Get all unique residue or patch names that have impropers.
-    templates = set(templates)
     residue_templates = sorted(
         (template for template in templates if template in residue_impropers), key=residue_order.get
     )

charmm/files/charmm36.yaml

@@ -77,6 +77,11 @@
       # that causes conflicting atom type names to be generated
       - DMPR
 
+      # Names are incompatible with improper lookup logic
+      - SB3-10
+      - SB3-12
+      - SB3-14
+
     exclude_patches:
       # Exclude CNH3 (CYTP to protonated cytosine) as it is incorrectly marked
       # as applying elsewhere
@@ -133,6 +138,14 @@
               - ./Residues/Residue/AllowPatch[@name='DEO5_0']
 
     fixes:
+      # Remove proline-specific patch from non-proline residues
+      - action: remove
+        target:
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_0']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_1']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_2']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_3']
+
       # Here we manually set nucleic acid external bonds because CHARMM doesn't
       # provide a reliable way to get at this information and ParmEd fails to
       # guess it correctly for nucleic acids; we also fix RNA-to-DNA patches

charmm/files/charmm36_protein.yaml

@@ -29,6 +29,14 @@
         output: charmm36_protein_d.xml
 
     fixes:
+      # Remove proline-specific patch from non-proline residues
+      - action: remove
+        target:
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_0']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_1']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_2']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_3']
+
       # Here we manually add in the DISU (cysteine disulfide) patch to the force
       # field since ParmEd doesn't convert it automatically
       - action: append

charmm/files/drude2023.yaml

@@ -21,15 +21,20 @@
       # Won't support single amino acids since the terminal patches are
       # different and cause ambiguity for assigning templates to polypeptides,
       # but this shouldn't be a significant use case
+      - CNES
+      - CT2S
       - CTES
+      - GNTS
       - NTES
     split:
       - input:
         - toppar/drude/drude_toppar_2023/toppar_drude_d_aminoacids_2023a.str
         output: charmm_polar_2023_d.xml
     fixes:
       - action: remove
-        target: ./Residues/Residue[@name!='PRO']/AllowPatch[@name='CTEP_0']
+        target:
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='CTEP_0']
+          - ./Residues/Residue[@name!='PRO']/AllowPatch[@name='PROP_0']
   References:
    protein:
      - >-