inital draft

bin/martini_sour

@@ -77,18 +77,25 @@ def main():
     parser_conv_itp.add_argument('-v', dest='verbosity', action='count',
                                 help='Enable debug logging output. Can be given '
                                 'multiple times.', default=0)
+    parser_conv_itp.add_argument('-ff', dest='forcefield',type=str,
+                                help='Force-Field for which to do conversion.', default="martini3_beta")
 
     file_group = parser_conv_itp.add_argument_group('Input and output options')
     file_group.add_argument('-f', dest='inpath', required=False, type=Path,
-                            help='Input file (ITP)', nargs="*")
+                            help='Input file (ITP)')
     file_group.add_argument('-o', dest='outpath', type=Path,
                             help='Output ITP (ITP)')
 
     titrate_group = parser_conv_itp.add_argument_group('Titratable Bead definitions.')
-    titrate_group.add_argument('-bases', dest='bases', type=str, nargs='+',
-                           help='An enumeration of residues to convert to bases.')
-    titrate_group.add_argument('-acids', dest='bases', type=str, nargs='+',
-                           help='An enumeration of residues to convert to acids.')
+    titrate_group.add_argument('-tres', dest='specs', action='append',
+                                type=lambda s: s.split(':'), default=[],
+                                help='Add a modification to a residue. Desired '
+                                'modification is specified as, e.g. A-ASP45:ASP0. '
+                                'The format is <chain>-<resname><resid>:<modification>.'
+                                ' Elements of the specification can be omitted as '
+                                'required.')
+    titrate_group.add_argument('-pkas', dest='pKas', type=lambda x: x.split(':'), action='append',
+                           help='pKa per residue')
     titrate_group.add_argument('-auto', dest='auto', type=bool, default=False,
                            help='Identify acids/bases automatically from known building blocks.')
 
@@ -161,13 +168,6 @@ def main():
 
     args = parser.parse_args()
 
-    if args.list_blocks:
-        force_field = load_library("libs", args.list_blocks, [])
-        msg = 'The following Blocks are present in the following libraries: {}:'
-        print(msg.format(args.list_blocks))
-        for block in force_field.blocks:
-            print(block)
-        parser.exit()
 
     loglevels = {0: logging.INFO, 1: logging.DEBUG, 2: 5}
     LOGGER.setLevel(loglevels[args.verbosity])

martini_sour/data/martini3_beta/citations.bib

@@ -0,0 +1,10 @@
+@article{Martini3,
+abstract={The coarse-grained Martini force field is widely used in biomolecular simulations. Here we present the refined model, Martini 3 (http://cgmartini.nl), with an improved interaction balance, new bead types and expanded ability to include specific interactions representing, for example, hydrogen bonding and electronic polarizability. The updated model allows more accurate predictions of molecular packing and interactions in general, which is exemplified with a vast and diverse set of applications, ranging from oil/water partitioning and miscibility data to complex molecular systems, involving protein–protein and protein–lipid interactions and material science applications as ionic liquids and aedamers.},
+author={Souza, Paulo C T and Alessandri, Riccardo and Barnoud, Jonathan and Thallmair, Sebastian and Faustino, Ignacio and Grünewald, Fabian and Patmanidis, Ilias and Abdizadeh, Haleh and Bruininks, Bart M H and Wassenaar, Tsjerk A and Kroon, Peter C and Melcr, Josef and Nieto, Vincent and Corradi, Valentina and Khan, Hanif M and Domański, Jan and Javanainen, Matti and Martinez-Seara, Hector and Reuter, Nathalie and Best, Robert B and Vattulainen, Ilpo and Monticelli, Luca and Periole, Xavier and Tieleman, D Peter and de Vries, Alex H and Marrink, Siewert J},
+doi={10.1038/s41592-021-01098-3},                                             
+issn={1548-7105},                                                             
+journal={Nature Methods},                                                     
+title={{Martini 3: a general purpose force field for coarse-grained molecular dynamics}},
+url={https://doi.org/10.1038/s41592-021-01098-3},                             
+year={2021}                                                                   
+}

martini_sour/data/martini3_beta/links.ff

@@ -0,0 +1,94 @@
+[ link ]
+[ atoms ]
+BAS {"degree": "0", "atype":"N5d", "replace": {"suffix": "F"}}
+D0 {"replace": {"atype": "DB1"}}
+DP {"replace": {"atype": "DB2"}}
+[ bonds ]
+BAS D0 1 0.187 4000
+[ constraints ]
+BAS DP 1 0.000
+
+[ link ]
+[ atoms ]
+BAS {"degree": "1", "atype":"N5d", "replace": {"suffix": "R"}}
+D0 {"replace": {"atype": "DB1"}}
+DP {"replace": {"atype": "DB2"}}
+[ bonds ]
+BAS D0 1 0.187 4000
+[ constraints ]
+BAS DP 1 0.000
+[ angles ]
+D0  BAS  R   2  0  500
+
+[ link ]
+[ atoms ]
+BAS {"degree": "2", "atype":"N5d", "replace": {"suffix": ""}}
+D0 {"replace": {"atype": "DB1"}}
+DP {"replace": {"atype": "DB2"}}
+[ bonds ]
+BAS D0 1 0.187 4000
+[ constraints ]
+BAS DP 1 0.000
+[ angles ]
+D0  BAS  R   2  0  500
+
+[ link ]
+[ atoms ]
+BAS {"degree": "0", "atype":"SN6d", "replace": {"suffix": "F"}}
+D0 {"replace": {"atype": "DB4"}}
+DP {"replace": {"atype": "DB3"}}
+[ bonds ]
+BAS D0 1 0.187 4000
+[ constraints ]
+BAS DP 1 0.000
+
+[ link ]
+[ atoms ]
+BAS {"degree": "1", "atype":"SN6d", "replace": {"suffix": "R"}}
+D0 {"replace": {"atype": "DB4"}}
+DP {"replace": {"atype": "DB3"}}
+[ bonds ]
+BAS D0 1 0.187 4000
+[ constraints ]
+BAS DP 1 0.000
+[ angles ]
+D0  BAS  R   2  0  500
+
+[ link ]
+[ atoms ]
+BAS {"degree": "2", "atype":"SN6d", "replace": {"suffix": ""}}
+D0 {"replace": {"atype": "DB4"}}
+DP {"replace": {"atype": "DB3"}}
+[ bonds ]
+BAS D0 1 0.187 4000
+[ constraints ]
+BAS DP 1 0.000
+[ angles ]
+D0  BAS  R   2  0  500
+
+[ link ]
+[ atoms ]
+ACD {"atype":"P2", "replace": {"suffix": ""}}
+DN {"atype": "DA1"}
+DP {"atype": "DA2"}
+[ constraints ]
+ACD DN 1 0.166
+ACD DP 1 0.200
+[ angles ]
+DN  ACD  DP   2   126 500
+[ edges ]
+ACD DN
+ACD DP
+
+[ link ]
+[ atoms ]
+ACD {"atype":"P2", "replace": {"suffix": ""}}
+DN {"atype": "DA1"}
+DP {"atype": "DA2"}
+R { }
+[ angles ]
+DN  ACD  R    2  0.0  500
+[ edges ]
+ACD DN
+ACD DP
+ACD R

martini_sour/data/martini3_beta/modifications.ff

@@ -0,0 +1,23 @@
+[ citations ]
+Martini3
+TitrMartini
+
+[ modification ]
+acid
+[ atoms ]
+ACD {"replace": {"charge": 0.0}}
+DN  {"add": {"atype": "DA1", "charge": -1.0, "idx": 1}}
+DP  {"add": {"atype": "DA2", "charge": 0.0, "idx": 2}}
+[ edges ]
+ACD DN
+ACD DP
+
+[ modification ]
+base
+[ atoms ]
+BAS {"replace": {"charge": -1.0}}
+D0  {"add": {"charge": 0.0, "idx": 1}}
+DP  {"add": {"charge": 1.0, "idx": 2}}
+[ edges ]
+BAS D0
+BAS DP

martini_sour/src/conv_itp.py

@@ -11,7 +11,191 @@
 # 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.
+import os
+import sys
+import networkx as nx
+import vermouth
+from vermouth.forcefield import ForceField
+from vermouth.gmx.itp_read import read_itp
+from vermouth.graph_utils import make_residue_graph
+from vermouth.log_helpers import StyleAdapter, get_logger
+from vermouth.file_writer import open, DeferredFileWriter
+from vermouth.citation_parser import citation_formatter
+from vermouth.processors.annotate_mut_mod import AnnotateMutMod
+from vermouth.processors.processor import Processor
+from vermouth.processors.do_links import DoLinks
+from martini_sour import DATA_PATH
+
+LOGGER = StyleAdapter(get_logger(__name__))
+
+class AnnotatepKas(Processor):
+
+    def __init__(self, pKas, res_graph=None):
+        self.pKas = {res: pKa for res, pKa in pKas}
+        self.res_graph = res_graph
+
+    def run_molecule(self, molecule):
+        if not self.res_graph:
+            self.res_graph = make_residue_graph(molecule)
+
+        for resid in self.res_graph:
+            graph = self.res_graph.nodes[resid]['graph']
+            for node in graph.nodes:
+                if "modification" in graph.nodes[node]\
+                and graph.nodes[node]["atomname"] in ["BAS", "ACD"]:
+                    resname = graph.nodes[node]["resname"]
+                    pKa = self.pKas[resname]
+                    atomtype = graph.nodes[node]["atype"]
+                    suffix = graph.nodes[node]["suffix"]
+                    new_atype = atomtype + "_" + pKa + suffix
+                    molecule.nodes[node]["atype"] = new_atype
+
+def _relable_molecule(molecule, mapping):
+    relabled_graph = nx.relabel_nodes(molecule, mapping=mapping, copy=True)
+    new_molecule = molecule.copy()
+    new_molecule.add_nodes_from(relabled_graph.nodes(data=True))
+    new_molecule.remove_edges_from(molecule.edges)
+    new_molecule.add_edges_from(relabled_graph.edges)
+    return new_molecule
+
+
+class AddTitratableBeads(Processor):
+
+    def __init__(self, force_field, res_graph=None):
+        """
+        """
+        self.res_graph = res_graph
+        self.force_field = force_field
+
+    def insert_node(self, molecule, attrs, insert_idx):
+        """
+        Insert a node with attributes `attrs` into an
+        existing `molecule` at the position `insert_idx`.
+        All interactions are relabeled as well.
+        """
+        # establish a correspondance between the old nodes
+        # and the new nodes
+        nodes = list(molecule.nodes)
+        correspondance = dict(zip(nodes, nodes))
+        temp_idx = len(nodes)
+        for from_node, to_node in correspondance.items():
+            if from_node >= insert_idx:
+                correspondance[from_node] = to_node + 1
+        correspondance[temp_idx] = insert_idx
+        # add the new node at the end of the graph
+        molecule.add_node(temp_idx, **attrs)
+
+        # also add it to the residue graph
+        resid = attrs['resid'] - 1
+        self.res_graph.nodes[resid]['graph'].add_node(temp_idx, **attrs)
+
+        # relabel the molecule with the new nodes
+        molecule = _relable_molecule(molecule, mapping=correspondance)
+
+        # relabel residue graph
+        for res in self.res_graph:
+            nx.relabel_nodes(self.res_graph.nodes[res]['graph'], mapping=correspondance)
+
+        # update the interactions lists
+        for inter_type in molecule.interactions:
+            for inter in molecule.interactions[inter_type]:
+                new_atoms = tuple([correspondance[atom] for atom in inter.atoms])
+                inter._replace(atoms=new_atoms)
+
+        return molecule, correspondance, temp_idx
+
+    def replace_node(self, molecule, node, attrs):
+        """
+        Replace the attributes of an existing node.
+        """
+        for target_node in molecule.nodes:
+            if molecule.nodes[target_node]["atomname"] == node:
+                molecule.nodes[target_node].update(attrs)
+                return target_node
+
+    def convert_titratable_beads(self, molecule):
+        for res_node in self.res_graph:
+            graph = self.res_graph.nodes[res_node]['graph']
+            node = list(graph.nodes)[0]
+            if "modification" in graph.nodes[node]:
+                mod_name = graph.nodes[node]["modification"][0]
+                modf = self.force_field.modifications[mod_name]
+                modf_to_mol = {}
+                for node in modf.nodes:
+                    if "replace" in modf.nodes[node]:
+                        attrs = modf.nodes[node]["replace"]
+                        mol_node = self.replace_node(molecule, node, attrs)
+                        modf_to_mol[node] = mol_node
+                    elif "add" in modf.nodes[node]:
+                        attrs = modf.nodes[node]["add"]
+                        idx = attrs.pop("idx")
+                        insert_idx = idx + min(graph.nodes)
+                        attrs["resid"] = graph.nodes[res_node]["resid"]
+                        attrs["resname"] = graph.nodes[res_node]["resname"]
+                        attrs["atomname"] = node
+                        attrs["charge_group"] = graph.nodes[res_node]["charge_group"]
+                        molecule, mapping, mol_node = self.insert_node(molecule, attrs, insert_idx)
+                        modf_to_mol[node] = mol_node
+                        for node, mol_node in modf_to_mol.items():
+                            modf_to_mol[node] = mapping[mol_node]
+
+                for idx, jdx in modf.edges:
+                    molecule.add_edge(modf_to_mol[idx], modf_to_mol[jdx])
+
+        return molecule
+
+    def run_molecule(self, molecule):
+        if not self.res_graph:
+            self.res_graph = make_residue_graph(molecule)
+        molecule = self.convert_titratable_beads(molecule)
+        return molecule
 
 def conv_itp(args):
-    print("Not implemented yet.")
-    return None
+    # load force-field
+    force_field = ForceField(os.path.join(DATA_PATH,args.forcefield))
+
+    # read itp and load into molecule
+    with open(args.inpath) as _file:
+        lines = _file.readlines()
+
+    read_itp(lines, force_field)
+    mol_name = list(force_field.blocks.keys())[0]
+    mol = force_field.blocks[mol_name].to_molecule()
+
+    # set mol labels and make edges
+    mol.make_edges_from_interaction_type("bonds")
+    degrees = {node: str(deg) for node, deg in mol.degree()}
+    nx.set_node_attributes(mol, degrees, "degree")
+
+    # set specifications for residues to
+    # convert to titratable
+    if args.auto:
+        resspecs = ["ASP", "GLU", "LYS"]
+    else:
+        resspecs = args.specs
+
+    # annotate all titratable residues with modifications
+    AnnotateMutMod(modifications=resspecs).run_molecule(mol)
+
+    # convert beads to titratable beads
+    mol = AddTitratableBeads(force_field).run_molecule(mol)
+
+    # apply links
+    DoLinks().run_molecule(mol)
+
+    # annotate pKas
+    AnnotatepKas(pKas=args.pKas).run_molecule(mol)
+
+    # write the converted molecule
+    with open(args.outpath, 'w') as outpath:
+        header = [' '.join(sys.argv) + "\n"]
+        header.append("Please cite the following papers:")
+        for citation in mol.citations:
+            cite_string = citation_formatter(mol.force_field.citations[citation])
+            LOGGER.info("Please cite: " + cite_string)
+            header.append(cite_string)
+
+        vermouth.gmx.itp.write_molecule_itp(mol, outpath,
+                                            moltype=mol_name,
+                                            header=header)
+    DeferredFileWriter().write()

setup.cfg

@@ -33,7 +33,7 @@ install-requires =  # ?? requires-dist?
     numpy
     decorator == 4.4.2
     networkx ~= 2.0
-    vermouth >= 0.7.1
+    vermouth >= 0.7.2
     MDAnalysis
 zip-safe = False