add draft for ff-output

vermouth/ffoutput.py

@@ -0,0 +1,140 @@
+# Copyright 2020 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.
+
+"""
+Write .ff files.
+
+The FF file format describes molecule components for a given force field.
+
+The format is built on top of a subset of the ITP format. Describing a block
+is done in the same way an ITP file describes a molecule.
+"""
+
+class ForceFieldDirectiveWriter():
+    """
+    Write force-field files according to the
+    vermouth force-field definition.
+    """
+    def __init__(self, forcefield, stream):
+        """
+        Parameters
+        ----------
+        forcefield: `:class:vermouth.forcefield.ForceField`
+            the force-field object to write
+
+        stream: ``
+            the stream to which to write; must have a write method
+        """
+        self.forcefield = forcefield
+        self.stream = stream
+        # these attributes have a specific order in the moleculetype section
+        self.normal_order_block_atoms = ["atype", "resid", "resname",
+                                         "atomname", "charge_group", "charge"] #, "mass"]
+
+    def write(self):
+        """
+        Write the forcefield to file.
+        """
+        for name, block in self.forcefield.blocks.items():
+            self.stream.write("[ moleculetype ]\n")
+            excl = str(block.nrexcl)
+            self.stream.write(f"{name} {excl}\n")
+            self.write_atoms_block(block.nodes(data=True))
+            self.write_interaction_dict(block.interactions)
+
+        for link in self.forcefield.links:
+            self.write_link_header()
+            self.write_atoms_link(link.nodes(data=True))
+            self.write_interaction_dict(link.interactions)
+            self.write_edges(link.edges)
+
+    def write_interaction_dict(self, inter_dict):
+        """
+        Writes interactions to `self.stream`, with a new
+        interaction directive per type. Meta attributes
+        are kept and written as json parasable dicts.
+
+        Parameters
+        ----------
+        inter_dict: `class:dict[list[vermouth.molecule.Interaction]]`
+            the interaction dict to write
+        """
+        for inter_type in inter_dict:
+            self.stream.write(f"[ {inter_type} ]\n")
+            for interaction in inter_dict[inter_type]:
+                atom_string = " ".join(interaction.atoms)
+                param_string = " ".join(interaction.parameters)
+                meta_string = "{" + " ,".join([f"\"{key}\": \"{value}\"" for key, value in interaction.meta.items()]) + "}"
+                line = atom_string + " " + param_string + " " + meta_string + "\n"
+                self.stream.write(line)
+
+    def write_edges(self, edges):
+        """
+        Writes edges to `self.stream` into the edges directive.
+
+        Parameters
+        ----------
+        edges: abc.iteratable
+            pair-wise iteratable edge list
+        """
+        self.stream.write("[ edges ]\n")
+        for idx, jdx in edges:
+            self.stream.write(f"{idx} {jdx}\n")
+
+    def write_atoms_block(self, nodes):
+        """
+        Writes the nodes/atoms of the block atomtype directive to `self.stream`.
+        All attributes are written following the GROMACS atomtype directive
+        style.
+
+        Parameters
+        ----------
+        edges: abc.iteratable
+            pair-wise iteratable edge list
+        """
+        self.stream.write("[ atoms ]\n")
+        for idx, (node, attrs) in enumerate(nodes):
+            idx += 1
+            attr_line = " ".join([str(attrs[attr]) for attr in self.normal_order_block_atoms ])
+            line = f"{idx} " + attr_line + "\n"
+            self.stream.write(line)
+
+    def write_atoms_link(self, nodes):
+        """
+        Writes the nodes/atoms of the link atomtype directive to `self.stream`.
+        All attributes are written as json style dicts.
+
+        Parameters:
+        -----------
+        nodes: abc.itertable[tuple(abc.hashable, dict)]
+            list of nodes in form of a list with hashable node-key and dict
+            of attributes. The format is the same as returned by networkx.nodes(data=True)
+        """
+        self.stream.write("[ atoms ]\n")
+        for node_key, attributes  in nodes:
+            attr_line = " {" + " ,".join([f"\"{key}\": \"{value}\"" for key, value in attributes.items()]) + "}"
+            line = str(node_key) + attr_line + "\n"
+            self.stream.write(line)
+
+    def write_link_header(self):
+        """
+        Write the link directive header, with the resnames written
+        in form readable to geenerate a `:class:vermouth.molecule.Choice`
+        object.
+
+        Prameters
+        ---------
+        resnames: `abc.itertable[str]`
+        """
+        self.stream.write("[ link ]\n")

vermouth/gmx/rtp.py

@@ -260,21 +260,44 @@ def _complete_block(block, bondedtypes):
     # have some control over these dihedral angles through the bondedtypes
     # section.
     all_dihedrals = []
-    for center, dihedrals in itertools.groupby(
-            sorted(block.guess_dihedrals(), key=_dihedral_sorted_center),
-            _dihedral_sorted_center):
+    had_atoms = []
+    for center, dihedrals in itertools.groupby(sorted(block.guess_dihedrals(), key=_dihedral_sorted_center), _dihedral_sorted_center):
+        #print(center)
         if _keep_dihedral(center, block, bondedtypes):
             # TODO: Also sort the dihedrals by index.
             # See src/gromacs/gmxpreprocess/gen_add.cpp::dcomp in the
             # Gromacs source code (see version 2016.3 for instance).
-            atoms = sorted(dihedrals, key=_count_hydrogens)[0]
-            all_dihedrals.append(Interaction(atoms=atoms, parameters=[], meta={}))
+     #       atoms = sorted(dihedrals, key=_count_hydrogens)[0]
+            for atoms in dihedrals:
+                if atoms[::-1] not in had_atoms:
+                    all_dihedrals.append(Interaction(atoms=atoms, parameters=[], meta={}))
+                    had_atoms.append(atoms)
     # TODO: Sort the dihedrals by index
     block.interactions['dihedrals'] = (
         block.interactions.get('dihedrals', []) + all_dihedrals
     )
 
-    # TODO: generate 1-4 interactions between pairs of hydrogen atoms
+    # Add all the angles
+    all_angles = []
+    had_angles = []
+    for atoms in block.guess_angles():
+        if frozenset(atoms) not in had_angles:
+            all_angles.append(Interaction(atoms=atoms, parameters=[], meta={}))
+            had_angles.append(frozenset(atoms))
+
+    # TODO: Sort the dihedrals by index
+    block.interactions['angles'] = (block.interactions.get('angles', []) + all_angles)
+
+    # TODO: generate 1-4 interactions between all atoms and keep pairs of hydrogen atoms
+    # if HH14 is set to 1
+    pairs = []
+    for node in block.nodes:
+        paths = nx.single_source_shortest_path(G=block, source=node, cutoff=4)
+        neighbours = [node for node, path in paths.items() if 4 == len(path)]
+        for ndx in neighbours:
+             if frozenset([ndx, node]) not in pairs:
+                block.interactions['pairs'].append(Interaction(atoms=(node, ndx), parameters=[], meta={}))
+                pairs.append(frozenset([ndx, node]))
 
     # Add function types to the interaction parameters. This is done as a
     # post processing step to cluster as much interaction specific code
@@ -285,12 +308,13 @@ def _complete_block(block, bondedtypes):
     # twice. Yet, none of the RTP files distributed with Gromacs 2016.3 causes
     # issue.
     functypes = {
-        'bonds': bondedtypes.bonds,
-        'angles': bondedtypes.angles,
-        'dihedrals': bondedtypes.dihedrals,
-        'impropers': bondedtypes.impropers,
-        'exclusions': 1,
-        'cmap': 1,
+        'bonds': str(bondedtypes.bonds),
+        'angles': str(bondedtypes.angles),
+        'dihedrals': str(bondedtypes.dihedrals),
+        'impropers': str(bondedtypes.impropers),
+        'exclusions': '1',
+        'cmap': '1',
+        'pairs': '1',
     }
     for name, interactions in block.interactions.items():
         for interaction in interactions:
@@ -466,6 +490,8 @@ def _dihedral_sorted_center(atoms):
     #return sorted(atoms[1:-1])
     return atoms[1:-1]
 
+def _angle_sorted_center(atoms):
+    return atoms[1]
 
 def read_rtp(lines, force_field):
     """
@@ -519,6 +545,5 @@ def read_rtp(lines, force_field):
     # inter-residues information. We need to split the pre-blocks into
     # blocks and links.
     blocks, links = _split_blocks_and_links(pre_blocks)
-
     force_field.blocks.update(blocks)
     force_field.links.extend(links)