From a6a0f7b1ed4621236d78b92a480fab9b40e33d9d Mon Sep 17 00:00:00 2001 From: Peter Kroon Date: Mon, 25 Nov 2024 14:43:21 +0100 Subject: [PATCH] Split up rtp links to interaction sized bites --- vermouth/gmx/rtp.py | 111 ++++++++++++++++++++++++-------------------- 1 file changed, 61 insertions(+), 50 deletions(-) diff --git a/vermouth/gmx/rtp.py b/vermouth/gmx/rtp.py index 08ec20c5c..5d3bc0f4b 100644 --- a/vermouth/gmx/rtp.py +++ b/vermouth/gmx/rtp.py @@ -21,7 +21,9 @@ import networkx as nx from ..molecule import Block, Link, Interaction -from .. import utils + +from ..log_helpers import StyleAdapter, get_logger +LOGGER = StyleAdapter(get_logger(__name__)) __all__ = ['read_rtp'] @@ -293,13 +295,15 @@ def _split_blocks_and_links(pre_blocks): Split an individual pre-block into a block and a link. """ blocks = {} - links = [] + all_links = [] for name, pre_block in pre_blocks.items(): - block, link = _split_block_and_link(pre_block) + block, links = _split_block_and_link(pre_block) blocks[name] = block - if link: - links.append(link) - return blocks, links + if links: + for link in links: + if link: + all_links.append(link) + return blocks, all_links def _split_block_and_link(pre_block): @@ -320,7 +324,7 @@ def _split_block_and_link(pre_block): ------- block: Block All the intra-residue information. - link: Link + links: List[Link] All the inter-residues information. """ block = Block(force_field=pre_block.force_field) @@ -339,19 +343,24 @@ def _split_block_and_link(pre_block): except AttributeError: pass - link_atoms = set() + links = [] - for _, interactions in pre_block.interactions.items(): + for interaction_type, interactions in pre_block.interactions.items(): for interaction in interactions: # Atoms that are not defined in the atoms section but only in interactions # (because they're part of the neighbours), don't have an atomname atomnames = [pre_block.nodes[n_idx].get('atomname') or n_idx for n_idx in interaction.atoms] if any(atomname[0] in '+-' for atomname in atomnames): - link_atoms.update(interaction.atoms) - # Here we make 1 link consisting of all the atoms in interactions that - # contain at least one atom from a neighbouring residue, which may be an - # issue since it constrains where the link applies. - link = Link(pre_block.subgraph(link_atoms)) + link = Link(pre_block.subgraph(interaction.atoms)) + if not nx.is_connected(link): + LOGGER.info('Discarding {} between atoms {} for link' + ' defined in residue {} in force field {}' + ' because it is disconnected.', + interaction_type, ' '.join(interaction.atoms), + pre_block.name, pre_block.force_field.name, + type='inconsistent-data') + continue + links.append(link) block = pre_block block.remove_nodes_from([n for n in pre_block if pre_block.nodes[n].get('atomname', n)[0] in '+-']) @@ -368,46 +377,48 @@ def _split_block_and_link(pre_block): # RTP files convey the order by prefixing the names with + or -. We need to # get rid of these prefixes. order = {'+': +1, '-': -1} - relabel_mapping = {} - # (Deep)copy interactions, since relabel_nodes removes/adds nodes, which - # wipes out all interactions otherwise. - interactions = copy.deepcopy(link.interactions) - remove_attrs = 'atype charge charge_group'.split() - for idx, node in enumerate(link.nodes()): - atomname = node - if node[0] in '+-': - link.nodes[node]['order'] = order[node[0]] - atomname = atomname[1:] - else: - link.nodes[node]['order'] = 0 - link.nodes[node]['resname'] = block.name - for attr in remove_attrs: - if attr in link.nodes[node]: - del link.nodes[node][attr] - link.nodes[node]['atomname'] = atomname - relabel_mapping[node] = idx - nx.relabel_nodes(link, relabel_mapping, copy=False) - - # By relabelling the nodes, we lost the relations between interactions and - # nodes, so we need to relabel the atoms in the interactions - new_interactions = collections.defaultdict(list) - for name, interactions in interactions.items(): - for interaction in interactions: - atoms = tuple(relabel_mapping[atom] for atom in interaction.atoms) - if not any(link.nodes[node]['order'] for node in atoms): - continue - new_interactions[name].append(Interaction( - atoms=atoms, - parameters=interaction.parameters, - meta=interaction.meta - )) - link.interactions = new_interactions + + for link in links: + relabel_mapping = {} + # (Deep)copy interactions, since relabel_nodes removes/adds nodes, which + # wipes out all interactions otherwise. + interactions = copy.deepcopy(link.interactions) + remove_attrs = 'atype charge charge_group'.split() + for idx, node in enumerate(link.nodes()): + atomname = node + if node[0] in '+-': + link.nodes[node]['order'] = order[node[0]] + atomname = atomname[1:] + else: + link.nodes[node]['order'] = 0 + link.nodes[node]['resname'] = block.name + for attr in remove_attrs: + if attr in link.nodes[node]: + del link.nodes[node][attr] + link.nodes[node]['atomname'] = atomname + relabel_mapping[node] = idx + nx.relabel_nodes(link, relabel_mapping, copy=False) + + # By relabelling the nodes, we lost the relations between interactions and + # nodes, so we need to relabel the atoms in the interactions + new_interactions = collections.defaultdict(list) + for name, interactions in interactions.items(): + for interaction in interactions: + atoms = tuple(relabel_mapping[atom] for atom in interaction.atoms) + if not any(link.nodes[node]['order'] for node in atoms): + continue + new_interactions[name].append(Interaction( + atoms=atoms, + parameters=interaction.parameters, + meta=interaction.meta + )) + link.interactions = new_interactions + link.interactions = dict(link.interactions) # Revert the interactions back to regular dicts to avoid creating # keys when querying them. block.interactions = dict(block.interactions) - link.interactions = dict(link.interactions) - return block, link + return block, links def _clean_lines(lines):