Implementation of the Martini3 Go-model

bin/martinize2

@@ -800,7 +800,7 @@ def entry():
         print(", ".join(known_force_fields[args.to_ff].modifications))
         parser.exit()
 
-    if args.elastic and args.govs_includes:
+    if args.elastic and args.go:
         parser.error(
             "A rubber band elastic network and GoMartini are not "
             "compatible. The -elastic and -govs-include flags cannot "
@@ -953,15 +953,13 @@ def entry():
         LOGGER.info("Reading Go model contact map.", type="step")
         go_map = read_go_map(args.go_map)
         LOGGER.info("Generating the Go model.", type="step")
-        GoPipeline(moltype=args.govs_moltype,
-                   contact_map=go_map,
-                   cutoff_short=args.go_low,
-                   cutoff_long=args.go_up,
-                   go_eps=args.go_eps,
-                   res_dist=args.go_res_dist,
-                   prefix=args.govs_moltype,
-                   res_graph=res_graph,
-                  ).run_system(system)
+        GoPipeline.run_system(system,
+                              moltype=args.govs_moltype,
+                              contact_map=go_map,
+                              cutoff_short=args.go_low,
+                              cutoff_long=args.go_up,
+                              go_eps=args.go_eps,
+                              res_dist=args.go_res_dist,)
 
         defines = ("GO_VIRT",)
         itp_paths = ["go_atomtypes.itp", "go_nbparams.itp"]

vermouth/gmx/topology.py

@@ -92,15 +92,18 @@
                     a1, a2 = nb_params.atoms
                 # self interaction
                 else:
                     a1 = nb_params.atoms[0]
                     a2 = nb_params.atoms[0]
 
                 if C6C12:
                     nb1, nb2 = convert_sigma_epsilon(nb_params.sigma, nb_params.epsilon)
                 else:
                     nb1, nb2 = nb_params.sigma, nb_params.epsilon
 
-                comments = ";" + " ".join(nb_params.meta['comment'])
+                if nb_params.meta.get('comment'):
+                    comments = ";" + " ".join(nb_params.meta['comment'])
+                else:
+                    comments = ""
                 itp_file.write(f"{a1} {a2} 1 {nb1:3.8F} {nb2:3.8F} {comments}\n")
 
 def write_gmx_topology(system, top_path, itp_paths=[], C6C12=False, defines=(), header=()):
@@ -119,21 +122,32 @@
     system: :class:vermouth.System
     top_path: pathlib.Path
         path for topology file
+    itp_paths: list[pathlib.Path]
+        list of paths for writing the topology parameters
+        like atomtypes. 
+    C6C12: bool
+        write non-bonded interaction parameters using LJ
+        C6C12 form
     defines: tuple(str)
         define statments to include in the topology
     header: tuple(str)
         any comment lines to include at the beginning
     """
     if not system.molecules:
         raise ValueError("No molecule in the system. Nothing to write.")
 
+    itp_paths = itp_paths[::-1]
     include_string = ""
-    # First we write the atomtypes and nonbondparams directive
-    if itp_paths:
-        write_atomtypes(system, itp_paths[0], C6C12)
-        write_nonbond_params(system, itp_paths[1], C6C12)
-        include_string = "\n".join('#include "{}"'.format(path) for path in itp_paths[:2])
-        include_string += "\n"
+    # First we write the atomtypes directive
+    if "atomtypes" in system.gmx_topology_params:
+        _path = itp_paths.pop()
+        write_atomtypes(system, _path, C6C12)
+        include_string += f'\n #include "{_path}"'
+    # Next we write the nonbond_params directive
+    if "nonbond_params" in system.gmx_topology_params:
+        _path = itp_paths.pop()
+        write_nonbond_params(system, _path, C6C12)
+        include_string += f'\n #include "{_path}"\n'
     # Write the ITP files for the molecule types, and prepare writing the
     # [ molecules ] section of the top file.
     # * We write one ITP file for each different moltype in the system, the
@@ -166,11 +180,11 @@
                 header[-1] = header[-1] + "\n"
                 header.append("Pleas cite the following papers:")
                 for citation in molecule.citations:
                     cite_string = citation_formatter(
                         molecule.force_field.citations[citation]
                     )
                     LOGGER.info("Please cite: " + cite_string)
                     header.append(cite_string)
                 vermouth.gmx.itp.write_molecule_itp(molecule, outfile, header=header)
             this_moltype_len = len(molecule.meta["moltype"])
             if this_moltype_len > max_name_length:

vermouth/processors/water_bias.py

@@ -14,15 +14,9 @@
 
 from .processor import Processor
 from ..graph_utils import make_residue_graph
-from ..rcsu.go_utils import get_go_type_from_attributes, _get_bead_size
+from ..rcsu.go_utils import get_go_type_from_attributes, _get_bead_size, _in_resid_region
 from ..gmx.topology import NonbondParam
 
-def _in_region(resid, regions):
-    for start, stop in regions:
-        if start <= resid <= stop:
-            return True
-    return False
-
 class ComputeWaterBias(Processor):
     """
     Processor which computes the water  bias for
@@ -91,13 +85,13 @@
             chain = res_graph.nodes[res_node]['chain']
             resname = res_graph.nodes[res_node]['resname']
 
-            if _in_region(resid, self.idr_regions):
+            if _in_resid_region(resid, self.idr_regions):
                 eps = self.water_bias.get('idr', 0.0)
             elif self.auto_bias:
                 sec_struc = res_graph.nodes[res_node]['cgsecstruct']
                 eps = self.water_bias.get(sec_struc, 0.0)
             else:
                 continue
 
             if eps == 0.0:
                 continue
@@ -126,13 +120,13 @@
         Assign water bias for a single molecule
         """
         if not self.system:
             raise IOError('This processor requires a system.')
 
         if not molecule.meta.get('moltype'):
             raise ValueError('The molecule does not have a moltype name.')
 
         if hasattr(molecule, 'res_graph'):
             res_graph = molecule.res_graph
         else:
             res_graph = make_residue_graph(molecule)
 

vermouth/rcsu/contact_map.py

@@ -42,8 +42,9 @@ def read_go_map(file_path):
             continue
 
         if tokens[0] == "R" and len(tokens) == 18:
-            contacts.append((int(tokens[1]), tokens[2], int(tokens[4]), tokens[5]))
+            contacts.append((int(tokens[5]), tokens[4], int(tokens[9]), tokens[8]))
 
     if len(contacts) == 0:
         raise IOError("You contact map is empty. Are you sure it has the right formatting?")
+
     return contacts

vermouth/rcsu/go_pipeline.py

@@ -25,35 +25,31 @@
     """
     Wrapping all processors for the go model.
     """
-    def __init__(self, processor_list, **kwargs):
+    def __init__(self, processor_list):
         self.processor_list = processor_list
-        self.kwargs = kwargs
+        self.kwargs = {}
 
-    def prepare_run(self, system):
+    def prepare_run(self, system, moltype):
         """
         Things to do before running the pipeline.
         """
         # merge all molecules in the system
         # this will eventually become deprecated
         # with the proper Go-model for multimers
         vermouth.MergeAllMolecules().run_system(system)
-        for molecule in system.molecules:
-            res_graph = vermouth.graph_utils.make_residue_graph(system.molecules[0])
-            molecule.residue_graph = res_graph
+        molecule = system.molecules[0]
+#        res_graph = vermouth.graph_utils.make_residue_graph(molecule)
+#        molecule.res_graph = res_graph
+        molecule.meta['moltype'] = moltype
 
-    def postprocess_run(self, system):
-        """
-        Do required post-processing.
-        """
-        pass
-
-    def run_system(self, system):
-        self.prepare_run(self, system)
+    def run_system(self, system, **kwargs):
+        self.kwargs = kwargs
+        self.prepare_run(system, moltype=kwargs['moltype'])
         for processor in self.processor_list:
             process_args = inspect.getfullargspec(processor).args
-            process_args_values = {self.kwargs[arg] for arg in process_args}
+            process_args_values = {arg:self.kwargs[arg] for arg in kwargs.keys() if arg in process_args}
             processor(**process_args_values).run_system(system)
         return system
 
 GoPipeline = GoProcessorPipline([SetMoleculeMeta,
                                  VirtualSideCreator,

vermouth/rcsu/go_structure_bias.py

@@ -16,9 +16,11 @@
 """
 import numpy as np
 import networkx as nx
+from ..graph_utils import make_residue_graph
 from ..molecule import Interaction
 from ..processors.processor import Processor
 from ..selectors import filter_minimal, select_backbone
+from ..gmx.topology import NonbondParam
 from .go_utils import get_go_type_from_attributes
 
 class ComputeStructuralGoBias(Processor):
@@ -48,8 +50,8 @@
                  cutoff_long,
                  go_eps,
                  res_dist,
-                 res_graph=None,
-                 prefix="VS"):
+                 moltype,
+                 res_graph=None):
         """
         Initialize the Processor with arguments required
         to setup the Go model structural bias.
@@ -71,21 +73,28 @@
             the residue graph they are ignored; this
             is similar to sequence distance but takes
             into account disulfide bridges for example
+        moltype: str
+            name of the molecule to treat
         res_graph: :class:vermouth.Molecule
             residue graph of the molecule; if None it
             get's generated automatically
-        prefix: str
-            prefix for all virtual-site atomtypes
+        system: :class:vermouth.System
+            the system
+        magic_number: float
+            magic number for Go contacts from the old
+            GoVirt script.
         """
         self.contact_map = contact_map
         self.cutoff_short = cutoff_short
         self.cutoff_long = cutoff_long
         self.go_eps = go_eps
         self.res_dist = res_dist
-        self.prefix = prefix
-        self.res_graph = res_graph
+        self.moltype = moltype
         # don't modify
+        self.res_graph = None
+        self.system = None
         self.__chain_id_to_resnode = {}
+        self.magic_number = 1.12246204830
 
     # do not overwrite when subclassing
     def _chain_id_to_resnode(self, chain, resid):
@@ -107,20 +116,20 @@
             a dict matching the chain,resid to the self.res_graph node
         """
         if self.__chain_id_to_resnode:
-            return self.__chain_id_to_resnode
+            return self.__chain_id_to_resnode[(chain, resid)]
 
         # for each residue collect the chain and residue in a dict
         # we use this later for identifying the residues from the
         # contact map
         for resnode in self.res_graph.nodes:
             chain = self.res_graph.nodes[resnode].get('chain', None)
             # in vermouth within a molecule all resid are unique
             # when merging multiple chains we store the old resid
             # the go model always references the input resid i.e.
             # the _old_resid
             resid = self.res_graph.nodes[resnode].get('_old_resid')
             self.__chain_id_to_resnode[(chain, resid)] = resnode
-        return self.__chain_id_to_resnode
+        return self.__chain_id_to_resnode[(chain, resid)]
 
     def contact_selector(self, molecule):
         """
@@ -140,47 +149,47 @@
             list of node keys and distance
         """
         # distance_matrix of elegible pairs as tuple(node, node, dist)
         contact_matrix = []
         # find all pairs of residues that are within bonded distance of
         # self.res_dist
         connected_pairs = dict(nx.all_pairs_shortest_path_length(self.res_graph,
                                                             cutoff=self.res_dist))
         for contact in self.contact_map:
             resIDA, chainA, resIDB, chainB = contact
             # identify the contact in the residue graph based on
             # chain ID and resid
-            resA = self.chain_id_to_resnode[(chainA, resIDA)]
-            resB = self.chain_id_to_resnode[(chainB, resIDB)]
+            resA = self._chain_id_to_resnode(chainA, resIDA)
+            resB = self._chain_id_to_resnode(chainB, resIDB)
             # make sure that both residues are not connected
             # note: contacts should be symmteric so we only
             # check against one
             if resB not in connected_pairs[resA]:
                 # now we lookup the backbone nodes within the residue contact
-                bb_node_A = next(filter_minimal(self.res_graph.nodes[resB], select_backbone))
-                bb_node_B = next(filter_minimal(self.res_graph.nodes[resB], select_backbone))
+                bb_node_A = next(filter_minimal(self.res_graph.nodes[resA]['graph'], select_backbone))
+                bb_node_B = next(filter_minimal(self.res_graph.nodes[resB]['graph'], select_backbone))
                 # compute the distance between bb-beads
                 dist = np.linalg.norm(molecule.nodes[bb_node_A]['position'] -
                                       molecule.nodes[bb_node_B]['position'])
                 # verify that the distance between BB-beads satisfies the
                 # cut-off criteria
-                if self.cut_off_large > dist > self.cutoff_short:
+                if self.cutoff_long > dist > self.cutoff_short:
                     # find the go virtual-sites for this residue
                     # probably can be done smarter but mehhhh
-                    atype_a = get_go_type_from_attributes(molecule,
-                                                          resid=resIDA,
-                                                          chain=chainA,
-                                                          prefix=self.prefix)
-                    atype_b = get_go_type_from_attributes(molecule,
-                                                          resid=resIDB,
-                                                          chain=chainB,
-                                                          prefix=self.prefix)
+                    atype_a = next(get_go_type_from_attributes(self.res_graph.nodes[resA]['graph'],
+                                                               resid=resIDA,
+                                                               chain=chainA,
+                                                               prefix=self.moltype))
+                    atype_b = next(get_go_type_from_attributes(self.res_graph.nodes[resB]['graph'],
+                                                               resid=resIDB,
+                                                               chain=chainB,
+                                                               prefix=self.moltype))
                     # generate backbone backbone exclusions
                     # perhaps one day can be it's own function
                     excl = Interaction(atoms=(bb_node_A, bb_node_B), parameters=[], meta={"group": "Go model exclusion"})
                     molecule.interactions['exclusions'].append(excl)
 
                     contact_matrix.append((atype_a, atype_b, dist))
         return contact_matrix
 
     def compute_go_interaction(self, contacts):
         """
@@ -198,15 +207,27 @@
         dict[frozenset(str, str): float]
             dict of interaction parameters indexed by atomtype
         """
         go_inters = {}
         for atype_a, atype_b, dist in contacts:
             # compute the LJ sigma paramter for this contact
             # 1.12246204830 is a magic number by Sebastian
-            sigma = dist / 1.12246204830
+            sigma = dist / self.magic_number
             # find the go virtual-sites for this residue
             # probably can be done smarter but mehhhh
-            go_inters[frozenset(atype_a, atype_b)] = (sigma, self.go_eps)
-        return go_inters
+            contact_bias = NonbondParam(atoms=(atype_a, atype_b),
+                                        sigma=sigma,
+                                        epsilon=self.go_eps,
+                                        meta={"comment": ["go bond {resid_a} {resid_b}"]})
+            self.system.gmx_topology_params["nonbond_params"].append(contact_bias)
+
+    def run_molecule(self, molecule):
+        self.res_graph = make_residue_graph(molecule)
+        # compute the contacts; this also creates
+        # the exclusions
+        contacts = self.contact_selector(molecule)
+        # compute the interaction parameters
+        self.compute_go_interaction(contacts)
+        return molecule
 
     def run_system(self, system):
         """
@@ -217,11 +238,6 @@
         system: vermouth.system.System
             The system to process. Is modified in-place.
         """
+        self.system = system
         for molecule in system.molecules:
-            # compute the contacts; this also creates
-            # the exclusions
-            contacts = self.contact_selector(molecule)
-            # compute the interaction parameters
-            inters = self.compute_go_interaction(contacts)
-            # update the contacts for the system
-            self.system.gmx_topology_params.upadte(inters)
+            self.run_molecule(molecule)

vermouth/rcsu/go_utils.py

@@ -44,10 +44,9 @@
         attrs = molecule.nodes[node]
         if attributes_match(attrs, kwargs) and attrs['atype'].startswith(prefix):
             yield attrs['atype']
-
     resid = kwargs['resid']
     chain = kwargs['chain']
     raise IOError(f"Could not find GoVs with resid {resid} in chain {chain}.")
 
 def _in_resid_region(resid, regions):
     """
@@ -64,7 +63,9 @@
     -------
     bool
     """
-    for low, up in regions:
+    for limits in regions:
+        # perhaps someone gives them as reversed
+        low, up = sorted(limits)
         if low <= resid <= up:
             return True
     return False

vermouth/rcsu/go_vs_includes.py

@@ -45,26 +45,26 @@
         Assign molecule type names to the molecules in a system.
     :func:`add_virtual_sites`
     """
-    def __init__(self, sections=('exclusions', )):
+    def __init__(self, sections=()):
         self.sections = sections
         self.system = None
 
     def run_molecule(self, molecule):
         moltype = molecule.meta.get('moltype')
         if not moltype:
             raise ValueError('The molecule does not have a moltype name.')
 
         if not self.system:
             raise IOError('This processor requires a system.')
 
         self.add_virtual_sites(molecule, prefix=moltype)
 
         includes = molecule.meta.get('post_section_lines', {})
         for section in self.sections:
             section_includes = includes.get(section, [])
             section_includes.append('#include "{moltype}_{section}_VirtGoSites.itp"'
                                     .format(moltype=moltype, section=section))
             includes[section] = section_includes
         molecule.meta['post_section_lines'] = includes
         return molecule
 
@@ -106,7 +106,7 @@
         # If there are no atoms, then there is nothing to do. We can exit early,
         # avoiding to deal with empty iterators.
         if not molecule.nodes:
             return
         virtual_site_nodes = []
         virtual_sites = []
         atomtypes = {}

vermouth/tests/gmx/test_itp.py

@@ -41,7 +41,11 @@ def dummy_molecule():
             'charge_group': 1, 'charge': 0, 'mass': 72
         }),
         (1, {
-            'atype': 'A', 'resid': 1, 'resname': 'X', 'atomname': 'A',
+            'atype': 'B', 'resid': 1, 'resname': 'X', 'atomname': 'B',
+            'charge_group': 1, 'charge': 0, 'mass': 72
+        }),
+        (2, {
+            'atype': 'C', 'resid': 1, 'resname': 'X', 'atomname': 'C',
             'charge_group': 1, 'charge': 0, 'mass': 72
         }),
     ))