diff --git a/vermouth/ffoutput.py b/vermouth/ffoutput.py new file mode 100644 index 000000000..48ab2a8fe --- /dev/null +++ b/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") diff --git a/vermouth/gmx/rtp.py b/vermouth/gmx/rtp.py index b4d4e92ed..6b8d80f21 100644 --- a/vermouth/gmx/rtp.py +++ b/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)