From 2e3b68f66af4325af8cbdf13dcec0b378cc243f2 Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Sat, 7 Mar 2020 22:18:40 +1100 Subject: [PATCH 1/6] added secondary_structure module --- .../analysis/secondary_structure.py | 451 ++++++++++++++++++ package/MDAnalysis/coordinates/PDB.py | 2 +- 2 files changed, 452 insertions(+), 1 deletion(-) create mode 100644 package/MDAnalysis/analysis/secondary_structure.py diff --git a/package/MDAnalysis/analysis/secondary_structure.py b/package/MDAnalysis/analysis/secondary_structure.py new file mode 100644 index 00000000000..c2458a80710 --- /dev/null +++ b/package/MDAnalysis/analysis/secondary_structure.py @@ -0,0 +1,451 @@ +# -*- 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 + +This module contains classes for computing secondary structure. MDAnalysis provides +a wrapper for the external DSSP_ and STRIDE_ programs to be run on an MD trajectory. + +Both DSSP_ and STRIDE_ need to be installed separately. + +.. _DSSP: https://swift.cmbi.umcn.nl/gv/dssp/ +.. _STRIDE: http://webclu.bio.wzw.tum.de/stride/ + + +Classes and functions +--------------------- + +.. autoclass:: DSSP +.. autoclass:: STRIDE + +""" + +import errno +import subprocess +import os +import numpy as np +import pandas as pd + +from ..core.topologyattrs import ResidueAttr +from ..lib import util +from .base import AnalysisBase + +class SecondaryStructure(ResidueAttr): + """Single letter code for secondary structure""" + attrname = 'secondary_structures' + singular = 'secondary_structure' + dtype = '2s}\n"), + "{tempFactor:6.2f} {segID:<4s}{element:>2s} \n"), 'REMARK': "REMARK {0}\n", 'COMPND': "COMPND {0}\n", 'HEADER': "HEADER {0}\n", From c5bf9bdababf0045748c86b5254b9492b3ec55fa Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Mon, 9 Mar 2020 18:54:22 +1100 Subject: [PATCH 2/6] added mdtraj dssp --- .travis.yml | 2 +- .../analysis/secondary_structure/__init__.py | 27 ++ .../analysis/secondary_structure/base.py | 234 ++++++++++++++++ .../analysis/secondary_structure/dssp.py | 251 ++++++++++++++++++ .../wrapper.py} | 232 +++------------- .../source/documentation_pages/references.rst | 32 +++ .../analysis/test_secondary_structure.py | 194 ++++++++++++++ .../MDAnalysisTests/data/adk_oplsaa_dssp.npy | Bin 0 -> 17248 bytes .../data/adk_oplsaa_dssp_simple.npy | Bin 0 -> 17248 bytes testsuite/MDAnalysisTests/datafiles.py | 6 + 10 files changed, 783 insertions(+), 195 deletions(-) create mode 100644 package/MDAnalysis/analysis/secondary_structure/__init__.py create mode 100644 package/MDAnalysis/analysis/secondary_structure/base.py create mode 100644 package/MDAnalysis/analysis/secondary_structure/dssp.py rename package/MDAnalysis/analysis/{secondary_structure.py => secondary_structure/wrapper.py} (58%) create mode 100644 testsuite/MDAnalysisTests/analysis/test_secondary_structure.py create mode 100644 testsuite/MDAnalysisTests/data/adk_oplsaa_dssp.npy create mode 100644 testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy diff --git a/.travis.yml b/.travis.yml index 2818163d465..e39f7cfe005 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ 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" + - CONDA_DEPENDENCIES="${CONDA_MIN_DEPENDENCIES} seaborn>=0.7.0 clustalw=2.1 netcdf4 scikit-learn joblib>=0.12 mdtraj" - CONDA_CHANNELS='biobuilds conda-forge' - CONDA_CHANNEL_PRIORITY=True - PIP_DEPENDENCIES="duecredit parmed" 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..a84562e3239 --- /dev/null +++ b/package/MDAnalysis/analysis/secondary_structure/base.py @@ -0,0 +1,234 @@ +# -*- 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 + + +""" + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +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='4Aq5QbsXaux0jk_!@1Au14m1qcaJP$7#LQ6PeBM1!~rZkQnaG_N$;=j@-z zPMTM<*35k0%sONP|9tb}>u-N}=i<+czph^1-v4xW_0i?k$1k2=U0+`P{OazHyPMyB ze0BHo_P_JbZ+^YMeSLoa%gyh%ujBWgeQAK6!W-e&wO_lF#;2S5~L}#8Vx;-izHMy(y+QQQg-4KlXb2>mJPC z_nl(-Q=g;vKp!Wc)44kPK=Zgq_Y>F7tA9E9uKrrx*6GVf@Ar2u5OYyW{=STaegG+wh zn^&EEVDFdaP`8vf-SbiMThha+Z_njVtbXhFj`^r=`&ze{oZ?CrPu{$A(TT&4T6{aPokPWy=EqdI$WUb@HDHKsRFpS||~So_$c zJ(wH&zEi9Zw%<|n^I1+ld)_+!aql_4diekEORTT@7V}wd-FNlvVQb&1p8x7Lp6a)r z(?|0@_3xE=u=^C7*YE07s>f5`fuAMT2cISVo{QbL-=yz*p}utAdgt_dFLsZub6#Ds zKIr_^{Xf;b&TCJ7*QHn;+Mk#X?zJ!d=-`sy-mB9O&4GLQp8C+`2kTcXAE&eD*$x!(`_JD2sjRNvp7Sl!xvKb1P~V_xhYr}uU)Up+cMY9IPf@7*`aYHN_S*kr?bH34+a8zR+r7F{pQGm2zvOF)t+Nj_ zk38pl&!aD$a+*UO?&Yicd-ZJ}zOlNETfa7@*MIKC-d)}j`JCR#xxVsC_iN{BU42(i zOqb77sUEApHs*WDhrEev>wQ1WWj|$g)irO5^;ajJ^8G$b`D^$6RO;;8y!M4l@9kcF zsn1dS(7)tsiLJAbQr=N>(4$jMbEw0;{qeDPU;X$-dK2ZT-;4TUf5*I-@T|Wp4+!|TlfFi>%CL=G`;V1 z-g>cp9<>Mk%gN_-uHHVdb2y*4bxuCY$#;Eh)mLAA=SR^T)opyVJ|Es)jQqGq=hHhm zZ@xNozjaQ2OFAW8K2K%)Z1t_4@5u*Ud8K=t`m|1b>i0@M+9O@@xt#R--&JD$?eA^n Vv}gOC-nr-U?G2qDwLkvp{0kE2f!hE8 literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy b/testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy new file mode 100644 index 0000000000000000000000000000000000000000..c1e18db1c63fc390b4e89cc23c3769e581c7b8c8 GIT binary patch literal 17248 zcmds!y=qiZ6oq5!Q%pAsDNL{uJ5g*Dn-o?um=Oy}WFj`=Q~1J0%-1YftbNWsXBh5n z7W-%YoPGbizk2uKd-!qK-F&(|T<`z&-|W8K z?eF_{Uw7a4_w|#DXXl&Ci|0?zH@`Q1-{V(j|3C2cL+y~gD!zTtn}sfi=4bC%IcE8w z%X<>CKN_8{^Ktw1^?y+_Pw0Gm`^k?vIjaZ{?9B%c9nj@q-ioT{`0Tm4>yP5~>Tt)k zqPshN?v{{!hqLdezA3NHmA|y^eXsM?ywtki>#Mw3y6Qu%ey(PJG&(=+J74zrFAmvD z=zMzn$&b0q?fUS*-hA-T0o9j&FIM%~b3ShF`c8*_mB*bF(b4|eYToIPJ0^5|`m^qa zzGBMDTzRRU^$v8MKY111K6<^Vy1hD8p8C1Es@J;L>gV?L|7h)a+HcM8_Ek64nm5&| z=W=LYd9GG}w^!{~_1LTGvUhAAzM$UuST+9OEMN2$xqB6rul1?s*PN+N?|`57+#J55 z`wzP3bhWR%R4adUo&1@P_C*hUMdpgSNA3O-efAQn-_knyF(+pg;ektY!bi{LVBU(V z=lJZox$BSO_3CiHwW7N_eeRZ!eTTE|=loJ$=8o_Bv)-Y)>Q7!p-J@2$)*2nfR-Nj; z`nbBP*Sgp0=l1pgXzh5~Z_V%aRX5d|H`S`=a%f+9u2z4ySM68z*mJsW?)pxLelHI` zsCQml%{vvjUlo;~>QQ%eKGodpIUVhf%I$R2kD8o{%t5Ey+D8v;<)?jgTQTV_ook+| z9z6X&TCRa7Qe#|v*s^NjvJZTR-=g+)~tvc*i^^djZ?uNgcucgmB6Y~EYCjP8@ zC100U`${T|Bu#HJ(sJsU-e^8^QIae*qYDnqu1&W5A5Z`L#K+fdFftmzLtLa58OXt z(wTKH=bQ2}cYN1B_B+TyUPb2M=zOzu$%%d|!XJgL{%&9YkJgT-{iW|xy=tv_Q?0tA zb31?aO?9gd`$_lI=1cE?>U?)^{#O;94*gkoKwq(%n?2?C)=S^%e95V(`%C@d&(cL- zQT<%4yiv9DYwhd*(YmU)^j%dS-P63OR{ql5@X>Yt%&XX{!+uqN^q%wQK8pBje{D7I QbjaNjIzIhbcS9fi3xC*Q>;M1& literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index f6873c7a7f0..97c892048c8 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -186,6 +186,8 @@ "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) ] from pkg_resources import resource_filename @@ -509,5 +511,9 @@ 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/adk_oplsaa_dssp.npy') +ADK_DSSP_SIMPLE = resource_filename(__name__, 'data/adk_oplsaa_dssp_simple.npy') + # This should be the last line: clean up namespace del resource_filename From a341b1468140d12b919eecc777438ab60242ce69 Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Tue, 10 Mar 2020 17:10:47 +1100 Subject: [PATCH 3/6] added tests for wrapper classes --- .travis.yml | 18 + maintainer/install_dssp.sh | 126 + maintainer/install_stride.sh | 94 + .../analysis/secondary_structure/base.py | 13 +- .../analysis/secondary_structure/wrapper.py | 66 +- package/MDAnalysis/converters/__init__.py | 0 package/MDAnalysis/converters/base.py | 2415 +++++++++++++++++ .../analysis/test_secondary_structure.py | 294 +- .../data/{ => dssp}/adk_oplsaa_dssp.npy | Bin .../{ => dssp}/adk_oplsaa_dssp_simple.npy | Bin .../data/dssp/adk_oplsaa_mkdssp.npy | Bin 0 -> 8688 bytes .../data/dssp/adk_oplsaa_stride.npy | Bin 0 -> 8688 bytes testsuite/MDAnalysisTests/datafiles.py | 8 +- 13 files changed, 2892 insertions(+), 142 deletions(-) create mode 100755 maintainer/install_dssp.sh create mode 100755 maintainer/install_stride.sh create mode 100644 package/MDAnalysis/converters/__init__.py create mode 100644 package/MDAnalysis/converters/base.py rename testsuite/MDAnalysisTests/data/{ => dssp}/adk_oplsaa_dssp.npy (100%) rename testsuite/MDAnalysisTests/data/{ => dssp}/adk_oplsaa_dssp_simple.npy (100%) create mode 100644 testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp.npy create mode 100644 testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride.npy diff --git a/.travis.yml b/.travis.yml index e39f7cfe005..53e95cfbb2a 100644 --- a/.travis.yml +++ b/.travis.yml @@ -34,6 +34,8 @@ env: - 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..36ad8585398 --- /dev/null +++ b/maintainer/install_dssp.sh @@ -0,0 +1,126 @@ +#!/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 + ;; + 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..99b35c62299 --- /dev/null +++ b/maintainer/install_stride.sh @@ -0,0 +1,94 @@ +#!/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/ +# +# +# +# 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/base.py b/package/MDAnalysis/analysis/secondary_structure/base.py index a84562e3239..09817aa94b7 100644 --- a/package/MDAnalysis/analysis/secondary_structure/base.py +++ b/package/MDAnalysis/analysis/secondary_structure/base.py @@ -34,8 +34,6 @@ """ import numpy as np -import pandas as pd -import matplotlib.pyplot as plt from ...core.topologyattrs import ResidueAttr from ..base import AnalysisBase @@ -129,7 +127,8 @@ def __init__(self, universe, select='backbone', verbose=False, super(SecondaryStructureBase, self).__init__(universe.universe.trajectory, verbose=verbose) self._universe = universe.universe - self.residues = universe.select_atoms(select).residues + self.atomgroup = universe.select_atoms(select) + self.residues = self.atomgroup.residues self.n_residues = len(self.residues) self._add_topology_attr = add_topology_attr @@ -210,6 +209,14 @@ def plot_content(self, ax=None, simple=False, figsize=(8, 6), ax: :class: `matplotlib.axes.Axes` """ + try: + import pandas as pd + except ImportError: + raise ImportError('pandas is required for this function.') + try: + import matplotlib.pyplot as plt + except ImportError: + raise ImportError('matplotlib is required for this function.') if not simple: data = self.ss_counts diff --git a/package/MDAnalysis/analysis/secondary_structure/wrapper.py b/package/MDAnalysis/analysis/secondary_structure/wrapper.py index 3b38bc52f42..9109285434b 100644 --- a/package/MDAnalysis/analysis/secondary_structure/wrapper.py +++ b/package/MDAnalysis/analysis/secondary_structure/wrapper.py @@ -51,6 +51,8 @@ import errno import tempfile import subprocess +import warnings +import numpy as np from ...due import due, Doi from ...lib import util @@ -71,6 +73,7 @@ path="MDAnalysis.analysis.secondary_structure.wrapper", cite_module=True) + class SecondaryStructureWrapper(SecondaryStructureBase): """Base class for secondary structure analysis that wraps an external program. @@ -97,15 +100,20 @@ def __init__(self, executable, universe, select='protein', verbose=False, select=select, verbose=verbose, add_topology_attr=add_topology_attr) - + def _execute(self): """Run the wrapped program.""" + # ignore PDB warnings + warnings.filterwarnings("ignore", + message="Found no information for attr") + fd, pdbfile = tempfile.mkstemp(suffix='.pdb') os.close(fd) try: - self.residues.atoms.write(pdbfile) + self.atomgroup.write(pdbfile) cmd_args = [self.exe] + self.cmd.format(pdb=pdbfile).split() - proc = subprocess.Popen(cmd_args, stdout=subprocess.PIPE, + print(' '.join(cmd_args)) + proc = subprocess.Popen(cmd_args, stdout=subprocess.PIPE, stderr=subprocess.PIPE) stdout, stderr = proc.communicate() finally: @@ -113,7 +121,7 @@ def _execute(self): os.unlink(pdbfile) except OSError: pass - return output + return stdout.decode('utf-8') def _prepare(self): super(SecondaryStructureWrapper, self)._prepare() @@ -121,11 +129,11 @@ def _prepare(self): self.phi = np.zeros((nf, self.n_residues), dtype=float) self.psi = np.zeros((nf, self.n_residues), dtype=float) self.sasa = np.zeros((nf, self.n_residues), dtype=float) - + def _compute_dssp(self): output = self._execute() self._process_output(output) - + class DSSPWrapper(SecondaryStructureWrapper): """Runs :program:`mkdssp` on a trajectory. @@ -186,35 +194,40 @@ class DSSPWrapper(SecondaryStructureWrapper): """ exe_name = 'mkdssp' - cmd = '-i {pdb} -o stdout' + cmd = '-i {pdb}' columns = { - 'ss': 17, + 'ss': 16, 'sasa': (35, 39), - 'phi': (105, 110), + 'phi': (104, 109), 'psi': (111, 116), } def __init__(self, universe, select='protein', executable='mkdssp', verbose=False, add_topology_attr=False): - super(DSSPWrapper, self).__init__(executable, universe, - select=select, verbose=verbose, - add_topology_attr=add_topology_attr) - + super(DSSPWrapper, self).__init__(executable, universe, + select=select, verbose=verbose, + add_topology_attr=add_topology_attr) def _process_output(self, output): frame = self._frame_index per_res = output.split(' # RESIDUE AA STRUCTURE BP1 BP2 ACC')[1] - lines = per_res.split()[1:] - for i, line in enumerate(lines): + i = 0 + for line in per_res.split('\n')[1:]: + try: + if line[13] == '!': + continue + except IndexError: + continue ss = line[self.columns['ss']] - if not ss: + + if ss == ' ': ss = 'C' self.ss_codes[frame][i] = ss for kw in ('phi', 'psi', 'sasa'): _i, _j = self.columns[kw] getattr(self, kw)[frame][i] = line[_i:_j] - + i += 1 class STRIDEWrapper(SecondaryStructureWrapper): @@ -241,7 +254,7 @@ class STRIDEWrapper(SecondaryStructureWrapper): attribute ``secondary_structure`` to your residues. verbose: bool, optional Turn on more logging and debugging. - + Attributes ---------- residues: :class:`~MDAnalysis.core.groups.ResidueGroup` @@ -279,17 +292,18 @@ class STRIDEWrapper(SecondaryStructureWrapper): def __init__(self, universe, select='protein', executable='stride', verbose=False, add_topology_attr=False): - super(STRIDEWrapper, self).__init__(executable, universe, - select=select, verbose=verbose, - add_topology_attr=add_topology_attr) + super(STRIDEWrapper, self).__init__(executable, universe, + select=select, verbose=verbose, + add_topology_attr=add_topology_attr) def _process_output(self, output): lines = output.split('\nASG ')[1:] lines[-1] = lines[-1].split('\n')[0] frame = self._frame_index for i, line in enumerate(lines): - resname, chain, resid, resnum, ss, ssname, phi, psi, area = line.split() - self.ss_codes[frame][i] = ss - self.phi[frame][ss] = phi - self.psi[frame][ss] = psi - self.sasa[frame][ss] = area \ No newline at end of file + # resname chain resid resnum ss ssname phi psi area ~~~~ + fields = line.split() + self.ss_codes[frame][i] = fields[4] + self.phi[frame][i] = fields[6] + self.psi[frame][i] = fields[7] + self.sasa[frame][i] = fields[8] diff --git a/package/MDAnalysis/converters/__init__.py b/package/MDAnalysis/converters/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/package/MDAnalysis/converters/base.py b/package/MDAnalysis/converters/base.py new file mode 100644 index 00000000000..b7e6e79ddee --- /dev/null +++ b/package/MDAnalysis/converters/base.py @@ -0,0 +1,2415 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# + + +"""\ +Base classes --- :mod:`MDAnalysis.coordinates.base` +=================================================== + +Derive other Timestep, FrameIterator, Reader and Writer classes from the classes +in this module. The derived classes must follow the :ref:`Trajectory API` +in :mod:`MDAnalysis.coordinates.__init__`. + +Timestep +-------- + +A :class:`Timestep` holds information for the current time frame in +the trajectory. It is one of the central data structures in +MDAnalysis. + +.. class:: Timestep + + .. automethod:: __init__ + .. automethod:: from_coordinates + .. automethod:: from_timestep + .. automethod:: _init_unitcell + .. autoattribute:: n_atoms + .. attribute::`frame` + + frame number (0-based) + + .. versionchanged:: 0.11.0 + Frames now 0-based; was 1-based + + .. autoattribute:: time + .. autoattribute:: dt + .. autoattribute:: positions + .. autoattribute:: velocities + .. autoattribute:: forces + .. autoattribute:: has_positions + .. autoattribute:: has_velocities + .. autoattribute:: has_forces + .. attribute:: _pos + + :class:`numpy.ndarray` of dtype :class:`~numpy.float32` of shape + (*n_atoms*, 3) and internal FORTRAN order, holding the raw + cartesian coordinates (in MDAnalysis units, i.e. Å). + + .. Note:: + + Normally one does not directly access :attr:`_pos` but uses + the :meth:`~MDAnalysis.core.groups.AtomGroup.coordinates` + method of an :class:`~MDAnalysis.core.groups.AtomGroup` but + sometimes it can be faster to directly use the raw + coordinates. Any changes to this array are immediately + reflected in atom positions. If the frame is written to a new + trajectory then the coordinates are changed. If a new + trajectory frame is loaded, then *all* contents of + :attr:`_pos` are overwritten. + + .. attribute:: _velocities + + :class:`numpy.ndarray` of dtype :class:`~numpy.float32`. of shape + (*n_atoms*, 3), holding the raw velocities (in MDAnalysis + units, i.e. typically Å/ps). + + .. Note:: + + Normally velocities are accessed through the + :attr:`velocities` or the + :meth:`~MDAnalysis.core.groups.AtomGroup.velocities` + method of an :class:`~MDAnalysis.core.groups.AtomGroup` + + :attr:`~Timestep._velocities` only exists if the :attr:`has_velocities` + flag is True + + .. versionadded:: 0.7.5 + + .. attribute:: _forces + + :class:`numpy.ndarray` of dtype :class:`~numpy.float32`. of shape + (*n_atoms*, 3), holding the forces + + :attr:`~Timestep._forces` only exists if :attr:`has_forces` + is True + + .. versionadded:: 0.11.0 + Added as optional to :class:`Timestep` + + .. autoattribute:: dimensions + .. autoattribute:: triclinic_dimensions + .. autoattribute:: volume + .. automethod:: __getitem__ + .. automethod:: __eq__ + .. automethod:: __iter__ + .. automethod:: copy + .. automethod:: copy_slice + + +FrameIterators +-------------- + +Iterator classes used by the by the :class:`ProtoReader`. + +.. autoclass:: FrameIteratorBase + +.. autoclass:: FrameIteratorSliced + +.. autoclass:: FrameIteratorAll + +.. autoclass:: FrameIteratorIndices + + +Readers +------- + +Readers know how to take trajectory data in a given format and present it in a +common API to the user in MDAnalysis. There are two types of readers: + +1. Readers for *multi frame trajectories*, i.e., file formats that typically + contain many frames. These readers are typically derived from + :class:`ReaderBase`. + +2. Readers for *single frame formats*: These file formats only contain a single + coordinate set. These readers are derived from + :class`:SingleFrameReaderBase`. + +The underlying low-level readers handle closing of files in different +ways. Typically, the MDAnalysis readers try to ensure that files are always +closed when a reader instance is garbage collected, which relies on +implementing a :meth:`~ReaderBase.__del__` method. However, in some cases, this +is not necessary (for instance, for the single frame formats) and then such a +method can lead to undesirable side effects (such as memory leaks). In this +case, :class:`ProtoReader` should be used. + + +.. autoclass:: ReaderBase + :members: + :inherited-members: + +.. autoclass:: SingleFrameReaderBase + :members: + :inherited-members: + +.. autoclass:: ProtoReader + :members: + + + +Writers +------- + +Writers know how to write information in a :class:`Timestep` to a trajectory +file. + +.. autoclass:: WriterBase + :members: + :inherited-members: + + +Converters +---------- + +Converters output information to other libraries. + +.. autoclass:: ConverterBase + :members: + :inherited-members: + + +Helper classes +-------------- + +The following classes contain basic functionality that all readers and +writers share. + +.. autoclass:: IOBase + :members: + +""" +from __future__ import absolute_import +import six +from six.moves import range + +import numpy as np +import numbers +import copy +import warnings +import weakref + +from . import core +from .. import NoDataError +from .. import ( + _READERS, + _SINGLEFRAME_WRITERS, + _MULTIFRAME_WRITERS, + _CONVERTERS +) +from .. import units +from ..auxiliary.base import AuxReader +from ..auxiliary.core import auxreader +from ..core import flags +from ..lib.util import asiterable, Namespace + + +class Timestep(object): + """Timestep data for one frame + + :Methods: + + ``ts = Timestep(n_atoms)`` + + create a timestep object with space for n_atoms + + .. versionchanged:: 0.11.0 + Added :meth:`from_timestep` and :meth:`from_coordinates` constructor + methods. + :class:`Timestep` init now only accepts integer creation. + :attr:`n_atoms` now a read only property. + :attr:`frame` now 0-based instead of 1-based. + Attributes `status` and `step` removed. + """ + order = 'F' + + def __init__(self, n_atoms, **kwargs): + """Create a Timestep, representing a frame of a trajectory + + Parameters + ---------- + n_atoms : int + The total number of atoms this Timestep describes + positions : bool, optional + Whether this Timestep has position information [``True``] + velocities : bool (optional) + Whether this Timestep has velocity information [``False``] + forces : bool (optional) + Whether this Timestep has force information [``False``] + reader : Reader (optional) + A weak reference to the owning Reader. Used for + when attributes require trajectory manipulation (e.g. dt) + dt : float (optional) + The time difference between frames (ps). If :attr:`time` + is set, then `dt` will be ignored. + time_offset : float (optional) + The starting time from which to calculate time (in ps) + + + .. versionchanged:: 0.11.0 + Added keywords for `positions`, `velocities` and `forces`. + Can add and remove position/velocity/force information by using + the ``has_*`` attribute. + """ + # readers call Reader._read_next_timestep() on init, incrementing + # self.frame to 0 + self.frame = -1 + self._n_atoms = n_atoms + + self.data = {} + + for att in ('dt', 'time_offset'): + try: + self.data[att] = kwargs[att] + except KeyError: + pass + try: + # do I have a hook back to the Reader? + self._reader = weakref.ref(kwargs['reader']) + except KeyError: + pass + + # Stupid hack to make it allocate first time round + # ie we have to go from not having, to having positions + # to make the Timestep allocate + self._has_positions = False + self._has_velocities = False + self._has_forces = False + + # These will allocate the arrays if the has flag + # gets set to True + self.has_positions = kwargs.get('positions', True) + self.has_velocities = kwargs.get('velocities', False) + self.has_forces = kwargs.get('forces', False) + + self._unitcell = self._init_unitcell() + + # set up aux namespace for adding auxiliary data + self.aux = Namespace() + + + @classmethod + def from_timestep(cls, other, **kwargs): + """Create a copy of another Timestep, in the format of this Timestep + + .. versionadded:: 0.11.0 + """ + ts = cls(other.n_atoms, + positions=other.has_positions, + velocities=other.has_velocities, + forces=other.has_forces, + **kwargs) + ts.frame = other.frame + ts.dimensions = other.dimensions + try: + ts.positions = other.positions.copy(order=cls.order) + except NoDataError: + pass + try: + ts.velocities = other.velocities.copy(order=cls.order) + except NoDataError: + pass + try: + ts.forces = other.forces.copy(order=cls.order) + except NoDataError: + pass + + # Optional attributes that don't live in .data + # should probably iron out these last kinks + for att in ('_frame',): + try: + setattr(ts, att, getattr(other, att)) + except AttributeError: + pass + + if hasattr(ts, '_reader'): + other._reader = weakref.ref(ts._reader()) + + ts.data = copy.deepcopy(other.data) + + return ts + + @classmethod + def from_coordinates(cls, + positions=None, + velocities=None, + forces=None, + **kwargs): + """Create an instance of this Timestep, from coordinate data + + Can pass position, velocity and force data to form a Timestep. + + .. versionadded:: 0.11.0 + """ + has_positions = positions is not None + has_velocities = velocities is not None + has_forces = forces is not None + + lens = [len(a) for a in [positions, velocities, forces] + if a is not None] + if not lens: + raise ValueError("Must specify at least one set of data") + n_atoms = max(lens) + # Check arrays are matched length? + if not all(val == n_atoms for val in lens): + raise ValueError("Lengths of input data mismatched") + + ts = cls(n_atoms, + positions=has_positions, + velocities=has_velocities, + forces=has_forces, + **kwargs) + if has_positions: + ts.positions = positions + if has_velocities: + ts.velocities = velocities + if has_forces: + ts.forces = forces + + return ts + + def _init_unitcell(self): + """Create custom datastructure for :attr:`_unitcell`.""" + # override for other Timesteps + return np.zeros((6), np.float32) + + def __eq__(self, other): + """Compare with another Timestep + + .. versionadded:: 0.11.0 + """ + if not isinstance(other, Timestep): + return False + + if not self.frame == other.frame: + return False + + if not self.n_atoms == other.n_atoms: + return False + + if not self.has_positions == other.has_positions: + return False + if self.has_positions: + if not (self.positions == other.positions).all(): + return False + + if not self.has_velocities == other.has_velocities: + return False + if self.has_velocities: + if not (self.velocities == other.velocities).all(): + return False + + if not self.has_forces == other.has_forces: + return False + if self.has_forces: + if not (self.forces == other.forces).all(): + return False + + return True + + def __ne__(self, other): + return not self == other + + def __getitem__(self, atoms): + """Get a selection of coordinates + + ``ts[i]`` + + return coordinates for the i'th atom (0-based) + + ``ts[start:stop:skip]`` + + return an array of coordinates, where start, stop and skip + correspond to atom indices, + :attr:`MDAnalysis.core.groups.Atom.index` (0-based) + """ + if isinstance(atoms, numbers.Integral): + return self._pos[atoms] + elif isinstance(atoms, (slice, np.ndarray)): + return self._pos[atoms] + else: + raise TypeError + + def __len__(self): + return self.n_atoms + + def __iter__(self): + """Iterate over coordinates + + ``for x in ts`` + + iterate of the coordinates, atom by atom + """ + for i in range(self.n_atoms): + yield self[i] + + def __repr__(self): + desc = "< Timestep {0}".format(self.frame) + try: + tail = " with unit cell dimensions {0} >".format(self.dimensions) + except NotImplementedError: + tail = " >" + return desc + tail + + def copy(self): + """Make an independent ("deep") copy of the whole :class:`Timestep`.""" + return self.__deepcopy__() + + def __deepcopy__(self): + return self.from_timestep(self) + + def copy_slice(self, sel): + """Make a new `Timestep` containing a subset of the original `Timestep`. + + Parameters + ---------- + sel : array_like or slice + The underlying position, velocity, and force arrays are sliced + using a :class:`list`, :class:`slice`, or any array-like. + + Returns + ------- + :class:`Timestep` + A `Timestep` object of the same type containing all header + information and all atom information relevant to the selection. + + Note + ---- + The selection must be a 0 based :class:`slice` or array of the atom indices + in this :class:`Timestep` + + Example + ------- + Using a Python :class:`slice` object:: + + new_ts = ts.copy_slice(slice(start, stop, step)) + + Using a list of indices:: + + new_ts = ts.copy_slice([0, 2, 10, 20, 23]) + + + .. versionadded:: 0.8 + .. versionchanged:: 0.11.0 + Reworked to follow new Timestep API. Now will strictly only + copy official attributes of the Timestep. + + """ + # Detect the size of the Timestep by doing a dummy slice + try: + pos = self.positions[sel, :] + except NoDataError: + # It's cool if there's no Data, we'll live + pos = None + except Exception: + six.raise_from( + TypeError( + "Selection type must be compatible with slicing" + " the coordinates" + ), + None) + try: + vel = self.velocities[sel, :] + except NoDataError: + vel = None + except Exception: + six.raise_from( + TypeError("Selection type must be compatible with slicing" + " the coordinates"), + None) + try: + force = self.forces[sel, :] + except NoDataError: + force = None + except Exception: + six.raise_from( + TypeError( + "Selection type must be compatible with slicing" + " the coordinates" + ), + None) + + new_TS = self.__class__.from_coordinates( + positions=pos, + velocities=vel, + forces=force) + + new_TS._unitcell = self._unitcell.copy() + + new_TS.frame = self.frame + + for att in ('_frame',): + try: + setattr(new_TS, att, getattr(self, att)) + except AttributeError: + pass + + if hasattr(self, '_reader'): + new_TS._reader = weakref.ref(self._reader()) + + new_TS.data = copy.deepcopy(self.data) + + return new_TS + + @property + def n_atoms(self): + """A read only view of the number of atoms this Timestep has + + .. versionchanged:: 0.11.0 + Changed to read only property + """ + # In future could do some magic here to make setting n_atoms + # resize the coordinate arrays, but + # - not sure if that is ever useful + # - not sure how to manage existing data upon extension + return self._n_atoms + + @property + def has_positions(self): + """A boolean of whether this Timestep has position data + + This can be changed to ``True`` or ``False`` to allocate space for + or remove the data. + + .. versionadded:: 0.11.0 + """ + return self._has_positions + + @has_positions.setter + def has_positions(self, val): + if val and not self._has_positions: + # Setting this will always reallocate position data + # ie + # True -> 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/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py b/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py index a6c7f551ed9..3ff5c47e206 100644 --- a/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py +++ b/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py @@ -30,140 +30,73 @@ import MDAnalysis as mda from MDAnalysis.analysis import secondary_structure as ss -from MDAnalysisTests.datafiles import TPR, XTC, ADK_DSSP, ADK_DSSP_SIMPLE - -XTC_modes = ['C', 'E', 'E', 'E', 'E', 'E', 'C', 'C', 'T', 'T', 'S', 'C', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'C', - 'C', 'E', 'E', 'E', 'T', 'T', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'T', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'T', 'S', 'G', 'G', 'G', 'S', - 'S', 'C', 'E', 'E', 'E', 'E', 'S', 'C', 'C', 'C', 'S', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'T', 'C', 'C', 'C', 'S', - 'E', 'E', 'E', 'E', 'E', 'E', 'C', 'C', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'T', 'E', 'E', 'E', 'S', 'T', 'T', 'T', 'S', - 'S', 'E', 'E', 'E', 'T', 'T', 'T', 'B', 'C', 'S', 'S', 'S', 'T', - 'T', 'B', 'C', 'T', 'T', 'T', 'C', 'C', 'B', 'C', 'S', 'S', 'S', - 'G', 'G', 'G', 'S', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'S', 'C', 'C', 'E', 'E', 'E', - 'E', 'E', 'C', 'S', 'S', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'C', 'C'] - -XTC_modes_simple = ['Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', - 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Coil', 'Coil', 'Coil', 'Strand', 'Strand', - 'Strand', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Coil', 'Strand', 'Strand', 'Strand', - 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Coil', 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Strand', 'Strand', 'Strand', 'Coil', - 'Coil', 'Coil', 'Strand', 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', - 'Coil', 'Strand', 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', - 'Strand', 'Coil', 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', - 'Helix', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Coil', 'Coil', 'Coil', 'Coil', 'Strand', - 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', 'Coil', - 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil'] - -XTC_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])} +from MDAnalysisTests import executable_not_found +from MDAnalysisTests.datafiles import (TPR, XTC, + ADK_DSSP, ADK_DSSP_SIMPLE, + ADK_MKDSSP, + ADK_STRIDE, + ) -XTC_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])} +def mdtraj_found(): + try: + import mdtraj + except ImportError: + return False + else: + return True -class TestDSSP(object): - - prot_err_msg = 'Should not have empty secondary structure in protein residues' +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 = "DSSP.plot_content() did not produce an Axes instance" + plot_err_msg = 'plot_content() did not produce an Axes instance' - @pytest.fixture() + n_prot = 214 + + @pytest.fixture(scope='class') def universe(self): return mda.Universe(TPR, XTC) - - @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 mdtraj_codes(self): - return np.load(ADK_DSSP) - - @pytest.fixture() - def mdtraj_simple(self, dssp): - mdtraj_simple = np.load(ADK_DSSP_SIMPLE).astype(object) - for code in ['H', 'E', 'C']: - mdtraj_simple[mdtraj_simple == - code] = dssp.ss_codes_to_simple[code] - return mdtraj_simple - - def test_correct_assignments(self, dssp, mdtraj_codes, mdtraj_simple): - assert_equal(dssp.ss_codes, mdtraj_codes) - assert_equal(dssp.ss_simple, mdtraj_simple) - assert_equal(dssp.ss_mode, XTC_modes) - assert_equal(dssp.simple_mode, XTC_modes_simple) + + 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], XTC_code_counts[k]) + 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], XTC_simple_counts[k]) + assert_almost_equal(dssp.simple_counts[k], self.simple_counts[k]) def test_add_topology_attr(self, universe): assert not hasattr(universe._topology, 'secondary_structures') - dssp = ss.DSSP(universe, select='backbone', o_name='O O1', - add_topology_attr=True) + dssp = self.analysis_class(universe, select='backbone', + add_topology_attr=True) dssp.run() assert hasattr(universe._topology, 'secondary_structures') secstruct = universe.residues.secondary_structures - assert not np.any(secstruct[:214] == ''), self.prot_err_msg - assert np.all(secstruct[214:] == ''), self.other_err_msg + 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_non_protein(self, universe, mdtraj_codes, mdtraj_simple): - n_prot = 214 - dssp = ss.DSSP(universe, select='all', o_name='O O1') + def test_non_protein(self, universe, all_codes, all_simple): + dssp = self.analysis_class(universe, select='all') dssp.run() - assert_equal(dssp.ss_codes[:, :n_prot], mdtraj_codes) - assert_equal(dssp.ss_simple[:, :n_prot], mdtraj_simple) - assert_equal(dssp.ss_mode[:n_prot], XTC_modes) - assert_equal(dssp.simple_mode[:n_prot], XTC_modes_simple) + assert_equal(dssp.ss_codes[:, :self.n_prot], all_codes) + assert_equal(dssp.ss_simple[:, :self.n_prot], all_simple) + assert_equal(dssp.ss_mode[:self.n_prot: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], XTC_code_counts[k]) + 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], XTC_simple_counts[k]) + assert_almost_equal(dssp.simple_counts[k], self.simple_counts[k]) - assert np.all(dssp.ss_codes[:, n_prot:] == ''), self.other_err_msg - assert np.all(dssp.ss_simple[:, n_prot:] == ''), self.other_err_msg - assert np.all(dssp.ss_mode[n_prot:] == ''), self.other_err_msg - assert np.all(dssp.simple_mode[n_prot:] == ''), self.other_err_msg + assert np.all(dssp.ss_codes[:, self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.ss_simple[:, self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.ss_mode[self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.simple_mode[self.n_prot:] == ''), self.other_err_msg def test_plot_all_lines(self, dssp): ax = dssp.plot_content(kind='line', simple=False) @@ -192,3 +125,142 @@ def test_plot_n_xticks(self, dssp): 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): + assert_equal(dssp.ss_codes, all_codes) + assert_equal(dssp.ss_simple, all_simple) + 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]) + +@pytest.mark.skipif(executable_not_found('mkdssp'), + reason="Tests skipped because mkdssp (DSSP) not found") +class TestDSSPWrapper(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', '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) + + +@pytest.mark.skipif(executable_not_found('stride'), + reason="Tests skipped because STRIDE not found") +class TestStrideWrapper(BaseTestSecondaryStructure): + + 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.DSSPWrapper(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/data/adk_oplsaa_dssp.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy similarity index 100% rename from testsuite/MDAnalysisTests/data/adk_oplsaa_dssp.npy rename to testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy diff --git a/testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy similarity index 100% rename from testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy rename to testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy 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 0000000000000000000000000000000000000000..61e1e9441cd0eed353aa6d71a2f856674194909c GIT binary patch literal 8688 zcmdUwJ&Ti35QSswuSge#EnKk@E5)xu5!Av;7B^xcu4KhV{1yIT^AryphBx;kY(QVQ zoS8c_XU*J!!dwCMc zH4fsVC694=^+CCe=_kU4zohw;cj$8t`sYryxzC&#J=F3Y;-&{N_vIe)&RDB2`5twH zBen993y1zFT4(=R4mC2ncD^t8awhED6LI`w^-^2;^_{85t#0e!A&eX`kXcAbw42Y?ge{B1NnR{Zes3BEi*mH%{`}fdBGkhIyk+& z8qeOxaHx^l^Y6iXICCa!_hdi(s%yzl&aAEc`sTE_6<%C+q|V24Wka!?!1BF7}fCQ{K^= zzG!mRit0B{y(7A)JJt%u+GUauObhqYqOY58hoc5QgvJ9f^5 z?Vjw1-|9h3?y8sC%C9efk6V80Kx&tIW=_VJ-j#9Z1@ov6%=-V~{}q~@1$$)2=8(_# E7k{~X2><{9 literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..142a2871c4f512f816b92572d157a0077b7b3010 GIT binary patch literal 8688 zcmd6p&1zIz5QW>RPtmhh$iRpraisV&Q5=|oBMBODAVv~#B0hyLY%lPGk7A#*_fFF4 zg=W>ys#U8t_uSlHuTEaRetT>4ee+}Y;q3DL#qQyD_vqw!ceLGoJiqvMar))m`NfB` z)SsPxzC7F4FF&1rJ=^DZj_)09?;k(7d$j$zz1lywpKLap{GK*`Z}zpd--~^}_e0-1 z68GlZN5i#|m^tIfoH?0YdxP}-c9omH_10qhB(}C=`fGah!eusb%OiU5%z6B022QKz zTHUZTnam()^YuxEK z;bstT!SH7Nume3~de*=0j(Uf+d%p46(mCnYGs%`^^&5zsgR}5LbP2`b+-t z0eX%ah*l6Dy&&HXyjIV(#NpQODz(a6@_(qh_j7AS?gVGOFjsYxx8A^>hXUUTsfKV?V#^jJXd{k`pmY_dY*Hg7(IB@gE#qIu?KAV zK4gu(L2j+++u_f6Vy5aQPtW>^Q=8EbmmNXxnbARGW~T1Gjz)b4=e{Yid$VtbM`P6| zhqLr1)raT!PT-&wghww3A0D}T<|hs}y!@C8*8gTM`9DjPl_zXQVfS^&Oo1ro`^ezL|XzSABB&qjxX|R-e>C_x&BYj)zh7 zEb*xSo8qxI$gLI4$!#;c{T5X>dHU8*Ji8-)T6?JV{K5l!{Mu*ESUJ5sv-KUE`|VoX z)OtJU+b6O09n)ua4a4=U{S&Wo_^#N=Opu;y@ao^H)c^J0E4+bz7vg%$+S}G1xd&DnVXpoO0bH>rxdizrDn-Y5$@#s$Y#;Q*aXH=Uxu--!*bbsl)GRHCAK>ywP T-~K-kkK7&fEsjUuI-0)$v%P#B literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 97c892048c8..8934ed578d6 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -188,6 +188,8 @@ "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) ] from pkg_resources import resource_filename @@ -512,8 +514,10 @@ NAMDBIN = resource_filename(__name__, 'data/adk_open.coor') -ADK_DSSP = resource_filename(__name__, 'data/adk_oplsaa_dssp.npy') -ADK_DSSP_SIMPLE = resource_filename(__name__, 'data/adk_oplsaa_dssp_simple.npy') +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') # This should be the last line: clean up namespace del resource_filename From d7804a5f8c8ae45555a787d68213bc5da3d8fde8 Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Sat, 7 Mar 2020 22:18:40 +1100 Subject: [PATCH 4/6] added secondary_structure module --- .../analysis/secondary_structure.py | 451 ++++++++++++++++++ package/MDAnalysis/coordinates/PDB.py | 2 +- 2 files changed, 452 insertions(+), 1 deletion(-) create mode 100644 package/MDAnalysis/analysis/secondary_structure.py diff --git a/package/MDAnalysis/analysis/secondary_structure.py b/package/MDAnalysis/analysis/secondary_structure.py new file mode 100644 index 00000000000..c2458a80710 --- /dev/null +++ b/package/MDAnalysis/analysis/secondary_structure.py @@ -0,0 +1,451 @@ +# -*- 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 + +This module contains classes for computing secondary structure. MDAnalysis provides +a wrapper for the external DSSP_ and STRIDE_ programs to be run on an MD trajectory. + +Both DSSP_ and STRIDE_ need to be installed separately. + +.. _DSSP: https://swift.cmbi.umcn.nl/gv/dssp/ +.. _STRIDE: http://webclu.bio.wzw.tum.de/stride/ + + +Classes and functions +--------------------- + +.. autoclass:: DSSP +.. autoclass:: STRIDE + +""" + +import errno +import subprocess +import os +import numpy as np +import pandas as pd + +from ..core.topologyattrs import ResidueAttr +from ..lib import util +from .base import AnalysisBase + +class SecondaryStructure(ResidueAttr): + """Single letter code for secondary structure""" + attrname = 'secondary_structures' + singular = 'secondary_structure' + dtype = '2s}\n"), + "{tempFactor:6.2f} {segID:<4s}{element:>2s} \n"), 'REMARK': "REMARK {0}\n", 'COMPND': "COMPND {0}\n", 'HEADER': "HEADER {0}\n", From 0ad1730407d853ca1272ef61b7e5a739693221dc Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Mon, 9 Mar 2020 18:54:22 +1100 Subject: [PATCH 5/6] added mdtraj dssp --- .travis.yml | 2 +- .../analysis/secondary_structure/__init__.py | 27 ++ .../analysis/secondary_structure/base.py | 234 ++++++++++++++++ .../analysis/secondary_structure/dssp.py | 251 ++++++++++++++++++ .../wrapper.py} | 232 +++------------- .../source/documentation_pages/references.rst | 32 +++ .../analysis/test_secondary_structure.py | 194 ++++++++++++++ .../MDAnalysisTests/data/adk_oplsaa_dssp.npy | Bin 0 -> 17248 bytes .../data/adk_oplsaa_dssp_simple.npy | Bin 0 -> 17248 bytes testsuite/MDAnalysisTests/datafiles.py | 6 + 10 files changed, 783 insertions(+), 195 deletions(-) create mode 100644 package/MDAnalysis/analysis/secondary_structure/__init__.py create mode 100644 package/MDAnalysis/analysis/secondary_structure/base.py create mode 100644 package/MDAnalysis/analysis/secondary_structure/dssp.py rename package/MDAnalysis/analysis/{secondary_structure.py => secondary_structure/wrapper.py} (58%) create mode 100644 testsuite/MDAnalysisTests/analysis/test_secondary_structure.py create mode 100644 testsuite/MDAnalysisTests/data/adk_oplsaa_dssp.npy create mode 100644 testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy diff --git a/.travis.yml b/.travis.yml index 1249fc9a70d..bb9b06f9886 100644 --- a/.travis.yml +++ b/.travis.yml @@ -28,7 +28,7 @@ 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" 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..a84562e3239 --- /dev/null +++ b/package/MDAnalysis/analysis/secondary_structure/base.py @@ -0,0 +1,234 @@ +# -*- 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 + + +""" + +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +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='4Aq5QbsXaux0jk_!@1Au14m1qcaJP$7#LQ6PeBM1!~rZkQnaG_N$;=j@-z zPMTM<*35k0%sONP|9tb}>u-N}=i<+czph^1-v4xW_0i?k$1k2=U0+`P{OazHyPMyB ze0BHo_P_JbZ+^YMeSLoa%gyh%ujBWgeQAK6!W-e&wO_lF#;2S5~L}#8Vx;-izHMy(y+QQQg-4KlXb2>mJPC z_nl(-Q=g;vKp!Wc)44kPK=Zgq_Y>F7tA9E9uKrrx*6GVf@Ar2u5OYyW{=STaegG+wh zn^&EEVDFdaP`8vf-SbiMThha+Z_njVtbXhFj`^r=`&ze{oZ?CrPu{$A(TT&4T6{aPokPWy=EqdI$WUb@HDHKsRFpS||~So_$c zJ(wH&zEi9Zw%<|n^I1+ld)_+!aql_4diekEORTT@7V}wd-FNlvVQb&1p8x7Lp6a)r z(?|0@_3xE=u=^C7*YE07s>f5`fuAMT2cISVo{QbL-=yz*p}utAdgt_dFLsZub6#Ds zKIr_^{Xf;b&TCJ7*QHn;+Mk#X?zJ!d=-`sy-mB9O&4GLQp8C+`2kTcXAE&eD*$x!(`_JD2sjRNvp7Sl!xvKb1P~V_xhYr}uU)Up+cMY9IPf@7*`aYHN_S*kr?bH34+a8zR+r7F{pQGm2zvOF)t+Nj_ zk38pl&!aD$a+*UO?&Yicd-ZJ}zOlNETfa7@*MIKC-d)}j`JCR#xxVsC_iN{BU42(i zOqb77sUEApHs*WDhrEev>wQ1WWj|$g)irO5^;ajJ^8G$b`D^$6RO;;8y!M4l@9kcF zsn1dS(7)tsiLJAbQr=N>(4$jMbEw0;{qeDPU;X$-dK2ZT-;4TUf5*I-@T|Wp4+!|TlfFi>%CL=G`;V1 z-g>cp9<>Mk%gN_-uHHVdb2y*4bxuCY$#;Eh)mLAA=SR^T)opyVJ|Es)jQqGq=hHhm zZ@xNozjaQ2OFAW8K2K%)Z1t_4@5u*Ud8K=t`m|1b>i0@M+9O@@xt#R--&JD$?eA^n Vv}gOC-nr-U?G2qDwLkvp{0kE2f!hE8 literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy b/testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy new file mode 100644 index 0000000000000000000000000000000000000000..c1e18db1c63fc390b4e89cc23c3769e581c7b8c8 GIT binary patch literal 17248 zcmds!y=qiZ6oq5!Q%pAsDNL{uJ5g*Dn-o?um=Oy}WFj`=Q~1J0%-1YftbNWsXBh5n z7W-%YoPGbizk2uKd-!qK-F&(|T<`z&-|W8K z?eF_{Uw7a4_w|#DXXl&Ci|0?zH@`Q1-{V(j|3C2cL+y~gD!zTtn}sfi=4bC%IcE8w z%X<>CKN_8{^Ktw1^?y+_Pw0Gm`^k?vIjaZ{?9B%c9nj@q-ioT{`0Tm4>yP5~>Tt)k zqPshN?v{{!hqLdezA3NHmA|y^eXsM?ywtki>#Mw3y6Qu%ey(PJG&(=+J74zrFAmvD z=zMzn$&b0q?fUS*-hA-T0o9j&FIM%~b3ShF`c8*_mB*bF(b4|eYToIPJ0^5|`m^qa zzGBMDTzRRU^$v8MKY111K6<^Vy1hD8p8C1Es@J;L>gV?L|7h)a+HcM8_Ek64nm5&| z=W=LYd9GG}w^!{~_1LTGvUhAAzM$UuST+9OEMN2$xqB6rul1?s*PN+N?|`57+#J55 z`wzP3bhWR%R4adUo&1@P_C*hUMdpgSNA3O-efAQn-_knyF(+pg;ektY!bi{LVBU(V z=lJZox$BSO_3CiHwW7N_eeRZ!eTTE|=loJ$=8o_Bv)-Y)>Q7!p-J@2$)*2nfR-Nj; z`nbBP*Sgp0=l1pgXzh5~Z_V%aRX5d|H`S`=a%f+9u2z4ySM68z*mJsW?)pxLelHI` zsCQml%{vvjUlo;~>QQ%eKGodpIUVhf%I$R2kD8o{%t5Ey+D8v;<)?jgTQTV_ook+| z9z6X&TCRa7Qe#|v*s^NjvJZTR-=g+)~tvc*i^^djZ?uNgcucgmB6Y~EYCjP8@ zC100U`${T|Bu#HJ(sJsU-e^8^QIae*qYDnqu1&W5A5Z`L#K+fdFftmzLtLa58OXt z(wTKH=bQ2}cYN1B_B+TyUPb2M=zOzu$%%d|!XJgL{%&9YkJgT-{iW|xy=tv_Q?0tA zb31?aO?9gd`$_lI=1cE?>U?)^{#O;94*gkoKwq(%n?2?C)=S^%e95V(`%C@d&(cL- zQT<%4yiv9DYwhd*(YmU)^j%dS-P63OR{ql5@X>Yt%&XX{!+uqN^q%wQK8pBje{D7I QbjaNjIzIhbcS9fi3xC*Q>;M1& literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index f6873c7a7f0..97c892048c8 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -186,6 +186,8 @@ "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) ] from pkg_resources import resource_filename @@ -509,5 +511,9 @@ 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/adk_oplsaa_dssp.npy') +ADK_DSSP_SIMPLE = resource_filename(__name__, 'data/adk_oplsaa_dssp_simple.npy') + # This should be the last line: clean up namespace del resource_filename From aded6deae5b00c7798d4280a58030655738cd375 Mon Sep 17 00:00:00 2001 From: Lily Wang Date: Tue, 10 Mar 2020 17:10:47 +1100 Subject: [PATCH 6/6] added tests for wrapper classes --- .travis.yml | 18 + maintainer/install_dssp.sh | 126 +++++++ maintainer/install_stride.sh | 117 +++++++ .../analysis/secondary_structure/base.py | 13 +- .../analysis/secondary_structure/wrapper.py | 67 ++-- .../source/documentation_pages/references.rst | 11 +- .../analysis/test_secondary_structure.py | 328 ++++++++++++------ .../MDAnalysisTests/coordinates/test_pdb.py | 43 ++- .../data/{ => dssp}/adk_oplsaa_dssp.npy | Bin .../{ => dssp}/adk_oplsaa_dssp_simple.npy | Bin .../data/dssp/adk_oplsaa_mkdssp.npy | Bin 0 -> 8688 bytes .../data/dssp/adk_oplsaa_stride.npy | Bin 0 -> 8688 bytes .../data/dssp/mkdssp_phi_psi_sasa.npy | Bin 0 -> 51488 bytes .../data/dssp/stride_phi_psi_sasa.npy | Bin 0 -> 51488 bytes testsuite/MDAnalysisTests/datafiles.py | 12 +- 15 files changed, 576 insertions(+), 159 deletions(-) create mode 100755 maintainer/install_dssp.sh create mode 100755 maintainer/install_stride.sh rename testsuite/MDAnalysisTests/data/{ => dssp}/adk_oplsaa_dssp.npy (100%) rename testsuite/MDAnalysisTests/data/{ => dssp}/adk_oplsaa_dssp_simple.npy (100%) create mode 100644 testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_mkdssp.npy create mode 100644 testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_stride.npy create mode 100755 testsuite/MDAnalysisTests/data/dssp/mkdssp_phi_psi_sasa.npy create mode 100755 testsuite/MDAnalysisTests/data/dssp/stride_phi_psi_sasa.npy diff --git a/.travis.yml b/.travis.yml index bb9b06f9886..ee68a9d0f5d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -34,6 +34,8 @@ env: - 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..36ad8585398 --- /dev/null +++ b/maintainer/install_dssp.sh @@ -0,0 +1,126 @@ +#!/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 + ;; + 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/base.py b/package/MDAnalysis/analysis/secondary_structure/base.py index a84562e3239..09817aa94b7 100644 --- a/package/MDAnalysis/analysis/secondary_structure/base.py +++ b/package/MDAnalysis/analysis/secondary_structure/base.py @@ -34,8 +34,6 @@ """ import numpy as np -import pandas as pd -import matplotlib.pyplot as plt from ...core.topologyattrs import ResidueAttr from ..base import AnalysisBase @@ -129,7 +127,8 @@ def __init__(self, universe, select='backbone', verbose=False, super(SecondaryStructureBase, self).__init__(universe.universe.trajectory, verbose=verbose) self._universe = universe.universe - self.residues = universe.select_atoms(select).residues + self.atomgroup = universe.select_atoms(select) + self.residues = self.atomgroup.residues self.n_residues = len(self.residues) self._add_topology_attr = add_topology_attr @@ -210,6 +209,14 @@ def plot_content(self, ax=None, simple=False, figsize=(8, 6), ax: :class: `matplotlib.axes.Axes` """ + try: + import pandas as pd + except ImportError: + raise ImportError('pandas is required for this function.') + try: + import matplotlib.pyplot as plt + except ImportError: + raise ImportError('matplotlib is required for this function.') if not simple: data = self.ss_counts diff --git a/package/MDAnalysis/analysis/secondary_structure/wrapper.py b/package/MDAnalysis/analysis/secondary_structure/wrapper.py index 3b38bc52f42..8b5ec9c812b 100644 --- a/package/MDAnalysis/analysis/secondary_structure/wrapper.py +++ b/package/MDAnalysis/analysis/secondary_structure/wrapper.py @@ -51,6 +51,8 @@ import errno import tempfile import subprocess +import warnings +import numpy as np from ...due import due, Doi from ...lib import util @@ -66,11 +68,12 @@ path="MDAnalysis.analysis.secondary_structure.wrapper", cite_module=True) -due.cite(Doi("10.1093/nar/gkh429"), +due.cite(Doi("10.1002/prot.340230412"), description="STRIDE", path="MDAnalysis.analysis.secondary_structure.wrapper", cite_module=True) + class SecondaryStructureWrapper(SecondaryStructureBase): """Base class for secondary structure analysis that wraps an external program. @@ -97,15 +100,19 @@ def __init__(self, executable, universe, select='protein', verbose=False, select=select, verbose=verbose, add_topology_attr=add_topology_attr) - + def _execute(self): """Run the wrapped program.""" + # ignore PDB warnings + warnings.filterwarnings("ignore", + message="Found no information for attr") + fd, pdbfile = tempfile.mkstemp(suffix='.pdb') os.close(fd) try: - self.residues.atoms.write(pdbfile) + self.atomgroup.write(pdbfile) cmd_args = [self.exe] + self.cmd.format(pdb=pdbfile).split() - proc = subprocess.Popen(cmd_args, stdout=subprocess.PIPE, + proc = subprocess.Popen(cmd_args, stdout=subprocess.PIPE, stderr=subprocess.PIPE) stdout, stderr = proc.communicate() finally: @@ -113,7 +120,7 @@ def _execute(self): os.unlink(pdbfile) except OSError: pass - return output + return stdout.decode('utf-8') def _prepare(self): super(SecondaryStructureWrapper, self)._prepare() @@ -121,11 +128,11 @@ def _prepare(self): self.phi = np.zeros((nf, self.n_residues), dtype=float) self.psi = np.zeros((nf, self.n_residues), dtype=float) self.sasa = np.zeros((nf, self.n_residues), dtype=float) - + def _compute_dssp(self): output = self._execute() self._process_output(output) - + class DSSPWrapper(SecondaryStructureWrapper): """Runs :program:`mkdssp` on a trajectory. @@ -186,35 +193,40 @@ class DSSPWrapper(SecondaryStructureWrapper): """ exe_name = 'mkdssp' - cmd = '-i {pdb} -o stdout' + cmd = '-i {pdb}' columns = { - 'ss': 17, + 'ss': 16, 'sasa': (35, 39), - 'phi': (105, 110), + 'phi': (104, 109), 'psi': (111, 116), } def __init__(self, universe, select='protein', executable='mkdssp', verbose=False, add_topology_attr=False): - super(DSSPWrapper, self).__init__(executable, universe, - select=select, verbose=verbose, - add_topology_attr=add_topology_attr) - + super(DSSPWrapper, self).__init__(executable, universe, + select=select, verbose=verbose, + add_topology_attr=add_topology_attr) def _process_output(self, output): frame = self._frame_index per_res = output.split(' # RESIDUE AA STRUCTURE BP1 BP2 ACC')[1] - lines = per_res.split()[1:] - for i, line in enumerate(lines): + i = 0 + for line in per_res.split('\n')[1:]: + try: + if line[13] == '!': + continue + except IndexError: + continue ss = line[self.columns['ss']] - if not ss: + + if ss == ' ': ss = 'C' self.ss_codes[frame][i] = ss for kw in ('phi', 'psi', 'sasa'): _i, _j = self.columns[kw] getattr(self, kw)[frame][i] = line[_i:_j] - + i += 1 class STRIDEWrapper(SecondaryStructureWrapper): @@ -241,7 +253,7 @@ class STRIDEWrapper(SecondaryStructureWrapper): attribute ``secondary_structure`` to your residues. verbose: bool, optional Turn on more logging and debugging. - + Attributes ---------- residues: :class:`~MDAnalysis.core.groups.ResidueGroup` @@ -279,17 +291,18 @@ class STRIDEWrapper(SecondaryStructureWrapper): def __init__(self, universe, select='protein', executable='stride', verbose=False, add_topology_attr=False): - super(STRIDEWrapper, self).__init__(executable, universe, - select=select, verbose=verbose, - add_topology_attr=add_topology_attr) + super(STRIDEWrapper, self).__init__(executable, universe, + select=select, verbose=verbose, + add_topology_attr=add_topology_attr) def _process_output(self, output): lines = output.split('\nASG ')[1:] lines[-1] = lines[-1].split('\n')[0] frame = self._frame_index for i, line in enumerate(lines): - resname, chain, resid, resnum, ss, ssname, phi, psi, area = line.split() - self.ss_codes[frame][i] = ss - self.phi[frame][ss] = phi - self.psi[frame][ss] = psi - self.sasa[frame][ss] = area \ No newline at end of file + # resname chain resid resnum ss ssname phi psi area ~~~~ + fields = line.split() + self.ss_codes[frame][i] = fields[4] + self.phi[frame][i] = fields[6] + self.psi[frame][i] = fields[7] + self.sasa[frame][i] = fields[8] diff --git a/package/doc/sphinx/source/documentation_pages/references.rst b/package/doc/sphinx/source/documentation_pages/references.rst index d6d103c4303..945d2fbc3a1 100644 --- a/package/doc/sphinx/source/documentation_pages/references.rst +++ b/package/doc/sphinx/source/documentation_pages/references.rst @@ -164,13 +164,14 @@ If you use the DSSP class in :mod:`MDAnalysis.analysis.secondary_structure.dssp` .. _`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 [Heinig2004]_ . +please cite [Frishman1995]_ . -.. [Heinig2004] Matthias Heinig and Dmitrij Frishman. - STRIDE: a web server for secondary structure assignment from known atomic coordinates of proteins. - *Nucleic Acids Research* **32 (Web Server issue)** (2004), W500-W502. doi: `10.1093/nar/gkh429`_ +.. [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.1093/nar/gkh429`: https://doi.org/10.1093/nar/gkh429 +.. _`10.1002/prot.340230412`:https://doi.org/10.1002/prot.340230412 .. _citations-using-duecredit: diff --git a/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py b/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py index a6c7f551ed9..57cbcaf79cc 100644 --- a/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py +++ b/testsuite/MDAnalysisTests/analysis/test_secondary_structure.py @@ -30,140 +30,74 @@ import MDAnalysis as mda from MDAnalysis.analysis import secondary_structure as ss -from MDAnalysisTests.datafiles import TPR, XTC, ADK_DSSP, ADK_DSSP_SIMPLE - -XTC_modes = ['C', 'E', 'E', 'E', 'E', 'E', 'C', 'C', 'T', 'T', 'S', 'C', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'C', 'C', - 'C', 'E', 'E', 'E', 'T', 'T', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'C', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'T', 'C', 'C', 'C', 'C', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'T', 'S', 'G', 'G', 'G', 'S', - 'S', 'C', 'E', 'E', 'E', 'E', 'S', 'C', 'C', 'C', 'S', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'T', 'C', 'C', 'C', 'S', - 'E', 'E', 'E', 'E', 'E', 'E', 'C', 'C', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'T', 'E', 'E', 'E', 'S', 'T', 'T', 'T', 'S', - 'S', 'E', 'E', 'E', 'T', 'T', 'T', 'B', 'C', 'S', 'S', 'S', 'T', - 'T', 'B', 'C', 'T', 'T', 'T', 'C', 'C', 'B', 'C', 'S', 'S', 'S', - 'G', 'G', 'G', 'S', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'H', 'H', 'T', 'S', 'C', 'C', 'E', 'E', 'E', - 'E', 'E', 'C', 'S', 'S', 'C', 'H', 'H', 'H', 'H', 'H', 'H', 'H', - 'H', 'H', 'H', 'H', 'C', 'C'] - -XTC_modes_simple = ['Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', - 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Coil', 'Coil', 'Coil', 'Strand', 'Strand', - 'Strand', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Coil', 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Coil', 'Strand', 'Strand', 'Strand', - 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Coil', 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', 'Coil', - 'Coil', 'Coil', 'Coil', 'Strand', 'Strand', 'Strand', 'Coil', - 'Coil', 'Coil', 'Strand', 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', - 'Coil', 'Strand', 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', 'Coil', - 'Strand', 'Coil', 'Coil', 'Coil', 'Coil', 'Helix', 'Helix', - 'Helix', 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Coil', 'Coil', 'Coil', 'Coil', 'Strand', - 'Strand', 'Strand', 'Strand', 'Strand', 'Coil', 'Coil', 'Coil', - 'Coil', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', - 'Helix', 'Helix', 'Helix', 'Helix', 'Helix', 'Coil', 'Coil'] - -XTC_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])} +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, + ) -XTC_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])} +def mdtraj_found(): + try: + import mdtraj + except ImportError: + return False + else: + return True -class TestDSSP(object): - prot_err_msg = 'Should not have empty secondary structure in protein residues' +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 = "DSSP.plot_content() did not produce an Axes instance" + plot_err_msg = 'plot_content() did not produce an Axes instance' - @pytest.fixture() + n_prot = 214 + + @pytest.fixture(scope='class') def universe(self): return mda.Universe(TPR, XTC) - @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 mdtraj_codes(self): - return np.load(ADK_DSSP) - - @pytest.fixture() - def mdtraj_simple(self, dssp): - mdtraj_simple = np.load(ADK_DSSP_SIMPLE).astype(object) - for code in ['H', 'E', 'C']: - mdtraj_simple[mdtraj_simple == - code] = dssp.ss_codes_to_simple[code] - return mdtraj_simple - - def test_correct_assignments(self, dssp, mdtraj_codes, mdtraj_simple): - assert_equal(dssp.ss_codes, mdtraj_codes) - assert_equal(dssp.ss_simple, mdtraj_simple) - assert_equal(dssp.ss_mode, XTC_modes) - assert_equal(dssp.simple_mode, XTC_modes_simple) + 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], XTC_code_counts[k]) + 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], XTC_simple_counts[k]) + assert_almost_equal(dssp.simple_counts[k], self.simple_counts[k]) def test_add_topology_attr(self, universe): assert not hasattr(universe._topology, 'secondary_structures') - dssp = ss.DSSP(universe, select='backbone', o_name='O O1', - add_topology_attr=True) + dssp = self.analysis_class(universe, select='backbone', + add_topology_attr=True) dssp.run() assert hasattr(universe._topology, 'secondary_structures') secstruct = universe.residues.secondary_structures - assert not np.any(secstruct[:214] == ''), self.prot_err_msg - assert np.all(secstruct[214:] == ''), self.other_err_msg + 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_non_protein(self, universe, mdtraj_codes, mdtraj_simple): - n_prot = 214 - dssp = ss.DSSP(universe, select='all', o_name='O O1') + def test_non_protein(self, universe, all_codes): + dssp = self.analysis_class(universe, select='all') dssp.run() - assert_equal(dssp.ss_codes[:, :n_prot], mdtraj_codes) - assert_equal(dssp.ss_simple[:, :n_prot], mdtraj_simple) - assert_equal(dssp.ss_mode[:n_prot], XTC_modes) - assert_equal(dssp.simple_mode[:n_prot], XTC_modes_simple) + assert_equal(dssp.ss_codes[:, :self.n_prot], all_codes) + assert_equal(dssp.ss_mode[:self.n_prot: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], XTC_code_counts[k]) + 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], XTC_simple_counts[k]) + assert_almost_equal(dssp.simple_counts[k], self.simple_counts[k]) - assert np.all(dssp.ss_codes[:, n_prot:] == ''), self.other_err_msg - assert np.all(dssp.ss_simple[:, n_prot:] == ''), self.other_err_msg - assert np.all(dssp.ss_mode[n_prot:] == ''), self.other_err_msg - assert np.all(dssp.simple_mode[n_prot:] == ''), self.other_err_msg + assert np.all(dssp.ss_codes[:, self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.ss_simple[:, self.n_prot:] + == ''), self.other_err_msg + assert np.all(dssp.ss_mode[self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.simple_mode[self.n_prot:] == ''), self.other_err_msg def test_plot_all_lines(self, dssp): ax = dssp.plot_content(kind='line', simple=False) @@ -192,3 +126,179 @@ def test_plot_n_xticks(self, dssp): 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, all_codes, all_simple): + dssp = self.analysis_class(universe, select='all') + dssp.run() + + assert_equal(dssp.ss_codes[:, :self.n_prot], all_codes) + assert_equal(dssp.ss_simple[:, :self.n_prot], all_simple) + assert_equal(dssp.ss_mode[:self.n_prot: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]) + + assert np.all(dssp.ss_codes[:, self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.ss_simple[:, self.n_prot:] + == ''), self.other_err_msg + assert np.all(dssp.ss_mode[self.n_prot:] == ''), self.other_err_msg + assert np.all(dssp.simple_mode[self.n_prot:] == ''), self.other_err_msg + + +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) + 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) + + +@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/adk_oplsaa_dssp.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy similarity index 100% rename from testsuite/MDAnalysisTests/data/adk_oplsaa_dssp.npy rename to testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp.npy diff --git a/testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy b/testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy similarity index 100% rename from testsuite/MDAnalysisTests/data/adk_oplsaa_dssp_simple.npy rename to testsuite/MDAnalysisTests/data/dssp/adk_oplsaa_dssp_simple.npy 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 0000000000000000000000000000000000000000..61e1e9441cd0eed353aa6d71a2f856674194909c GIT binary patch literal 8688 zcmdUwJ&Ti35QSswuSge#EnKk@E5)xu5!Av;7B^xcu4KhV{1yIT^AryphBx;kY(QVQ zoS8c_XU*J!!dwCMc zH4fsVC694=^+CCe=_kU4zohw;cj$8t`sYryxzC&#J=F3Y;-&{N_vIe)&RDB2`5twH zBen993y1zFT4(=R4mC2ncD^t8awhED6LI`w^-^2;^_{85t#0e!A&eX`kXcAbw42Y?ge{B1NnR{Zes3BEi*mH%{`}fdBGkhIyk+& z8qeOxaHx^l^Y6iXICCa!_hdi(s%yzl&aAEc`sTE_6<%C+q|V24Wka!?!1BF7}fCQ{K^= zzG!mRit0B{y(7A)JJt%u+GUauObhqYqOY58hoc5QgvJ9f^5 z?Vjw1-|9h3?y8sC%C9efk6V80Kx&tIW=_VJ-j#9Z1@ov6%=-V~{}q~@1$$)2=8(_# E7k{~X2><{9 literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..142a2871c4f512f816b92572d157a0077b7b3010 GIT binary patch literal 8688 zcmd6p&1zIz5QW>RPtmhh$iRpraisV&Q5=|oBMBODAVv~#B0hyLY%lPGk7A#*_fFF4 zg=W>ys#U8t_uSlHuTEaRetT>4ee+}Y;q3DL#qQyD_vqw!ceLGoJiqvMar))m`NfB` z)SsPxzC7F4FF&1rJ=^DZj_)09?;k(7d$j$zz1lywpKLap{GK*`Z}zpd--~^}_e0-1 z68GlZN5i#|m^tIfoH?0YdxP}-c9omH_10qhB(}C=`fGah!eusb%OiU5%z6B022QKz zTHUZTnam()^YuxEK z;bstT!SH7Nume3~de*=0j(Uf+d%p46(mCnYGs%`^^&5zsgR}5LbP2`b+-t z0eX%ah*l6Dy&&HXyjIV(#NpQODz(a6@_(qh_j7AS?gVGOFjsYxx8A^>hXUUTsfKV?V#^jJXd{k`pmY_dY*Hg7(IB@gE#qIu?KAV zK4gu(L2j+++u_f6Vy5aQPtW>^Q=8EbmmNXxnbARGW~T1Gjz)b4=e{Yid$VtbM`P6| zhqLr1)raT!PT-&wghww3A0D}T<|hs}y!@C8*8gTM`9DjPl_zXQVfS^&Oo1ro`^ezL|XzSABB&qjxX|R-e>C_x&BYj)zh7 zEb*xSo8qxI$gLI4$!#;c{T5X>dHU8*Ji8-)T6?JV{K5l!{Mu*ESUJ5sv-KUE`|VoX z)OtJU+b6O09n)ua4a4=U{S&Wo_^#N=Opu;y@ao^H)c^J0E4+bz7vg%$+S}G1xd&DnVXpoO0bH>rxdizrDn-Y5$@#s$Y#;Q*aXH=Uxu--!*bbsl)GRHCAK>ywP T-~K-kkK7&fEsjUuI-0)$v%P#B literal 0 HcmV?d00001 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 0000000000000000000000000000000000000000..bcb84129cc46f5e1e440621a3688c4ee28aa5372 GIT binary patch literal 51488 zcmbWAU5I7Zb?0~DnI30SA;PU>eEpGJmE_Wg%W_LsmerDM9hF+8uPdpglGI(cY?n&X zwVm{g6<5a4Oa=}F^i&T64Fh+`LxG_n&<}=!frfw&3%ac%;0eJ|K#}o zr)N*xzdXKw_LHCf=}+JLcR%>aPe1tZ&h^#zu55i6*SCK3-bWwC`2&;X@gwJ!$B#bp z*n`XC|2Y2t`2Fj0=hy$NtH<#Vy8kkd-{|u2KlbA{eJ79CH@eB4-`VeXvpc^t?{)S5 z_Qj!&pTC*M_3!2VPu}4CJNeu-#*NnVeCySr&U5r`p09EL>YGDdXX_Q_o#t`$dXAs1 z<@xP5hAMCWy_CCqr%Smva{bk{9RKkp=J_u2Zsz#mg}lG}TCR6|F^^|2vhNr2x%KDs z`s9@yfArc=>hB`+uQC3WyuR^vm-pY!>-*00W{zL_PRe<(mgCDSIc}MKPmt?qn)f@_ zU1I-7Z|3-07xVtk<*wZM?R>YZdb{>^t}}ZT`H?g1k@lRjU*x`m96j+F9!#6u)sC}wc(eO zgPSAPAD{Lwo(o<~J?c03(fT`i{Rmudy`1ay*g5c(I7eTB$7=rraC#OwjyUeVik|V~ zkf-83_@vzV#k}?WI({;R=l=Y?uHGMhkowvNFJ&|Bxr>}L?0W|NU;TdC`8@tUxxoC& zwUYAA(QDW}{Zsv(+Gp4!J_9f8F|NZuw)M03(*6_p z?HC@Nvi}8n@Qtpnu<`a#k7Mvwuy^1zue0B!)uEmjc}AQP-_Bgh`#b8x{lM43W9i@c z%Psu(2>v@{UlD(j$KaER^E2;FaRXk?uYn_;o5C9d_6U5WUEV=2m-2j@cz3hu1H6LA zm^b}^{m))TF6~CV4&JYN8D7i#w~%Yj{%^yNH;{9Ud3T9B(--pj{!9qMEtG#6mQ`t zRgRm)rJKk(g|An?&-fQpjwN`Y@tpjDyi?-Q9DJ?6l=|tgf8e98xA&c#_lo>Rp7hu6 zr5}y3^OrtM`S$RiWBk8^-zVe|yYJ?6x3SxK?Qb`>as+|tm@78L{HUA#_I^m}}{u*-sCGjuvy=tGS@dkhU zoADG zJbn_V(~jm(;L&)xZ9e-{>TNFm!OJRo8qAl;(*nPVYy9cprL@x-^l(MHm}indhTqkF z4%7qtFXcRukJR}N%sY5K@|i0Cn*4eqPaVIG{lee#`3pII3p@m0q@2W=8}Q@*LQkvc zY3B#2*Z;iPRrw0Mzd?QxxUBZP1&>BNspr?lFLD*~;k}oUN1lRT0FyaPE$ z*lCVkgAc3xSK-eT?j&9n_&w~N`Y@gbwWkgW7M-J|jNV`y=x_c%(zG zXNWHg@bCfocg&wSryh?uTjd)2Jo)fnyq|dc@n+ibEWEk~4zH5GpU2OZ@Xy=CsdqQK zI&c3}>Tkk+C+Ii!pZkOV<)Prm%!B0xaMl|ysQ;oasL$UvPMG&F|37ek7QC+#M|bsi z;{RVE|Hs__#5zyDG9O|8*R&7uw2*fU*zYssjn}MWiI=hOD(^b?w_Z>Em>>4?u6Dpr zq7JF{TEOpa5bqA~`@g_GTg1)C%M-`Qvqro+gMUuGnfJ~Aj>+3U)}M*%7t!CA`470+ z#*fC}YyT4SS?5q+U30(i?=Pqy%6*?+KlD`hw~fEv#GkiUhU5kP3jYawR6Rx>n0f(6 zfrE;hPmHfG zUN_7`>~o1L`0|>03HzuI$3ND~YwRofdUf43`fY_g^Dn@`n0Qst&pz>ULtKe7_P45i z^keNNKalsV%bWHK%*%~;@(%G)J@3N{Q4i;q(DxE?_7MG?C!abnUaIexz#Dx%^UewR zP{Cd|PWU$CAM->%BK3ve59Im!W#sJ@_kow>J@&nU|6Gwbz}W;G>*w;(5jc)Gn(G@U z;KS?kGy2}KzIg;3k+0mq{^5^R?%6$pRR)=`xa3br$3R8L>^K1a}|6{@rNn# z;+l0id3D6$YM()U%2OY9b>4mOHg+EQ7CU`r)>Sz_v0hRhc=G51_BUU`kE3sr>$2a- zOX@sFt6XRNd2t~hkiVHXM}N1_>l1`LN?J;GgjCy5A-6Hv(U4taFonN1rIi$$!{wLO9eV3}&;KTadD*43-`9h!dzUv>on|fL>Zvr>lxI$z{_b^l}LxdxxK_Zam5cgbV7jHlY!z9Bd{gI^Ewnf)a28}iomcHqhApH#Vm z#}gOCjlZH#@!!D52J`(G`+q`z_7?FV`r^4?@Vm~sQ3uueCmi?Sq0-{7OC&%FvS zf6Q?PPsTCz1rFEjo58F0`Rb(cPv9x-1%4L9x%Nlyq#^3tlO zDR$eDC$WFXT}9})9(u0lVUN@U_Pfsf_8FGtWAt`|zT!vZ@eA1b=gP_XP4ri&)1$7c z@-4H@GV3gHf7kdgKZ2ue^Dp&@p3j<(fWwF{?jv!-zK+>{#0PO*<+x=Ycbcah2R|=Z zFTnKK=5){~+yf1^w;WFCi{0Urv0kzu#4SMEtM&SRk*8cu}wS z7uY9xUhr{_V_oyx5qj8`m#|;N(<)c+ca?KSeG&0J{R7-}_QjENfx7Br^Em3ktCw@0 zGvKA*hcDT`U?1DO-?)KZZ=3hqH#4rmgExE+VqU3z@V}0Dbg%K}`P6UZ6Lnv^=EM5M zNnEM=S_XgNPgNfi{Qg(utDSkYc{BJqAa9TNtcfG^CBH=cNS?5-iJlHP2JTY-_C=BR z6Y&WS&Y!|j^eyVV@m@Wj7e8OaZs6@E{2hIfI&MvU_+Ep!HOCGSC#ycTc~4_rIk3N| zv+DhQ__)t@xk{I)WK8VM?atTYyH(rUHSq1u*&+^U9u<9mttULiGaKMFtEv#(&h#O?z) z9MKO7+}7vg{c=4<9#HKX^3`?1zw7nryVw1+ej|Q(1NkHG%lodcUm));^2PTXb=*FF zaYlYJpX2$s>Z>hojlKSl4+=dbJi)@I(1?-T0z z`IEdM*I{4B%=6dqdc*%}%3~fD>uAS{tI5hx*YE8Y%NKm#6!}%kiM&^k^9udeZTP6; zxsiPf@Ozy&5&lxwv(LRlKJ;JN*MsOka8czsgeR`S53>()|8x9i8{XaF{FeUtc8-hp z2f42|@_e0f(Lbtk*w>DDSNHcWax|W!?hJmc_qQ*m{Kw$#s_};WKkAn{Z{+J$zen=r zh1^f{532oc=-2qmi1S-!#m4nC;k4=&_!oqT3E=>t@~ zbov#37!FW&oRbe7 zf_wE8`AzDdeDN&nY{*OCAmUe*KlGRQMt}A>7g#6kR>y@tv{OFEzv3a@D=Tl>6+fGJ z1pl|dZSdl!>^3ugsVDG$1m8zLEawFmbM;L>YS%iz{Kx#jdIUWk5KrddDC&x;&-3QN zZ^TlatKddV?A+r0M3Fhs(vDVRXxt|zxbZ7t{;8W)E|6iex#qs z13mbPI9>G-_RV?48TE10_jMin(lNfSb5*_oho6{-DJS}le7(-Qzmmr*;tsuktQ_EG z#{MG@t?L|Ehd7UZf&ZUnT;Qw9dB(WJd$5R$b$^lHrC!L(BOk8ke+HhOHZRs5j;DTT z?%itjTQuTWD=W?IuaZh{mJ}mq`<@s4x$M4$z1ZTIA`_GB5|5IGsFL}31`(8$$ zA26T1kNm->u2YZZ;rrzId*~(N+NZ|tyxrCF;K#bo@ul3)EpWVUoFKnjx_~^)7w@m@ zzLt!e9K#-U-T^-ung65rwNrSHzDpf%erLa^Lr)QBt6k1BU&vL*clK48*Y^v9{PZOD zqb`a1w>}sAP>&&J+VAB1g?imQG~z@(AB%Io2MZjhzVvJQ1QF-!_`CSN>F30!>-J&D zALIR0U3URKIRsCi7zeP&mVAXiqrXx0FqgOBYx}#=$Ek8%L#{L6_6qBSJ?j0ad#jv1 zyb|9j)^q;}+*Z4u0e83Tb7S|xddNOAynj`^khl8YEBHD6{+*mR>Xy{s^Qr#@^IZ4u zdl_&t1@B?^#2fs$#(4Sv8uDJ}yY8@KU1yKDHIg^5=Pve`eJjV^BwvU=WL+oXQ|iI| z0ltlXSG~S!eg*HDXT^F|l|A+7KhPfzen|V_-y6&m@u$u=qt5D(``CI?KK1>oc(#uN zFI{^lakY<~Z|XPXe<63popsTS{_+5CH9zFJb@L?qSIBz=zm1u10iM5&UT64U#PO>C z;MIhddA@x+^mz@xJ|MsQbNYeH#H;JzF5bUXJ$+1l)opfl-VQnAy=A>0-&t2Z2OjHv z<5Bow>X~@7VE!v_i!1lbzT}ZKWeaiCvd_o$cZzPEgjN<1p=pZUl1|8L_* z#?vFezeoOpej@)$ee!-a-iuUwUbp_`{lbBH195&EUib`obbODP`v4cm_+jM#_-P(* zTtc7rx4`oi>=fTqB!5!Jjj#%T|?dH*^_uxTwMo3Ov;L2I9tu{_RzH zoc&DTkLaV+@exPr`cdc9^Bwb3aI{XJaRK|CXT69kiCg;w_)}+`M!%u=)C+aUu6Z2$ zh&r;aKiB`LTP{)u|9A5b>qXzEXm5@;nJ4-ub-jpNbw3^TN%Y<8^(E?#RpWwr2DlBL zO8lVDsJ~J#@ZA;lhd$TX=RWJz!SRoIR)@ds(5JkC{-U2#?HB#dy576QZz6BG z?)wSWSztdm&{xEvy6?bY-H-Ws;l0OP|2LoT`$+Qb=no|R^$YlL9XY4)XxOQmqUGE6 z)8$WaAMv5icZ|PYCx84zJAjYL$^0q5@9W`bxmjAMxxThn4@;34`Bc@BPmP2VN#o9iKe)O%GQhvtuOcWEbhX8qOF*O=!5Pbm+2GoM~1K1O|2 z=ig@CLvcJ!Irr&P?!iaaH5>3w#LG{w-@)tl6RbDPbNwC*-i_}I>-x*!-S^3Xqm-L< z;_oNxz6;-%3_gb)g7@n9cn^>`VPB2E^dHu_CEu98THjbN5WnKReD$BJ;+%Xd`ao5l zL*kY1bq4a8crTmsm{+66o7i~`e=r}KAm1_k5`0#lkN0hL{dMi)dwh6wnSDf@tK;Ik zlzJQyA8v!gea2Z2q%v!~%5~7E{rVl`+=K7u#GhSwKkB;r+%@wj`w8&LocCOXI$;4E zTmdh>e;gSP_RZ}DAM>g9>V z#~Jdw_zofQ3;*tekL&Qlg9kFI{|b>n38XYzUW8-L%Me84>M_Y^sg zag^&@@X?m}uKgeOv233bTtxq|t`q!H?Hb?h)Or28oh5i)zm9%Uy?>qj`!n=0CgA#p z_=f+7F6Zy`8|uRN`-Hlle~!)p@VidDG#=an7k|cg2kztN`1v*UkDo-o zRM(9@V8vH_zfsTch1Wxls*lir;sKlnj+4Jm>Z+6yC_sL72#tz{>Ro>`F)MLb-dK`$4zhGbCFX>P0XIcM-r{g_F zy&v(X9@}>+*z1t-!Oy9eSGzol592#{GTuhtGWCR98@$(CP#$n_7JlAf{iO>j&l&4} z@E-l&x^D1C9?_fs-Gn3XKf!+@-qvx!(^Zb})6@@qy`}!mC!R~cl4n+tGwQs=J$$ia zyl3B$2iJAZv;U~$>i7lZ3H#Oc_ry6oao)ZWeWu67A-p{GyC?Fo=ojQZ><=3MSuf&N zuFJlo?yUPQ_9c+OG|tNXtJ?@yoUQjQCW zr`dyf|KMU?pFYp?oGEw4){_T|n8ual;cLsTfWA6cX3@g zSg${aJ}1oob@cKe`@5Le*U{tABaHhE#w{}L4Cn6S4c2`yuP>mtWAkhFc_5yi%W<*a zitFifxz53A&fo0!81sr-<+UCFCxOevgZgBiPxOPQa^Bca^5;0mkG{$L=(iJZ7ucV8 zcm}y%V!ZpsZ$^)C|5fHcJER{^xUL;u%jbLi==iNX?_LHUPcc7upWVmxXLFv}!^rLT zh2O~W12~EIMadJ-=5uQo^LPyYZ>=%@1@`wuo{ujx|Kd>Y{08HZe)lO1Rf1}k}scRee@i0Ca>RveO}6A@Je2Pysay)s(u#0(;oW>JC!@X z{>7aCSUbQ^jraB6U;nbsv#)-LPb2(!3Xjfye@MR}j@`mNd+4#A4#e@vxAMGTkKm2egZ$%l@{$?&-(BIkhjRauUpA z^)o}_Pkk#lc#Qa$xH8U=H;>5UBK{|C*yk92ia1s7{0@!t+UsoMX!a=kBmdiz$A}*V zJ+H%~@=nCPx?b43>SvvNbVPn}MBJO2|B_GaJ)ZKEr_;W>@NneO6$e|a6Y(c;aVigJ z`~aRJ?@YWsjULHIBL3XjcIn8=Bmcj1TFPBrJ(O=@w>NSh3+y-iCFe8WmB$aBQUCS+ z>`VFld~vAPuV2i0d-i$EzKur<#>>ZZ{@)AzfmK8`2l@@`Mc$VllyBqr}gzpQ-yT)N`?Z>gn$?zD=pGMr+9b_0&&~ z+(+au@qF$Vy^Qg<{kw+5JMrdY=Oh1#dNlQIo}<2rFQKQpuPx&8o_WR1U#32PFXdf^ zcekx~tn=ju{plR?tM6~HU+wYt8TYo}&voO5a{1iXa=!=WGq0zdd(2-RYm!#dkFjW0mhs^HJ>r43kzWoLJ1mG?5{Iugr$`N^e%JCL*^V|_~ zPvEVnXX?0{mstP*q})gPwS5!vCH;9{UjE;Sn+4;$cqZOid~|_C(4jKEzQEjwAl&{&_CqPu<3tnlZOxB zV}xHty_nCTzleXSZ+M{fC&aNhdYb+JULHs1@)&%X&%?h*#M54!g1hL0<+{wXZXSJC z?z5NYVY_FJ$q!{W}OO7j`Zh+d@lGRaRi^tzik>koBV2@ z`mN~ylyWU;tT5tRMMlmHWtf z?6ay@LIC=HKD|vlto=+TJ*KgqQ=x-!0@Q3i%I-h;gxpCO~ z*!%`N8HWdP_LZD({uT6llFug29|rIE$;dq5k-Q##i7L;?I)nOeq@O>@ezh;W)IW-y z(0|m!6~57Dt@b;XH>g{Wy7W7Pr_&z# z$LDxX{{SD?!O;!-D&T3${*HYG?fUs4akrA|MP8ZflCO=>*SdWd{Z8In%;)2N`Xl*8 z;Hu6$HxD$9+n<+L%zGXGb@mB;Ke9ykD zm(xhm_p!fX-_m-*{tP@D{pjRh{BzrQat_>mHIGZ`*?y66z3VlU=UhN)zrkr2i zz)ylF5?|I$+S&IE#9QNiLC)D-;O8OqO&?<}e^Srw$?K=(RqxcNGvaOZKdS12r}AiD z=J1)sEB%_eykmZczeJxc_i26y{~NF8@`m}6c^22BKbL-TH~O|8W`BY_W=Z=nKH^e7 zXFuQki#QVbeC`LmTi1@rBSz%uk-yY&k>}@p|6-`e=!4{S-`|n1&pwyeo4&AHtE==A z^c(%O)Q|c;rrt3wk#F_rWd?q?(8C@$i~OU`6ZvnIZyVnECv}o~IqBP^9sl;1_4(N) z_(3_9gZOb|{GOl>aCYNzK6lpV(et`>Gk$VF-8!{@XujsW#&i65c9eaZwh!|r`-$*G zM?CC|bC0Ba3$7=hnBV>0kbGw!48J%}eHs0P)Cc*C^>O6I^|^SzUGFd1UzE4aJA7V! zd~vAbqrR#~dAA_vlzAdQN_%{HsQ8Ndr5+3Mt(Px}-_hT!{JclJp6RdF9eY|aksNB zj(nm%ca{D~Z_EKc(kuRjSnaJQ?Q{fD&lvrIDu3h;iIcD7x$k9tpBe8*>i)w{Ic_QOU>`T$SLQs_ zWzkot>ji(-=c2wz{kg919qpsai{etdSjX_5=mz=a+*jmdCwjAgK4+b>Fw}iVpDFXa`lOWE0XvCv z;zq1r=Zk(`l{fey?Fs%`eB}BrytqVsJJ$ZJ+v7LUZ%Mllcjorz^ndcfc>kX3+b6t_ z{oTWnylDPgsh@a%QRfZ0tGo;FoPIz5;*k2Vo^_w(*Gu$+CcmBIg6|S{`WJQXmhql^ z*1R(KMm?n**gNvIs<+v1BDcJPe$Lo$H;(?nkoKa^*(WdY{q{f|(!IgF)_34F^2NH( z=&z-ojc@v${Rr!`x7g=Wjz4eSD!#ype!mYcBYviS%sb7$f1Q2V$UAf2+V5*QKKi9; zSMz@HXno{+3+rt3G~+$iF8B()*7c`;r=tGMJHh>deLnM6^0Ynl6ZuS)CvchaT93h- z@f}6#6TKVXqfe3gz46O>-~4U%26}0z%ixuGpH%Nh-BfXPOuP+z)$7H0yo`O! zgPa$AEf_bMFV@R>tylF=`U10ihV%#gGWuePYvs2OA#TmP%%AKJS_k4sdyMxz*50W; z!7b~cE3OO-QJr#oOtnlz3aZ^&k#=z*!Pxw zAmZX4`AodmO?l|IM87zl!`KpVMCEhxUCQX21W6;}=&xMbIe^L+l`GWl^>p<5v&MY&YIxF(*lv6(@PaR)E4*j2Fd^c3ri|-Ee zxi950`pbzY-lsQyNAK}HLh{kqa^CpfA@L7SY^cYZy%U-w+zG|#gi$oEp= zAN9H7Ir>@LiFM}W8S!0A&gVXf58L+Jj2E}^`@rDQy1({5fV?c^ zsr!uY>+1dZPP@+gkG_|v`^U4IXDbI zNjt3M`r)^!Z}Tt4b@XXxyl0F&JMBfBYxxv@8U4V-f&D-HPd|vhW7Wfg@yhpP)D!VO zOZ#iTJeQw;m3f=*q2qgs)T4Ej`N{7NdEI_G^G9DkufuP}zL9_LK-`M{SJlrkei`q( za$nTd@!fywPyRlGeXTc7`d_Ip`waS({dDypuD+1-EbG6F6ZfI-a_84uhiVsb;F0%B z|ESMX<%st-`8;`F=SU-vmD;j(MLIcF1{+t6xq1M*W}9|2FH{ zFGK#Qhx57nxo`hj{w40j_kB4|K7g0QpVE$NDbE5tyi5N$@LZpZ`nKHp9iq=U^=N#D zP_G-Wq8{`-pEIBMO3oAcYP}!%VeVf(J2@s^>{Ef`xp-ur=tEYy?Z+<2qvUg;m)w_l zwlAo^$_M6a)E`IYPx7~XV*TVezq=&QX?<|+?;5Y|E4q$;1dlF>TX5s|1aZCUA-=!O zecU_L!nluspy}i{_XQ<_b&{2eTn~(^!jG#Q@qKCf8TTjfZ{$mLe>3aB zyN0wM_03uHG}ei}W7@@j5B6NQ?_mExK4(7>5A(T`dN%Q{{N@=aI7|JZ-*}H(&s+cL ze(YqQ1^o5Yr(5<-h&L1aM7$3Tyyrf@f}OBeyno5-`VD;7@(}zX_LJk(vwSq+I(g8R zc@*z2n_cfl?|j!cp?-|_RcR0TM?PF;-b>v79Q)_Hi<#d^TJK}uCG+Rc4LR?4sCI~W zRFC>Ye798Z@1Eo#X|FHy9QAwaR~vUv`@Q~eUSvF^K8t=qu1`J@{pN~;s7n){_GQV_ zkHJIq!4k*DDdNS*zJdO4AJF%8)~nhw?M|_`H*S!ue_UEyPnP)+eN`?9*uzVX!j z93Hi=)SJJVuh92vb<=&UXPlM)zRdpb&GnTd@|fHo^G)raVE5iW4gIu{n%&P`$6kzzHc}FTfdN>F3=x{`Xu+Se$<1&Z>4-& z;C)0r7~d%+&dn3_(-Hd-Kj?G8d`f=8Poket_Y-~NI$zXpDIaku`uz2Lk98*Wmx_LI z!26ir#lG?q@wa`Sa6j^bqv%K0_3V$A&+@$cV%=l%iun7%9@hCC^lH4Y4}jix&D+d-@5y!JJyq(9x-I^WA@v0hUwbp}?+`EJ z{bS;pd~`O>`@MOr?+5ME$^-Q0M&^mu7l+x;iNA|ZIj~E7FIDB(($40!#sT~4EBRdX zpIvjvsLJ{U~zu z;C|=Bygz(5kJddAUrPNga^$_y0B#H4R}|kT8NcOO`}$*^_q(G_&K<9u_BE9EtsJKu z@tsDw^E>xZ9t*r#FUCJj{&%?V4fN8#%XRS9G5_!#uCuSedwp_vkL&0s`fBxhd7C5l zIdmL$YkE1yIG=0w&+APd^@$ye@4>j$neb)z&NdQ_m}mS1zwR z-g!F4f#1P>-^KNkNAaM28cr|g{a$;@GuA!#@jTy8F*B#%npXNWxgPzn+<1cX4A%F90>mlO#hWf$IjYq|U*TuDVeKzImJLFvC zoa;O9bNxw<@RU5M{fhloaey5=@htB@kn?otw=6QRa%z9Z%ga6f){dw7NomFPkG1<_=F@N0>*w=&zrqpy7I{v+J;d|kLcE;g zy7m*#@8$L8A6z$H%8QYwmOHpPp=2dCx>gZlhY z-gln^e%0*CbK;=!s5lY-@&|GoXPO@?KkIcI$NxUpmrv!14xAYOig8eUl>9Y6V!s1^ zCmzLp`9ZEzI*!Dv9iD4?gC7g?wCjcE%!ja7d_P<24;)2ZT~D=}-}RRjuA^V|ZM>4F zcF;!C`U6X$RE)nCM~^L9CJ{P#9W{e3Bq9rwq`AK!1}ef>Jm zO5=O&?fmMG{Ls9z`48(9{IJDMpSMnWk#!&EsN9}^Hm~ciWsz}QZ~P1%8qU=#&;N{ z{>(4E4j+l59e5%ByURJAI2Zpv703U!E9WtvZ1olSob|HLH~$77 z@__hnawA{-9Z~K`+#sj+w!RYI@_&zg+k9i@0q?!>`02dfuW%n+7xf@so`9&UMYt-#$OxnGk+0x@&`O+yfUvV-q&xsPw&h9ijP6Q)1Khke5c`?yrKhd9rEn_ zB9MOuMQq@oVbto0_D!HU1z?Q z2N|ESTQP5fZ=1h^1LNcl_G>s|KmK>v^>1<1#!Zk%-Vjg7+4ulH>Z}*EkNzlMsW0_p ze!#w(f8mGn-=Xmedlz}2YAgD9o}2%1y&cB6kN;=9l9!xE{}hkfAG{ak)^FgIPQNhj zfEVQ(uzRmw%nMHG(>_&EZ|dE6buRbQujElZHGW}z^CI!1y~O?aNBMlmI`%uw--iv> zF;8bd%`e;@SJmfl6+JELBaor)Oc4X)V>xgITxn@t}O3CN7Q(0j? z{G@Z z9sbt;G}mqO%cI00>xyIr%rxe%y|d>AAK~t zz)!}j_8j&s=(pv6h2N7oZ@-eq_&b|Ye-GqQTpn)bdGp_Ib00s~e(n8_>n-veeA@2A z8%6(APW$58pZz+2<1KmOt(KI8GX@;UKRu=^N3?#+YDtDOgZHaVYV9{t#O{V;NY zgO+#rJaM9v|4-u=yW4MQxRb}>N#j*nK@RLc#{aEL)U$qSzVG`F{j2D|@(9NbaCbUh zH9do`ZWB4ML$i~4gnHK=CC4{?yFc~KKIJ=cR9HtpF8tKzAa8l&^if=a$L9Z^?Q$e~ zY<9z+jO%TF!N4dJv~Yxy7^og8UtPz5}lnZ9#LKUP2D^Tk|uFwWq=13WHHtb2O&W?s~I zL;LCf;>CSXPgwW6Uhlf-x&Kkh(-Eh}#8LBs#*gwaxYhr5z>PfB;tV{Z{|@-uaE1Nj zM>`*32kpr|Hjt|akDY!eKZCc1GjJ_l#djX#4YFAzw_O7xO4D z?%;>UYu~FeP0gps>%imv3hS6}fd6hY zueW`h=ZnL_bMQpN$0FDDAN`4OO&1By|ut#saBA%Lun8!9em`})?_-ipv z7JlHn`NI>e2Or7{=HvDel+$?b`&aWG@zHHCpZ@X)=f(l;r#;jsJlFaL<~7u1J#t$= z^qVwpX3pF4 z1bG5Iw0r@7Z+UpbU&){Nmv1b6CcK7$H|)=co4tNr@efy zP~=JEE-T2tnMdmI*XAMLt0kl(;> z2ju1^0oSZk^8J8 z{+eCQH^A8j{?+Pt^wQSVzulj>lb@f?_1gT5Q_t=F%>y~E@t5<+@A&m0{?zS4WByib#F z!Kw0Ze2~}eugj;#fet?{i;RN@dgFKD=X$OW{;Y?a-%uA8>XkwNb{+MgoEKQvy!CvZ zH@?DO#7nOpuxE?o)Df+pjozF`{weYfd5JjDe_Fg^KK;NrED!EnNjV4fYd=su57?!} z4{-)w4ma5+eplermOpSlX8Z>075L)f7VCqP*4O-DUU%I>yfVKtpHaR+y*ZC|CJwZ? zuV0{7{io?=k^Nxz=FiyKzFfiHEnmBv=U(K<{@VR6ueW%mT;Qzr113BV-o(GSF5*kv znU6h^&voFU)9 z^Muw%xi7EV4=nOM{@eVBIPZGagZgjtYj~mI7C-FuvlZ5JqR;1f!v}H}$H^1&wf6D6 z4lvBgQQi<`z<@R9ywo+ci|%g$z3XXxPr z``F6U>UZT?%<=6#-!kVlvr*lNm4duaz#B1jilY)nj|Vp`VI(uns`Z9{*Eb`F6k#jSo5RLq>HVF;kA5$%TYZn7wV(cKo!Rg@VSgQZ{g89)(DEeZ z1~=NR*ItYHocodg+kNd1F2rfeKbhZn){BRNT1e{FR*`)&F5i|mX2+7BuCVG*y|{fQhuZ19|VFZt{6M|qOVrO9i2z;D!-{?&;? z*W2vs6l3)x4wMTV>Hm%Y9jATeS8<78Ew6MwaMI!^^AF&7jDLz-`*GbY<&-xYAG6NM zA;-0P&wR(Y`q{j{fqwO0_1F9jJBXvEH~ocq3OLz7Zu6{Gw>!>wg?@VRb&B8IulQ>D zyzhO$eUleEihtv8n`ewZRyZ1u&vEYh%Bg+iC-v6lc;)^ld}Y0e9gFc_oVah}iuhJ< zuBToeWn`|ibIXzIc@1mYTI_n; z*$r!3?0iDnl$06Mr=-lBkv?%s%9SbqAAbX5NB-I`jy68?`Go2->JmFfuKT3ul!`?4 zkqZT@6JJcQ|ND=69f#K@?AIM5`!z)2_YH-KiAVnD_Tt3%hi%W_S6LH<-=8be^VOA6 z@sR;7pV^@-hE)-D*BU!r~duw3)rQI+_)UEX~|qs|}C*L)AI z)OkgsuTa}!6FJ8bH5-T%o-t>^G6-M_9xzpq@L2z0D3i^BT*mL-{8wblBg-;zXeTw0u4HL4npece&OxWtGPHuvq8& z3$?zxDzxv->~B+<#y@w3)>~DmeLq{I@mCdSod-&|uOt!pb9J%KKVGJF4KLI4$LIC( zA<*@#BGGGCf##oDtovJ5i$5oE?7u|oyrNjw^SJ(2Df5?$ei_Ke&PMI;(hBCQ(0%P? z;`8)!$y<88uCFcC_19{&p5vRiZ!!E^-n+lRzgLP9N7#{nA^W_%M06g%T-S>V#fN0% zZ%i5TS}pkwF{^aFGFNmSSCBX&_G_-Bq6wfOS$8p+Y$vL!d4RWe_d*72=I(J8YOd7Bx< zj{FrCi63`VYCo4R(K=44io$+s;^i zzf|{+DiohzD8de4Kg=!$dAqv+In0CJ*wgpRMAx=un&(RN+w5AfA8AFp|FmM+v)#zW z&?4#0`$gLK))g9mS&8IobfM@qy;$;b1@iMijrhN^P}l#xLH9pXr}GWf;@3YLH2)U# zwY5_7^}}9VUMqgjVILVeqStoT`$?Ac?@sJpX_3y8R*LVBuhR3j&)`=cU8UcjC`6uC zYF$4lmHh9DBR8vc-~8%ceg^tpU9IQ()a&@jY8@A^7C&yR7yl|6wEwBg#Fx*Ck(0V8 z$ltf!4i-gW+=eneKQbRW{7>+&FIQ^3X)CZd#iG}L73ny$82!U9*>k}!78dBa zr{L4%V(n}2>L}QwndMr?4cM7Ck*6OVmABxhGE2mlZ1#WdGV!HrC30M#`7dMq1bqLv zSbVrW9tA)DkJY+wXNlyk1b=-$bYGc|9IO-_E-lwQR}@MgPF}@)g`&?O?4kKZ;S2U# z2QE1kJJ+*H@^n3ZCI`RpH1ph!T;EzC`*LrI_&TUW^B=B@0zJ=emK+T!)%828B+oxW z-Y%_;0zP`9N$Yt7y(@@AzgpI@M0%XMO#CoA7x?@K?8<$s#J5#tqQmVa8n1Mz_PwxJ z*XJxppDLrkr*quC6(N@mqIWWM9MK~B{t~++U}C+`2dtZmiIFX7_^K|1HenZjQLpEp&DTB(s30D)Bpaq4a-9nfTHQ z9{#F8&mB`21-!YkR_CT4L2uWr6rES1kB!j(=Ot0VP2-ES&o9y2@nw?V;SD-w9QD&- z{hdpYd*pgPe(^`;;?LdlkURX;x50U)mqG5H01woo2PYuMcOmC*FV+70l}g^u#lO7) z9{FU2^wsn(*u$qQB(Elg{A(FFYnkg8__B%UBQo_BQOx=}wnY?eC>h?f#KF$@JEp6K2>`D zB=j*}t9j>FAWwDTYgeu0a1`r)xmfo3A2m_1Z*MSvd0iB6%>?A`L~v0?vFPwwspvVY zjCHJ(e|#D{c?Q~Pnu*;&Tf>y9g02|BZrq&MX@7)bBJTcV^0s&>i()m$>oN5$P;wSD3M&A$NKl6 zzj^3iJ9>8aDy^#sIlIdJMP(#AUa4`v0QaYmZ?N&h@88tx`5|@UcQ$@>i|af1@QDV| zZyo;XMV?R1h5jqGkNqBJ7V7?fOC`r&=ZUY+luK^=G>Fc(fM;{6_4`cZFW!KiUKIuZ zWcV=X{pD_V!9fxHGdmFCtAXg)uQlIM*jiNN=iD(ri?^n4e5H9rvK>1O1ivR?ar zxIyxLU`+4);v;zLdg8X%z|9-8(4Q5ee*y8?)%Xp=fx#ankq7z*e6c)Iz%!TS>$xf9 zYf6ZhZU={&eGKz&1s{HizxyT&J%L~QR%jioR%zUq$~3OUAF^lIt?QU~kK+L1hEXed z9(;cb`H84H3Vgq#9J+fRs$TQ#sn9&{;J=>1-fgKCKiga%IPY7j`wdqrP9SfSQlNEy z2JX0me9z15!|11c4)lDYTzq|CrPgt9X)oV}ugIevz`p(SsGJ5rQ(7&4SJa5#h2XR}t@B zv{ZaKu2B28bq9KsuN2?UUk3lk$1Yi^`3y${dvj@}_*jS?U0SGle!N=pcs=%Ro9A81 z^xO{oRa=YfL066D8ws5Xkhiv~DD1aytLX9^c;@HC)61JQUw4l5r4RBSUj|*kA`Lg%Ep912hA>^5t)nX5j^Bv@k_MxvM zD)ju^8qIe&7rVMh`kJ{?`x}n^dIbLeV5R8wB6-;J@%Q&Pf;(1dUqv;V|3}!xhyEqb z7z2Nj$>$Z3-+B|gc@6$_+%nmL33a0DEco(iocUHsFBX!ATU@RC^6-zRAn&vBhvo-^ z-+I4V&krX4ZLf@iAKh3hx(vh~q^y$Mc9v@2ZzJD7sfdC;e78(`(Fh*@ioC+_8^rHT z_2S=SwUXPPVXu~!NUjHz>-V!(;xAT+Zhu;$-)Aq?{@;fFzhobqvvvJ`a7G^a)d|eM zi@5b~)G_V^_gP#aJua8LSsg?2M!n#E^4n{%Cx%}GT`nLGaA}_O*Ya|~KRnAmW|#Kr zX&CPv=KmS|zXbnKR1*dJa6_5Kd6GQE@4#6NfPNQ8 z!Jb;aJ@~Dpdgi6hFfJ~B{uTTvLy1;>+uV}$4_@K?+Gik@3X6-fCuhfE&hB5Tt5aq-+^9Q9WAWG?6mSN zp4S5RnthbNC!aimIPx~nkHCj1@Y&*~u+IwOlNs2PdFat&k@PB!c&`u~{ulN)Wtrr& zu~ByH^JdZY{A%n}h4l4R{M?iZ;r_}dU0+0gKB-3bX=tPLX^jN$W7H#eBX^7I^;}hj==N>&;+k6gWtHUX67p|-$#Whi zo?qG^KVIp2kSBZpYscU4?}zN8U!L~&KJ}2tSIdrnSc9J;AActK@9b>RcQpE)l8YWv zKf9+u^QB_H%wNTh{Fz?~c5e!H&-5wCS4JH;t^vPPAU^dtZiX&LuWwX%nzC(bL-^M(_IXFU&9 z9m&2@r}>O}&nV=0UWNGjDfP#X;9CT~SPSl)RH^a5#7|yAKK)n3%U@GB{$FM{7Z4Lo-#c&=fE_;ee%-Rk1P!Q{g-u`jPq75dE~t>s=_`~=NXTl z9jc3hzq5K&s6RXm9_m{ih4Z!0<-+PH-~q#{Ay54ZezKPM;90Mipf3+oCmQ7Sh9#2s zlqK5da<4~|@2#nnyzVQHg8cn=h3I7QoZ`tU(bwpzyhL2T8$J*F{Z6&UJ&!ufxydmWIOGB3=_lj`%)SIUe33fj->8GGZWdqe zB@b3wqw#O251_C-3j5pGD7pTCac`&;Jqzk&M>ghShgj#E#B+V9KUf?U@&X+$clnyf z`dET~Uvo4c73ArAm`Fs6QP~oYP(}I?Ro0+zI5%2VxJt0JpwE{_A(x zMeAz`c5s^GhBEE*1?W5{OZ>hXJo;<57d2Y9^<{+T|5mN}&&3aYhxJ%}G05dY>ITUL z;)BHrl2>r=tyR+B!NlKpRf#W;GR~w%$(hMTc<$Xg^o)FIYqQ3+{;S}pZt}h$;AW0&@M@X%)r0>sKN{q;x?cY33gYVX;O9rwQGQey2@m4mF2&B> zutM`*1WtUEyyOQyPLAmJPvq)1;40ILFz*QX@>lfS>I1TK-fux(ViW%2KIS(+7UXIR z@%<0VM6bKSb-Qy!&vNACad4E?J%WBOcD?fYzSn`k2d}S|oLK%T@P8xg8?r)j9;a^h zRb>?7hwdaqOJs|H3?fK3Xpi^f7!B?*Cy9 zc8c}XVehJm1CCiKxv>5P*#Yt}-SlTAL6=KbW2eF8AL7?8`* z4Dh}y^3Zo6U+2Z6V23iXC+~Rv27Gcqar<=gOO|&G_U|_OPhMUjfBG)*^Ud(}1@`rp z*VmTF9$4Ko%xCqlV5dK3yz-hT_`z?J2hJ*&99-z}myd%#9O!vY^m#e_sUT0i2fB?X z{z~$`A8_1XDn;+-z|H?tMqCX3dW`+No`oHA|3SU&-?1x}XAQWs1;1te3}GLC0&kk$ z3U3rko^M3|P6S_?T?qOx06Q}VKd_@%^82Ref$-DT-xuuNRrQhVC~?4Rap}XW;Hbeh zlJD1Qqp&{f>k8`_iQT=9y!FG4lAmjq%TH9$x3Y%$$B{>|xSxupl8dD^@|%XAfjsZ>x`n=hrL~9jzZL%=bs~W)G96?sj`dzRBVf<*UJ$*OS+O1brI1RQoVL z5cXroK;L_-qJTq=Lp~3IE3c`Se7s5j-GY_i6!!5ZIAcRI^5=OE>ML83Gn*&WQ4Y~> z`}_JRtn=64fz3zn3;g+6nc{+d^nHC+F8%r|cCFm=4y(kE(K*tGH{;NqeE(SV!16a? zzk9*6Pou}l(AWBugTE`LpQj5tOhW!gg3k)c7g?S^*xd^5o6C4{^k)J6w&T#Nbgv)b zhZA!px9@=`3;6v?=-}(TnZAQRa6RF9Ch+D6?3MMqhWL0J{?q%OTj(ROcqY*Ggl6n* zqxiL*IM(v1(&I|W+41DV{zAX%vlWuxG~&)S??a{^;_@8n-|hIX3G`bHr$6K*k0Z#J z^(Edqk2+g|d`}sAeieT4oJy_h+B)&ipu_tO`vbj)k+VO8lu3b zOX2@!^8eNJvm1^Hzkjt_*z;!ikEP_1rsZheecayG$PS!H{qJV)-=?l`b+hO+2|KfmJjWy0$#vyX zki&@nr7Y~iOUUtkHInB`h?_?8{=@@br@-$fd7mQq;CSlwMWyU#rQ~Ebboe8Er`NKV4t^c6T^&RZo$Gq2Yj@Mtww_5+8Fee6vEK+TK8AYwLDp+^g@6w}&DH%cW8V(17X2SDj{<)mTrN7i zK)=#==)1Chil9f=^1RuHU>`ToFS?2PM;f?h5dP8XEn)v=XTtp7I8LjMf`6DBGXpzf{hh*9;LB?2i+3~biQu1jP9m(+@M5sHbIBW9 z{21c0ChA{*%8~wLxu5oUia6zd$6Mg6>%l9(2H)RS6NT}Ha=xlr`uS45{MRDpYXdh9 zt`$BitJV5mL+}4V-11U|@JbYyomveJp1f3a%|-w2z_0Zbh)=)AK3E(W=u-=w?(%#W z`S5G$Ln}l-US{5Ju^yMBg}g7anfTuF4$7Yp*NsIls~RPD=X-v^?G*lMIQh;8Y9;^Y zy}Cl-wVK{GCsHkx%{U zeQ?2P^q;k%KM&JiyMuX7AkXp`aijSm)t{Ju8FC+63SEd(tsf@Xm9aI#<8L<#$N#iI z&kymwM*4o@j(@3FUI3oE!t4IvJ?kS3{%H<=?+|g?=jiKp{Nr=b%iezq@;}{iAaZPV z{E)wR4ZSgc5$w*%$oI?iDa^p04dVHS>4!UjowT}^a0u%i!Fy+ep+g_?5jXSRhrJgP z{AEcsewO`MeLB$f!bZ_+4t1%YVecENWC!Q>9)(BFZ|v2tAeYIN@Snc468ehoWgRof zuQz+2wBud;bq(*CtS5g{0bVJ=ZXTfiy5HZUAP>|||Fzkjpm%psCoQ61+4@ogA2wqL z-UC1G1@BuvF~oTbu(zk7--oEf{0DqAyBXr0ua=08ZRArT^k)_J>q~IoM(-a2Z(QYd zC;E_I^LqcqQCOepy&r>oOs|TB7sd_sXKczh6`;o~u})b)@jzK>Yi3&qsi#^W(bjfal?fhyLREl!7SO z{g=s8KEixy6nxl+_pF|D{p7vS$CrpczslD-tZo#>?X_S2z6JLEQ^yCUw_J1X= z+KfNxYCf|dsZ(AeC$WxW;=GPki7YSfOk&7Z=OpY*!us(5A5Gf=%3-dkTs`F_CXx&%SNl!Co6Z zgS>ej*!r%5y{)W@TuV3 zJhfc*-tdFsACJ41N>2Yy{W8PjWa{YlK4;*60r~RBu)lk#t5}^P_>nKE2cANH{A~2- zZ{V&t`h6|+=1shJGMxVP`CbPjAG4YK+CKU{p2k0Jr|;!`_XqgjGac`Ff9wj$eRf>* z{UQGU7~=RplV{nFT^qVW_u2c!!Jhxt{X6nx_ABu3YVWTi|Jgvlj@3nio>?6!$lEY~ z&pIyLHa;2KCLe*(8cQN zvWLVGk2mPKbkC!+uI0_LuNP3)vV4K`qf&l!5P6cB^fh1S{bB6;Wb)dF$zMLeI^MwF zjQ9KpdGzmNx2z97?4t{OY;jFk|1I!&CiUKJp2r{#`c1jwlymE&U>~n(hyuSRRSS3R zZ;k@JUT=tkJ-oVEzYnaC9F06d_VEVl!*?%|oZUZ4bbEce@btqOiC_orPlhp>z&c4>zf**K>xku zf%YwcPWTOrzde3aUW7W}wk6`%T{)89LkqQDdoM2RV?q}CFjo70yGiSLg?_^A<7Lk> z(lpLO`hdpgX`GH6$(!+6b=+KyXX^=a_)fOgW#1JDxOPsI&_1%FpkG%uiT(@g(OcfH z-PR~My?&9#d15i^oTGhjSgiG2#rHh!ibT)*$Ui)@RO>mA&Gn_)|EcuB#Lh(j2S>p! z<~M0QUoX`9+NaB}EK3vpw@iay6UCQDnzfE$O`_YSHKN;%ad3B)YOlDEt1B`l{y5HVI4g7ddBq z&aah;uICUBUfQVZo1tH7qxSPsh1UBra(7vc=rW*6@!|8eQQ&WTgZMS5S$b7Kz5iAE z1h+0CZlsUy8}cPzr)$4M(zUPg-p4*v>wPs%bXwY^^_d+H^VK)%`e&u$LpJZ>ZJ+nQ z?^`=o_WG~&?1T4Fo@>xJg}it9bcEcK_Z_lWa!^_r1;70@^IzI3obr6W{6Sis_B*12 zeUFwrA3$FRG)4gz=FSlx3`Yh(lf7Kyw4Sc@{{y)$Zjd}&+ok+L$lKU^Zoyu@fgNkeWPfud|A!*2?^eE(^ZTqs?8sk)owD!I z1RUR=ez^x0L$547e|e_pd^&wem(SDh$1N6L^CPYE62?v8{fq0TqHnpPXIqZWZ^^>0 z@E%O=0@;GjM;`H7(ok@CdI!7%^sLnM8>3HyG~`ysQm-zAH&d-SuX7wY4GLb8Cuum`C9+{Ea~gINPHMMPjr8j^`8e{uU#npE&BUY!9PFIq;>qF zMRfm8LlpeZd%Wjk^zeq;f4;s>!Ox$xylw{xlD z-=+r9wHtgqY@YZweSSjumWA5aUGVqsya(Edzj`=J^7Opd#p%1r!0%M!&$i4JefRsm z8^pH^^8D5p6ZCydj`V)ZEYbg`*x8ZT$NssAAg_sZ?f25zTHjUVp?dI-DLKLkH!c=` za_8!N``?}lcA>3FcJ1Lb>C;fYFO{DL-%rx@p6TM}MJ*@F{X}6 z1bmdqeoC68XH$^-_9pq0&zv7w%(p;%K8)Nvh<@CeDS2C-BRRd_^KF^37x!gpeVsl(;IF$#QrjlRAYR^un~wC@RXW#8W}ki8i& zSNDI&`znvs>-V+I$Qj?ydVK=?r7kgfp4L6JNPgp0b?QZqTjQ%vNvxuicdGTX#M7=gWX!&BENfx z_klMK72nUBDSBPhpyz+Zdqa1Tmpix^y+?n4l`Z{wFBAFq{$%?4$1IY*c1OrnzW88$ z+Uh@>g*?sEejng{gU9nk$F`-i3v2j3#ep21CoN7WJ_gsu(-VQ83FP>l`H5hc(|KQI z2Kt=0L~=hOOLTnZ1kwGN1=`2)EqX2o95VxXtipeP&HF(o$7N4vfV0+)5H9N&Eql2e ze6ZK?BKtP{7wqowywA3KmiA-%8u0Z3=y4JIEP%e*9EI zPV_6lF5fj%auDh;ZThC)NY^_0csw#*dU$NJ>|3tO?F{_xBF$s+ zrhFjZCtT^@*~t|hF389JfREl>EIFT)Eq!@~@0^Xyl>VOFBs=jzz36>^z2q$3Bs+Ot zru1SLdHTy;epYJSm;Jp}`s?1N&*v-V8?;#Kn#TLLzh5kz*74aQ#e-?`Q*Sp&4yH`L6Lgj&O^W7d`jyduKA`G%kKT9QR{gM z9NLdQ^b?w;Pq7-=g~mG3|3dWQ5%^LO*Sc%5&nJU#r|{jb+)T;;>t|}+?;jMYoxyTGy|c|NG72+iaHyaLAjye=_VF$w~is zVLgv8lU;pdzT|J?p{D}gT#tVL2syu-{KCjq*`Ml0>BT>rb$%yt>CH2w&ujb0zg*lP zevRh+pa$aEyXNTkt;Ag>M zKUn1Tg@9+CnJ4@H$ztpS@;W9>bhY>8#9!>}68D>}Na+`$<0KjArTWEA%tGg1=9mEWYm^ zD?Mz(|L>Tm>*Zrbr#r~QrhyO3=+co&B3mde~ zyE#q-4?pI5OZ>ioJXFy!;;X$cqP%IA#(QV3=zaGB*=h4b!5+MqA-x>=ZSd1P@oO*m z%HCfO_ONg<^0S!zd%r0D>IMA4H{i-3Evn0nA7@m%Z`oo{KX8>!SWKpZ})4G zoNTL=TrTl^bDjMDb2H?R4JXBp{B5q0+!l}@IhK4}ChOntc7r(LL$@dSiST?%p7_$2 z_fzb>k&xfNeqJKvHS(v4&cjlqUze=Z{*G&m0-mjg&iPrA_r2I#t9Pj%0sX(NTJ-o9 zcJ249KGLF3LO?-cF zHTdP-7Rk{RaM~ZqGoL+6^zNJ{eNH#OuBJg29uOrf5Q4xv1w^NsRb&34I38O`~ zN`BvQrudaiykvRBAV)7`OTIH^B?5o$Pt$P%`NcPo&-WoBr-PI4A?|)0|5^ghpS)0VzU4&bC9m;jv-Wig^gh?|6S(BVTH*Z#eBU#LKAf)k z;@4O6_55D!{&mgb+cNBTlqr7oAwSsU(8_owN(EyVp6e};Kt!?nK~sN4Mp{d|%7 z$#wXv<&pSrc>?JX`LOelyW8oj-RpROxc6@OF=8Hejedz5LK&@SWA;LLH#>f1V0Qs}lviA58qZj=W

IC8h3%?v~2tl^t+inYb+OeJcH3@eC@5^|1Zf$+4m1rH%$kZAXg21S0U5CD>XxO9@Z`Y^h3VS z@yVH@^Xt>1fJ-l@uCoz;cr)Kgs-S<>zK5jwz%Axi6kiP&U*1A)&LjT2vjzUt!Vhrf zO6uX2IilkqcwhH=aD|N<RoNw=ODxO~^IWl=tJ!g^Db*NhN8-58m<~i~){{&Y(L>{}KRq|_bw{RWxtzl0#)bIYj5FCLLEo)D zBzr$s`uitv+tqw;aAH*y@|A&bgrU08Wl`X*bCVcpcR*HI_< z0r`Rh)MZ;o%5RpUhx^9s{CCAsuuBg$ihn6B;+y5!f?mAOdq898GjF5s=jTn@??1gi zXr}gM_9^h)zKb37?a9Y}xZ2pDt40 zY@O_2J@PT|OyrXOrVFrBmM>TTMVkD;+zip}Q~wSl?{S~EC=v2i+mNq^h|8?b5!U^@ z;}EX@4SQ^P`XEmaIc~+iJUK`Ewthv~?OF0O=S{=@%$6MNbARmlK=8w<^dENYd@9(5 zEc7CsyvF;@vhx|>FthVv9j}(lu3yUcV8^;0$$+2q$$Zr$eYlyvyGq_4eJ+w;H~kCq zHnOkRsQX<>{AB-*La>Kc-&B4hP4r9RJ>47lzCcQe_+;`B|j!@ zrNB>`p}nN6lt7-35pv=5?_p&sD0c+UhEF!elhWI z58wanNXO1#U+<>=Y5#6Tke^}XWBa9|cjN(@sK;-plDy6Zw_Dywc`Lq$eMKbPvuD2K z@DTaa@yP3|;HUS<&p(6R?#NKy#rhn@U%v0QZi(!{p=QxH;dRzo!Zk%xw9Y}~`=;Zs zp3T!f2j&o8`FFT;u$%DtbmaPkNt%D{D9wK#dU*`-%^L7(C-rEn7l(Cj#Q&@#KmR=S zhkD{}dw(GKfyem1`Wy7EJ#;kB5cF{keGS$(BRi5N`79v5yWjgvsKbmLuYIp{{hq3M zu5CiDk*C}m>5qMfG05*BzH9Lu`q~EmHa!b=cIpI;I}4m*|DKTY(8P(I;M4CBKUzE$ z=r9exWcm~ANPqO^Y_FS+#Xijz{(QDv{JfOBg7rly4kK@HHaPu7kE`NQkk4OLMu9$G zS8Ba?5#JksglnmTy^nocHB0<|m;Rz#sYfUAURl`^?e`VWkGr3lB0Ifhs`P3!`BAH* ztG}mF>rY`H9{1Q$xC{Lol^#i+D@A|HF9*B+3HjPCzNc6|L4LZ7{B^^biP(|9lw#fg zed3&Htx>pc^8`E-)rc-{tj3Pal^o1&?7c6nw`-dAaZHi?QyzLT@^jTaEglbe`g8g& z-zDB{VX%fCSphatWGD~@A-M~nB~(#{@?6vkh2@Jjs>I(*As_Jd9Pz)_@jiV~wce*Y zMZaG~zG(vfvB2xx;N9-|l80x>!|h1IzKxQ9*h~FtJNft9iSz%RJi#97oG*BvMNJgw z)jmsnTf0#B;!CdwdA@^ua@PXk0;`9F=k^kpRL|G`w^I+?F<gV&l!@>TNY~_ zH)qQ}w9-fO$y{&)ap1^BvJ^6R@~OyI5_%S2Bs=u^6n@XuF`0kg?*`(+lgVR#vX8#N zIU3jc3PN7V>Kd|x)9|Y&AP@9&T0cRM<0VbHuZ8;gmGqB&m%iyC)soZLhl}JF(GRP4 z1v~O7dT~m#5#OrxmUj#8-FIZ)<%EVgFvYxtcnm;pu=E%F~5Y?#&f{ z?+0&~pH)6%B=&MCdOjIBnj`vGQjfT#I?{N=ac{>X*Slh2KA-FMe-%q>+T0igJ1_%$`2YCfF#h*jHE-PM$omxe<23$P#E-U* zFwa2xi?+|veloqT>-|Q&Pu%7D>iH_~f6Wlz-lrelyriA&m{6Gk1mRYkH%=e2b(oc z8Fi{N$zvP@=X^#!{u}shb;5w_T2B?-VyU8kjQ5-NQTMVwG{te$6ZTQBpW%Ii#39M} zr^Ab7Z}(-$K8)dgfjf|^{@DM~g^|lcz)$mvB&WY6-+m73znOo>VH171Gn*u5n~v5q z13zve&oJ;y>I|%}e7^kLhPk5SNS1}?~CoAZ{6fK$khq>tFQdKiNuxG{~dn+8^6zS ze?@-qtSBM;Iau;M>{RK)alDUohF@i{>)+^I7P#pxaP>y=@K25xU6RiP=kpz#Ytic({MDJwl0VCz1$j#0 zJB(M-pE_)u_FGE+p^xM1aoSfz{`a@kr4s0ay%!Pe^I`8#13!17_n*(xxZeX$+xrW# zBY*94WOvpLQ@l{czl*t?_Xc)?4_d&D`-zt}xqmH}-5y8&c|CsPWcnB~;ZrVg+IZ-^ zdkAt4?%VBg6?(7-J$i5ndY*|qoG!djG)4A&5`A{d$p_l|hex(`s->#7#w!EG2mDgRI|MXpa<$6mUI2Am0x__Sndlk)){@VBD!u)CUefFCVp7S`4 zyrg}%GSEMrxPK$^YJKXW5A~XS8#~uJbUV+^7_UZ7QqoSFj z|H!ejiw)H0Q^}A07ja}g|IW!@$!~r|UCRF5`+%#K<)Qb~`L4+!Zl8rb@$U`T`wz;G zc$~0E&tE@R_x-#~`+XaGX77&$J-CuM|6TImm;1jT!uKLFSjWfo@7+bde3-dv#Z1MdrLP#m>A z1%6;hw>L?TF5~^VdE|e$@xI18x#G*fCBiX3;5}#i_b-Ee7&cP;-p_jiFB8w7JX7m^ znD+%U(uo&JYw>KX(J~Z5NQI=!-pmgSe%E{Oz5d zSK|GOcZqxVqpuT?lSlo%xv`qpz6%!YMA6&oXWutS_S@pOfJZiYJVd_iUEaSwm3;lp zIig1?ad**i@-Jsm9~?kE-QGtJ`|AhZOT#Yud;N={ zZ~M1pcdQN`>}dgZvW>ozdmOjT$M1|%TsQh?JQn!*JL+qfwrKw^65kb4XTRl4;hG)v zL$s#ozQOoA`;L0>52JXWbp!R=YyG`m-V@z@D)w~mQ^J+}yRI?%v$oUEc_IEgpFW8n z%+UO|N6;I*Xa6p3@aMbdNDr34mm$2rYWQ8h(+^b2`%``K-&4u2-%S2*40*{B^b1Yq zJ>&^Vk@5-jJD=7fxrq~>-c3Dl3OIXjt@Q4w<8 zbj@M>J$iqqt&nxK=(r)yJXsnqHHUSL(R0b^{9eLy7A! zjkBvn<7E1G5{h)+hC4w0&effW8S5VJ5Ki}`MM|Z z`*My8xGxqHjy=3cdNY%Kj@5B1>x&K6d3S-H-?K``xXVS6&aF>fIQ@9dUr?-d?uxVi zS&Z-Vk7XSz81F*e-?O@R9m31SI`3b>{pjQV8eN}tuEw8Iq~BBAPo1Ijdh(;|xqmhF zf`ZlhJu21F`uxNChqh|mJpa9w_0QHiFM}^#hsD<^oEv||f8?qypLJE~Sb$u1&(Qe> z=wo^LFz*>5Ryovq(Dou{MGGw>^wdQSUXtv`0Gjwx~7m&SV(Yuv6R>Gy=M zpMIMT?9o>EVDWTVPiG56gxjK$@%xBzej_}v|mBTn4XXv?txw^lpg>~of`*8g}1^;>YEbd2tTiw3m z-`=_aJ~&`H(9(pKo#i9B?XZ%QEdLphFVWFBy;`}@n>&$-O84TTAepl+bt@%9fesP7y>B8<=pRn-qDxTjZdCG1S-RqIt zt`ixjzplre4*j{FtLyE|)A4P_8^SuFOTiN6$L=RP{us>l<1}8k-@l0QN9btplLtBL z1~2Yetnv3C*RiRZw~hN+C-6M}tlj;E`>SJF=M2_~zuU09cRhhV8w!};c?otZL0#j=jhgqeGX1`ZV`r1*>3=qSXP-;mU(M0|LmgirZ@XA;7kZG&cxm9> zJoi)JR+A6e3H(v~Y|(lDSsK6ZMbO>H8LsQ*#{(axfD?A{ed7)A-|&LQh3*rHKQ^;} z!yV!IZTM5OU*Z0(k?0Y}-frSKv-=tk{M>wIekNmODNjh%G)-h=b{2j-*OSS$q z;>=jO&N~-JVU+q7_E)U)Xt|CF>V|2_tlQ;+c%vOZJ`}oJ92&-p9Csd-FVPLYC6JS$ z=*QL)&0m6?+WYK5j+0%!7isT-PL11^~2?^iGCrBhgs)!RgG?2P#@&5yt8 zaGXdr9#9OW4+mu4a|GVVP?1TIOc&^>? z8o1HkqYUfr#9t+r>Hhurqc)eb3$&l23-x<}=a0ba{m~D%N6-6vnr-)8Q zM_8ZZ_Fnz0*Z!@KF3>T)O6T?DKRhmIJ&EiQr>Eyr+zt%XJeEHa{~V{fU1S|*XSFZI{XP5dOjvbZ3OmrP%c=XKh_ zeLYA0U$}1%__!}|YBxC3@=sdVx3v!QSCZ>9b-r~O>su&#o4wKfUf-U=^>JEXC-M=) zf5v!zJ^Fna{;)T1=zi~Hp2to-56&sT-Y3W#b)BU3C9&?E-_`nT{IG7%JFFp3v547)8?I&{k3jK#-r}u(Sb~64u_%X{b2D)TIzbV8$owIfQ;8Fi9 zJB7cqe3JNzzuJHu_dH~)=O?f;9njPBP_N;y;+_x0&umyKy1d1@b~--f-#hN{`~mrY z`)+aA|9(S$bXKDQh3z^6Bq8acN zJ>TGdd?@#MTK2{JF&O%Wg35IqkgY)jnRBv*!>>(mh9sV7G2}er@JxI z8sq4~j4TZwbbKLt7JBrjBr|7inPq~XW5!rzX` zT6c%zYtMf*==x6VaewMa*SP&YPxJIRo+sYwnk_n8{a1E$8t2S!@n5hDJK;|&bw#s# zVLu7aH<35Cd|SA_nf!IqTJF0Lev+5Be1EvFzDCEr$mwRjm)GSudXS!5kKSFzcVpsi zca!z|^FA+jse6Ie?KpRn=M&stPzP$qU!@|Kme&e)(&Q=ZI~lx|nx%Edz?YtSxFRCo7nO7WbR7-E(Ndkc%6f|DBw zaC__|jn_I_*L(A#3%TBno;m)@*L^YaZ~Z-vAb)M~V(iGD%U_A}b*Si`L>zw@d)tP- z7~jNKkKczfkN1JNeL#>uw{Di*IF8$p-+J-^X|qJv4zJ7N*An1Q!~K$5^5Y%gAj>;S&d$|53D3uo zN3gg!{O)mB*JQ>ynr{rhC*cPz{tokOLN9y$3vt3ZXKS6mI7{PhTf@BQ=Vj#0U7rt< z7qUEBphJZHweOq;`ferOaXT5|zsz2Sbqrm@ec0W4{FcSDqQ_{Bd(+vvZ$0tcdi-xH ze%f$Rn9ukTo_o*f=JlN!y557G-G;um6Q5b$OwWT)FLQn%mljur-}AsNo;Ob8{g18) z{u3`7{t0yNLNEFs^;hC2dSdU%hI6Z<1iiJqv-ZtCTj6IXc4hy4UZxp6q#K_v_RF&3}dQW8iP=C((I|j_t%NR`(3oqa4;vy=@KlEDbxg z^=Lj&`}25!c+L87!Z`LH{eCj^Bd+SrqaQ0eC4)1pUnbn&>3KEmmErtA&yL~ZTYvIS ztsakjB)c^QzUMV+{(9uc>WlJM(7Ox#?d!DsPmqV!QLK}ESF+QU_{Hq8)|bwGqxIac zJ^$$Zayi4FKHIAGn4OLt`Lp-?B_H_ZFObtX{IU9NkR!vrVO`zC2llAw>3V`bSpQwP z-|Fe%{`J_KIR3kfeBd_x!+g&NrfVGs!Mhfph4IJKLyz&%k>$p$uO*D{eL$AS2$v0y z=sIzh)uF{F&$oLX41C<}_zU@IZ)DvrU*M6HMvZUz6wPz=zZ(A~@Fo7!>PO;>%MtxB zU88i}{GRBIei#lE9bGTbFUu>7-o8%iPbM$nd5a4)4(skBzvp(p)%_2+rq%O-@F6Ly zcR!lHFZ7$raiYdIzodN<*RH`&tZ{tre%$v#o!Z`?3FEe*XM2#}jx_Dp`Xe=Oveset zGIr$8@SpgIfAzS{`bon&eu2DLUNF$d`mMrs%iD$LH@khE#rW>GvHRAC5$-<+JDtKl z%nyg>EiVzqeYZ{bT}S?DB>A~+_H)qZqi(eZJeTHuAC4z;G@jK-!nn_~Kkuh)Vn0b9 zH=@7SFrVjd4QB_stO2i=pob}5r_0s)4L=3CZgTn%k0sEnR`Oh)PhIbQDGNj|@0)Jt zyNPX{KN|+$j`}O{?_}nkqV-uHYoKrPJg$$`zFWx?UITvb$<7DI@8T@T^?QoyeZ}}QTomXL&(-hdXZ&+t>2x8rn<{q%bh{@?Nm`W?U8MLy2H_Y?NJ-ti21Qo|cs_gB=_(vb)3 zeCrt2HBa+coE*mK=03xdl9ySG1O9!jQRkND4EMKE_po{7gkL`6i@p=(*YW644?+WB|7wh+U;$d)QU*ZGr@3y#C`~?R+;qjIG zY0n?ypN$@(`#6r|ee1!Uh9^WP@>glZA=WPt)@OD%{N4uL+<$wX<6``G$G_~G?H`G+ zBu6fH3z^U3cJ#RGIL&8%Iy}0UxXSbO9`Ci|@2&4bdV(K|umjcyqWzLrT@Nlz!Y*39 zTl+mt^A)(9xIPe1ZuPj`H!@hyS$#&&lTR;jz4v?|ajMm;Bwz4<1N!gvTl1@;3;hNi zZofUBmd)>{YCe}w>wok4!^*4=BAoB?{tNPQmd_U*vPDm`$2xcY%3%D_`h72UY8!I# zJa*jT!`P9(dhlE)eY-o?=y}8aVScOohW&2MijLeKvwm^OkJsy=!xQL9486B}V0hm0 zYvI_5zv=OQAmX19%=ZNPG8BKA?0As)Vhj0vtM7*Stv(g@-HR{3rR&@NCH`u4+%ZJI zdml-w$6X%3J6)kutDle3zReE=I;?jbHCW?WfA0~p#1dY=K+i0%EI9!OSbQZq%;P$J zLfcxPFL~co{FUdO)6jG4GYRu-M*giH7tRx&cN@U?9uJT&SPx#>!*g4~jhX2GP_IK_ ze_YQFg2Q^i1MQ4weTb58aPhYq*^m1(@-vnn(0++$>^-coUi**e(x2xBYPTYTkursHg#d%s2Nl-~KmyeS@E(N|-AYTD1>-ogXd|r>o zUcz62Q&K!XG*aV4S?q&+ZaeXB%yDxX`~5=WG~oxV9uVZ%>PcZe3C^uASNq2fCcr7y z4;}7H!C!4X>aW86TiNGU_-Z$cK5q9Nm&}1bo(~{?_qezXdFgJ{eC_l_SwEc(X@BfL zh5Smp`zx0Z$7|%JTjye@#=>84QG)!ydhhpee>D&Of~zC!?0j(BL2$OMMD~yC@uP8N zSZ43<-sj=(9XY>d!I$x@Ybg9+J)OjRj&t_7KlHflAo$Mns9BNb#qW9F(jN50>X))x zjN3X$^RM@O${6k2;-jEHmcP;dykBGp=V|;tLFZm~h>g^_;ac$>`#%M{Yxx8pPwRAi z=yZ4e$k6zmZm&+!z6@VUufcuEIm|;mGNM`cMU9&0Ap0_06`tD!9^3$5V+(YDjC~ls z4eQ%Ozn0^eP2@}C$erugUhii=OXCkcPwTPwU}R^I3wsYJcI0m}cJ4az0nax<58~P0 zxX9~a=L>F{+`$K_-draO`}e1 z@0&;-X0bo)lf{)`UkT^mIIY9#R^ht&3(3J?_DOxMojPQJuZK828Tu#VkF9wx8dC6Soe+6^WC2RI9B`dzOJqby3gWhJ%^lnJRb8rh3B!zTG)=wcl09PKw?>mm{A)nmaujcs~`s1zcA^A#WoY9*&l1qX>I-2%$u9T!nKb-P zH}<^M=gH8xsjSD|U(r17ue{!Ss>VGByxm27(}jIcWnRn22R?5G$F>pQx4K^jADVv) z_t}3mANC>D?JC!8y@6i6b$4*H*Y)#|uWddLa<-SehxZkGU+OO8C)52RdBZOE(?{cF z(UEvF4ZScq2?2-vbHqT;TXT2Ax~E4j%LIoj7)seD=|J zy*Mz3bI;!=Aur66)r3oH)M zQYY!y8%H>fBj0Ork?2YO+4~fdop1O*^RJ>C`Vc|q-nv*-gCJ=JOF>*?p25e=i-rw>!@CIGH?FEB46nTv)f~>yxM(b#lFfJd@dx zaG%wKB_|$V_h+BTZEyW{9&$(>()@am8_R!&V<-Kyox}$toDT58{bZ;6MdZff;xK;) zcs2$O+vagOxXbpV`%V{KtPUH_58|&5;+OW|pY46vK)*HMYwK4J{4(4tx>9fUdWYpd X!anT%wJ@&R-2(5^oT%q5|MY(W^v;zg literal 0 HcmV?d00001 diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 97c892048c8..92cec23ed56 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -188,6 +188,10 @@ "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 @@ -512,8 +516,12 @@ NAMDBIN = resource_filename(__name__, 'data/adk_open.coor') -ADK_DSSP = resource_filename(__name__, 'data/adk_oplsaa_dssp.npy') -ADK_DSSP_SIMPLE = resource_filename(__name__, 'data/adk_oplsaa_dssp_simple.npy') +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