Source code for MDAnalysis.coordinates.PQR

# -*- 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
#


"""
PQR file format --- :mod:`MDAnalysis.coordinates.PQR`
=====================================================

Read atoms with charges from a PQR_ file (as written by PDB2PQR_). The
following is adopted from the description of the PQR_ format as used by APBS_:

*MDAnalysis* reads very loosely-formatted PQR files: all fields are
**whitespace-delimited** rather than the strict column formatting mandated
by the PDB_ format. This more liberal formatting allows coordinates
which are larger/smaller than ±999 Å.

MDAnalysis reads data on a per-line basis from PQR files using the following format::

   recordName serial atomName residueName chainID residueNumber X Y Z charge radius

If this fails it is assumed that the *chainID* was omitted and the shorter
format is read::

   recordName serial atomName residueName residueNumber X Y Z charge radius

Anything else will raise a :exc:`ValueError`.

The whitespace is the most important feature of this format: fields
*must* be separated by at least one space or tab character. The fields
are:

*recordName*
    A string which specifies the type of PQR entry and should either be ATOM or
    HETATM.
*serial*
    An integer which provides the atom index (but note that MDAnalysis renumbers
    atoms so one cannot rely on the *serial*)
*atomName*
    A string which provides the atom name.
*residueName*
    A string which provides the residue name.
*chainID*
    An optional string which provides the chain ID of the atom.
*residueNumber*
    An integer which provides the residue index.
*X Y Z*
    Three floats which provide the atomic coordiantes.
*charge*
    A float which provides the atomic charge (in electrons).
*radius*
    A float which provides the atomic radius (in Å).

Clearly, this format can deviate wildly from PDB_ due to the use of whitespaces
rather than specific column widths and alignments. This deviation can be
particularly significant when large coordinate values are used.

Output should look like this (although the only real requirement is
*whitespace* separation between *all* entries). The chainID is optional
and can be omitted::

  ATOM      1  N    MET     1     -11.921   26.307   10.410 -0.3000 1.8500
  ATOM     36  NH1  ARG     2      -6.545   25.499    3.854 -0.8000 1.8500
  ATOM     37 HH11  ARG     2      -6.042   25.480    4.723  0.4600 0.2245


.. Warning:: Fields *must be white-space separated* or data are read
             incorrectly. PDB formatted files are *not* guaranteed to be
             white-space separated so extra care should be taken when quickly
             converting PDB files into PQR files using simple scripts.

For example, PQR files created with PDB2PQR_ and the `--whitespace`
option are guaranteed to conform to the above format::

   pdb2pqr --ff=charmm --whitespace 4ake.pdb 4ake.pqr

.. _PQR:     https://apbs-pdb2pqr.readthedocs.io/en/latest/formats/pqr.html
.. _APBS:    https://apbs-pdb2pqr.readthedocs.io/en/latest/apbs/index.html
.. _PDB2PQR: https://apbs-pdb2pqr.readthedocs.io/en/latest/pdb2pqr/index.html
.. _PDB:     http://www.wwpdb.org/documentation/file-format
"""
from __future__ import absolute_import
from six.moves import zip

import itertools
import numpy as np
import warnings

from ..lib import util
from . import base


[docs]class PQRReader(base.SingleFrameReaderBase): """Read a PQR_ file into MDAnalysis. .. _PQR: https://apbs-pdb2pqr.readthedocs.io/en/latest/formats/pqr.html .. versionchanged:: 0.11.0 Frames now 0-based instead of 1-based """ format = 'PQR' units = {'time': None, 'length': 'Angstrom'} # how to slice fields[x:y] to grab coordinates _SLICE_INDICES = { 'ORIGINAL': (-5, -2), 'NO_CHAINID': (-5, -2), 'GROMACS': (-6, -3), } def _read_first_frame(self): from ..topology.PQRParser import PQRParser flavour = None coords = [] unitcell = np.zeros(6, dtype=np.float32) with util.openany(self.filename) as pqrfile: for line in pqrfile: if line.startswith(('ATOM', 'HETATM')): if flavour is None: flavour = PQRParser.guess_flavour(line) x, y = self._SLICE_INDICES[flavour] fields = line.split() # convert all entries at the end once for optimal speed coords.append(fields[x:y]) self.n_atoms = len(coords) self.ts = self._Timestep.from_coordinates( coords, **self._ts_kwargs) self.ts._unitcell[:] = unitcell self.ts.frame = 0 # 0-based frame number if self.convert_units: # in-place ! self.convert_pos_from_native(self.ts._pos) self.convert_pos_from_native(self.ts._unitcell[:3])
[docs] def Writer(self, filename, **kwargs): """Returns a PQRWriter for *filename*. Parameters ---------- filename : str filename of the output PQR file Returns ------- :class:`PQRWriter` """ return PQRWriter(filename, **kwargs)
[docs]class PQRWriter(base.WriterBase): """Write a single coordinate frame in whitespace-separated PQR format. Charges ("Q") are taken from the :attr:`MDAnalysis.core.groups.Atom.charge` attribute while radii are obtaine from the :attr:`MDAnalysis.core.groups.Atom.radius` attribute. * If the segid is 'SYSTEM' then it will be set to the empty string. Otherwise the first letter will be used as the chain ID. * The serial number always starts at 1 and increments sequentially for the atoms. The output format is similar to ``pdb2pqr --whitespace``. .. versionadded:: 0.9.0 """ format = 'PQR' units = {'time': None, 'length': 'Angstrom'} # serial, atomName, residueName, chainID, residueNumber, XYZ, charge, radius fmt_ATOM = ("ATOM {serial:6d} {name:<4} {resname:<3} {chainid:1.1}" " {resid:4d} {pos[0]:-8.3f} {pos[1]:-8.3f}" " {pos[2]:-8.3f} {charge:-7.4f} {radius:6.4f}\n") fmt_remark = "REMARK {0} {1}\n" def __init__(self, filename, convert_units=True, **kwargs): """Set up a PQRWriter with full whitespace separation. Parameters ---------- filename : str output filename convert_units: bool (optional) units are converted to the MDAnalysis base format; [``True``] remarks : str (optional) remark lines (list of strings) or single string to be added to the top of the PQR file """ self.filename = util.filename(filename, ext='pqr') self.convert_units = convert_units # convert length and time to base units self.remarks = kwargs.pop('remarks', "PQR file written by MDAnalysis")
[docs] def write(self, selection, frame=None): """Write selection at current trajectory frame to file. Parameters ---------- selection : AtomGroup or Universe MDAnalysis AtomGroup or Universe frame : int (optional) optionally move to frame index `frame`; by default, write the current frame .. versionchanged:: 0.11.0 Frames now 0-based instead of 1-based """ # write() method that complies with the Trajectory API u = selection.universe if frame is not None: u.trajectory[frame] # advance to frame else: try: frame = u.trajectory.ts.frame except AttributeError: frame = 0 # should catch cases when we are analyzing a single frame(?) atoms = selection.atoms # make sure to use atoms (Issue 46) coordinates = atoms.positions # can write from selection == Universe (Issue 49) if self.convert_units: self.convert_pos_to_native(coordinates) # inplace because coordinates is already a copy # Check atom attributes # required: # - name # - resname # - chainid # - resid # - position # - charge # - radius attrs = {} missing_topology = [] for attr, dflt in ( ('names', itertools.cycle(('X',))), ('resnames', itertools.cycle(('UNK',))), ('resids', itertools.cycle((1,))), ('charges', itertools.cycle((0.0,))), ('radii', itertools.cycle((1.0,))), ): try: attrs[attr] = getattr(atoms, attr) except AttributeError: attrs[attr] = dflt missing_topology.append(attr) # chainids require special handling # try chainids, then segids # if neither, use ' ' # if 'SYSTEM', use ' ' try: attrs['chainids'] = atoms.chainids except AttributeError: try: attrs['chainids'] = atoms.segids except AttributeError: pass if not 'chainids' in attrs or all(attrs['chainids'] == 'SYSTEM'): attrs['chainids'] = itertools.cycle((' ',)) if 'charges' in missing_topology: total_charge = 0.0 else: total_charge = atoms.total_charge() if missing_topology: warnings.warn( "Supplied AtomGroup was missing the following attributes: " "{miss}. These will be written with default values. " "".format(miss=', '.join(missing_topology))) with util.openany(self.filename, 'w') as pqrfile: # Header / Remarks # The *remarknumber* is typically 1 but :program:`pdb2pgr` # also uses 6 for the total charge and 5 for warnings. for rem in util.asiterable(self.remarks): pqrfile.write(self.fmt_remark.format(rem, 1)) pqrfile.write(self.fmt_remark.format( "Input: frame {0} of {1}".format(frame, u.trajectory.filename), 5)) pqrfile.write(self.fmt_remark.format( "total charge: {0:+8.4f} e".format(total_charge), 6)) # Atom descriptions and coords for atom_index, (pos, name, resname, chainid, resid, charge, radius) in enumerate(zip( coordinates, attrs['names'], attrs['resnames'], attrs['chainids'], attrs['resids'], attrs['charges'], attrs['radii']), start=1): # pad so that only 4-letter atoms are left-aligned name = " " + name if len(name) < 4 else name pqrfile.write(self.fmt_ATOM.format( serial=atom_index, name=name, resname=resname, chainid=chainid, resid=resid, pos=pos, charge=charge, radius=radius))