# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
#
# MDAnalysis --- https://www.mdanalysis.org
# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
# (see the file AUTHORS for the full list of names)
#
# Released under the GNU Public Licence, v2 or any higher version
#
# Please cite your use of MDAnalysis in published work:
#
# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
# doi: 10.25080/majora-629e541a-00e
#
# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
#
"""
AMBER PRMTOP topology parser
============================
Reads an AMBER top file to build the system.
Amber keywords are turned into the following attributes:
+----------------------------+----------------------+
| AMBER flag | MDAnalysis attribute |
+----------------------------+----------------------+
| ATOM_NAME | names |
+----------------------------+----------------------+
| CHARGE | charges |
+----------------------------+----------------------+
| ATOMIC_NUMBER | elements |
+----------------------------+----------------------+
| MASS | masses |
+----------------------------+----------------------+
| BONDS_INC_HYDROGEN | bonds |
| BONDS_WITHOUT_HYDROGEN | |
+----------------------------+----------------------+
| ANGLES_INC_HYDROGEN | angles |
| ANGLES_WITHOUT_HYDROGEN | |
+----------------------------+----------------------+
| DIHEDRALS_INC_HYDROGEN | dihedrals / improper |
| DIHEDRALS_WITHOUT_HYDROGEN | |
+----------------------------+----------------------+
| ATOM_TYPE_INDEX | type_indices |
+----------------------------+----------------------+
| AMBER_ATOM_TYPE | types |
+----------------------------+----------------------+
| RESIDUE_LABEL | resnames |
+----------------------------+----------------------+
| RESIDUE_POINTER | residues |
+----------------------------+----------------------+
TODO:
Add support for Chamber-style topologies
More stringent tests
.. Note::
The Amber charge is converted to electron charges as used in
MDAnalysis and other packages. To get back Amber charges, multiply
by 18.2223.
Chamber-style Amber topologies (i.e. topologies generated via parmed
conversion of a CHARMM topology to an AMBER one) are not currently supported.
Support will likely be added in future MDAnalysis releases.
.. _`PARM parameter/topology file specification`:
http://ambermd.org/formats.html#topology
Classes
-------
.. autoclass:: TOPParser
:members:
:inherited-members:
"""
from __future__ import absolute_import, division
from six.moves import range, zip
from six import raise_from
import numpy as np
import functools
from math import ceil
import itertools
from . import guessers
from .tables import Z2SYMB
from ..lib.util import openany, FORTRANReader
from .base import TopologyReaderBase
from ..core.topology import Topology
from ..core.topologyattrs import (
Atomnames,
Atomtypes,
Atomids,
Charges,
Elements,
Masses,
Resnames,
Resids,
Resnums,
Segids,
AtomAttr,
Bonds,
Angles,
Dihedrals,
Impropers
)
import warnings
import logging
logger = logging.getLogger('MDAnalysis.topology.TOPParser')
class TypeIndices(AtomAttr):
"""Numerical type of each Atom"""
attrname = 'type_indices'
singular = 'type_index'
level = 'atom'
[docs]class TOPParser(TopologyReaderBase):
"""Reads topology information from an AMBER top file.
Reads the following Attributes if in topology:
- Atomnames
- Charges
- Masses
- Elements
- Atomtypes
- Resnames
- Type_indices
- Bonds
- Angles
- Dihedrals (inc. impropers)
Guesses the following attributes:
- Elements (if not included in topology)
The format is defined in `PARM parameter/topology file
specification`_. The reader tries to detect if it is a newer
(AMBER 12?) file format by looking for the flag "ATOMIC_NUMBER".
.. _`PARM parameter/topology file specification`:
http://ambermd.org/formats.html#topology
.. versionchanged:: 0.7.6
parses both amber10 and amber12 formats
.. versionchanged:: 0.19.0
parses bonds, angles, dihedrals, and impropers
.. versionchanged:: 1.0.0
warns users that chamber-style topologies are not current supported
"""
format = ['TOP', 'PRMTOP', 'PARM7']
[docs] def parse(self, **kwargs):
"""Parse Amber PRMTOP topology file *filename*.
Returns
-------
A MDAnalysis Topology object
"""
# Sections that we grab as we parse the file
sections = {
"ATOM_NAME": (1, 20, self.parse_names, "name", 0),
"CHARGE": (1, 5, self.parse_charges, "charge", 0),
"ATOMIC_NUMBER": (1, 10, self.parse_elements, "elements", 0),
"MASS": (1, 5, self.parse_masses, "mass", 0),
"ATOM_TYPE_INDEX": (1, 10, self.parse_type_indices, "type_indices", 0),
"AMBER_ATOM_TYPE": (1, 20, self.parse_types, "types", 0),
"RESIDUE_LABEL": (1, 20, self.parse_resnames, "resname", 11),
"RESIDUE_POINTER": (1, 10, self.parse_residx, "respoint", 11),
"BONDS_INC_HYDROGEN": (3, 10, self.parse_bonded, "bondh", 2),
"BONDS_WITHOUT_HYDROGEN": (3, 10, self.parse_bonded, "bonda", 3),
"ANGLES_INC_HYDROGEN": (4, 10, self.parse_bonded, "angh", 4),
"ANGLES_WITHOUT_HYDROGEN": (4, 10, self.parse_bonded, "anga", 5),
"DIHEDRALS_INC_HYDROGEN": (5, 10, self.parse_bonded, "dihh", 6),
"DIHEDRALS_WITHOUT_HYDROGEN": (5, 10, self.parse_bonded, "diha", 7)
}
attrs = {} # empty dict for attrs that we'll fill
# Open and check top validity
# Reading header info POINTERS
with openany(self.filename) as self.topfile:
header = next(self.topfile)
if not header.startswith("%VE"):
raise ValueError(
"{0} is not a valid TOP file. %VE Missing in header"
"".format(self.filename))
title = next(self.topfile).split()
if not (title[1] == "TITLE"):
# Raise a separate warning if Chamber-style TOP is detected
if title[1] == "CTITLE":
emsg = ("{0} is detected as a Chamber-style TOP file. "
"At this time MDAnalysis does not support such "
"topologies".format(self.filename))
else:
emsg = ("{0} is not a valid TOP file. "
"'TITLE' missing in header".format(self.filename))
raise ValueError(emsg)
while not header.startswith('%FLAG POINTERS'):
header = next(self.topfile)
next(self.topfile)
topremarks = [next(self.topfile).strip() for i in range(4)]
sys_info = [int(k) for i in topremarks for k in i.split()]
header = next(self.topfile)
# grab the next section title
next_section = header.split("%FLAG")[1].strip()
while next_section is not None:
try:
(num_per_record, per_line,
func, name, sect_num) = sections[next_section]
except KeyError:
def next_getter():
return self.skipper()
else:
num = sys_info[sect_num] * num_per_record
numlines = (num // per_line)
if num % per_line != 0:
numlines += 1
attrs[name] = func(num_per_record, numlines)
def next_getter():
return next(self.topfile)
try:
line = next_getter()
# Capture case where section is empty w/ 1 empty line
if numlines == 0 and not line.strip():
line = next_getter()
except StopIteration:
next_section = None
else:
try:
next_section = line.split("%FLAG")[1].strip()
except IndexError:
raise_from(
IndexError((
"%FLAG section not found, formatting error "
"for PARM7 file {0} ").format(self.filename)),
None)
# strip out a few values to play with them
n_atoms = len(attrs['name'])
resptrs = attrs.pop('respoint')
resptrs.append(n_atoms)
residx = np.zeros(n_atoms, dtype=np.int32)
for i, (x, y) in enumerate(zip(resptrs[:-1], resptrs[1:])):
residx[x:y] = i
n_res = len(attrs['resname'])
# Deal with recreating bonds and angle records here
attrs['bonds'] = Bonds([i for i in itertools.chain(
attrs.pop('bonda'), attrs.pop('bondh'))])
attrs['angles'] = Angles([i for i in itertools.chain(
attrs.pop('anga'), attrs.pop('angh'))])
attrs['dihedrals'], attrs['impropers'] = self.parse_dihedrals(
attrs.pop('diha'), attrs.pop('dihh'))
# Guess elements if not in topology or if any dummy atoms are found
if not 'elements' in attrs:
msg = ("ATOMIC_NUMBER record not found, guessing atom elements "
"based on their atom types")
logger.warning(msg)
warnings.warn(msg)
attrs['elements'] = Elements(
guessers.guess_types(attrs['types'].values),
guessed=True)
elif np.any(attrs['elements'].values == "DUMMY"):
# This approach assumes np.any is much faster than np.where
attrs['elements']._guessed = True
for dummy_idx in np.where(attrs['elements'].values == 'DUMMY')[0]:
atom_type = attrs['types'].values[dummy_idx]
atom_ele = guessers.guess_atom_element(atom_type)
msg = ("Unknown ATOMIC_NUMBER value found, guessing atom "
"element from type: {0} assigned to {1}".format(
atom_type, atom_ele))
logger.warning(msg)
warnings.warn(msg)
attrs['elements'].values[dummy_idx] = atom_ele
# atom ids are mandatory
attrs['atomids'] = Atomids(np.arange(n_atoms) + 1)
attrs['resids'] = Resids(np.arange(n_res) + 1)
attrs['resnums'] = Resnums(np.arange(n_res) + 1)
attrs['segids'] = Segids(np.array(['SYSTEM'], dtype=object))
top = Topology(n_atoms, n_res, 1,
attrs=list(attrs.values()),
atom_resindex=residx,
residue_segindex=None)
return top
[docs] def skipper(self):
"""TOPParser :class: helper function, skips lines of input parm7 file
until we find the next %FLAG entry and return that
Returns
-------
line : string
String containing the current line of the parm7 file
"""
line = next(self.topfile)
while not line.startswith("%FLAG"):
line = next(self.topfile)
return line
[docs] def parse_names(self, num_per_record, numlines):
"""Extracts atoms names from parm7 file
Parameters
----------
num_per_record : int
The number of entries for each record in the section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
attr : :class:`Atomnames`
A :class:`Atomnames` instance containing the names of each atom as
defined in the parm7 file
"""
vals = self.parsesection_mapper(numlines, lambda x: x)
attr = Atomnames(np.array(vals, dtype=object))
return attr
[docs] def parse_resnames(self, num_per_record, numlines):
"""Extracts the names of each residue
Parameters
----------
num_per_record : int
The number of entries for each recrod in section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
attr : :class:`Resnames`
A :class:`Resnames` instance containing the names of each residue
as defined in the parm7 file
"""
vals = self.parsesection_mapper(numlines, lambda x: x)
attr = Resnames(np.array(vals, dtype=object))
return attr
[docs] def parse_charges(self, num_per_record, numlines):
"""Extracts the partial charges for each atom
Parameters
----------
num_per_record : int
The number of entries for each record in section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
attr : :class:`Charges`
A :class:`Charges` instance containing the partial charges of each
atom as defined in the parm7 file
"""
vals = self.parsesection_mapper(numlines, lambda x: float(x))
charges = np.array(vals, dtype=np.float32)
charges /= 18.2223 # to electron charge units
attr = Charges(charges)
return attr
[docs] def parse_masses(self, num_per_record, numlines):
"""Extracts the mass of each atom
Parameters
----------
num_per_record : int
The number of entries for each record in section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
attr : :class:`Masses`
A :class:`Masses` instance containing the mass of each atom as
defined in the parm7 file
"""
vals = self.parsesection_mapper(numlines, lambda x: float(x))
attr = Masses(vals)
return attr
[docs] def parse_elements(self, num_per_record, numlines):
"""Extracts the atomic numbers of each atom and converts to element type
Parameters
----------
num_per_record : int
The number of entries for each record in section(unused input)
numlines : int
The number of lines to be pasred in current section
Returns
-------
attr : :class:`Elements`
A :class:`Elements` instance containing the element of each atom
as defined in the parm7 file
Note
----
If the record contains atomic numbers <= 0, these will be treated as
dummy elements and an attempt will be made to guess the element based
on atom type. See issue #2306 for more details.
"""
vals = self.parsesection_mapper(
numlines,
lambda x: Z2SYMB[int(x)] if int(x) > 0 else "DUMMY")
attr = Elements(np.array(vals, dtype=object))
return attr
[docs] def parse_types(self, num_per_record, numlines):
"""Extracts the force field atom types of each atom
Parameters
----------
num_per_record : int
The number of entries for each record in section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
attr : :class:`Atomtypes`
A :class:`Atomtypes` instance containing the atom types for each
atom as defined in the parm7 file
"""
vals = self.parsesection_mapper(numlines, lambda x: x)
attr = Atomtypes(np.array(vals, dtype=object))
return attr
[docs] def parse_type_indices(self, num_per_record, numlines):
"""Extracts the index of atom types of the each atom involved in Lennard
Jones (6-12) interactions.
Parameters
----------
num_per_record : int
The number of entries for each record in section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
attr :class:`TypeIndices`
A :class:`TypeIndices` instance containing the LJ 6-12 atom type
index for each atom
"""
vals = self.parsesection_mapper(numlines, lambda x: int(x))
attr = TypeIndices(np.array(vals, dtype=np.int32))
return attr
[docs] def parse_residx(self, num_per_record, numlines):
"""Extracts the residue pointers for each atom
Parameters
----------
num_per_record : int
The number of entries for each record in section (unused input)
numlines : int
The number of lines to be parsed in current section
Returns
-------
vals : list of int
A list of zero-formatted residue pointers for each atom
"""
vals = self.parsesection_mapper(numlines, lambda x: int(x) - 1)
return vals
[docs] def parse_chunks(self, data, chunksize):
"""Helper function to parse AMBER PRMTOP bonds/angles.
Parameters
----------
data : list of int
Input list of the parm7 bond/angle section, zero-indexed
num_per_record : int
The number of entries for each record in the input list
Returns
-------
vals : list of int tuples
A list of tuples containing the atoms involved in a given bonded
interaction
Note
----
In the parm7 format this information is structured in the following
format: [ atoms 1:n, internal index ]
Where 1:n represent the ids of the n atoms involved in the bond/angle
and the internal index links to a given set of FF parameters.
Therefore, to extract the required information, we split out the list
into chunks of size num_per_record, and only extract the atom ids.
"""
vals = [tuple(data[x:x+chunksize-1])
for x in range(0, len(data), chunksize)]
return vals
[docs] def parse_bonded(self, num_per_record, numlines):
"""Extracts bond information from PARM7 format files
Parameters
----------
num_per_record : int
The number of entries for each record in section
numlines : int
The number of lines to be parsed for this section
Note
----
For the bond/angle sections of parm7 files, the atom numbers are set to
coordinate array index values. As detailed in
http://ambermd.org/formats.html to recover the actual atom number, one
should divide the values by 3 and add 1. Here, since we want to satisfy
zero-indexing, we only divide by 3.
"""
fields = self.parsesection_mapper(numlines, lambda x: int(x) // 3)
section = self.parse_chunks(fields, num_per_record)
return section
[docs] def parsesection_mapper(self, numlines, mapper):
"""Parses FORTRAN formatted section, and returns a list of all entries
in each line
Parameters
----------
numlines : int
The number of lines to be parsed in this section
mapper : lambda operator
Operator to format entries in current section
Returns
-------
section : list
A list of all entries in a given parm7 section
"""
section = []
y = next(self.topfile).strip("%FORMAT(")
y.strip(")")
x = FORTRANReader(y)
for i in range(numlines):
l = next(self.topfile)
for j in range(len(x.entries)):
val = l[x.entries[j].start:x.entries[j].stop].strip()
if val:
section.append(mapper(val))
return section
[docs] def parse_dihedrals(self, diha, dihh):
"""Combines hydrogen and non-hydrogen containing AMBER dihedral lists
and extracts sublists for conventional dihedrals and improper angles
Parameters
----------
diha : list of tuples
The atom ids of dihedrals not involving hydrogens
dihh : list of tuples
The atom ids of dihedrals involving hydrogens
Returns
-------
dihedrals : :class:`Dihedrals`
A :class:`Dihedrals` instance containing a list of all unique
dihedrals as defined by the parm7 file
impropers : :class:`Impropers`
A :class:`Impropers` instance containing a list of all unique
improper dihedrals as defined by the parm7 file
Note
----
As detailed in http://ambermd.org/formats.html, the dihedral sections
of parm7 files contain information about both conventional dihedrals
and impropers. The following must be accounted for:
1) If the fourth atom in a dihedral entry is given a negative value,
this indicates that it is an improper.
2) If the third atom in a dihedral entry is given a negative value,
this indicates that it 1-4 NB interactions are ignored for this
dihedrals. This could be due to the dihedral within a ring, or if it is
part of a multi-term dihedral definition or if it is an improper.
"""
improp = []
dihed = []
for i in itertools.chain(diha, dihh):
if i[3] < 0:
improp.append(i[:2]+(abs(i[2]),)+(abs(i[3]),))
elif i[2] < 0:
vals = i[:2] + (abs(i[2]),) + i[3:]
dihed.append(vals)
else:
dihed.append(i)
dihed = sorted(set(dihed))
dihedrals = Dihedrals(dihed)
impropers = Impropers(improp)
return dihedrals, impropers