(protein) termini patching

.gitignore

@@ -24,3 +24,5 @@ tests/test_data/gen_seq/output/*.itp
 doc/build
 doc/source/api
 doc/source/.doctrees
+
+/polyply/data/martini3007/*

bin/polyply

@@ -77,6 +77,10 @@ def main(): # pylint: disable=too-many-locals,too-many-statements
     dna_group = parser_gen_itp.add_argument_group('DNA specifc options')
     dna_group.add_argument('-dsdna', dest='dsdna', action='store_true',
                            help='complement single sequence to dsDNA sequence')
+    protein_termini_group = parser_gen_itp.add_argument_group('Protein termini group')
+    protein_termini_group.add_argument('-ter', dest='ter', action='store_true',
+                                       help='patch termini with correct parameters'
+                                       )
 
     parser_gen_itp.set_defaults(func=gen_itp)
 

polyply/data/martini3/aminoacids.ff

@@ -802,64 +802,3 @@ resname $protein_resnames
 CA  BB
 [ edges ]
 CA BB
-
-[ link ]
-resname $protein_resnames
-[ atoms ]
-CA { }
-BB { }
-+BB { }
-[ exclusions ]
-#meta {"comment": "CA-BB"}
-CA +BB
-[ edges ]
-BB +BB
-
-;; Protein termini. These links should be applied last.
-[ link ]
-resname $protein_resnames
-[ atoms ]
-BB {"replace": {"atype": "Q5", "charge": 1}}
-[ non-edges ]
-BB -BB
-
-[ link ]
-resname $protein_resnames
-[ atoms ]
-BB {"replace": {"atype": "Q5", "charge": -1}}
-[ non-edges ]
-BB +BB
-
-[ link ]
-resname $protein_resnames
-[ atoms ]
-BB {"replace": {"atype": "P6", "charge": 0}}
-[ non-edges ]
-BB +BB
-BB -BB
-
-[ link ]
-resname $protein_resnames
-[ molmeta ]
-neutral_termini true
-[ atoms ]
-BB {"replace": {"atype": "P5", "charge": 0}}
-[ non-edges ]
-BB -BB
-
-[ link ]
-resname $protein_resnames
-[ molmeta ]
-neutral_termini true
-[ atoms ]
-BB {"replace": {"atype": "P6", "charge": 0}}
-[ non-edges ]
-BB +BB
-
-;; Cystein bridge
-[ link ]
-resname "CYS"
-[ constraints ]
-SC1 >SC1 1 0.24 {"comment": "Disulfide bridge"}
-[ features ]
-disulfide

polyply/data/martini3/modifications.ff

@@ -0,0 +1,37 @@
+[ citations ]
+Martini3
+
+[ modification ]
+C-ter
+[ atoms ]
+BB {"replace": {"atype": "Q5", "charge": -1.0}}
+
+[ modification ]
+N-ter
+[ atoms ]
+BB {"replace": {"atype": "Q5", "charge": 1.0}}
+
+[ modification ]
+zwitter
+[ atoms ]
+BB {"replace": {"atype": "Q5"}}
+
+[ modification ]
+COOH-ter
+[ atoms ]
+BB {"replace": {"atype": "P6", "charge": 0.0}}
+
+[ modification ]
+NH2-ter
+[ atoms ]
+BB {"replace": {"atype": "P6", "charge": 0.0}}
+
+[ modification ]
+CCAP-ter
+[ atoms ]
+BB {"replace": {"atype": "P1", "charge": 0.0}}
+
+[ modification ]
+NCAP-ter
+[ atoms ]
+BB {"replace": {"atype": "P2", "charge": 0.0}}

polyply/src/apply_links.py

@@ -230,16 +230,16 @@ def match_link_and_residue_atoms(meta_molecule, link, link_to_resid):
         # relative resid has been asserted before so we can
         # exclude it here
         ignore = ['order', 'charge_group', 'replace', 'resid']
-        matchs = list(find_atoms(block, ignore=ignore, **attrs))
+        matches = list(find_atoms(block, ignore=ignore, **attrs))
 
-        if len(matchs) == 1:
-            link_to_mol[node] = matchs[0]
-        elif len(matchs) == 0:
-            msg = "Found no matchs for node {} in resiue {}. Cannot apply link."
+        if len(matches) == 1:
+            link_to_mol[node] = matches[0]
+        elif len(matches) == 0:
+            msg = "Found no matches for node {} in resiue {}. Cannot apply link."
             raise MatchError(msg.format(node, resid))
         else:
             msg = "Found {} matches for node {} in resiue {}. Cannot apply link."
-            raise MatchError(msg.format(len(matchs), node, resid))
+            raise MatchError(msg.format(len(matches), node, resid))
 
     return link_to_mol
 
@@ -475,8 +475,8 @@ def run_molecule(self, meta_molecule):
                                              res_link,
                                              node_match=_res_match,
                                              edge_match=_linktype_match)
-            raw_matchs = GM.subgraph_isomorphisms_iter()
-            for match in raw_matchs:
+            raw_matches = GM.subgraph_isomorphisms_iter()
+            for match in raw_matches:
                 nodes = match.keys()
                 resids =[meta_molecule.nodes[node]["resid"] for node in nodes]
                 orders = [ res_link.nodes[match[node]]["order"] for node in nodes]

polyply/src/apply_termini.py

@@ -0,0 +1,86 @@
+# Copyright 2024 University of Groningen
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#    http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+from vermouth.log_helpers import StyleAdapter, get_logger
+LOGGER = StyleAdapter(get_logger(__name__))
+
+def apply_terminal_mod(meta_molecule, term, res):
+    """
+
+    Apply a terminal modification.
+    Note this assumes that the modification is additive, ie. no atoms are
+    removed
+
+    Parameters
+    ----------
+    meta_molecule: :class:`polyply.src.meta_molecule.MetaMolecule`
+    term: str
+        name modification in forcefield to apply
+    res: int
+        resid of where to apply the modification
+
+    Returns
+    ----------
+    meta_molecule
+    """
+
+    molecule = meta_molecule.molecule
+
+    mod_atoms = {}
+    for i in molecule.force_field.modifications[term].atoms:
+        if 'replace' in i:
+            mod_atoms[i['atomname']] = i['replace']
+        else:
+            mod_atoms[i['atomname']] = {}
+
+
+    anum_dict = {}
+    for atomid, node in enumerate(meta_molecule.molecule.nodes):
+        if meta_molecule.molecule.nodes[node]['resid'] == res:
+            # Ensure unique resnames for unique blocks, necessary for coordinate building
+            # This does mean that the resnames are not automatically recognised as protein
+            # ones by gromacs, but what's fun if not making your own index file?
+            # Raise some info about it so people know
+            resname = meta_molecule.molecule.nodes[node]['resname']
+            meta_molecule.molecule.nodes[node]['resname'] = resname + term[0]
+
+            # add the modification from the modifications dict
+            aname = meta_molecule.molecule.nodes[node]['atomname']
+            if aname in mod_atoms.keys():
+                anum_dict[aname] = atomid
+                for key, value in mod_atoms[aname].items():
+                    meta_molecule.molecule.nodes[node][key] = value
+
+    mod_interactions = molecule.force_field.modifications[term].interactions
+    for i in mod_interactions:
+        for j in molecule.force_field.modifications[term].interactions[i]:
+            molecule.add_interaction(i,
+                                     (anum_dict[j.atoms[0]],
+                                      anum_dict[j.atoms[1]]),
+                                     j.parameters,
+                                     meta=j.meta)
+    return meta_molecule
+
+
+def prot_termini(meta_molecule):
+
+    max_res = max([meta_molecule.molecule.nodes[i]['resid'] for i in meta_molecule.molecule.nodes])
+    modifications = [('N-ter', 1),
+                     ('C-ter', max_res)]
+
+    LOGGER.info("Changing terminal resnames by their patch. Check your index file.")
+    for mod, res in modifications:
+        meta_molecule = apply_terminal_mod(meta_molecule, mod, res)
+
+    return meta_molecule

polyply/src/gen_itp.py

@@ -34,6 +34,7 @@
 from polyply.src.graph_utils import find_missing_edges
 from .load_library import load_ff_library
 from .gen_dna import complement_dsDNA
+from .apply_termini import prot_termini
 
 LOGGER = StyleAdapter(get_logger(__name__))
 
@@ -60,7 +61,10 @@ def split_seq_string(sequence):
         monomers.append(Monomer(resname=resname, n_blocks=n_blocks))
     return monomers
 
-def gen_params(name="polymer", outpath=Path("polymer.itp"), inpath=[], lib=None, seq=None, seq_file=None, dsdna=False):
+def gen_params(name="polymer", outpath=Path("polymer.itp"), inpath=[],
+               lib=None, seq=None, seq_file=None,
+               dsdna=False, ter=False):
+
     """
     Top level function for running the polyply parameter generation.
     Parameters seq and seq_file are mutually exclusive. Set the other
@@ -107,6 +111,9 @@ def gen_params(name="polymer", outpath=Path("polymer.itp"), inpath=[], lib=None,
     LOGGER.info("applying links between residues",  type="step")
     meta_molecule = ApplyLinks().run_molecule(meta_molecule)
 
+    if ter:
+        prot_termini(meta_molecule)
+
     # Raise warning if molecule is disconnected
     msg = "Missing a link between residue {idxA} {resA} and residue {idxB} {resB}."
     for missing in find_missing_edges(meta_molecule, meta_molecule.molecule):