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
+#
+# https://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)