From 04c3bebb09a4f60a5099043ac4cc1ceec836750e Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Wed, 6 Jul 2022 11:27:09 +0200 Subject: [PATCH 1/4] add draft for ff-output --- vermouth/ffoutput.py | 140 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 vermouth/ffoutput.py diff --git a/vermouth/ffoutput.py b/vermouth/ffoutput.py new file mode 100644 index 00000000..12009b8f --- /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(link) + 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") From 0ca9cc64fe3e7843c06845efbbe83c5cecf28ba1 Mon Sep 17 00:00:00 2001 From: Fabian Grunewald <32294573+fgrunewald@users.noreply.github.com> Date: Wed, 31 Aug 2022 11:59:13 +0200 Subject: [PATCH 2/4] Skip EN generation for no-coordiantes If the selector returns no coordinates the EN code break. --- vermouth/processors/apply_rubber_band.py | 62 ++++++++++++------------ 1 file changed, 32 insertions(+), 30 deletions(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 4bf7fc4a..3f5233b9 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -303,36 +303,38 @@ def apply_rubber_band(molecule, selector, raise ValueError('All atoms from the selection must have coordinates. ' 'The following atoms do not have some: {}.' .format(' '.join(missing))) - coordinates = np.stack(coordinates) - distance_matrix = self_distance_matrix(coordinates) - constants = compute_force_constants(distance_matrix, lower_bound, - upper_bound, decay_factor, decay_power, - base_constant, minimum_force) - - connected = build_connectivity_matrix(molecule, res_min_dist, node_to_idx, - selected_nodes=selection) - - same_domain = build_pair_matrix(molecule, domain_criterion, idx_to_node, - selected_nodes=selection) - - can_be_linked = (~connected) & same_domain - # Multiply the force constant by 0 if the nodes cannot be linked. - constants *= can_be_linked - distance_matrix = distance_matrix.round(5) # For compatibility with legacy - for from_idx, to_idx in zip(*np.triu_indices_from(constants)): - # note the indices in the matrix are not anymore the idx of - # the full molecule but the subset of nodes in selection - from_key = idx_to_node[selection[from_idx]] - to_key = idx_to_node[selection[to_idx]] - force_constant = constants[from_idx, to_idx] - length = distance_matrix[from_idx, to_idx] - if force_constant > minimum_force: - molecule.add_interaction( - type_='bonds', - atoms=(from_key, to_key), - parameters=[bond_type, length, force_constant], - meta={'group': 'Rubber band'}, - ) + + if coordinates: + coordinates = np.stack(coordinates) + distance_matrix = self_distance_matrix(coordinates) + constants = compute_force_constants(distance_matrix, lower_bound, + upper_bound, decay_factor, decay_power, + base_constant, minimum_force) + + connected = build_connectivity_matrix(molecule, res_min_dist, node_to_idx, + selected_nodes=selection) + + same_domain = build_pair_matrix(molecule, domain_criterion, idx_to_node, + selected_nodes=selection) + + can_be_linked = (~connected) & same_domain + # Multiply the force constant by 0 if the nodes cannot be linked. + constants *= can_be_linked + distance_matrix = distance_matrix.round(5) # For compatibility with legacy + for from_idx, to_idx in zip(*np.triu_indices_from(constants)): + # note the indices in the matrix are not anymore the idx of + # the full molecule but the subset of nodes in selection + from_key = idx_to_node[selection[from_idx]] + to_key = idx_to_node[selection[to_idx]] + force_constant = constants[from_idx, to_idx] + length = distance_matrix[from_idx, to_idx] + if force_constant > minimum_force: + molecule.add_interaction( + type_='bonds', + atoms=(from_key, to_key), + parameters=[bond_type, length, force_constant], + meta={'group': 'Rubber band'}, + ) def always_true(*args, **kwargs): # pylint: disable=unused-argument From 527f0219443d7a0064b6d0fdf4c0c8fac724fae3 Mon Sep 17 00:00:00 2001 From: fgrunewald Date: Wed, 31 Aug 2022 14:55:20 +0200 Subject: [PATCH 3/4] skip EN if no coordinates qualify; raise warning for nans --- vermouth/processors/apply_rubber_band.py | 75 ++++++++++++++---------- 1 file changed, 43 insertions(+), 32 deletions(-) diff --git a/vermouth/processors/apply_rubber_band.py b/vermouth/processors/apply_rubber_band.py index 3f5233b9..c5ad1636 100644 --- a/vermouth/processors/apply_rubber_band.py +++ b/vermouth/processors/apply_rubber_band.py @@ -23,7 +23,9 @@ from .processor import Processor from .. import selectors from ..graph_utils import make_residue_graph +from ..log_helpers import StyleAdapter, get_logger +LOGGER = StyleAdapter(get_logger(__name__)) # the bond type of the RB DEFAULT_BOND_TYPE = 6 @@ -303,38 +305,47 @@ def apply_rubber_band(molecule, selector, raise ValueError('All atoms from the selection must have coordinates. ' 'The following atoms do not have some: {}.' .format(' '.join(missing))) - - if coordinates: - coordinates = np.stack(coordinates) - distance_matrix = self_distance_matrix(coordinates) - constants = compute_force_constants(distance_matrix, lower_bound, - upper_bound, decay_factor, decay_power, - base_constant, minimum_force) - - connected = build_connectivity_matrix(molecule, res_min_dist, node_to_idx, - selected_nodes=selection) - - same_domain = build_pair_matrix(molecule, domain_criterion, idx_to_node, - selected_nodes=selection) - - can_be_linked = (~connected) & same_domain - # Multiply the force constant by 0 if the nodes cannot be linked. - constants *= can_be_linked - distance_matrix = distance_matrix.round(5) # For compatibility with legacy - for from_idx, to_idx in zip(*np.triu_indices_from(constants)): - # note the indices in the matrix are not anymore the idx of - # the full molecule but the subset of nodes in selection - from_key = idx_to_node[selection[from_idx]] - to_key = idx_to_node[selection[to_idx]] - force_constant = constants[from_idx, to_idx] - length = distance_matrix[from_idx, to_idx] - if force_constant > minimum_force: - molecule.add_interaction( - type_='bonds', - atoms=(from_key, to_key), - parameters=[bond_type, length, force_constant], - meta={'group': 'Rubber band'}, - ) + + if not coordinates: + return + + coordinates = np.stack(coordinates) + if np.any(np.isnan(coordinates)): + LOGGER.warning("Found nan coordinates in molecule {}. " + "Will not generate an EN for it. ", + molecule.moltype, + type='unmapped-atom') + return + + distance_matrix = self_distance_matrix(coordinates) + constants = compute_force_constants(distance_matrix, lower_bound, + upper_bound, decay_factor, decay_power, + base_constant, minimum_force) + + connected = build_connectivity_matrix(molecule, res_min_dist, node_to_idx, + selected_nodes=selection) + + same_domain = build_pair_matrix(molecule, domain_criterion, idx_to_node, + selected_nodes=selection) + + can_be_linked = (~connected) & same_domain + # Multiply the force constant by 0 if the nodes cannot be linked. + constants *= can_be_linked + distance_matrix = distance_matrix.round(5) # For compatibility with legacy + for from_idx, to_idx in zip(*np.triu_indices_from(constants)): + # note the indices in the matrix are not anymore the idx of + # the full molecule but the subset of nodes in selection + from_key = idx_to_node[selection[from_idx]] + to_key = idx_to_node[selection[to_idx]] + force_constant = constants[from_idx, to_idx] + length = distance_matrix[from_idx, to_idx] + if force_constant > minimum_force: + molecule.add_interaction( + type_='bonds', + atoms=(from_key, to_key), + parameters=[bond_type, length, force_constant], + meta={'group': 'Rubber band'}, + ) def always_true(*args, **kwargs): # pylint: disable=unused-argument From 6c595dea107a790fbb36b68d30d6019f1321662e Mon Sep 17 00:00:00 2001 From: "f.grunewald" Date: Sun, 30 Apr 2023 14:24:59 +0200 Subject: [PATCH 4/4] patch rtp and ffoutput --- vermouth/ffoutput.py | 4 ++-- vermouth/gmx/rtp.py | 51 +++++++++++++++++++++++++++++++++----------- 2 files changed, 40 insertions(+), 15 deletions(-) diff --git a/vermouth/ffoutput.py b/vermouth/ffoutput.py index 12009b8f..48ab2a8f 100644 --- a/vermouth/ffoutput.py +++ b/vermouth/ffoutput.py @@ -40,7 +40,7 @@ def __init__(self, forcefield, stream): 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"] + "atomname", "charge_group", "charge"] #, "mass"] def write(self): """ @@ -54,7 +54,7 @@ def write(self): self.write_interaction_dict(block.interactions) for link in self.forcefield.links: - self.write_link_header(link) + self.write_link_header() self.write_atoms_link(link.nodes(data=True)) self.write_interaction_dict(link.interactions) self.write_edges(link.edges) diff --git a/vermouth/gmx/rtp.py b/vermouth/gmx/rtp.py index b4d4e92e..6b8d80f2 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)