Source code for

# Copyright 2018 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
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.

Handle the RTP format from Gromacs.

import collections
import itertools
import networkx as nx

from ..molecule import Block, Link, Interaction
from .. import utils

__all__ = ['read_rtp']

# Name of the subsections in RTP files.
# Names starting with a '_' are for internal use.
RTP_SUBSECTIONS = ('atoms', 'bonds', 'angles', 'dihedrals',
                   'impropers', 'cmap', 'exclusions',

_BondedTypes = collections.namedtuple(
    'bonds angles dihedrals impropers all_dihedrals nrexcl HH14 remove_dih'

class _IterRTPSubsectionLines:
    Iterate over the lines of an RTP file within a subsection.

    def __init__(self, parent):
        self.parent = parent
        self.lines = parent.lines
        self.running = True

    def __next__(self):
        if not self.running:
            raise StopIteration
        line = next(self.lines)
        if line.strip().startswith('['):
            self.running = False
            raise StopIteration
        return line

    def __iter__(self):
        return self

    def flush(self):
        Move the iterator after the last line of the subsection.
        for _ in self:

class _IterRTPSubsections:
    Iterate over the subsection of a RTP file within a section.

    For each subsections, yields its name and  an iterator over its lines.

    def __init__(self, parent):
        self.parent = parent
        self.lines = parent.lines
        self.buffer = collections.deque()
        self.current_subsection = None
        self.running = True

    def __next__(self):
        if not self.running:
            raise StopIteration
        if self.current_subsection is not None:
        if self.buffer:
            line = self.buffer.popleft()
            line = next(self.lines)
        stripped = line.strip()
        if stripped.startswith('['):
            # A section header looks like "[ name ]". It matches the following
            # regexp: r"^\s*\[\s*(?P<name>)\s*\]\s*$". Once stripped, the
            # trailing spaces at the beginning and at the end are removed so
            # the string starts with '[' and ends with ']'. The slicing remove
            # these brackets. The final call to strip remove the potential
            # white characters between the brackets and the section name.
            name = stripped[1:-1].strip()
            if name in RTP_SUBSECTIONS:
                subsection = _IterRTPSubsectionLines(self)
                self.current_subsection = subsection
                return name, subsection
            self.running = False
            raise StopIteration
        raise IOError('I am almost sure I should not be here...')

    def __iter__(self):
        return self

    def flush(self):
        Move the iterator after the last subsection of the section.
        for _ in self:

class _IterRTPSections:
    Iterate over the sections of a RTP file.

    For each section, yields the name of the sections and an iterator over its

    def __init__(self, lines):
        self.lines = lines
        self.buffer = collections.deque()
        self.current_section = None

    def __next__(self):
        if self.current_section is not None:
        if self.buffer:
            line = self.buffer.popleft()
            line = next(self.lines)
        stripped = line.strip()
        if stripped.startswith('['):
            name = stripped[1:-1].strip()
            section = _IterRTPSubsections(self)
            self.current_section = section
            # The "[ bondedtypes ]" is special in the sense that it does
            # not have subsection, but instead have its content directly
            # under it. This breaks the neat hierarchy the rest of the file
            # has. Here we restore the hierarchy by faking that the file
            # contains:
            #    [ bondedtypes ]
            #     [ _bondedtypes ]
            # For now, I do that on the basis of the section name. If other
            # sections are special that way, I'll detect them with a look
            # ahead.
            if name == 'bondedtypes':
                section.buffer.append(' [ _bondedtypes ]')
            return name, section
        raise IOError('Hum... There is a bug in the RTP reader.')

    def __iter__(self):
        return self

def _atoms(subsection, block):
    for line in subsection:
        name, atype, charge, charge_group = line.split()
        atom = {
            'atomname': name,
            'atype': atype,
            'charge': float(charge),
            'charge_group': int(charge_group),

def _base_rtp_parser(interaction_name, natoms):
    def wrapped(subsection, block):
        Parse the lines from a RTP subsection and populate the block.
        interactions = []
        for line in subsection:
            splitted = line.strip().split()
            atoms = splitted[:natoms]
            parameters = splitted[natoms:]
        block.interactions[interaction_name] = interactions
    return wrapped

def _parse_bondedtypes(section):
    # Default taken from
    # 'src/gromacs/gmxpreprocess/resall.cpp::read_resall' in the Gromacs
    # source code.
    defaults = _BondedTypes(bonds=1, angles=1, dihedrals=1,
                            impropers=1, all_dihedrals=0,
                            nrexcl=3, HH14=1, remove_dih=1)

    # The 'bondedtypes' section contains its line directly under it. In
    # order to match the hierarchy model of the rest of the file, the
    # iterator actually yields a subsection named '_bondedtypes'. We need
    # to read the fist line of that first virtual subsection.
    _, lines = next(section)
    line = next(lines)
    read = [int(x) for x in line.split()]

    # Fill with the defaults. The file gives the values in order so we
    # need to append the missing values from the default at the end.
    bondedtypes = _BondedTypes(*(read + list(defaults[len(read):])))

    # Make sure there is no unexpected lines in the section.
    # Come on Jonathan! There must be a more compact way of doing it.
    except StopIteration:
        raise IOError('"[ bondedtypes ]" section is missformated.')
    except StopIteration:
        raise IOError('"[ bondedtypes ]" section is missformated.')
    return bondedtypes

def _count_hydrogens(names):
    return len([name for name in names if utils.first_alpha(name) == 'H'])

def _keep_dihedral(center, block, bondedtypes):
    if (not bondedtypes.all_dihedrals) and block.has_dihedral_around(center):
        return False
    if bondedtypes.remove_dih and block.has_improper_around(center):
        return False
    return True

def _complete_block(block, bondedtypes):
    Add information from the bondedtypes section to a block.

    Generate implicit dihedral angles, and add function types to the

    # Generate missing dihedrals
    # As pdb2gmx generates all the possible dihedral angles by default,
    # RTP files are written assuming they will be generated. A RTP file
    # 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),
        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={}))
    # 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 function types to the interaction parameters. This is done as a
    # post processing step to cluster as much interaction specific code
    # into this method.
    # I am not sure the function type can be set explicitly in the RTP
    # file except through the bondedtypes section. If it is possible, then the
    # following code can break and atoms can have the function type written
    # 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,
    for name, interactions in block.interactions.items():
        for interaction in interactions:
            interaction.parameters.insert(0, functypes[name])

    # Set the nrexcl to the block.
    block.nrexcl = bondedtypes.nrexcl

def _split_blocks_and_links(pre_blocks):
    Split all the pre-blocks from `pre_block` into blocks and links.

    pre_blocks: dict
        A dict with residue names as keys and instances of :class:`Block`
        as values.

    blocks: dict
        A dict like `pre_block` with all inter-residues information
        stripped out.
    links: list
        A list of instances of :class:`Link` containing the inter-residues
        information from `pre_blocks`.

    See Also
        Split an individual pre-block into a block and a link.
    blocks = {}
    links = []
    for name, pre_block in pre_blocks.items():
        block, link = _split_block_and_link(pre_block)
        blocks[name] = block
    return blocks, links

def _split_block_and_link(pre_block):
    Split `pre_block` into a block and a link.

    A pre-block is a block as read from an RTP file. It contains both
    intra-residue and inter-residues information. This method split this
    information so that the intra-residue interactions are put in a block
    and the inter-residues ones are put in a link.

    pre_block: Block
        The block to split.

    block: Block
        All the intra-residue information.
    link: Link
        All the inter-residues information.
    block = Block(force_field=pre_block.force_field)
    link = Link()

    # It is easier to fill the interactions using a defaultdict,
    # yet defaultdicts are more annoying when reading as querying them
    # creates the keys. So the interactions are revert back to regular
    # dict at the end.
    block.interactions = collections.defaultdict(list)
    link.interactions = collections.defaultdict(list) =
        block.nrexcl = pre_block.nrexcl
    except AttributeError:

    # Filter the particles from neighboring residues out of the block.
    for atom in pre_block.atoms:
        if not atom['atomname'].startswith('+-'):
            atom['resname'] =

    # Create the edges of the link and block based on the edges in the pre-block.
    # This will create too many edges in the link, but the useless ones will be
    # pruned latter.
    block.add_edges_from(edge for edge in pre_block.edges
                         if not any(node[0] in '+-' for node in edge))

    # Split the interactions from the pre-block between the block (for
    # intra-residue interactions) and the link (for inter-residues ones).
    # The "relevant_atoms" set keeps track of what particles are
    # involved in the link. This will allow to prune the link without
    # iterating again through its interactions.
    relevant_atoms = set()
    for name, interactions in pre_block.interactions.items():
        for interaction in interactions:
            for_link = any(atom[0] in '+-' for atom in interaction.atoms)
            if for_link:

    # Prune the link to keep only the edges and particles that are
    # relevant.
    nodes = set(link.nodes())
    link.remove_nodes_from(nodes - relevant_atoms)
    # Some interactions do not generate nodes (impropers for instance). If a
    # node is only described in one of such interactions, then the node never
    # appears in the link. Here we make sure these nodes exists even if they
    # are note connected.

    # Atoms from a links are matched against a molecule based on its node
    # attributes. The name is a primary criterion, but other criteria can be
    # the residue name (`resname` key) or the residue order in the chain
    # (`order` key). The residue order is 0 when refering to the current
    # residue, +int to refer to residues after the current one in the sequence,
    # -int to refer to a previous residue in the sequence, and '*' for any
    # residue but the current one.
    # RTP files convey the order by prefixing the names with + or -. We need to
    # get rid of these prefixes.
    order = {'+': +1, '-': -1}
    relabel_mapping = {}
    for idx, node in enumerate(link.nodes()):
        atomname = node
        if node[0] in '+-':
            link.nodes[node]['order'] = order[node[0]]
            atomname = atomname[1:]
            link.nodes[node]['order'] = 0
            link.nodes[node]['resname'] =
        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 link.interactions.items():
        for interaction in interactions:
            atoms = tuple(relabel_mapping[atom] for atom in interaction.atoms)
    link.interactions = new_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

def _clean_lines(lines):
    # TODO: merge continuation lines
    for line in lines:
        splitted = line.split(';', 1)
        if splitted[0].strip():
            yield splitted[0]

def _dihedral_sorted_center(atoms):
    #return sorted(atoms[1:-1])
    return atoms[1:-1]

[docs] def read_rtp(lines, force_field): """ Read blocks and links from a Gromacs RTP file to populate a force field Parameters ---------- lines: An iterator over the lines of a RTP file (e.g. a file handle, or a list of string). force_field: vermouth.forcefield.ForceField The force field to populate in place. Raises ------ IOError Something in the file could not be parsed. """ _subsection_parsers = { 'atoms': _atoms, 'bonds': _base_rtp_parser('bonds', natoms=2), 'angles': _base_rtp_parser('angles', natoms=3), 'impropers': _base_rtp_parser('impropers', natoms=4), 'cmap': _base_rtp_parser('cmap', natoms=5), } # An RTP file contains both the blocks and the links. # We first read everything in "pre-blocks"; then we separate the # blocks from the links. pre_blocks = {} bondedtypes = None cleaned = _clean_lines(lines) for section_name, section in _IterRTPSections(cleaned): if section_name == 'bondedtypes': bondedtypes = _parse_bondedtypes(section) continue block = Block(force_field=force_field) pre_blocks[section_name] = block = section_name for subsection_name, subsection in section: if subsection_name in _subsection_parsers: _subsection_parsers[subsection_name](subsection, block) # Pre-blocks only contain the interactions that are explicitly # written in the file. Some are incomplete (missing implicit defaults) # or must be built from the "bondedtypes" rules. for pre_block in pre_blocks.values(): _complete_block(pre_block, bondedtypes) # At this point, the pre-blocks contain both the intra- and # 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)