# -*- 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 Lesser GNU Public Licence, v2.1 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. As of version 2.0.0, elements are no longer guessed if ATOMIC_NUMBER records are missing. In those scenarios, if elements are necessary, users will have to invoke the element guessers after parsing the topology file. Please see :mod:`MDAnalysis.guessers` for more details. .. _`PARM parameter/topology file specification`: https://ambermd.org/FileFormats.php#topo.cntrl Classes ------- .. autoclass:: TOPParser :members: :inherited-members: """ import numpy as np import itertools from ..guesser.tables import Z2SYMB from ..lib.util import openany, FORTRANReader from .base import TopologyReaderBase, change_squash from ..core.topology import Topology from ..core.topologyattrs import ( Atomnames, Atomtypes, Atomids, Charges, Elements, Masses, Resnames, Resids, Resnums, Segids, ChainIDs, 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" 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) - ChainIDs (from %RESIDUE_CHAINID) - Segids (from %RESIDUE_CHAINID) 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`: https://ambermd.org/FileFormats.php#topo.cntrl Additionally, the RESIDUE_CHAINID non-standard flag is supported. This can be added with the `addPDB`_ command from parmed: .. _`addPDB`: https://parmed.github.io/ParmEd/html/parmed.html#addpdb Notes ----- Elements are obtained from the atomic numbers (if present). If a given input atomic number does not belong to an element (usually either -1 or 0), the element will be assigned an empty record. .. 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 currently supported .. versionchanged:: 2.0.0 no longer guesses elements if missing .. versionchanged:: 2.7.0 gets Segments and chainIDs from flag RESIDUE_CHAINID, when present """ format = ["TOP", "PRMTOP", "PARM7"] 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, ), "RESIDUE_CHAINID": (1, 20, self.parse_chainids, "segids", 11), } 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: errmsg = ( f"%FLAG section not found, formatting error " f"for PARM7 file {self.filename} " ) raise IndexError(errmsg) from 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") ) # Warn user if elements not in topology if "elements" not in attrs: msg = ( "ATOMIC_NUMBER record not found, elements attribute will " "not be populated. If needed these can be guessed using " "universe.guess_TopologyAttrs(to_guess=['elements'])." ) logger.warning(msg) warnings.warn(msg) elif np.any(attrs["elements"].values == ""): # only send out one warning that some elements are unknown msg = ( "Unknown ATOMIC_NUMBER value found for some atoms, these " "have been given an empty element record. If needed these " "can be guessed using " "universe.guess_TopologyAttrs(to_guess=['elements'])." ) logger.warning(msg) warnings.warn(msg) # 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) # Amber's 'RESIDUE_CHAINID' is a by-residue attribute, turn it into # a by-atom attribute when present. See PR #4007. if "segids" in attrs and len(attrs["segids"]) == n_res: segidx, (segids,) = change_squash( (attrs["segids"],), (attrs["segids"],) ) chainids = [attrs["segids"][r] for r in residx] attrs["segids"] = Segids(segids) attrs["ChainIDs"] = ChainIDs(chainids) n_segs = len(segids) else: if "segids" in attrs: msg = ( f"Number of residues ({n_res}) does not match number of " f"%RESIDUE_CHAINID ({len(attrs['segids'])}). Skipping section." ) logger.warning(msg) warnings.warn(msg) attrs["segids"] = Segids(np.array(["SYSTEM"], dtype=object)) segidx = None n_segs = 1 top = Topology( n_atoms, n_res, n_segs, attrs=list(attrs.values()), atom_resindex=residx, residue_segindex=segidx, ) return top 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 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 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 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 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 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 unknown atomic numbers (e.g. <= 0), these will be treated as unknown elements and assigned an empty string value. See issues #2306 and #2651 for more details. .. versionchanged:: 2.0.0 Unrecognised elements will now return a empty string. The parser will no longer attempt to guess the element by default. """ vals = self.parsesection_mapper( numlines, lambda x: Z2SYMB[int(x)] if int(x) > 0 else "" ) attr = Elements(np.array(vals, dtype=object)) return attr 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 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 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 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 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 `the specification`_, 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. .. _`the specification`: https://ambermd.org/FileFormats.php#topo.cntrl """ fields = self.parsesection_mapper(numlines, lambda x: int(x) // 3) section = self.parse_chunks(fields, num_per_record) return section 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 = [] def get_fmt(file): """Skips '%COMMENT' lines until it gets the FORMAT specification for the section.""" line = next(file) if line[:7] == "%FORMAT": return line[8:].split(")")[0] elif line[:8] == "%COMMENT": return get_fmt(file) else: raise ValueError( "Invalid header line. Does not begin with either %FLAG, %FORMAT " f"nor %COMMENT:\n{line}" ) # There may be %COMMENT lines between %FLAG and %FORMAT statements. Skip them. fmt = get_fmt(self.topfile) x = FORTRANReader(fmt) 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 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 `the specification`_, 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. .. _`the specification`: https://ambermd.org/FileFormats.php#topo.cntrl """ 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 def parse_chainids(self, num_per_record: int, numlines: int): """Extracts the chainID of each residue 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 : numpy array A numpy array containing the chainID of each residue as defined in the parm7 file """ vals = self.parsesection_mapper(numlines, lambda x: x) attr = np.array(vals) return attr