Source code for MDAnalysis.coordinates.PDBQT
# -*- 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
#
"""
PDBQT structure files in MDAnalysis --- :mod:`MDAnalysis.coordinates.PDBQT`
===========================================================================
MDAnalysis reads coordinates from PDBQT_ files and additional optional
data such as B-factors, partial charge and AutoDock_ atom types. It
is also possible to substitute a PDBQT file for a PSF file in order to
define the list of atoms (but no connectivity information will be
available in this case).
.. _PDBQT:
http://autodock.scripps.edu/faqs-help/faq/what-is-the-format-of-a-pdbqt-file
.. _AutoDock:
http://autodock.scripps.edu/
"""
import os
import errno
import itertools
import numpy as np
import warnings
from ..lib import util
from . import base
[docs]class PDBQTReader(base.SingleFrameReaderBase):
"""PDBQTReader that reads a PDBQT-formatted file, no frills.
Records read:
- CRYST1 for unitcell A,B,C, alpha,beta,gamm
- ATOM. HETATM for x,y,z
Original `PDB format documentation`_ with `AutoDOCK extensions`_
.. _PDB format documentation:
http://www.wwpdb.org/documentation/file-format-content/format32/v3.2.html
.. _AutoDOCK extensions:
http://autodock.scripps.edu/faqs-help/faq/what-is-the-format-of-a-pdbqt-file
============= ============ =========== =============================================
COLUMNS DATA TYPE FIELD DEFINITION
============= ============ =========== =============================================
1 - 6 Record name "CRYST1"
7 - 15 Real(9.3) a a (Angstroms).
16 - 24 Real(9.3) b b (Angstroms).
25 - 33 Real(9.3) c c (Angstroms).
34 - 40 Real(7.2) alpha alpha (degrees).
41 - 47 Real(7.2) beta beta (degrees).
48 - 54 Real(7.2) gamma gamma (degrees).
1 - 6 Record name "ATOM "
7 - 11 Integer serial Atom serial number.
13 - 16 Atom name Atom name.
17 Character altLoc Alternate location indicator. IGNORED
18 - 21 Residue name resName Residue name.
22 Character chainID Chain identifier.
23 - 26 Integer resSeq Residue sequence number.
27 AChar iCode Code for insertion of residues. IGNORED
31 - 38 Real(8.3) x Orthogonal coordinates for X in Angstroms.
39 - 46 Real(8.3) y Orthogonal coordinates for Y in Angstroms.
47 - 54 Real(8.3) z Orthogonal coordinates for Z in Angstroms.
55 - 60 Real(6.2) occupancy Occupancy.
61 - 66 Real(6.2) tempFactor Temperature factor.
67 - 76 Real(10.4) partialChrg Gasteiger PEOE partial charge *q*.
79 - 80 LString(2) atomType AutoDOCK atom type *t*.
============= ============ =========== =============================================
We ignore torsion notation and just pull the partial charge and atom type columns::
COMPND NSC7810
REMARK 3 active torsions:
REMARK status: ('A' for Active; 'I' for Inactive)
REMARK 1 A between atoms: A7_7 and C22_23
REMARK 2 A between atoms: A9_9 and A11_11
REMARK 3 A between atoms: A17_17 and C21_21
ROOT
123456789.123456789.123456789.123456789.123456789.123456789.123456789.123456789. (column reference)
ATOM 1 A1 INH I 1.054 3.021 1.101 0.00 0.00 0.002 A
ATOM 2 A2 INH I 1.150 1.704 0.764 0.00 0.00 0.012 A
ATOM 3 A3 INH I -0.006 0.975 0.431 0.00 0.00 -0.024 A
ATOM 4 A4 INH I 0.070 -0.385 0.081 0.00 0.00 0.012 A
ATOM 5 A5 INH I -1.062 -1.073 -0.238 0.00 0.00 0.002 A
ATOM 6 A6 INH I -2.306 -0.456 -0.226 0.00 0.00 0.019 A
ATOM 7 A7 INH I -2.426 0.885 0.114 0.00 0.00 0.052 A
ATOM 8 A8 INH I -1.265 1.621 0.449 0.00 0.00 0.002 A
ATOM 9 A9 INH I -1.339 2.986 0.801 0.00 0.00 -0.013 A
ATOM 10 A10 INH I -0.176 3.667 1.128 0.00 0.00 0.013 A
ENDROOT
BRANCH 9 11
ATOM 11 A11 INH I -2.644 3.682 0.827 0.00 0.00 -0.013 A
ATOM 12 A16 INH I -3.007 4.557 -0.220 0.00 0.00 0.002 A
ATOM 13 A12 INH I -3.522 3.485 1.882 0.00 0.00 0.013 A
ATOM 14 A15 INH I -4.262 5.209 -0.177 0.00 0.00 -0.024 A
ATOM 15 A17 INH I -2.144 4.784 -1.319 0.00 0.00 0.052 A
ATOM 16 A14 INH I -5.122 4.981 0.910 0.00 0.00 0.012 A
ATOM 17 A20 INH I -4.627 6.077 -1.222 0.00 0.00 0.012 A
ATOM 18 A13 INH I -4.749 4.135 1.912 0.00 0.00 0.002 A
ATOM 19 A19 INH I -3.777 6.285 -2.267 0.00 0.00 0.002 A
ATOM 20 A18 INH I -2.543 5.650 -2.328 0.00 0.00 0.019 A
BRANCH 15 21
ATOM 21 C21 INH I -0.834 4.113 -1.388 0.00 0.00 0.210 C
ATOM 22 O1 INH I -0.774 2.915 -1.581 0.00 0.00 -0.644 OA
ATOM 23 O3 INH I 0.298 4.828 -1.237 0.00 0.00 -0.644 OA
ENDBRANCH 15 21
ENDBRANCH 9 11
BRANCH 7 24
ATOM 24 C22 INH I -3.749 1.535 0.125 0.00 0.00 0.210 C
ATOM 25 O2 INH I -4.019 2.378 -0.708 0.00 0.00 -0.644 OA
ATOM 26 O4 INH I -4.659 1.196 1.059 0.00 0.00 -0.644 OA
ENDBRANCH 7 24
TORSDOF 3
123456789.123456789.123456789.123456789.123456789.123456789.123456789.123456789. (column reference)
.. versionchanged:: 0.11.0
Frames now 0-based instead of 1-based
"""
format = 'PDBQT'
units = {'time': None, 'length': 'Angstrom'}
def _read_first_frame(self):
coords = []
unitcell = np.zeros(6, dtype=np.float32)
with util.openany(self.filename) as pdbfile:
for line in pdbfile:
# Should only break at the 'END' of a model definition
# and prevent premature exit for a torsion termination
# , eg, ENDBRANCH
if line.startswith('END\n'):
break
if line.startswith('CRYST1'):
# lengths
x, y, z = np.float32((line[6:15], line[15:24], line[24:33]))
# angles
A, B, G = np.float32((line[33:40], line[40:47], line[47:54]))
unitcell[:] = x, y, z, A, B, G
if line.startswith(('ATOM', 'HETATM')):
# convert all entries at the end once for optimal speed
coords.append([line[30:38], line[38:46], line[46:54]])
self.n_atoms = len(coords)
self.ts = self._Timestep.from_coordinates(
coords,
**self._ts_kwargs)
self.ts.dimensions = unitcell
self.ts.frame = 0 # 0-based frame number
if self.convert_units:
# in-place !
self.convert_pos_from_native(self.ts._pos)
if self.ts.dimensions is not None:
self.convert_pos_from_native(self.ts.dimensions[:3])
[docs] def Writer(self, filename, **kwargs):
"""Returns a permissive (simple) PDBQTWriter for *filename*.
Parameters
----------
filename : str
filename of the output PDBQT file
Returns
-------
:class:`PDBQTWriter`
"""
return PDBQTWriter(filename, **kwargs)
[docs]class PDBQTWriter(base.WriterBase):
"""PDBQT writer that implements a subset of the PDB_ 3.2 standard and the PDBQT_ spec.
.. _PDB: http://www.wwpdb.org/documentation/file-format-content/format32/v3.2.html
.. _PDBQT: http://autodock.scripps.edu/faqs-help/faq/what-is-the-format-of-a-pdbqt-file
"""
fmt = {
'ATOM': ("ATOM {serial:5d} {name:<4.4s} {resName:<4.4s}"
"{chainID:1.1s}{resSeq:4d}{iCode:1.1s}"
" {pos[0]:8.3f}{pos[1]:8.3f}{pos[2]:8.3f}{occupancy:6.2f}"
"{tempFactor:6.2f} {charge:< 1.3f} {element:<2.2s}\n"),
'REMARK': "REMARK {0}\n",
'TITLE': "TITLE {0}\n",
'CRYST1': ("CRYST1{box[0]:9.3f}{box[1]:9.3f}{box[2]:9.3f}"
"{ang[0]:7.2f}{ang[1]:7.2f}{ang[2]:7.2f} "
"{spacegroup:<11s}{zvalue:4d}\n"),
}
format = 'PDBQT'
units = {'time': None, 'length': 'Angstrom'}
pdb_coor_limits = {"min": -999.9995, "max": 9999.9995}
def __init__(self, filename, **kwargs):
self.filename = util.filename(filename, ext='pdbqt')
self.pdb = util.anyopen(self.filename, 'w')
[docs] def write(self, selection, frame=None):
"""Write selection at current trajectory frame to file.
Parameters
----------
selection : AtomGroup
The selection to be written
frame : int (optional)
optionally move to frame index `frame` before writing; the default
is to write the current trajectory frame
Note
----
The first letter of the
:attr:`~MDAnalysis.core.groups.Atom.segid` is used as the PDB
chainID.
.. versionchanged:: 0.11.0
Frames now 0-based instead of 1-based
"""
try:
u = selection.universe
except AttributeError:
errmsg = "Input obj is neither an AtomGroup or Universe"
raise TypeError(errmsg) from None
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 PDB (?)
atoms = selection.atoms # make sure to use atoms (Issue 46)
coor = atoms.positions # can write from selection == Universe (Issue 49)
# Check attributes
attrs = {}
missing_topology = []
for attr, dflt in (
('altLocs', ' '),
('charges', 0.0),
('icodes', ' '),
('names', 'X'),
('occupancies', 1.0),
('resids', 1),
('resnames', 'UNK'),
('tempfactors', 0.0),
('types', ' '),
):
try:
attrs[attr] = getattr(atoms, attr)
except AttributeError:
attrs[attr] = itertools.cycle((dflt,))
missing_topology.append(attr)
# Order of preference: chainids -> segids -> blank string
try:
attrs['chainids'] = atoms.chainids
except AttributeError:
try:
attrs['chainids'] = atoms.segids
except AttributeError:
attrs['chainids'] = itertools.cycle((' ',))
missing_topology.append('chainids')
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)))
# check if any coordinates are illegal (coordinates are already
# in Angstroem per package default)
if not self.has_valid_coordinates(self.pdb_coor_limits, coor):
self.close()
try:
os.remove(self.filename)
except OSError as err:
if err.errno == errno.ENOENT:
pass
raise ValueError(
"PDB files must have coordinate values between {0:.3f}"
" and {1:.3f} Angstroem: No file was written."
"".format(self.pdb_coor_limits["min"],
self.pdb_coor_limits["max"]))
# Write title record
# http://www.wwpdb.org/documentation/file-format-content/format32/sect2.html
line = "FRAME " + str(frame) + " FROM " + str(u.trajectory.filename)
self.pdb.write(self.fmt['TITLE'].format(line))
# Write CRYST1 record
# http://www.wwpdb.org/documentation/file-format-content/format32/sect8.html
box = self.convert_dimensions_to_unitcell(u.trajectory.ts)
self.pdb.write(self.fmt['CRYST1'].format(box=box[:3], ang=box[3:],
spacegroup='P 1', zvalue=1))
# Write atom records
# http://www.wwpdb.org/documentation/file-format-content/format32/sect9.html
for serial, (pos, name, resname, chainid, resid, icode,
occupancy, tempfactor, charge, element) in enumerate(
zip(coor, attrs['names'], attrs['resnames'], attrs['chainids'],
attrs['resids'], attrs['icodes'], attrs['occupancies'],
attrs['tempfactors'], attrs['charges'], attrs['types']),
start=1):
serial = util.ltruncate_int(serial, 5) # check for overflow here?
resid = util.ltruncate_int(resid, 4)
name = name[:4]
if len(name) < 4:
name = " " + name # customary to start in column 14
chainid = chainid.strip()[-1:] # take the last character
self.pdb.write(self.fmt['ATOM'].format(
serial=serial,
name=name,
resName=resname,
chainID=chainid,
resSeq=resid,
iCode=icode,
pos=pos,
occupancy=occupancy,
tempFactor=tempfactor,
charge=charge,
element=element,
))
self.close()