diff --git a/.travis.yml b/.travis.yml index 1249fc9a70d..ee68a9d0f5d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,12 +28,14 @@ env: - SETUP_CMD="${PYTEST_FLAGS}" - BUILD_CMD="pip install -e package/ && (cd testsuite/ && python setup.py build)" - CONDA_MIN_DEPENDENCIES="mmtf-python mock six biopython networkx cython matplotlib scipy griddataformats hypothesis gsd codecov" - - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles" + - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 chemfiles mdtraj" - CONDA_CHANNELS='biobuilds conda-forge' - CONDA_CHANNEL_PRIORITY=True - PIP_DEPENDENCIES="duecredit parmed" - NUMPY_VERSION=stable - INSTALL_HOLE="true" + - INSTALL_DSSP="true" + - INSTALL_STRIDE="true" - CYTHON_TRACE_NOGIL=1 matrix: @@ -55,6 +57,8 @@ matrix: BUILD_DOCS=true BUILD_CMD="cd ${TRAVIS_BUILD_DIR}/package && python setup.py build_ext --inplace" INSTALL_HOLE="false" + INSTALL_DSSP="false" + INSTALL_STRIDE="false" PIP_DEPENDENCIES="${PIP_DEPENDENCIES} sphinx==1.8.5 sphinx-sitemap sphinx_rtd_theme" - env: NAME="Lint" @@ -64,6 +68,8 @@ matrix: BUILD_CMD="" CONDA_DEPENDENCIES="" INSTALL_HOLE="false" + INSTALL_DSSP="false" + INSTALL_STRIDE="false" - env: NAME="python 2.7" PYTHON_VERSION=2.7 @@ -97,6 +103,8 @@ matrix: PIP_DEPENDENCIES="" CONDA_DEPENDENCIES=${CONDA_MIN_DEPENDENCIES} INSTALL_HOLE="false" + INSTALL_DSSP="false" + INSTALL_STRIDE="false" CODECOV="true" SETUP_CMD="${PYTEST_FLAGS} --cov=MDAnalysis" allow_failures: @@ -120,6 +128,16 @@ install: HOLE_BINDIR="${HOME}/hole2/exe"; \ export PATH=${PATH}:${HOLE_BINDIR}; \ fi + if [[ $INSTALL_DSSP == 'true' ]]; then \ + bash ./maintainer/install_dssp.sh $TRAVIS_OS_NAME "${HOME}"; \ + DSSP_BINDIR="${HOME}/dssp-3.1.4"; \ + export PATH=${PATH}:${DSSP_BINDIR}; \ + fi + if [[ $INSTALL_STRIDE == 'true' ]]; then \ + bash ./maintainer/install_stride.sh $TRAVIS_OS_NAME "${HOME}"; \ + STRIDE_BINDIR="${HOME}/stride"; \ + export PATH=${PATH}:${STRIDE_BINDIR}; \ + fi - git clone git://github.com/astropy/ci-helpers.git - source ci-helpers/travis/setup_conda.sh - eval $BUILD_CMD diff --git a/maintainer/install_dssp.sh b/maintainer/install_dssp.sh new file mode 100755 index 00000000000..4a9433ca6ad --- /dev/null +++ b/maintainer/install_dssp.sh @@ -0,0 +1,127 @@ +#!/usr/bin/env bash +# +# Written by Lily Wang, 2020 +# +# This script is placed into the Public Domain using the CC0 1.0 Public Domain +# Dedication https://creativecommons.org/publicdomain/zero/1.0/ +# +# +# NOTE: IF YOU RUN THIS SCRIPT YOU MUST BE ABLE TO COMPLY WITH THE DSSP +# LICENSE BELOW. +# +# Boost Software License - Version 1.0 - August 17th, 2003 +# Permission is hereby granted, free of charge, to any person or organization +# obtaining a copy of the software and accompanying documentation covered by +# this license (the "Software") to use, reproduce, display, distribute, execute, +# and transmit the Software, and to prepare derivative works of the Software, +# and to permit third-parties to whom the Software is furnished to do so, all +# subject to the following: + +# The copyright notices in the Software and this entire statement, including +# the above license grant, this restriction and the following disclaimer, +# must be included in all copies of the Software, in whole or in part, and all +# derivative works of the Software, unless such copies or derivative works are +# solely in the form of machine-executable object code generated by a source +# language processor. + +# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +# FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +# SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE FOR +# ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING +# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS +# IN THE SOFTWARE. +# +# +# Install DSSP from https://swift.cmbi.umcn.nl/gv/dssp/ +# +# arguments +# OSNAME : linux | darwin +# PREFIX : "path/to/installdir" +# +# PREFIX can be relative to CWD or an absolute path; it will be created +# and DSSP unpacked under PREFIX (the tar file contains "dssp-xx.xx.xx/...") +# +# + + +set -o errexit -o nounset + +SCRIPT=$(basename $0) + + +DOWNLOAD_URL='https://github.com/cmbi/dssp/archive/3.1.4.tar.gz' +TARFILE='3.1.4.tar.gz' + + +# path to dir with executables in the current HOLE distribution +DSSP_EXE_DIR=dssp-3.1.4 + +function die () { + local msg="$1" err=$2 + echo "[${SCRIPT}] ERROR: $msg [$err]" + exit $err +} + +function is_native_executable () { + local filename="$1" OSNAME="$2" + file "${filename}" | grep -qi ${OFORMAT} + return $? +} + +OSNAME="$1" +case "$OSNAME" in + Linux|linux) + OSNAME=Linux + OFORMAT=Linux + sudo apt-get -y install libboost-all-dev autoconf automake autotools-dev + ;; + OSX|osx|Darwin|darwin) + OSNAME=Darwin + OFORMAT="Mach-O" + ;; + *) + die "OS '${OSNAME}' not supported." 10;; +esac + + +PREFIX="$2" +test -n "$PREFIX" || die "missing argument PREFIX (installation directory)" 10 + +#------------------------------------------------------------ +# start installation +#------------------------------------------------------------ + +mkdir -p "$PREFIX" && cd "$PREFIX" || die "Failed to create and cd to $PREFIX" 1 +if ! test -f ${TARFILE}; then + echo "Downloading ${TARFILE} from ${DOWNLOAD_URL}..." + # fixing curl on travis/anaconda, see http://stackoverflow.com/a/31060428/334357 + export CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt + curl -L ${DOWNLOAD_URL} -o ${TARFILE} || \ + die "Failed to download ${TARFILE} from ${DOWNLOAD_URL}" 1 +fi + +# only install the executables in DSSP_EXE_DIR +tar xvf ${TARFILE} ${DSSP_EXE_DIR} || die "Failed 'tar xvf ${TARFILE} ${DSSP_EXE_DIR}'" 1 + +MDA_DSSP_BINDIR="${PWD}/${DSSP_EXE_DIR}" +DSSP_EXE="${MDA_DSSP_BINDIR}/mkdssp" + +cd $MDA_DSSP_BINDIR + +echo "Configuring" +./autogen.sh +./configure +echo "Making" +make + +test -f "${DSSP_EXE}" || die "mkdssp executable ${DSSP_EXE} not installed" 2 +is_native_executable ${DSSP_EXE} ${OFORMAT} || \ + { file ${DSSP_EXE}; \ + die "${DSSP_EXE} will not run on ${OSNAME} (object format ${OFORMAT})" 3; } + +echo "mkdssp was installed into ${MDA_DSSP_BINDIR}" +echo ">>> export PATH=\${PATH}:${MDA_DSSP_BINDIR}" + +# repeat this line in .travis.yml +export PATH=${PATH}:${MDA_DSSP_BINDIR} diff --git a/maintainer/install_stride.sh b/maintainer/install_stride.sh new file mode 100755 index 00000000000..be955b3ec5d --- /dev/null +++ b/maintainer/install_stride.sh @@ -0,0 +1,117 @@ +#!/usr/bin/env bash +# +# Written by Lily Wang, 2020 +# +# This script is placed into the Public Domain using the CC0 1.0 Public Domain +# Dedication https://creativecommons.org/publicdomain/zero/1.0/ +# +# NOTE: IF YOU RUN THIS SCRIPT YOU MUST BE ABLE TO COMPLY WITH THE DSSP +# LICENSE BELOW. +# +# All rights reserved, whether the whole or part of the program is +# concerned. Permission to use, copy, and modify this software and its +# documentation is granted for academic use, provided that: +# +# i. this copyright notice appears in all copies of the software and +# related documentation; +# +# ii. the reference given below (Frishman and Argos, 1995) must be +# cited in any publication of scientific results based in part or +# completely on the use of the program; +# +# iii. bugs will be reported to the authors. +# +# The use of the software in commercial activities is not allowed +# without a prior written commercial license agreement. +# +# IF YOU USE STRIDE IN PUBLISHED WORK, PLEASE CITE: +# +# Frishman, D and Argos, P. (1995) Knowledge-based secondary structure +# assignment. Proteins: structure, function and genetics, 23, +# 566-579. +# +# The source code is available at http://webclu.bio.wzw.tum.de/stride/ +# +# arguments +# OSNAME : linux | darwin +# PREFIX : "path/to/installdir" +# +# PREFIX can be relative to CWD or an absolute path; it will be created +# and STRIDE unpacked under PREFIX/stride. +# +# + + +set -o errexit -o nounset + +SCRIPT=$(basename $0) + + +DOWNLOAD_URL='http://webclu.bio.wzw.tum.de/stride/stride.tar.gz' +TARFILE='stride.tar.gz' + + +# path to dir with executables in the current HOLE distribution +STRIDE_EXE_DIR=stride + +function die () { + local msg="$1" err=$2 + echo "[${SCRIPT}] ERROR: $msg [$err]" + exit $err +} + +function is_native_executable () { + local filename="$1" OSNAME="$2" + file "${filename}" | grep -qi ${OFORMAT} + return $? +} + +OSNAME="$1" +case "$OSNAME" in + Linux|linux) + OSNAME=Linux + OFORMAT=Linux + ;; + OSX|osx|Darwin|darwin) + OSNAME=Darwin + OFORMAT="Mach-O" + ;; + *) + die "OS '${OSNAME}' not supported." 10;; +esac + + +PREFIX="$2" +test -n "$PREFIX" || die "missing argument PREFIX (installation directory)" 10 +PREFIX="${PREFIX}/${STRIDE_EXE_DIR}" + +#------------------------------------------------------------ +# start installation +#------------------------------------------------------------ + +mkdir -p "$PREFIX" && cd "$PREFIX" || die "Failed to create and cd to $PREFIX" 1 +if ! test -f ${TARFILE}; then + echo "Downloading ${TARFILE} from ${DOWNLOAD_URL}..." + # fixing curl on travis/anaconda, see http://stackoverflow.com/a/31060428/334357 + export CURL_CA_BUNDLE=/etc/ssl/certs/ca-certificates.crt + curl -L ${DOWNLOAD_URL} -o ${TARFILE} || \ + die "Failed to download ${TARFILE} from ${DOWNLOAD_URL}" 1 +fi + +# only install the executables in STRIDE_EXE_DIR +tar xvf ${TARFILE} || die "Failed 'tar xvf ${TARFILE}'" 1 +STRIDE_EXE="${PWD}/stride" + +echo "Making" +make + +test -f "${STRIDE_EXE}" || die "stride executable ${STRIDE_EXE} not installed" 2 +is_native_executable ${STRIDE_EXE} ${OFORMAT} || \ + { file ${STRIDE_EXE}; \ + die "${STRIDE_EXE} will not run on ${OSNAME} (object format ${OFORMAT})" 3; } + +echo "stride was installed into ${PWD}" +echo ">>> export PATH=\${PATH}:${PWD}" + +# repeat this line in .travis.yml +export PATH=${PATH}:${PWD} diff --git a/package/MDAnalysis/analysis/secondary_structure/__init__.py b/package/MDAnalysis/analysis/secondary_structure/__init__.py new file mode 100644 index 00000000000..3a2acaac689 --- /dev/null +++ b/package/MDAnalysis/analysis/secondary_structure/__init__.py @@ -0,0 +1,27 @@ +# -*- 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-2020 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 +# + +from __future__ import absolute_import + +from .dssp import DSSP +from .wrapper import DSSPWrapper, STRIDEWrapper \ No newline at end of file diff --git a/package/MDAnalysis/analysis/secondary_structure/base.py b/package/MDAnalysis/analysis/secondary_structure/base.py new file mode 100644 index 00000000000..c097055edd9 --- /dev/null +++ b/package/MDAnalysis/analysis/secondary_structure/base.py @@ -0,0 +1,243 @@ +# -*- 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-2020 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 +# + +"""Secondary structure analysis --- :mod:`MDAnalysis.analysis.secondary_structure` +================================================================================== + +:Authors: Lily Wang +:Year: 2020 +:Copyright: GNU Public License v3 + +.. versionadded:: 1.0.0 + + +""" + +from __future__ import absolute_import, division + +import numpy as np + +from ...core.topologyattrs import ResidueAttr +from ..base import AnalysisBase + + +class SecondaryStructure(ResidueAttr): + """Single letter code for secondary structure""" + attrname = 'secondary_structures' + singular = 'secondary_structure' + dtype = '`. Please cite + [Kabsch1983]_ , [Touw2015]_ and [McGibbon2015]_ if you use this in + published work. + + + Parameters + ---------- + universe: Universe or AtomGroup + The Universe or AtomGroup to apply the analysis to. As secondary + structure is a residue property, the analysis is applied to every + residue in your chosen atoms. + select: string, optional + The selection string for selecting atoms from ``universe``. The + analysis is applied to the residues in this subset of atoms. + c_name: string, optional + Atom name for amide carbon in protein backbone. If there are + multiple possibilities, separate the names with a space, e.g. + "C C2 Cz". This gets passed into :meth:`Universe.select_atoms`. + Names earlier in the string will get preferenced if multiple atoms + fitting the selection are found in one residue. If multiple atoms + have the same name, the earlier one (i.e. with lower index) is chosen. + o_name: string, optional + Atom name for amide oxygen in protein backbone. Analogous to + ``c_name``. + n_name: string, optional + Atom name for amide nitrogen in protein backbone. Analogous to + ``c_name``. + ca_name: string, optional + Atom name for alpha-carbon in protein residues. Analogous to + ``c_name``. + pro_name: string, optional + Residue name for prolines. Analogous to ``c_name``. + add_topology_attr: bool, optional + Whether to add the most common secondary structure as a topology + attribute ``secondary_structure`` to your residues. + verbose: bool, optional + Turn on more logging and debugging. + + + Attributes + ---------- + residues: :class:`~MDAnalysis.core.groups.ResidueGroup` + The residues to which the analysis is applied. + ss_codes: :class:`numpy.ndarray` of shape (n_frames, n_residues) + Single-letter secondary structure codes + ss_names: :class:`numpy.ndarray` of shape (n_frames, n_residues) + Secondary structure names + ss_simple: :class:`numpy.ndarray` of shape (n_frames, n_residues) + Simplified secondary structure names + ss_counts: dict of {code: :class:`numpy.ndarray` of shape (n_frames,)} + Dictionary that counts number of residues with each secondary + structure code, for each frame + simple_counts: dict of {name: :class:`numpy.ndarray` of shape (n_frames,)} + Dictionary that counts number of residues with each simplified + secondary structure, for each frame + ss_mode: :class:`numpy.ndarray` of shape (n_residues,) + The most common secondary structure for each residue + simple_mode: :class:`numpy.ndarray` of shape (n_residues,) + The most common simple secondary structure for each residue + ss_codes_to_names: dict of {code: name} + Dictionary converting each single-letter code to the full name of + the secondary structure + ss_codes_to_simple: dict of {code: name} + Dictionary converting each single-letter code to simplified + secondary structures + atomgroup: :class:`~MDAnalysis.core.groups.AtomGroup` + An ordered AtomGroup of backbone atoms + nco_indices: :class:`numpy.ndarray` of shape (n_residues, 3) + An array of the atom indices of N, C, O atoms, *relative* to the ``atomgroup``. + An index is -1 if it is not found in the residue. + ca_indices: :class:`numpy.ndarray` of shape (n_residues,) + An array of the atom indices of CA atoms, *relative* to the ``atomgroup``. + An index is -1 if it is not found in the residue. + is_proline: :class:`numpy.ndarray` of shape (n_residues,) + An integer array indicating whether a residue is a proline. 1 is True, + 0 is False. + not_protein: :class:`numpy.ndarray` of shape (n_residues,) + A boolean array indicating whether a residue has all of the N, C, O, CA atoms. + chain_ids: :class:`numpy.ndarray` of shape (n_residues,) + An integer array representing the index of the chainID or segment each residue + belongs to. + + + Example + ------- + + :: + + import MDAnalysis as mda + from MDAnalysis.tests.datafiles import PSF, DCD + from MDAnalysis.analysis import secondary_structure as ss + + u = mda.Universe(PSF, DCD) + dssp = ss.DSSP(u, select='backbone', o_name='O O1', + add_topology_attr=True) + dssp.run() + ss0 = u.residues[0].secondary_structure + print('The first residue has secondary structure {}'.format(ss0)) + + The full secondary structure codes can be accessed at :attr`DSSP.ss_codes`. + + + """ + + def __init__(self, universe, select='backbone', verbose=False, + add_topology_attr=False, c_name='C', o_name='O O1', + n_name='N', ca_name='CA', pro_name='PRO'): + # TODO: implement this on its own w/o mdtraj? + try: + from mdtraj.geometry._geometry import _dssp + except ImportError: + raise ValueError('DSSP requires mdtraj to be installed') + else: + self._mdtraj_dssp = _dssp + + super(DSSP, self).__init__(universe, select=select, verbose=verbose, + add_topology_attr=add_topology_attr) + self.names = {'C': c_name, + 'O': o_name, + 'N': n_name, + 'CA': ca_name, + 'PRO': pro_name} + + self._setup_atomgroups() + + def _setup_atomgroups(self): + """ + Set up atomgroup and indices for MDTraj. + """ + ag = self.residues.atoms + ridx = self.residues.resindices + + # grab n, c, o + ncoca = ag[[]] + ncoca_indices = np.array([[-1]*self.n_residues]*4, dtype=np.int32) + for i, atom in enumerate(['N', 'C', 'O', 'CA']): + name = self.names[atom] + selections = ['name {}'.format(n) for n in name.split()] + atoms = ag.select_atoms(*selections) + + # index of atoms, within atomgroup + residx, idx, counts = np.unique(atoms.resindices, + return_index=True, + return_counts=True) + + if np.count_nonzero(counts > 1): + # some residues have atoms with the same name + # take the first one? That's what MDTraj does + atoms = atoms[idx] + + # -1 means missing + found_idx = np.isin(ridx, residx, assume_unique=True) + ncoca_indices[i][found_idx] = np.arange(len(atoms)) + len(ncoca) + ncoca += atoms + + self.atomgroup = ncoca + self.nco_indices = np.ascontiguousarray(ncoca_indices[:3].T) + self.ca_indices = np.ascontiguousarray(ncoca_indices[-1]) + is_protein = (ncoca_indices != -1).sum(axis=0) == 4 + self.not_protein = ~is_protein + + pro = ag.select_atoms('resname {}'.format(self.names['PRO'])) + pro_idx = pro.residues.resindices + is_proline = np.isin(ridx, pro_idx, assume_unique=True) + self.is_proline_indices = is_proline.astype(np.int32) + + # chains need to be ints for mdtraj + if hasattr(self._universe._topology, 'chainIDs'): + chains = [r.atoms.chainIDs[0] for r in self.residues] + else: + chains = self.residues.segindices + cids, cidx = np.unique(chains, return_inverse=True) + self.chain_ids = np.arange(len(cids), dtype=np.int32)[cidx] + + def _compute_dssp(self): + # mdtraj uses angstrom + xyz = self.atomgroup.positions.reshape((1, -1, 3))/10 + xyz = np.ascontiguousarray(xyz.astype(np.float32)) + value = self._mdtraj_dssp(xyz, self.nco_indices, self.ca_indices, + self.is_proline_indices, + self.chain_ids) + arr = np.fromiter(value, dtype=' False -> True will wipe data from first True state + self._pos = np.zeros((self.n_atoms, 3), dtype=np.float32, + order=self.order) + self._has_positions = True + elif not val: + # Unsetting val won't delete the numpy array + self._has_positions = False + + @property + def positions(self): + """A record of the positions of all atoms in this Timestep + + Setting this attribute will add positions to the Timestep if they + weren't originally present. + + Returns + ------- + positions : numpy.ndarray with dtype numpy.float32 + position data of shape ``(n_atoms, 3)`` for all atoms + + Raises + ------ + :exc:`MDAnalysis.exceptions.NoDataError` + if the Timestep has no position data + + + .. versionchanged:: 0.11.0 + Now can raise :exc:`NoDataError` when no position data present + """ + if self.has_positions: + return self._pos + else: + raise NoDataError("This Timestep has no positions") + + @positions.setter + def positions(self, new): + self.has_positions = True + self._pos[:] = new + + @property + def _x(self): + """A view onto the x dimension of position data + + .. versionchanged:: 0.11.0 + Now read only + """ + return self.positions[:, 0] + + @property + def _y(self): + """A view onto the y dimension of position data + + .. versionchanged:: 0.11.0 + Now read only + """ + return self.positions[:, 1] + + @property + def _z(self): + """A view onto the z dimension of position data + + .. versionchanged:: 0.11.0 + Now read only + """ + return self.positions[:, 2] + + @property + def has_velocities(self): + """A boolean of whether this Timestep has velocity data + + This can be changed to ``True`` or ``False`` to allocate space for + or remove the data. + + .. versionadded:: 0.11.0 + """ + return self._has_velocities + + @has_velocities.setter + def has_velocities(self, val): + if val and not self._has_velocities: + self._velocities = np.zeros((self.n_atoms, 3), dtype=np.float32, + order=self.order) + self._has_velocities = True + elif not val: + self._has_velocities = False + + @property + def velocities(self): + """A record of the velocities of all atoms in this Timestep + + Setting this attribute will add velocities to the Timestep if they + weren't originally present. + + Returns + ------- + velocities : numpy.ndarray with dtype numpy.float32 + velocity data of shape ``(n_atoms, 3)`` for all atoms + + Raises + ------ + :exc:`MDAnalysis.exceptions.NoDataError` + if the Timestep has no velocity data + + + .. versionadded:: 0.11.0 + """ + if self.has_velocities: + return self._velocities + else: + raise NoDataError("This Timestep has no velocities") + + @velocities.setter + def velocities(self, new): + self.has_velocities = True + self._velocities[:] = new + + @property + def has_forces(self): + """A boolean of whether this Timestep has force data + + This can be changed to ``True`` or ``False`` to allocate space for + or remove the data. + + .. versionadded:: 0.11.0 + """ + return self._has_forces + + @has_forces.setter + def has_forces(self, val): + if val and not self._has_forces: + self._forces = np.zeros((self.n_atoms, 3), dtype=np.float32, + order=self.order) + self._has_forces = True + elif not val: + self._has_forces = False + + @property + def forces(self): + """A record of the forces of all atoms in this Timestep + + Setting this attribute will add forces to the Timestep if they + weren't originally present. + + Returns + ------- + forces : numpy.ndarray with dtype numpy.float32 + force data of shape ``(n_atoms, 3)`` for all atoms + + Raises + ------ + :exc:`MDAnalysis.exceptions.NoDataError` + if the Timestep has no force data + + + .. versionadded:: 0.11.0 + """ + if self.has_forces: + return self._forces + else: + raise NoDataError("This Timestep has no forces") + + @forces.setter + def forces(self, new): + self.has_forces = True + self._forces[:] = new + + @property + def dimensions(self): + """unitcell dimensions (*A*, *B*, *C*, *alpha*, *beta*, *gamma*) + + lengths *a*, *b*, *c* are in the MDAnalysis length unit (Å), and + angles are in degrees. + + Setting dimensions will populate the underlying native format + description of box dimensions + """ + # The actual Timestep._unitcell depends on the underlying + # trajectory format. It can be e.g. six floats representing + # the box edges and angles or the 6 unique components of the + # box matrix or the full box matrix. + return self._unitcell + + @dimensions.setter + def dimensions(self, box): + self._unitcell[:] = box + + @property + def volume(self): + """volume of the unitcell""" + return core.box_volume(self.dimensions) + + @property + def triclinic_dimensions(self): + """The unitcell dimensions represented as triclinic vectors + + Returns + ------- + numpy.ndarray + A (3, 3) numpy.ndarray of unit cell vectors + + Examples + -------- + The unitcell for a given system can be queried as either three + vectors lengths followed by their respective angle, or as three + triclinic vectors. + + >>> ts.dimensions + array([ 13., 14., 15., 90., 90., 90.], dtype=float32) + >>> ts.triclinic_dimensions + array([[ 13., 0., 0.], + [ 0., 14., 0.], + [ 0., 0., 15.]], dtype=float32) + + Setting the attribute also works:: + + >>> ts.triclinic_dimensions = [[15, 0, 0], [5, 15, 0], [5, 5, 15]] + >>> ts.dimensions + array([ 15. , 15.81138802, 16.58312416, 67.58049774, + 72.45159912, 71.56504822], dtype=float32) + + See Also + -------- + :func:`MDAnalysis.lib.mdamath.triclinic_vectors` + + + .. versionadded:: 0.11.0 + """ + return core.triclinic_vectors(self.dimensions) + + @triclinic_dimensions.setter + def triclinic_dimensions(self, new): + """Set the unitcell for this Timestep as defined by triclinic vectors + + .. versionadded:: 0.11.0 + """ + self.dimensions = core.triclinic_box(*new) + + @property + def dt(self): + """The time difference in ps between timesteps + + Note + ---- + This defaults to 1.0 ps in the absence of time data + + + .. versionadded:: 0.11.0 + """ + try: + return self.data['dt'] + except KeyError: + pass + try: + dt = self.data['dt'] = self._reader()._get_dt() + return dt + except AttributeError: + pass + warnings.warn("Reader has no dt information, set to 1.0 ps") + return 1.0 + + @dt.setter + def dt(self, new): + self.data['dt'] = new + + @dt.deleter + def dt(self): + del self.data['dt'] + + @property + def time(self): + """The time in ps of this timestep + + This is calculated as:: + + time = ts.data['time_offset'] + ts.time + + Or, if the trajectory doesn't provide time information:: + + time = ts.data['time_offset'] + ts.frame * ts.dt + + .. versionadded:: 0.11.0 + """ + offset = self.data.get('time_offset', 0) + try: + return self.data['time'] + offset + except KeyError: + return self.dt * self.frame + offset + + @time.setter + def time(self, new): + self.data['time'] = new + + @time.deleter + def time(self): + del self.data['time'] + + +class FrameIteratorBase(object): + """ + Base iterable over the frames of a trajectory. + + A frame iterable has a length that can be accessed with the :func:`len` + function, and can be indexed similarly to a full trajectory. When indexed, + indices are resolved relative to the iterable and not relative to the + trajectory. + + .. versionadded:: 0.19.0 + + """ + def __init__(self, trajectory): + self._trajectory = trajectory + + def __len__(self): + raise NotImplementedError() + + @staticmethod + def _avoid_bool_list(frames): + if isinstance(frames, list) and frames and isinstance(frames[0], bool): + return np.array(frames, dtype=bool) + return frames + + @property + def trajectory(self): + return self._trajectory + + +class FrameIteratorSliced(FrameIteratorBase): + """ + Iterable over the frames of a trajectory on the basis of a slice. + + Parameters + ---------- + trajectory: ProtoReader + The trajectory over which to iterate. + frames: slice + A slice to select the frames of interest. + + See Also + -------- + FrameIteratorBase + + .. versionadded:: 0.19.0 + + """ + def __init__(self, trajectory, frames): + # It would be easier to store directly a range object, as it would + # store its parameters in a single place, calculate its length, and + # take care of most the indexing. Though, doing so is not compatible + # with python 2 where xrange (or range with six) is only an iterator. + super(FrameIteratorSliced, self).__init__(trajectory) + self._start, self._stop, self._step = trajectory.check_slice_indices( + frames.start, frames.stop, frames.step, + ) + + def __len__(self): + return range_length(self.start, self.stop, self.step) + + def __iter__(self): + for i in range(self.start, self.stop, self.step): + yield self.trajectory[i] + self.trajectory.rewind() + + def __getitem__(self, frame): + if isinstance(frame, numbers.Integral): + length = len(self) + if not -length < frame < length: + raise IndexError('Index {} is out of range of the range of length {}.' + .format(frame, length)) + if frame < 0: + frame = len(self) + frame + frame = self.start + frame * self.step + return self.trajectory._read_frame_with_aux(frame) + elif isinstance(frame, slice): + step = (frame.step or 1) * self.step + if frame.start is None: + if frame.step is None or frame.step > 0: + start = self.start + else: + start = self.start + (len(self) - 1) * self.step + else: + start = self.start + (frame.start or 0) * self.step + if frame.stop is None: + if frame.step is None or frame.step > 0: + last = start + (range_length(start, self.stop, step) - 1) * step + else: + last = self.start + stop = last + np.sign(step) + else: + stop = self.start + (frame.stop or 0) * self.step + + new_slice = slice(start, stop, step) + frame_iterator = FrameIteratorSliced(self.trajectory, new_slice) + # The __init__ of FrameIteratorSliced does some conversion between + # the way indices are handled in slices and the way they are + # handled by range. We need to overwrite this conversion as we + # already use the logic for range. + frame_iterator._start = start + frame_iterator._stop = stop + frame_iterator._step = step + return frame_iterator + else: + # Indexing with a lists of bools does not behave the same in all + # version of numpy. + frame = self._avoid_bool_list(frame) + frames = np.array(list(range(self.start, self.stop, self.step)))[frame] + return FrameIteratorIndices(self.trajectory, frames) + + @property + def start(self): + return self._start + + @property + def stop(self): + return self._stop + + @property + def step(self): + return self._step + + +class FrameIteratorAll(FrameIteratorBase): + """ + Iterable over all the frames of a trajectory. + + Parameters + ---------- + trajectory: ProtoReader + The trajectory over which to iterate. + + See Also + -------- + FrameIteratorBase + + .. versionadded:: 0.19.0 + + """ + def __init__(self, trajectory): + super(FrameIteratorAll, self).__init__(trajectory) + + def __len__(self): + return self.trajectory.n_frames + + def __iter__(self): + return iter(self.trajectory) + + def __getitem__(self, frame): + return self.trajectory[frame] + + +class FrameIteratorIndices(FrameIteratorBase): + """ + Iterable over the frames of a trajectory listed in a sequence of indices. + + Parameters + ---------- + trajectory: ProtoReader + The trajectory over which to iterate. + frames: sequence + A sequence of indices. + + See Also + -------- + FrameIteratorBase + """ + def __init__(self, trajectory, frames): + super(FrameIteratorIndices, self).__init__(trajectory) + self._frames = [] + for frame in frames: + if not isinstance(frame, numbers.Integral): + raise TypeError("Frames indices must be integers.") + frame = trajectory._apply_limits(frame) + self._frames.append(frame) + self._frames = tuple(self._frames) + + def __len__(self): + return len(self.frames) + + def __iter__(self): + for frame in self.frames: + yield self.trajectory._read_frame_with_aux(frame) + + def __getitem__(self, frame): + if isinstance(frame, numbers.Integral): + frame = self.frames[frame] + return self.trajectory._read_frame_with_aux(frame) + else: + frame = self._avoid_bool_list(frame) + frames = np.array(self.frames)[frame] + return FrameIteratorIndices(self.trajectory, frames) + + @property + def frames(self): + return self._frames + + +class IOBase(object): + """Base class bundling common functionality for trajectory I/O. + + .. versionchanged:: 0.8 + Added context manager protocol. + """ + + #: dict with units of of *time* and *length* (and *velocity*, *force*, + #: ... for formats that support it) + units = {'time': None, 'length': None, 'velocity': None} + + def convert_pos_from_native(self, x, inplace=True): + """Conversion of coordinate array x from native units to base units. + + Parameters + ---------- + x : array_like + Positions to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input `x` is modified in place and also returned. + In-place operations improve performance because allocating new arrays + is avoided. + + + .. versionchanged:: 0.7.5 + Keyword `inplace` can be set to ``False`` so that a + modified copy is returned *unless* no conversion takes + place, in which case the reference to the unmodified `x` is + returned. + + """ + f = units.get_conversion_factor( + 'length', self.units['length'], flags['length_unit']) + if f == 1.: + return x + if not inplace: + return f * x + x *= f + return x + + def convert_velocities_from_native(self, v, inplace=True): + """Conversion of velocities array *v* from native to base units + + Parameters + ---------- + v : array_like + Velocities to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input *v* is modified in place and also returned. + In-place operations improve performance because allocating new arrays + is avoided. + + + .. versionadded:: 0.7.5 + """ + f = units.get_conversion_factor( + 'speed', self.units['velocity'], flags['speed_unit']) + if f == 1.: + return v + if not inplace: + return f * v + v *= f + return v + + def convert_forces_from_native(self, force, inplace=True): + """Conversion of forces array *force* from native to base units + + Parameters + ---------- + force : array_like + Forces to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input *force* is modified in place and also returned. + In-place operations improve performance because allocating new arrays + is avoided. + + .. versionadded:: 0.7.7 + """ + f = units.get_conversion_factor( + 'force', self.units['force'], flags['force_unit']) + if f == 1.: + return force + if not inplace: + return f * force + force *= f + return force + + def convert_time_from_native(self, t, inplace=True): + """Convert time *t* from native units to base units. + + Parameters + ---------- + t : array_like + Time values to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input `t` is modified in place and also returned + (although note that scalar values `t` are passed by value in Python and + hence an in-place modification has no effect on the caller.) In-place + operations improve performance because allocating new arrays is + avoided. + + + .. versionchanged:: 0.7.5 + Keyword `inplace` can be set to ``False`` so that a + modified copy is returned *unless* no conversion takes + place, in which case the reference to the unmodified `x` is + returned. + + """ + f = units.get_conversion_factor( + 'time', self.units['time'], flags['time_unit']) + if f == 1.: + return t + if not inplace: + return f * t + t *= f + return t + + def convert_pos_to_native(self, x, inplace=True): + """Conversion of coordinate array `x` from base units to native units. + + Parameters + ---------- + x : array_like + Positions to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input `x` is modified in place and also returned. + In-place operations improve performance because allocating new arrays + is avoided. + + + .. versionchanged:: 0.7.5 + Keyword `inplace` can be set to ``False`` so that a + modified copy is returned *unless* no conversion takes + place, in which case the reference to the unmodified `x` is + returned. + + """ + f = units.get_conversion_factor( + 'length', flags['length_unit'], self.units['length']) + if f == 1.: + return x + if not inplace: + return f * x + x *= f + return x + + def convert_velocities_to_native(self, v, inplace=True): + """Conversion of coordinate array `v` from base to native units + + Parameters + ---------- + v : array_like + Velocities to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input `v` is modified in place and also returned. + In-place operations improve performance because allocating new arrays + is avoided. + + + .. versionadded:: 0.7.5 + """ + f = units.get_conversion_factor( + 'speed', flags['speed_unit'], self.units['velocity']) + if f == 1.: + return v + if not inplace: + return f * v + v *= f + return v + + def convert_forces_to_native(self, force, inplace=True): + """Conversion of force array *force* from base to native units. + + Parameters + ---------- + force : array_like + Forces to transform + inplace : bool (optional) + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input `force` is modified in place and also returned. + In-place operations improve performance because allocating new arrays + is avoided. + + + .. versionadded:: 0.7.7 + """ + f = units.get_conversion_factor( + 'force', flags['force_unit'], self.units['force']) + if f == 1.: + return force + if not inplace: + return f * force + force *= f + return force + + def convert_time_to_native(self, t, inplace=True): + """Convert time *t* from base units to native units. + + Parameters + ---------- + t : array_like + Time values to transform + inplace : bool, optional + Whether to modify the array inplace, overwriting previous data + + Note + ---- + By default, the input *t* is modified in place and also + returned. (Also note that scalar values *t* are passed by + value in Python and hence an in-place modification has no + effect on the caller.) + + .. versionchanged:: 0.7.5 + Keyword *inplace* can be set to ``False`` so that a + modified copy is returned *unless* no conversion takes + place, in which case the reference to the unmodified *x* is + returned. + + """ + f = units.get_conversion_factor( + 'time', flags['time_unit'], self.units['time']) + if f == 1.: + return t + if not inplace: + return f * t + t *= f + return t + + def close(self): + """Close the trajectory file.""" + pass # pylint: disable=unnecessary-pass + + def __enter__(self): + return self + + def __exit__(self, exc_type, exc_val, exc_tb): + # see http://docs.python.org/2/library/stdtypes.html#typecontextmanager + self.close() + return False # do not suppress exceptions + + +class _Readermeta(type): + # Auto register upon class creation + def __init__(cls, name, bases, classdict): + type.__init__(type, name, bases, classdict) + try: + fmt = asiterable(classdict['format']) + except KeyError: + pass + else: + for f in fmt: + f = f.upper() + _READERS[f] = cls + + +class ProtoReader(six.with_metaclass(_Readermeta, IOBase)): + """Base class for Readers, without a :meth:`__del__` method. + + Extends :class:`IOBase` with most attributes and methods of a generic + Reader, with the exception of a :meth:`__del__` method. It should be used + as base for Readers that do not need :meth:`__del__`, especially since + having even an empty :meth:`__del__` might lead to memory leaks. + + See the :ref:`Trajectory API` definition in + :mod:`MDAnalysis.coordinates.__init__` for the required attributes and + methods. + + See Also + -------- + :class:`ReaderBase` + + + .. versionchanged:: 0.11.0 + Frames now 0-based instead of 1-based + """ + + #: The appropriate Timestep class, e.g. + #: :class:`MDAnalysis.coordinates.xdrfile.XTC.Timestep` for XTC. + _Timestep = Timestep + + def __init__(self): + # initialise list to store added auxiliary readers in + # subclasses should now call super + self._auxs = {} + self._transformations=[] + + def __len__(self): + return self.n_frames + + @classmethod + def parse_n_atoms(cls, filename, **kwargs): + """Read the coordinate file and deduce the number of atoms + + Returns + ------- + n_atoms : int + the number of atoms in the coordinate file + + Raises + ------ + NotImplementedError + when the number of atoms can't be deduced + """ + raise NotImplementedError("{} cannot deduce the number of atoms" + "".format(cls.__name__)) + + def next(self): + """Forward one step to next frame.""" + try: + ts = self._read_next_timestep() + except (EOFError, IOError): + self.rewind() + six.raise_from(StopIteration, None) + else: + for auxname in self.aux_list: + ts = self._auxs[auxname].update_ts(ts) + + ts = self._apply_transformations(ts) + + return ts + + def __next__(self): + """Forward one step to next frame when using the `next` builtin.""" + return self.next() + + def rewind(self): + """Position at beginning of trajectory""" + self._reopen() + self.next() + + @property + def dt(self): + """Time between two trajectory frames in picoseconds.""" + return self.ts.dt + + @property + def totaltime(self): + """Total length of the trajectory + + The time is calculated as ``(n_frames - 1) * dt``, i.e., we assume that + the first frame no time as elapsed. Thus, a trajectory with two frames will + be considered to have a length of a single time step `dt` and a + "trajectory" with a single frame will be reported as length 0. + + """ + return (self.n_frames - 1) * self.dt + + @property + def frame(self): + """Frame number of the current time step. + + This is a simple short cut to :attr:`Timestep.frame`. + """ + return self.ts.frame + + @property + def time(self): + """Time of the current frame in MDAnalysis time units (typically ps). + + This is either read straight from the Timestep, or calculated as + time = :attr:`Timestep.frame` * :attr:`Timestep.dt` + """ + return self.ts.time + + @property + def trajectory(self): + # Makes a reader effectively commpatible with a FrameIteratorBase + return self + + def Writer(self, filename, **kwargs): + """A trajectory writer with the same properties as this trajectory.""" + raise NotImplementedError( + "Sorry, there is no Writer for this format in MDAnalysis. " + "Please file an enhancement request at " + "https://github.com/MDAnalysis/mdanalysis/issues") + + def OtherWriter(self, filename, **kwargs): + """Returns a writer appropriate for *filename*. + + Sets the default keywords *start*, *step* and *dt* (if + available). *n_atoms* is always set from :attr:`Reader.n_atoms`. + + + See Also + -------- + :meth:`Reader.Writer` and :func:`MDAnalysis.Writer` + + """ + kwargs['n_atoms'] = self.n_atoms # essential + kwargs.setdefault('start', self.frame) + try: + kwargs.setdefault('dt', self.dt) + except KeyError: + pass + return core.writer(filename, **kwargs) + + def _read_next_timestep(self, ts=None): # pragma: no cover + # Example from DCDReader: + # if ts is None: + # ts = self.ts + # ts.frame = self._read_next_frame(etc) + # return ts + raise NotImplementedError( + "BUG: Override _read_next_timestep() in the trajectory reader!") + + def __iter__(self): + """ Iterate over trajectory frames. """ + self._reopen() + return self + + def _reopen(self): + """Should position Reader to just before first frame + + Calling next after this should return the first frame + """ + pass # pylint: disable=unnecessary-pass + + def _apply_limits(self, frame): + if frame < 0: + frame += len(self) + if frame < 0 or frame >= len(self): + raise IndexError("Index {} exceeds length of trajectory ({})." + "".format(frame, len(self))) + return frame + + def __getitem__(self, frame): + """Return the Timestep corresponding to *frame*. + + If `frame` is a integer then the corresponding frame is + returned. Negative numbers are counted from the end. + + If frame is a :class:`slice` then an iterator is returned that + allows iteration over that part of the trajectory. + + Note + ---- + *frame* is a 0-based frame index. + """ + if isinstance(frame, numbers.Integral): + frame = self._apply_limits(frame) + return self._read_frame_with_aux(frame) + elif isinstance(frame, (list, np.ndarray)): + if len(frame) != 0 and isinstance(frame[0], (bool, np.bool_)): + # Avoid having list of bools + frame = np.asarray(frame, dtype=np.bool) + # Convert bool array to int array + frame = np.arange(len(self))[frame] + return FrameIteratorIndices(self, frame) + elif isinstance(frame, slice): + start, stop, step = self.check_slice_indices( + frame.start, frame.stop, frame.step) + if start == 0 and stop == len(self) and step == 1: + return FrameIteratorAll(self) + else: + return FrameIteratorSliced(self, frame) + else: + raise TypeError("Trajectories must be an indexed using an integer," + " slice or list of indices") + + def _read_frame(self, frame): + """Move to *frame* and fill timestep with data.""" + raise TypeError("{0} does not support direct frame indexing." + "".format(self.__class__.__name__)) + # Example implementation in the DCDReader: + # self._jump_to_frame(frame) + # ts = self.ts + # ts.frame = self._read_next_frame(ts._x, ts._y, ts._z, + # ts._unitcell, 1) + # return ts + + def _read_frame_with_aux(self, frame): + """Move to *frame*, updating ts with trajectory, transformations and auxiliary data.""" + ts = self._read_frame(frame) # pylint: disable=assignment-from-no-return + for aux in self.aux_list: + ts = self._auxs[aux].update_ts(ts) + + ts = self._apply_transformations(ts) + + return ts + + def _sliced_iter(self, start, stop, step): + """Generator for slicing a trajectory. + + *start* *stop* and *step* are 3 integers describing the slice. + Error checking is not done past this point. + + A :exc:`NotImplementedError` is raised if random access to + frames is not implemented. + """ + # override with an appropriate implementation e.g. using self[i] might + # be much slower than skipping steps in a next() loop + try: + for i in range(start, stop, step): + yield self._read_frame_with_aux(i) + self.rewind() + except TypeError: # if _read_frame not implemented + six.raise_from( + TypeError( + "{0} does not support slicing." + "".format(self.__class__.__name__)), + None) + + def check_slice_indices(self, start, stop, step): + """Check frame indices are valid and clip to fit trajectory. + + The usage follows standard Python conventions for :func:`range` but see + the warning below. + + Parameters + ---------- + start : int or None + Starting frame index (inclusive). ``None`` corresponds to the default + of 0, i.e., the initial frame. + stop : int or None + Last frame index (exclusive). ``None`` corresponds to the default + of n_frames, i.e., it includes the last frame of the trajectory. + step : int or None + step size of the slice, ``None`` corresponds to the default of 1, i.e, + include every frame in the range `start`, `stop`. + + Returns + ------- + start, stop, step : tuple (int, int, int) + Integers representing the slice + + Warning + ------- + The returned values `start`, `stop` and `step` give the expected result + when passed in :func:`range` but gives unexpected behavior when passed + in a :class:`slice` when ``stop=None`` and ``step=-1`` + + This can be a problem for downstream processing of the output from this + method. For example, slicing of trajectories is implemented by passing + the values returned by :meth:`check_slice_indices` to :func:`range` :: + + range(start, stop, step) + + and using them as the indices to randomly seek to. On the other hand, + in :class:`MDAnalysis.analysis.base.AnalysisBase` the values returned + by :meth:`check_slice_indices` are used to splice the trajectory by + creating a :class:`slice` instance :: + + slice(start, stop, step) + + This creates a discrepancy because these two lines are not equivalent:: + + range(10, -1, -1) # [10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0] + range(10)[slice(10, -1, -1)] # [] + + """ + + slice_dict = {'start': start, 'stop': stop, 'step': step} + for varname, var in slice_dict.items(): + if isinstance(var, numbers.Integral): + slice_dict[varname] = int(var) + elif (var is None): + pass + else: + raise TypeError("{0} is not an integer".format(varname)) + + start = slice_dict['start'] + stop = slice_dict['stop'] + step = slice_dict['step'] + + if step == 0: + raise ValueError("Step size is zero") + + nframes = len(self) + step = step or 1 + + if start is None: + start = 0 if step > 0 else nframes - 1 + elif start < 0: + start += nframes + if start < 0: + start = 0 + + if step < 0 and start >= nframes: + start = nframes - 1 + + if stop is None: + stop = nframes if step > 0 else -1 + elif stop < 0: + stop += nframes + + if step > 0 and stop > nframes: + stop = nframes + + return start, stop, step + + def __repr__(self): + return ("<{cls} {fname} with {nframes} frames of {natoms} atoms>" + "".format( + cls=self.__class__.__name__, + fname=self.filename, + nframes=self.n_frames, + natoms=self.n_atoms + )) + + def add_auxiliary(self, auxname, auxdata, format=None, **kwargs): + """Add auxiliary data to be read alongside trajectory. + + Auxiliary data may be any data timeseries from the trajectory additional + to that read in by the trajectory reader. *auxdata* can be an + :class:`~MDAnalysis.auxiliary.base.AuxReader` instance, or the data + itself as e.g. a filename; in the latter case an appropriate + :class:`~MDAnalysis.auxiliary.base.AuxReader` is guessed from the + data/file format. An appropriate `format` may also be directly provided + as a key word argument. + + On adding, the AuxReader is initially matched to the current timestep + of the trajectory, and will be updated when the trajectory timestep + changes (through a call to :meth:`next()` or jumping timesteps with + ``trajectory[i]``). + + The representative value(s) of the auxiliary data for each timestep (as + calculated by the :class:`~MDAnalysis.auxiliary.base.AuxReader`) are + stored in the current timestep in the ``ts.aux`` namespace under *auxname*; + e.g. to add additional pull force data stored in pull-force.xvg:: + + u = MDAnalysis.Universe(PDB, XTC) + u.trajectory.add_auxiliary('pull', 'pull-force.xvg') + + The representative value for the current timestep may then be accessed + as ``u.trajectory.ts.aux.pull`` or ``u.trajectory.ts.aux['pull']``. + + See Also + -------- + :meth:`remove_auxiliary` + + Note + ---- + Auxiliary data is assumed to be time-ordered, with no duplicates. See + the :ref:`Auxiliary API`. + """ + if auxname in self.aux_list: + raise ValueError("Auxiliary data with name {name} already " + "exists".format(name=auxname)) + if isinstance(auxdata, AuxReader): + aux = auxdata + aux.auxname = auxname + else: + aux = auxreader(auxdata, format=format, auxname=auxname, **kwargs) + self._auxs[auxname] = aux + self.ts = aux.update_ts(self.ts) + + def remove_auxiliary(self, auxname): + """Clear data and close the :class:`~MDAnalysis.auxiliary.base.AuxReader` + for the auxiliary *auxname*. + + See Also + -------- + :meth:`add_auxiliary` + """ + aux = self._check_for_aux(auxname) + aux.close() + del aux + delattr(self.ts.aux, auxname) + + @property + def aux_list(self): + """ Lists the names of added auxiliary data. """ + return self._auxs.keys() + + def _check_for_aux(self, auxname): + """ Check for the existance of an auxiliary *auxname*. If present, + return the AuxReader; if not, raise ValueError + """ + if auxname in self.aux_list: + return self._auxs[auxname] + else: + raise ValueError("No auxiliary named {name}".format(name=auxname)) + + def next_as_aux(self, auxname): + """ Move to the next timestep for which there is at least one step from + the auxiliary *auxname* within the cutoff specified in *auxname*. + + This allows progression through the trajectory without encountering + ``NaN`` representative values (unless these are specifically part of the + auxiliary data). + + If the auxiliary cutoff is not set, where auxiliary steps are less frequent + (``auxiliary.dt > trajectory.dt``), this allows progression at the + auxiliary pace (rounded to nearest timestep); while if the auxiliary + steps are more frequent, this will work the same as calling + :meth:`next()`. + + See the :ref:`Auxiliary API`. + + See Also + -------- + :meth:`iter_as_aux` + """ + + aux = self._check_for_aux(auxname) + ts = self.ts + # catch up auxiliary if it starts earlier than trajectory + while aux.step_to_frame(aux.step + 1, ts) is None: + next(aux) + # find the next frame that'll have a representative value + next_frame = aux.next_nonempty_frame(ts) + if next_frame is None: + # no more frames with corresponding auxiliary values; stop iteration + raise StopIteration + # some readers set self._frame to -1, rather than self.frame, on + # _reopen; catch here or doesn't read first frame + while self.frame != next_frame or getattr(self, '_frame', 0) == -1: + # iterate trajectory until frame is reached + ts = self.next() + return ts + + def iter_as_aux(self, auxname): + """Iterate through timesteps for which there is at least one assigned + step from the auxiliary *auxname* within the cutoff specified in *auxname*. + + See Also + -------- + :meth:`next_as_aux` + :meth:`iter_auxiliary` + """ + aux = self._check_for_aux(auxname) + self._reopen() + aux._restart() + while True: + try: + yield self.next_as_aux(auxname) + except StopIteration: + return + + def iter_auxiliary(self, auxname, start=None, stop=None, step=None, + selected=None): + """ Iterate through the auxiliary *auxname* independently of the trajectory. + + Will iterate over the specified steps of the auxiliary (defaults to all + steps). Allows to access all values in an auxiliary, including those out + of the time range of the trajectory, without having to also iterate + through the trajectory. + + After interation, the auxiliary will be repositioned at the current step. + + Parameters + ---------- + auxname : str + Name of the auxiliary to iterate over. + (start, stop, step) : optional + Options for iterating over a slice of the auxiliary. + selected : lst | ndarray, optional + List of steps to iterate over. + + Yields + ------ + :class:`~MDAnalysis.auxiliary.base.AuxStep` object + + See Also + -------- + :meth:`iter_as_aux` + """ + aux = self._check_for_aux(auxname) + if selected is not None: + selection = selected + else: + selection = slice(start, stop, step) + for i in aux[selection]: + yield i + aux.read_ts(self.ts) + + def get_aux_attribute(self, auxname, attrname): + """Get the value of *attrname* from the auxiliary *auxname* + + Parameters + ---------- + auxname : str + Name of the auxiliary to get value for + attrname : str + Name of gettable attribute in the auxiliary reader + + See Also + -------- + :meth:`set_aux_attribute` + """ + aux = self._check_for_aux(auxname) + return getattr(aux, attrname) + + def set_aux_attribute(self, auxname, attrname, new): + """ Set the value of *attrname* in the auxiliary *auxname*. + + Parameters + ---------- + auxname : str + Name of the auxiliary to alter + attrname : str + Name of settable attribute in the auxiliary reader + new + New value to try set *attrname* to + + See Also + -------- + :meth:`get_aux_attribute` + :meth:`rename_aux` - to change the *auxname* attribute + """ + aux = self._check_for_aux(auxname) + if attrname == 'auxname': + self.rename_aux(auxname, new) + else: + setattr(aux, attrname, new) + + def rename_aux(self, auxname, new): + """ Change the name of the auxiliary *auxname* to *new*. + + Provided there is not already an auxiliary named *new*, the auxiliary + name will be changed in ts.aux namespace, the trajectory's + list of added auxiliaries, and in the auxiliary reader itself. + + Parameters + ---------- + auxname : str + Name of the auxiliary to rename + new : str + New name to try set + + Raises + ------ + ValueError + If the name *new* is already in use by an existing auxiliary. + """ + aux = self._check_for_aux(auxname) + if new in self.aux_list: + raise ValueError("Auxiliary data with name {name} already " + "exists".format(name=new)) + aux.auxname = new + self._auxs[new] = self._auxs.pop(auxname) + setattr(self.ts.aux, new, self.ts.aux[auxname]) + delattr(self.ts.aux, auxname) + + def get_aux_descriptions(self, auxnames=None): + """Get descriptions to allow reloading the specified auxiliaries. + + If no auxnames are provided, defaults to the full list of added + auxiliaries. + + Passing the resultant description to ``add_auxiliary()`` will allow + recreation of the auxiliary. e.g., to duplicate all auxiliaries into a + second trajectory:: + + descriptions = trajectory_1.get_aux_descriptions() + for aux in descriptions: + trajectory_2.add_auxiliary(**aux) + + + Returns + ------- + list + List of dictionaries of the args/kwargs describing each auxiliary. + + See Also + -------- + :meth:`MDAnalysis.auxiliary.base.AuxReader.get_description` + """ + if not auxnames: + auxnames = self.aux_list + descriptions = [self._auxs[aux].get_description() for aux in auxnames] + return descriptions + + @property + def transformations(self): + """ Returns the list of transformations""" + return self._transformations + + @transformations.setter + def transformations(self, transformations): + if not self._transformations: + self._transformations = transformations + else: + raise ValueError("Transformations are already set") + + def add_transformations(self, *transformations): + """Add all transformations to be applied to the trajectory. + + This function take as list of transformations as an argument. These + transformations are functions that will be called by the Reader and given + a :class:`Timestep` object as argument, which will be transformed and returned + to the Reader. + The transformations can be part of the :mod:`~MDAnalysis.transformations` + module, or created by the user, and are stored as a list `transformations`. + This list can only be modified once, and further calls of this function will + raise an exception. + + .. code-block:: python + + u = MDAnalysis.Universe(topology, coordinates) + workflow = [some_transform, another_transform, this_transform] + u.trajectory.add_transformations(*workflow) + + The transformations are applied in the order given in the list + `transformations`, i.e., the first transformation is the first + or innermost one to be applied to the :class:`Timestep`. The + example above would be equivalent to + + .. code-block:: python + + for ts in u.trajectory: + ts = this_transform(another_transform(some_transform(ts))) + + + Parameters + ---------- + transform_list : list + list of all the transformations that will be applied to the coordinates + in the order given in the list + + See Also + -------- + :mod:`MDAnalysis.transformations` + + """ + + try: + self.transformations = transformations + except ValueError: + six.raise_from( + ValueError( + "Can't add transformations again. " + "Please create new Universe object"), + None) + else: + self.ts = self._apply_transformations(self.ts) + + + # call reader here to apply the newly added transformation on the + # current loaded frame? + + def _apply_transformations(self, ts): + """Applies all the transformations given by the user """ + + for transform in self.transformations: + ts = transform(ts) + + return ts + + + +class ReaderBase(ProtoReader): + """Base class for trajectory readers that extends :class:`ProtoReader` with a + :meth:`__del__` method. + + New Readers should subclass :class:`ReaderBase` and properly implement a + :meth:`close` method, to ensure proper release of resources (mainly file + handles). Readers that are inherently safe in this regard should subclass + :class:`ProtoReader` instead. + + See the :ref:`Trajectory API` definition in + :mod:`MDAnalysis.coordinates.__init__` for the required attributes and + methods. + + See Also + -------- + :class:`ProtoReader` + + + .. versionchanged:: 0.11.0 + Most of the base Reader class definitions were offloaded to + :class:`ProtoReader` so as to allow the subclassing of ReaderBases without a + :meth:`__del__` method. Created init method to create common + functionality, all ReaderBase subclasses must now :func:`super` through this + class. Added attribute :attr:`_ts_kwargs`, which is created in init. + Provides kwargs to be passed to :class:`Timestep` + + """ + + def __init__(self, filename, convert_units=None, **kwargs): + super(ReaderBase, self).__init__() + + self.filename = filename + + if convert_units is None: + convert_units = flags['convert_lengths'] + self.convert_units = convert_units + + ts_kwargs = {} + for att in ('dt', 'time_offset'): + try: + val = kwargs[att] + except KeyError: + pass + else: + ts_kwargs[att] = val + + self._ts_kwargs = ts_kwargs + + def copy(self): + """Return independent copy of this Reader. + + New Reader will have its own file handle and can seek/iterate + independently of the original. + + Will also copy the current state of the Timestep held in + the original Reader + """ + new = self.__class__(self.filename, + n_atoms=self.n_atoms) + if self.transformations: + new.add_transformations(*self.transformations) + # seek the new reader to the same frame we started with + new[self.ts.frame] + # then copy over the current Timestep in case it has + # been modified since initial load + new.ts = self.ts.copy() + for auxname, auxread in self._auxs.items(): + new.add_auxiliary(auxname, auxread.copy()) + return new + + def __del__(self): + for aux in self.aux_list: + self._auxs[aux].close() + self.close() + + +class _Writermeta(type): + # Auto register this format upon class creation + def __init__(cls, name, bases, classdict): + type.__init__(type, name, bases, classdict) + try: + # grab the string which describes this format + # could be either 'PDB' or ['PDB', 'ENT'] for multiple formats + fmt = asiterable(classdict['format']) + except KeyError: + # not required however + pass + else: + # does the Writer support single and multiframe writing? + single = classdict.get('singleframe', True) + multi = classdict.get('multiframe', False) + + if single: + for f in fmt: + f = f.upper() + _SINGLEFRAME_WRITERS[f] = cls + if multi: + for f in fmt: + f = f.upper() + _MULTIFRAME_WRITERS[f] = cls + + +class WriterBase(six.with_metaclass(_Writermeta, IOBase)): + """Base class for trajectory writers. + + See Trajectory API definition in :mod:`MDAnalysis.coordinates.__init__` for + the required attributes and methods. + """ + + def convert_dimensions_to_unitcell(self, ts, inplace=True): + """Read dimensions from timestep *ts* and return appropriate unitcell. + + The default is to return ``[A,B,C,alpha,beta,gamma]``; if this + is not appropriate then this method has to be overriden. + """ + # override if the native trajectory format does NOT use + # [A,B,C,alpha,beta,gamma] + lengths, angles = ts.dimensions[:3], ts.dimensions[3:] + if not inplace: + lengths = lengths.copy() + lengths = self.convert_pos_to_native(lengths) + return np.concatenate([lengths, angles]) + + def write(self, obj): + """Write current timestep, using the supplied `obj`. + + Parameters + ---------- + obj : :class:`~MDAnalysis.core.groups.AtomGroup` or :class:`~MDAnalysis.core.universe.Universe` or a :class:`Timestep` + write coordinate information associate with `obj` + + Note + ---- + The size of the `obj` must be the same as the number of atoms provided + when setting up the trajectory. + """ + if isinstance(obj, Timestep): + ts = obj + else: + try: + ts = obj.ts + except AttributeError: + try: + # special case: can supply a Universe, too... + ts = obj.trajectory.ts + except AttributeError: + six.raise_from(TypeError("No Timestep found in obj argument"), None) + return self.write_next_timestep(ts) + + def __del__(self): + self.close() + + def __repr__(self): + try: + return "< {0!s} {1!r} for {2:d} atoms >".format(self.__class__.__name__, self.filename, self.n_atoms) + except (TypeError, AttributeError): + # no trajectory loaded yet or a Writer that does not need e.g. + # self.n_atoms + return "< {0!s} {1!r} >".format(self.__class__.__name__, self.filename) + + def has_valid_coordinates(self, criteria, x): + """Returns ``True`` if all values are within limit values of their formats. + + Due to rounding, the test is asymmetric (and *min* is supposed to be negative): + + min < x <= max + + Parameters + ---------- + criteria : dict + dictionary containing the *max* and *min* values in native units + x : numpy.ndarray + ``(x, y, z)`` coordinates of atoms selected to be written out + + Returns + ------- + bool + """ + x = np.ravel(x) + return np.all(criteria["min"] < x) and np.all(x <= criteria["max"]) + + # def write_next_timestep(self, ts=None) + + +class SingleFrameReaderBase(ProtoReader): + """Base class for Readers that only have one frame. + + To use this base class, define the method :meth:`_read_first_frame` to + read from file `self.filename`. This should populate the attribute + `self.ts` with a :class:`Timestep` object. + + .. versionadded:: 0.10.0 + .. versionchanged:: 0.11.0 + Added attribute "_ts_kwargs" for subclasses + Keywords "dt" and "time_offset" read into _ts_kwargs + """ + _err = "{0} only contains a single frame" + + def __init__(self, filename, convert_units=None, n_atoms=None, **kwargs): + super(SingleFrameReaderBase, self).__init__() + + self.filename = filename + if convert_units is None: + convert_units = flags['convert_lengths'] + self.convert_units = convert_units + + self.n_frames = 1 + self.n_atom = n_atoms + + ts_kwargs = {} + for att in ('dt', 'time_offset'): + try: + val = kwargs[att] + except KeyError: + pass + else: + ts_kwargs[att] = val + + self._ts_kwargs = ts_kwargs + self._read_first_frame() + + def copy(self): + """Return independent copy of this Reader. + + New Reader will have its own file handle and can seek/iterate + independently of the original. + + Will also copy the current state of the Timestep held in + the original Reader + """ + new = self.__class__(self.filename, + n_atoms=self.n_atoms) + new.ts = self.ts.copy() + for auxname, auxread in self._auxs.items(): + new.add_auxiliary(auxname, auxread.copy()) + # since the transformations have already been applied to the frame + # simply copy the property + new.transformations = self.transformations + + return new + + def _read_first_frame(self): # pragma: no cover + # Override this in subclasses to create and fill a Timestep + pass + + def rewind(self): + self._read_first_frame() + for auxname, auxread in self._auxs.items(): + self.ts = auxread.update_ts(self.ts) + super(SingleFrameReaderBase, self)._apply_transformations(self.ts) + + def _reopen(self): + pass + + def next(self): + raise StopIteration(self._err.format(self.__class__.__name__)) + + def __iter__(self): + yield self.ts + return + + def _read_frame(self, frame): + if frame != 0: + raise IndexError(self._err.format(self.__class__.__name__)) + + return self.ts + + def close(self): + # all single frame readers should use context managers to access + # self.filename. Explicitly setting it to the null action in case + # the IOBase.close method is ever changed from that. + pass + + def add_transformations(self, *transformations): + """ Add all transformations to be applied to the trajectory. + + This function take as list of transformations as an argument. These + transformations are functions that will be called by the Reader and given + a :class:`Timestep` object as argument, which will be transformed and returned + to the Reader. + The transformations can be part of the :mod:`~MDAnalysis.transformations` + module, or created by the user, and are stored as a list `transformations`. + This list can only be modified once, and further calls of this function will + raise an exception. + + .. code-block:: python + + u = MDAnalysis.Universe(topology, coordinates) + workflow = [some_transform, another_transform, this_transform] + u.trajectory.add_transformations(*workflow) + + Parameters + ---------- + transform_list : list + list of all the transformations that will be applied to the coordinates + + See Also + -------- + :mod:`MDAnalysis.transformations` + """ + #Overrides :meth:`~MDAnalysis.coordinates.base.ProtoReader.add_transformations` + #to avoid unintended behaviour where the coordinates of each frame are transformed + #multiple times when iterating over the trajectory. + #In this method, the trajectory is modified all at once and once only. + + super(SingleFrameReaderBase, self).add_transformations(*transformations) + for transform in self.transformations: + self.ts = transform(self.ts) + + def _apply_transformations(self, ts): + """ Applies the transformations to the timestep.""" + # Overrides :meth:`~MDAnalysis.coordinates.base.ProtoReader.add_transformations` + # to avoid applying the same transformations multiple times on each frame + + return ts + + +def range_length(start, stop, step): + if (step > 0 and start < stop): + # We go from a lesser number to a larger one. + return int(1 + (stop - 1 - start) // step) + elif (step < 0 and start > stop): + # We count backward from a larger number to a lesser one. + return int(1 + (start - 1 - stop) // (-step)) + else: + # The range is empty. + return 0 + + +class _Convertermeta(type): + # Auto register upon class creation + def __init__(cls, name, bases, classdict): + type.__init__(type, name, bases, classdict) + try: + fmt = asiterable(classdict['lib']) + except KeyError: + pass + else: + for f in fmt: + f = f.upper() + _CONVERTERS[f] = cls + +class ConverterBase(six.with_metaclass(_Convertermeta, IOBase)): + """Base class for converting to other libraries. + """ + + def __repr__(self): + return "<{cls}>".format(cls=self.__class__.__name__) + + def convert(self, obj): + raise NotImplementedError diff --git a/package/MDAnalysis/coordinates/PDB.py b/package/MDAnalysis/coordinates/PDB.py index bbc34e5529c..87b3beacad6 100644 --- a/package/MDAnalysis/coordinates/PDB.py +++ b/package/MDAnalysis/coordinates/PDB.py @@ -489,7 +489,7 @@ class PDBWriter(base.WriterBase): "ATOM {serial:5d} {name:<4s}{altLoc:<1s}{resName:<4s}" "{chainID:1s}{resSeq:4d}{iCode:1s}" " {pos[0]:8.3f}{pos[1]:8.3f}{pos[2]:8.3f}{occupancy:6.2f}" - "{tempFactor:6.2f} {segID:<4s}{element:>2s}\n"), + "{tempFactor:6.2f} {segID:<4s}{element:>2s} \n"), 'REMARK': "REMARK {0}\n", 'COMPND': "COMPND {0}\n", 'HEADER': "HEADER {0}\n", diff --git a/package/doc/sphinx/source/documentation_pages/references.rst b/package/doc/sphinx/source/documentation_pages/references.rst index 7c735b3d1b4..fdb300280d6 100644 --- a/package/doc/sphinx/source/documentation_pages/references.rst +++ b/package/doc/sphinx/source/documentation_pages/references.rst @@ -140,6 +140,39 @@ If you use the hydrogen bond analysis code in .. _`10.1039/C9CP01532A`: http://dx.doi.org/10.1039/C9CP01532A +If you use the DSSP code in :mod:`MDAnalysis.analysis.secondary_structure` please cite [Kabsch1983]_ and [Touw2015]_. + +.. [Kabsch1983] W. Kabsch and C. Sander. + Dictionary of protein secondary structure: Pattern recognition of hydrogen‐bonded and geometrical features. + *Biopolymers* **22** (1983), 2577-2637. doi: `10.1002/bip.360221211`_ + +.. _`10.1002/bip.360221211`: https://doi.org/10.1002/bip.360221211 + + +.. [Touw2015] Wouter G Touw, Coos Baakman, Jon Black, Tim AH te Beek, E Krieger, Robbie P Joosten, Gert Vriend. + A series of PDB related databases for everyday needs. + *Nucleic Acids Research* **43 (Database issue)** (2015), D364-D368. doi: `10.1093/nar/gku1028`_ + +.. _`10.1093/nar/gku1028`: https://doi.org/10.1093/nar/gku1028 + +If you use the DSSP class in :mod:`MDAnalysis.analysis.secondary_structure.dssp`, please *also* cite [McGibbon2015]_. + +.. [McGibbon2015] Robert T. McGibbon, Kyle A. Beauchamp, Matthew P. Harrigan, Christoph Klein, Jason M. Swails, Carlos X. Hernández, Christian R. Schwantes, Lee-Ping Wang, Thomas J. Lane, and Vijay S. Pande. + MDTraj: A Modern Open Library for the Analysis of Molecular Dynamics Trajectories. + *Biophysical Journal* **109 (8)** (2015), 1528–1532. doi: `10.1016/j.bpj.2015.08.015`_ + +.. _`10.1016/j.bpj.2015.08.015`: https://doi.org/10.1016/j.bpj.2015.08.015 + +If you use the STRIDE code in :mod:`MDAnalysis.analysis.secondary_structure`, +please cite [Frishman1995]_ . + +.. [Frishman1995] Dmitrij Frishman and Patrick Argos. + Knowledge-Based Protein Secondary Structure Assignment. + *Proteins: Structure, Function, and Bioinformatics* **23 (4)** (1995), 566–579. + doi: `10.1002/prot.340230412`_ + +.. _`10.1002/prot.340230412`: https://doi.org/10.1002/prot.340230412 + .. _citations-using-duecredit: Citations using Duecredit diff --git a/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py b/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py new file mode 100644 index 00000000000..4f7da0e764f --- /dev/null +++ b/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py @@ -0,0 +1,301 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# 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 +# +from __future__ import absolute_import + +import pytest +import matplotlib +import numpy as np + +from numpy.testing import assert_equal, assert_almost_equal + +import MDAnalysis as mda +from MDAnalysis.analysis import secondary_structure as ss +from MDAnalysisTests import executable_not_found +from MDAnalysisTests.datafiles import (TPR, XTC, + ADK_DSSP, ADK_DSSP_SIMPLE, + ADK_MKDSSP, MKDSSP_phi_psi_sasa, + ADK_STRIDE, STRIDE_phi_psi_sasa, + ) + + +def mdtraj_found(): + try: + import mdtraj + except ImportError: + return False + else: + return True + + +class BaseTestSecondaryStructure(object): + prot_err_msg = 'Protein residues do not match secondary structure found in {}' + other_err_msg = 'Should not have secondary structure code in non-protein residues' + plot_err_msg = 'plot_content() did not produce an Axes instance' + + n_prot = 214 + + @pytest.fixture(scope='class') + def universe(self): + return mda.Universe(TPR, XTC) + + @pytest.fixture() + def dssp_all(self, universe): + return self.analysis_class(universe, select='all').run() + + + def test_correct_assignments(self, dssp, all_codes): + assert_equal(dssp.ss_codes, all_codes) + assert_equal(dssp.ss_mode[::3], self.modes_3) + assert_equal(dssp.simple_mode[:10], self.modes_simple_10) + for k in dssp.ss_counts.keys(): + assert_almost_equal(dssp.ss_counts[k], self.code_counts[k]) + for k in dssp.simple_counts.keys(): + assert_almost_equal(dssp.simple_counts[k], self.simple_counts[k]) + + def test_non_protein(self, universe, dssp_all): + assert np.all(dssp_all.ss_codes[:, self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp_all.ss_simple[:, self.n_prot:] + == ''), self.other_err_msg + assert np.all(dssp_all.ss_mode[self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp_all.simple_mode[self.n_prot:] == ''), self.other_err_msg + + def test_add_topology_attr(self, universe): + assert not hasattr(universe._topology, 'secondary_structures') + dssp = self.analysis_class(universe, select='backbone', + add_topology_attr=True) + dssp.run() + assert hasattr(universe._topology, 'secondary_structures') + + secstruct = universe.residues.secondary_structures + modes = dssp.ss_mode[:self.n_prot] + prot_msg = self.prot_err_msg.format(self.analysis_class.__name__) + assert np.all(secstruct[:self.n_prot] == modes), prot_msg + assert np.all(secstruct[self.n_prot:] == ''), self.other_err_msg + + def test_plot_all_lines(self, dssp): + ax = dssp.plot_content(kind='line', simple=False) + assert isinstance(ax, matplotlib.axes.Axes), self.plot_err_msg + lines = ax.get_lines()[:] + assert len(lines) == 8 + + def test_plot_simple_lines(self, dssp): + ax = dssp.plot_content(kind='line', simple=True) + assert isinstance(ax, matplotlib.axes.Axes), self.plot_err_msg + lines = ax.get_lines()[:] + assert len(lines) == 3 + + def test_plot_all_bars(self, dssp): + ax = dssp.plot_content(kind='bar', simple=False) + assert isinstance(ax, matplotlib.axes.Axes), self.plot_err_msg + rects = [r for r in ax.get_children() if isinstance( + r, matplotlib.patches.Rectangle)] + assert len(rects) == (8 * dssp.n_frames) + 1 + + def test_plot_n_xticks(self, dssp): + ax = dssp.plot_content(n_xticks=3) + locs = ax.get_xticks() + assert len(locs) == 3 + + ax2 = dssp.plot_content(n_xticks=100) + locs2 = ax2.get_xticks() + assert len(locs2) == 10 + +@pytest.mark.skipif(not mdtraj_found(), + reason="Tests skipped because mdtraj not installed") +class TestDSSP(BaseTestSecondaryStructure): + + modes_3 = ['C', 'E', 'C', 'T', 'H', 'H', 'H', 'H', 'C', 'E', 'T', 'H', 'H', + 'H', 'C', 'H', 'H', 'H', 'H', 'C', 'H', 'H', 'H', 'H', 'T', 'G', + 'S', 'E', 'S', 'C', 'H', 'H', 'H', 'T', 'C', 'E', 'E', 'C', 'H', + 'H', 'H', 'E', 'T', 'S', 'E', 'T', 'C', 'S', 'B', 'T', 'C', 'S', + 'G', 'S', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'S', 'E', + 'E', 'S', 'H', 'H', 'H', 'H', 'C'] + + modes_simple_10 = ['Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', + 'Coil', 'Coil', 'Coil'] + + code_counts = {'G': np.array([7, 6, 3, 3, 3, 8, 3, 6, 0, 3]), + 'H': np.array([100, 96, 94, 97, 89, 101, 97, 96, 97, 96]), + 'I': np.array([0, 0, 0, 0, 5, 0, 0, 0, 0, 0]), + 'E': np.array([24, 28, 28, 27, 28, 29, 28, 31, 30, 30]), + 'B': np.array([5, 4, 3, 5, 4, 1, 3, 3, 1, 2]), + 'T': np.array([23, 29, 34, 28, 29, 19, 32, 27, 25, 27]), + 'S': np.array([20, 20, 22, 22, 23, 27, 22, 21, 31, 21]), + 'C': np.array([35, 31, 30, 32, 33, 29, 29, 30, 30, 35])} + + simple_counts = {'Helix': np.array([107, 102, 97, 100, 97, 109, 100, 102, 97, 99]), + 'Coil': np.array([78, 80, 86, 82, 85, 75, 83, 78, 86, 83]), + 'Strand': np.array([29, 32, 31, 32, 32, 30, 31, 34, 31, 32])} + + analysis_class = ss.DSSP + + @pytest.fixture() + def dssp(self, universe): + pytest.importorskip('mdtraj') + dssp = ss.DSSP(universe, select='backbone', o_name='O O1') + dssp.run() + return dssp + + @pytest.fixture() + def all_codes(self): + return np.load(ADK_DSSP) + + @pytest.fixture() + def all_simple(self, dssp): + mdtraj_simple = np.load(ADK_DSSP_SIMPLE).astype(object) + for code in ['H', 'E', 'C']: + matches = mdtraj_simple == code + mdtraj_simple[matches] = dssp.ss_codes_to_simple[code] + return mdtraj_simple + + def test_correct_assignments(self, dssp, all_codes, all_simple): + super(TestDSSP, self).test_correct_assignments(dssp, all_codes) + assert_equal(dssp.ss_simple, all_simple) + + def test_non_protein(self, universe, dssp_all, all_codes, all_simple): + super(TestDSSP, self).test_non_protein(universe, dssp_all) + + assert_equal(dssp_all.ss_codes[:, :self.n_prot], all_codes) + assert_equal(dssp_all.ss_simple[:, :self.n_prot], all_simple) + assert_equal(dssp_all.ss_mode[:self.n_prot:3], self.modes_3) + assert_equal(dssp_all.simple_mode[:10], self.modes_simple_10) + + for k in dssp_all.ss_counts.keys(): + assert_almost_equal(dssp_all.ss_counts[k], self.code_counts[k]) + for k in dssp_all.simple_counts.keys(): + assert_almost_equal(dssp_all.simple_counts[k], self.simple_counts[k]) + + +class BaseTestSecondaryStructureWrapper(BaseTestSecondaryStructure): + + phi_psi_sasa_file = None + + def test_executable_not_found(self, universe): + with pytest.raises(OSError) as exc: + dssp = self.analysis_class(universe, executable='foo') + assert 'executable not found' in str(exc.value) + + def test_correct_assignments(self, dssp, all_codes): + super(BaseTestSecondaryStructureWrapper, self).test_correct_assignments(dssp, all_codes) + phi, psi, sasa = np.load(self.phi_psi_sasa_file) + print(self.phi_psi_sasa_file) + assert_almost_equal(dssp.phi, phi) + assert_almost_equal(dssp.psi, psi) + assert_almost_equal(dssp.sasa, sasa) + + + +@pytest.mark.skipif(executable_not_found('mkdssp'), + reason="Tests skipped because mkdssp (DSSP) not found") +class TestDSSPWrapper(BaseTestSecondaryStructureWrapper): + + phi_psi_sasa_file = MKDSSP_phi_psi_sasa + + modes_3 = ['C', 'E', 'C', 'T', 'H', 'H', 'H', 'H', 'C', 'E', 'T', 'H', 'H', + 'H', 'C', 'H', 'H', 'H', 'H', 'C', 'H', 'H', 'H', 'H', 'T', 'G', + 'S', 'E', 'S', 'C', 'H', 'H', 'H', 'T', 'C', 'E', 'E', 'C', 'H', + 'H', 'H', 'E', 'T', 'C', 'E', 'T', 'C', 'S', 'B', 'T', 'C', 'C', + 'G', 'S', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'S', 'E', + 'E', 'S', 'H', 'H', 'H', 'H', ''] + + modes_simple_10 = ['Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', + 'Coil', 'Coil', 'Coil'] + + code_counts = {'G': np.array([7, 6, 3, 3, 3, 8, 3, 6, 0, 3]), + 'H': np.array([99, 96, 94, 97, 89, 101, 97, 96, 97, 96]), + 'I': np.array([0, 0, 0, 0, 5, 0, 0, 0, 0, 0]), + 'E': np.array([20, 28, 28, 23, 28, 29, 28, 31, 30, 30]), + 'B': np.array([5, 4, 3, 5, 4, 1, 3, 3, 1, 2]), + 'T': np.array([23, 28, 33, 28, 28, 19, 25, 27, 22, 26]), + 'S': np.array([13, 17, 19, 17, 17, 23, 19, 10, 25, 18]), + 'C': np.array([46, 34, 33, 40, 39, 32, 38, 40, 38, 38])} + + simple_counts = {'Helix': np.array([106, 102, 97, 100, 97, 109, 100, 102, 97, 99]), + 'Coil': np.array([82, 79, 85, 85, 84, 74, 82, 77, 85, 82]), + 'Strand': np.array([25, 32, 31, 28, 32, 30, 31, 34, 31, 32])} + analysis_class = ss.DSSPWrapper + + @pytest.fixture(scope='class') + def dssp(self, universe): + dssp = ss.DSSPWrapper(universe, executable='mkdssp', + select='backbone') + dssp.run() + return dssp + + @pytest.fixture() + def all_codes(self): + return np.load(ADK_MKDSSP) + + def test_non_protein(self, universe, dssp_all, all_codes): + super(TestDSSPWrapper, self).test_non_protein(universe, dssp_all) + + assert_equal(dssp_all.ss_codes[:, :self.n_prot], all_codes) + assert_equal(dssp_all.ss_mode[:self.n_prot:3], self.modes_3) + assert_equal(dssp_all.simple_mode[:10], self.modes_simple_10) + + for k in dssp_all.ss_counts.keys(): + assert_almost_equal(dssp_all.ss_counts[k], self.code_counts[k]) + for k in dssp_all.simple_counts.keys(): + assert_almost_equal(dssp_all.simple_counts[k], self.simple_counts[k]) + + +@pytest.mark.skipif(executable_not_found('stride'), + reason="Tests skipped because STRIDE not found") +class TestStrideWrapper(BaseTestSecondaryStructureWrapper): + + phi_psi_sasa_file = STRIDE_phi_psi_sasa + + modes_3 = ['C', 'E', 'C', 'T', 'H', 'H', 'H', 'H', 'C', 'E', 'C', 'H', 'H', + 'H', 'C', 'H', 'H', 'H', 'H', 'C', 'H', 'H', 'H', 'H', 'C', 'T', + 'T', 'E', 'T', 'T', 'H', 'H', 'H', 'C', 'C', 'E', 'E', 'C', 'H', + 'H', 'H', 'E', 'T', 'C', 'E', 'T', 'T', 'T', 'B', 'T', 'C', 'C', + 'G', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'E', + 'E', 'T', 'H', 'H', 'H', 'H', 'C'] + + modes_simple_10 = ['Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', + 'Coil', 'Coil', 'Coil'] + + code_counts = {'G': np.array([3, 3, 0, 3, 3, 0, 3, 3, 0, 3]), + 'H': np.array([101, 99, 100, 93, 102, 106, 102, 105, 101, 102]), + 'I': np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), + 'E': np.array([26, 29, 28, 28, 28, 32, 27, 30, 30, 30]), + 'B': np.array([3, 4, 3, 3, 6, 1, 3, 3, 1, 2]), + 'T': np.array([38, 41, 43, 48, 39, 38, 34, 34, 41, 38]), + 'S': np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0]), + 'C': np.array([43, 38, 40, 39, 36, 37, 45, 39, 41, 39])} + + simple_counts = {'Helix': np.array([104, 102, 100, 96, 105, 106, 105, 108, 101, 105]), + 'Coil': np.array([81, 79, 83, 87, 75, 75, 79, 73, 82, 77]), + 'Strand': np.array([29, 33, 31, 31, 34, 33, 30, 33, 31, 32])} + analysis_class = ss.STRIDEWrapper + + @pytest.fixture(scope='class') + def dssp(self, universe): + dssp = ss.STRIDEWrapper(universe, executable='stride', + select='backbone') + dssp.run() + return dssp + + @pytest.fixture() + def all_codes(self): + return np.load(ADK_STRIDE) diff --git a/testsuite/MDAnalysisTests/coordinates/test_pdb.py b/testsuite/MDAnalysisTests/coordinates/test_pdb.py index afde00786c3..c36a05ff538 100644 --- a/testsuite/MDAnalysisTests/coordinates/test_pdb.py +++ b/testsuite/MDAnalysisTests/coordinates/test_pdb.py @@ -46,6 +46,7 @@ class TestPDBReader(_SingleFrameReader): __test__ = True + def setUp(self): # can lead to race conditions when testing in parallel self.universe = mda.Universe(RefAdKSmall.filename) @@ -56,7 +57,8 @@ def setUp(self): def test_uses_PDBReader(self): from MDAnalysis.coordinates.PDB import PDBReader - assert isinstance(self.universe.trajectory, PDBReader), "failed to choose PDBReader" + assert isinstance(self.universe.trajectory, + PDBReader), "failed to choose PDBReader" def test_dimensions(self): assert_almost_equal( @@ -67,7 +69,8 @@ def test_dimensions(self): def test_ENT(self): from MDAnalysis.coordinates.PDB import PDBReader self.universe = mda.Universe(ENT) - assert isinstance(self.universe.trajectory, PDBReader), "failed to choose PDBReader" + assert isinstance(self.universe.trajectory, + PDBReader), "failed to choose PDBReader" class TestPDBMetadata(object): @@ -155,6 +158,7 @@ def test_REMARK(self, universe): class TestExtendedPDBReader(_SingleFrameReader): __test__ = True + def setUp(self): self.universe = mda.Universe(PDB_small, topology_format="XPDB", @@ -205,7 +209,6 @@ def u_no_resids(self): def u_no_names(self): return make_Universe(['resids', 'resnames'], trajectory=True) - def test_writer(self, universe, outfile): "Test writing from a single frame PDB file to a PDB file." "" universe.atoms.write(outfile) @@ -268,7 +271,7 @@ def test_write_single_frame_Writer(self, universe2, outfile): MDAnalysis.Writer (Issue 105)""" u = universe2 u.trajectory[50] - with mda.Writer(outfile) as W: + with mda.Writer(outfile) as W: W.write(u.select_atoms('all')) u2 = mda.Universe(outfile) assert_equal(u2.trajectory.n_frames, @@ -512,6 +515,7 @@ def helper(atoms, bonds): "the test reference; len(actual) is %d, len(desired) " "is %d" % (len(u._topology.bonds.values), len(desired))) + def test_conect_bonds_all(tmpdir): conect = mda.Universe(CONECT, guess_bonds=True) @@ -527,6 +531,7 @@ def test_conect_bonds_all(tmpdir): # assert_equal(len([b for b in conect.bonds if not b.is_guessed]), 1922) + def test_write_bonds_partial(tmpdir): u = mda.Universe(CONECT) # grab all atoms with bonds @@ -778,6 +783,7 @@ def test_serials(self, u): class TestPSF_CRDReader(_SingleFrameReader): __test__ = True + def setUp(self): self.universe = mda.Universe(PSF, CRD) self.prec = 5 # precision in CRD (at least we are writing %9.5f) @@ -793,7 +799,8 @@ def setUp(self): def test_uses_PDBReader(self): from MDAnalysis.coordinates.PDB import PDBReader - assert isinstance(self.universe.trajectory, PDBReader), "failed to choose PDBReader" + assert isinstance(self.universe.trajectory, + PDBReader), "failed to choose PDBReader" def test_write_occupancies(tmpdir): @@ -830,17 +837,26 @@ def test_atomname_alignment(self, writtenstuff): def test_atomtype_alignment(self, writtenstuff): result_line = ("ATOM 1 H5T GUA A 1 7.974 6.430 9.561" - " 1.00 0.00 RNAA H\n") + " 1.00 0.00 RNAA H \n") assert_equal(writtenstuff[3], result_line) @pytest.mark.parametrize('atom, refname', ((mda.coordinates.PDB.Pair('ASP', 'CA'), ' CA '), # Regular protein carbon alpha - (mda.coordinates.PDB.Pair('GLU', 'OE1'), ' OE1'), - (mda.coordinates.PDB.Pair('MSE', 'SE'), 'SE '), # Selenium like in 4D3L - (mda.coordinates.PDB.Pair('CA', 'CA'), 'CA '), # Calcium like in 4D3L - (mda.coordinates.PDB.Pair('HDD', 'FE'), 'FE '), # Iron from a heme like in 1GGE - (mda.coordinates.PDB.Pair('PLC', 'P'), ' P '), # Lipid phosphorus (1EIN) -)) + (mda.coordinates.PDB.Pair( + 'GLU', 'OE1'), ' OE1'), + # Selenium like in 4D3L + (mda.coordinates.PDB.Pair( + 'MSE', 'SE'), 'SE '), + # Calcium like in 4D3L + (mda.coordinates.PDB.Pair( + 'CA', 'CA'), 'CA '), + # Iron from a heme like in 1GGE + (mda.coordinates.PDB.Pair( + 'HDD', 'FE'), 'FE '), + # Lipid phosphorus (1EIN) + (mda.coordinates.PDB.Pair( + 'PLC', 'P'), ' P '), + )) def test_deduce_PDB_atom_name(atom, refname): # The Pair named tuple is used to mock atoms as we only need them to have a # ``resname`` and a ``name`` attribute. @@ -929,7 +945,8 @@ def test_write_pdb_zero_atoms(): def test_atom_not_match(tmpdir): # issue 1998 - outfile = str(tmpdir.mkdir("PDBReader").join('test_atom_not_match' + ".pdb")) + outfile = str(tmpdir.mkdir("PDBReader").join( + 'test_atom_not_match' + ".pdb")) u = mda.Universe(PSF, DCD) # select two groups of atoms protein = u.select_atoms("protein and name CA") diff --git a/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy new file mode 100644 index 00000000000..a06e781c2d2 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy new file mode 100644 index 00000000000..c1e18db1c63 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp.npy new file mode 100644 index 00000000000..61e1e9441cd Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp_simple.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp_simple.npy new file mode 100644 index 00000000000..00a2938e98e Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp_simple.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride.npy new file mode 100644 index 00000000000..142a2871c4f Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride_simple.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride_simple.npy new file mode 100644 index 00000000000..eef1b64cbf6 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride_simple.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/mkdssp_phi_psi_sasa.npy b/testsuite/MDAnalysisTests/data/dssp/mkdssp_phi_psi_sasa.npy new file mode 100755 index 00000000000..574bdbfad16 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/mkdssp_phi_psi_sasa.npy differ diff --git a/testsuite/MDAnalysisTests/data/dssp/stride_phi_psi_sasa.npy b/testsuite/MDAnalysisTests/data/dssp/stride_phi_psi_sasa.npy new file mode 100755 index 00000000000..01fec37eb2b Binary files /dev/null and b/testsuite/MDAnalysisTests/data/dssp/stride_phi_psi_sasa.npy differ diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index f6873c7a7f0..92cec23ed56 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -186,6 +186,12 @@ "GMX_DIR", # GROMACS directory "GMX_TOP_BAD", # file with an #include that doesn't exist "ITP_no_endif", # file missing an #endif + "ADK_DSSP", # DSSP of ADK, assigned by mdtraj.compute_dssp(simplified=False) (protein only) + "ADK_DSSP_SIMPLE", # DSSP of ADK, assigned by mdtraj.compute_dssp(simplified=True) (protein only) + "ADK_MKDSSP", # DSSP of ADK, assigned by mkdssp (protein only) + "ADK_STRIDE", # DSSP of ADK, assigned by stride (protein only) + "MKDSSP_phi_psi_sasa", # phi/psi angles and solvent-accessible surface area (mkdssp) + "STRIDE_phi_psi_sasa", # phi/psi angles and solvent-accessible surface area (stride) ] from pkg_resources import resource_filename @@ -509,5 +515,13 @@ ITP_no_endif = resource_filename(__name__, 'data/no_endif_spc.itp') NAMDBIN = resource_filename(__name__, 'data/adk_open.coor') + +ADK_DSSP = resource_filename(__name__, 'data/dssp/adk_oplsaa_dssp.npy') +ADK_DSSP_SIMPLE = resource_filename(__name__, 'data/dssp/adk_oplsaa_dssp_simple.npy') +ADK_MKDSSP = resource_filename(__name__, 'data/dssp/adk_oplsaa_mkdssp.npy') +ADK_STRIDE = resource_filename(__name__, 'data/dssp/adk_oplsaa_stride.npy') +MKDSSP_phi_psi_sasa = resource_filename(__name__, 'data/dssp/mkdssp_phi_psi_sasa.npy') +STRIDE_phi_psi_sasa = resource_filename(__name__, 'data/dssp/stride_phi_psi_sasa.npy') + # This should be the last line: clean up namespace del resource_filename