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 ITP file format from Gromacs.

import copy
import itertools

__all__ = ['write_molecule_itp', ]

def _attr_has_not_none_attr(obj, attr):
    Raise a Value error is 'obj' does not have an attribute 'attr' or if its
    value is ``None``.
        value = getattr(obj, attr)
    except AttributeError:
        value = None
    if value is None:
        raise ValueError('{} has no attribute "{}".'.format(obj, attr))

def _interaction_sorting_key(interaction):
    ifdef = interaction.meta.get('ifdef')
    ifndef = interaction.meta.get('ifndef')
    if ifdef is not None and ifndef is not None:
        raise ValueError('An interaction cannot have both an "ifdef" '
                         'and an "ifndef" meta attribute.')
    if ifdef is not None:
        conditional = (ifdef, True)
    elif ifndef is not None:
        conditional = (ifndef, False)
        conditional = ()

    group = interaction.meta.get('group')
    if group is None:
        group = ''

    return (conditional, group)

[docs] def write_molecule_itp(molecule, outfile, header=(), moltype=None, post_section_lines=None, pre_section_lines=None): """ Write a molecule in ITP format. The molecule must have a `nrexcl` attribute. Each atom in the molecule must have at least the following keys: `atype`, `resid`, `resname`, `atomname`, and `charge_group`. Atoms can also have a `charge` and a `mass` key. If the `moltype` argument is not provided, then the molecule must have a "moltype" meta attribute. Parameters ---------- molecule: Molecule The molecule to write. See above for the minimal information the molecule must contain. outfile: io.TextIOBase The file in which to write. header:[str] List of lines to write as comment at the beginning of the file. The comment character and the new line should not be included as they will be added in the function. moltype: str, optional The molecule type. If set to `None` (default), the molecule type is read from the "moltype" key of `molecule.meta`. post_section_lines: dict[str,[str]], optional List of lines to write at the end of some sections of the file. The argument is passed as a dict with the keys being the name of the sections, and the values being the lists of lines. If the argument is set to `None`, the lines will be read from the "post_section_lines" key of `molecule.meta`. pre_section_lines: dict[str,[str]], optional List of lines to write at the beginning of some sections, just after the section header. The argument is formatted in the same way as `post_section_lines`. If the argument is set to `None`, the lines will be read from the "post_section_lines" key of `molecule.meta`. Raises ------ ValueError The molecule is missing required information. """ # Make sure the molecule contains the information required to write the # header. if moltype is None: moltype = molecule.meta.get('moltype') if moltype is None: raise ValueError('A molecule must have a moltype to write an ' 'ITP, provide it with the moltype argument, or ' 'with the "moltype" meta attribute of the molecule.') _attr_has_not_none_attr(molecule, 'nrexcl') # Make sure the molecule contains the information required to write the # [atoms] section. The charge and mass can be ommited, if so gromacs take # them from the [atomtypes] section of the ITP file. for attribute in ('atype', 'resid', 'resname', 'atomname', 'charge_group'): if not all([attribute in atom for _, atom in molecule.atoms]): raise ValueError('Not all atom have a {}.'.format(attribute)) # Get the maximum length of each atom field so we can align the fields. # Atom indexes are written as a consecutive series starting from 1. # The maximum index of a 0-based series is `len(x) - 1`; because the # series starts at 1, the maximum value is `len(x). max_length = {'idx': len(str(len(molecule)))} for attribute in ('atype', 'resid', 'resname', 'atomname', 'charge_group', 'charge', 'mass'): max_length[attribute] = max(len(str(atom.get(attribute, ''))) for _, atom in molecule.atoms) # Write the header. # We want to follow the header with an empty line, only if there is a # header. The `has_header` variable is needed in case `header` is a # generator, in which case we cannot know before hand if it contains lines. has_header = False for line in header: outfile.write('; {}\n'.format(line)) has_header = True if has_header: outfile.write('\n') outfile.write('[ moleculetype ]\n') outfile.write('{} {}\n\n'.format(moltype, molecule.nrexcl)) # Get the post- and pre- section lines. These lines are will be written at # the end or at the beginning of the relevant sections. if post_section_lines is None: post_section_lines = molecule.meta.get('post_section_lines', {}) if pre_section_lines is None: pre_section_lines = molecule.meta.get('pre_section_lines', {}) seen_sections = set() # The atoms in the [atoms] section must be consecutively numbered, yet # there is no guarantee that the molecule fulfill that constrain. # Therefore we renumber the atoms. The `correspondence` dict allows to # keep track of the correspondence between the original and the new # numbering so we can apply the renumbering to the interactions. # The resid and charge_group should also be consecutive, though this is # left as the user responsibility. Make sure residues and charge groups are # correctly numbered. correspondence = {} outfile.write('[ atoms ]\n') seen_sections.add('atoms') for line in pre_section_lines.get('atoms', []): outfile.write(line + '\n') # for idx, (original_idx, atom) in enumerate(molecule.atoms, start=1): for idx, original_idx in enumerate(molecule.sorted_nodes, start=1): atom = molecule.nodes[original_idx] correspondence[original_idx] = idx new_atom = copy.copy(atom) # The charge and the mass can be blank and read from the [atomtypes] # section of the ITP file. new_atom['charge'] = new_atom.get('charge', '') new_atom['mass'] = new_atom.get('mass', '') outfile.write('{idx:>{max_length[idx]}} ' '{atype:<{max_length[atype]}} ' '{resid:>{max_length[resid]}} ' '{resname:<{max_length[resname]}} ' '{atomname:<{max_length[atomname]}} ' '{charge_group:>{max_length[charge_group]}} ' '{charge:>{max_length[charge]}} ' '{mass:>{max_length[mass]}}\n' .format(idx=idx, max_length=max_length, **new_atom)) for line in post_section_lines.get('atoms', []): outfile.write(line + '\n') outfile.write('\n') # Write the interactions conditional_keys = {True: '#ifdef', False: '#ifndef'} for name in molecule.sort_interactions(molecule.interactions): interactions = molecule.interactions[name] # Improper dihedral angles have their own section in the Molecule # object to distinguish them from the proper dihedrals. Yet, they # should be written under the [ dihedrals ] section of the ITP file. if name == 'impropers': name = 'dihedrals' outfile.write('[ {} ]\n'.format(name)) seen_sections.add(name) for line in pre_section_lines.get(name, []): outfile.write(line + '\n') interactions_group_sorted = sorted( interactions, key=_interaction_sorting_key ) interaction_grouped = itertools.groupby( interactions_group_sorted, key=_interaction_sorting_key ) for (conditional, group), interactions_in_group in interaction_grouped: if conditional: conditional_key = conditional_keys[conditional[1]] outfile.write('{} {}\n'.format(conditional_key, conditional[0])) if group: outfile.write('; {}\n'.format(group)) for interaction in interactions_in_group: atoms = ['{atom_idx:>{max_length[idx]}}' .format(atom_idx=correspondence[x], max_length=max_length) for x in interaction.atoms] parameters = ' '.join(str(x) for x in interaction.parameters) comment = '' if 'comment' in interaction.meta: comment = ' ; ' + interaction.meta['comment'] if name == 'virtual_sitesn': to_join = [atoms[0], parameters] + atoms[1:] else: to_join = atoms + [parameters] outfile.write(' '.join(to_join) + comment + '\n') if conditional: outfile.write('#endif\n') for line in post_section_lines.get(name, []): outfile.write(line + '\n') outfile.write('\n') # Some sections may have pre or post lines, but no other content. I that # case, we need to write the sections separately. remaining_sections = set(pre_section_lines) | set(post_section_lines) remaining_sections -= seen_sections for name in remaining_sections: outfile.write('[ {} ]\n'.format(name)) for line in pre_section_lines.get(name, []): outfile.write(line + '\n') for line in post_section_lines.get(name, []): outfile.write(line + '\n') outfile.write('\n')