diff --git a/package/MDAnalysis/analysis/__init__.py b/package/MDAnalysis/analysis/__init__.py index 1a0883882ea..4c6f52d904f 100644 --- a/package/MDAnalysis/analysis/__init__.py +++ b/package/MDAnalysis/analysis/__init__.py @@ -52,10 +52,6 @@ contains the often-used :func:`~MDAnalysis.analysis.distances.distance_array` function. -:mod:`~MDAnalysis.analysis.hbonds` - Analyze hydrogen bonds, including both the per frame results as well - as the dynamic properties and lifetimes. - :mod:`~MDAnalysis.analysis.helanal` Analysis of helices with the HELANAL_ algorithm. @@ -108,6 +104,9 @@ :mod:`~MDAnalysis.analysis.legacy.x3dna` was moved to the :mod:`MDAnalysis.analysis.legacy` package +.. versionchanged:: 2.0.0 + Removed deprecated :mod:`hbonds.hbond_analysis` + """ __all__ = [ @@ -117,7 +116,6 @@ 'density', 'distances', 'gnm', - 'hbonds', 'hydrogenbonds', 'helanal', 'hole', diff --git a/package/MDAnalysis/analysis/hbonds/__init__.py b/package/MDAnalysis/analysis/hbonds/__init__.py index 193b9a76242..a294bc0b572 100644 --- a/package/MDAnalysis/analysis/hbonds/__init__.py +++ b/package/MDAnalysis/analysis/hbonds/__init__.py @@ -20,14 +20,19 @@ # MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # -from __future__ import absolute_import + +import warning + +warnings.warn( + "This module was moved. " + "Please use MDAnalysis.analysis.hydrogenbonds instead.", + category=DeprecationWarning + ) __all__ = [ - 'HydrogenBondAnalysis', - 'HydrogenBondAutoCorrel', 'find_hydrogen_donors', + 'HydrogenBondAutoCorrel', 'WaterBridgeAnalysis', ] -from .hbond_analysis import HydrogenBondAnalysis -from .hbond_autocorrel import HydrogenBondAutoCorrel, find_hydrogen_donors -from .wbridge_analysis import WaterBridgeAnalysis +from ..hydrogenbonds import HydrogenBondAutoCorrel, find_hydrogen_donors +from ..hydrogenbonds import WaterBridgeAnalysis diff --git a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hbonds/hbond_analysis.py deleted file mode 100644 index e0d34ca721d..00000000000 --- a/package/MDAnalysis/analysis/hbonds/hbond_analysis.py +++ /dev/null @@ -1,1466 +0,0 @@ -# -*- 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 -# - -# Hydrogen Bonding Analysis -r"""Hydrogen Bond analysis (Deprecated) --- :mod:`MDAnalysis.analysis.hbonds.hbond_analysis` -============================================================================================ - -:Author: David Caplan, Lukas Grossar, Oliver Beckstein -:Year: 2010-2017 -:Copyright: GNU Public License v3 - -..Warning: - This module is deprecated and will be removed in version 2.0. - Please use :mod:`MDAnalysis.analysis.hydrogenbonds.hbond_analysis` instead. - -Given a :class:`~MDAnalysis.core.universe.Universe` (simulation -trajectory with 1 or more frames) measure all hydrogen bonds for each -frame between selections 1 and 2. - -The :class:`HydrogenBondAnalysis` class is modeled after the `VMD -HBONDS plugin`_. - -.. _`VMD HBONDS plugin`: http://www.ks.uiuc.edu/Research/vmd/plugins/hbonds/ - -Options: - - *update_selections* (``True``): update selections at each frame? - - *selection_1_type* ("both"): selection 1 is the: "donor", "acceptor", "both" - - donor-acceptor *distance* (Å): 3.0 - - Angle *cutoff* (degrees): 120.0 - - *forcefield* to switch between default values for different force fields - - *donors* and *acceptors* atom types (to add additional atom names) - -.. _Analysis Output: - -Output ------- - -The results are - - the **identities** of donor and acceptor heavy-atoms, - - the **distance** between the heavy atom acceptor atom and the hydrogen atom - that is bonded to the heavy atom donor, - - the **angle** donor-hydrogen-acceptor angle (180º is linear). - -Hydrogen bond data are returned per frame, which is stored in -:attr:`HydrogenBondAnalysis.timeseries` (In the following description, ``#`` -indicates comments that are not part of the output.):: - - results = [ - [ # frame 1 - [ # hbond 1 - , - , , , - , - ], - [ # hbond 2 - , - , , , - , - ], - .... - ], - [ # frame 2 - [ ... ], [ ... ], ... - ], - ... - ] - -Using the :meth:`HydrogenBondAnalysis.generate_table` method one can reformat -the results as a flat "normalised" table that is easier to import into a -database or dataframe for further processing. The table itself is a -:class:`numpy.recarray`. - -.. _Detection-of-hydrogen-bonds: - -Detection of hydrogen bonds ---------------------------- - -Hydrogen bonds are recorded based on a geometric criterion: - -1. The distance between acceptor and hydrogen is less than or equal to - `distance` (default is 3 Å). - -2. The angle between donor-hydrogen-acceptor is greater than or equal to - `angle` (default is 120º). - -The cut-off values `angle` and `distance` can be set as keywords to -:class:`HydrogenBondAnalysis`. - -Donor and acceptor heavy atoms are detected from atom names. The current -defaults are appropriate for the CHARMM27 and GLYCAM06 force fields as defined -in Table `Default atom names for hydrogen bonding analysis`_. - -Hydrogen atoms bonded to a donor are searched with one of two algorithms, -selected with the `detect_hydrogens` keyword. - -"distance" - Searches for all hydrogens (name "H*" or name "[123]H" or type "H") in the - same residue as the donor atom within a cut-off distance of 1.2 Å. - -"heuristic" - Looks at the next three atoms in the list of atoms following the donor and - selects any atom whose name matches (name "H*" or name "[123]H"). For - -The "distance" search is more rigorous but slower and is set as the -default. Until release 0.7.6, only the heuristic search was implemented. - -.. versionchanged:: 0.7.6 - Distance search added (see - :meth:`HydrogenBondAnalysis._get_bonded_hydrogens_dist`) and heuristic - search improved (:meth:`HydrogenBondAnalysis._get_bonded_hydrogens_list`) - -.. _Default atom names for hydrogen bonding analysis: - -.. table:: Default heavy atom names for CHARMM27 force field. - - =========== ============== =========== ==================================== - group donor acceptor comments - =========== ============== =========== ==================================== - main chain N O, OC1, OC2 OC1, OC2 from amber99sb-ildn - (Gromacs) - water OH2, OW OH2, OW SPC, TIP3P, TIP4P (CHARMM27,Gromacs) - - ARG NE, NH1, NH2 - ASN ND2 OD1 - ASP OD1, OD2 - CYS SG - CYH SG possible false positives for CYS - GLN NE2 OE1 - GLU OE1, OE2 - HIS ND1, NE2 ND1, NE2 presence of H determines if donor - HSD ND1 NE2 - HSE NE2 ND1 - HSP ND1, NE2 - LYS NZ - MET SD see e.g. [Gregoret1991]_ - SER OG OG - THR OG1 OG1 - TRP NE1 - TYR OH OH - =========== ============== =========== ==================================== - -.. table:: Heavy atom types for GLYCAM06 force field. - - =========== =========== ================== - element donor acceptor - =========== =========== ================== - N N,NT,N3 N,NT - O OH,OW O,O2,OH,OS,OW,OY - S SM - =========== =========== ================== - -Donor and acceptor names for the CHARMM27 force field will also work for e.g. -OPLS/AA (tested in Gromacs). Residue names in the table are for information -only and are not taken into account when determining acceptors and donors. -This can potentially lead to some ambiguity in the assignment of -donors/acceptors for residues such as histidine or cytosine. - -For more information about the naming convention in GLYCAM06 have a look at the -`Carbohydrate Naming Convention in Glycam`_. - -.. _`Carbohydrate Naming Convention in Glycam`: - http://glycam.ccrc.uga.edu/documents/FutureNomenclature.htm - -The lists of donor and acceptor names can be extended by providing lists of -atom names in the `donors` and `acceptors` keywords to -:class:`HydrogenBondAnalysis`. If the lists are entirely inappropriate -(e.g. when analysing simulations done with a force field that uses very -different atom names) then one should either use the value "other" for `forcefield` -to set no default values, or derive a new class and set the default list oneself:: - - class HydrogenBondAnalysis_OtherFF(HydrogenBondAnalysis): - DEFAULT_DONORS = {"OtherFF": tuple(set([...]))} - DEFAULT_ACCEPTORS = {"OtherFF": tuple(set([...]))} - -Then simply use the new class instead of the parent class and call it with -`forcefield` = "OtherFF". Please also consider to contribute the list of heavy -atom names to MDAnalysis. - -.. rubric:: References - -.. [Gregoret1991] L.M. Gregoret, S.D. Rader, R.J. Fletterick, and - F.E. Cohen. Hydrogen bonds involving sulfur atoms in proteins. Proteins, - 9(2):99–107, 1991. `10.1002/prot.340090204`_. - -.. _`10.1002/prot.340090204`: http://dx.doi.org/10.1002/prot.340090204 - - -How to perform HydrogenBondAnalysis ------------------------------------ - -All protein-water hydrogen bonds can be analysed with :: - - import MDAnalysis - import MDAnalysis.analysis.hbonds - - u = MDAnalysis.Universe('topology', 'trajectory') - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'protein', 'resname HOH', distance=3.0, angle=120.0) - h.run() - -.. Note:: - - Due to the way :class:`HydrogenBondAnalysis` is implemented, it is - more efficient to have the second selection (`selection2`) be the - *larger* group, e.g. the water when looking at water-protein - H-bonds or the whole protein when looking at ligand-protein - interactions. - - -The results are stored as the attribute -:attr:`HydrogenBondAnalysis.timeseries`; see :ref:`Analysis Output` for the -format and further options. - -A number of convenience functions are provided to process the -:attr:`~HydrogenBondAnalysis.timeseries` according to varying criteria: - -:meth:`~HydrogenBondAnalysis.count_by_time` - time series of the number of hydrogen bonds per time step -:meth:`~HydrogenBondAnalysis.count_by_type` - data structure with the frequency of each observed hydrogen bond -:meth:`~HydrogenBondAnalysis.timesteps_by_type` - data structure with a list of time steps for each observed hydrogen bond - -For further data analysis it is convenient to process the -:attr:`~HydrogenBondAnalysis.timeseries` data into a normalized table with the -:meth:`~HydrogenBondAnalysis.generate_table` method, which creates a new data -structure :attr:`HydrogenBondAnalysis.table` that contains one row for each -observation of a hydrogen bond:: - - h.generate_table() - -This table can then be easily turned into, e.g., a `pandas.DataFrame`_, and -further analyzed:: - - import pandas as pd - df = pd.DataFrame.from_records(h.table) - -For example, plotting a histogram of the hydrogen bond angles and lengths is as -simple as :: - - df.hist(column=["angle", "distance"]) - -.. _pandas.DataFrame: http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.html - - -.. TODO: notes on selection updating - - -Classes -------- - -.. autoclass:: HydrogenBondAnalysis - :members: - - .. attribute:: timesteps - - List of the times of each timestep. This can be used together with - :attr:`~HydrogenBondAnalysis.timeseries` to find the specific time point - of a hydrogen bond existence, or see :attr:`~HydrogenBondAnalysis.table`. - - .. attribute:: table - - A normalised table of the data in - :attr:`HydrogenBondAnalysis.timeseries`, generated by - :meth:`HydrogenBondAnalysis.generate_table`. It is a - :class:`numpy.recarray` with the following columns: - - 0. "time" - 1. "donor_index" - 2. "acceptor_index" - 3. "donor_resnm" - 4. "donor_resid" - 5. "donor_atom" - 6. "acceptor_resnm" - 7. "acceptor_resid" - 8. "acceptor_atom" - 9. "distance" - 10. "angle" - - It takes up more space than :attr:`~HydrogenBondAnalysis.timeseries` but - it is easier to analyze and to import into databases or dataframes. - - - .. rubric:: Example - - For example, to create a `pandas.DataFrame`_ from ``h.table``:: - - import pandas as pd - df = pd.DataFrame.from_records(h.table) - - - .. versionchanged:: 0.17.0 - The 1-based donor and acceptor indices (*donor_idx* and - *acceptor_idx*) are deprecated in favor of 0-based indices. - - .. automethod:: _get_bonded_hydrogens - - .. automethod:: _get_bonded_hydrogens_dist - - .. automethod:: _get_bonded_hydrogens_list - -""" -from __future__ import division, absolute_import -import six -from six.moves import range, zip - -import warnings -import logging -from collections import defaultdict - -import numpy as np - -from MDAnalysis import MissingDataWarning, NoDataError, SelectionError, SelectionWarning -from .. import base -from MDAnalysis.lib.log import ProgressBar -from MDAnalysis.lib.NeighborSearch import AtomNeighborSearch -from MDAnalysis.lib import distances -from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency - - -logger = logging.getLogger('MDAnalysis.analysis.hbonds') - -warnings.warn( - "This module is deprecated as of MDAnalysis version 1.0." - "It will be removed in MDAnalysis version 2.0" - "Please use MDAnalysis.analysis.hydrogenbonds.hbond_analysis instead.", - category=DeprecationWarning - ) - - -class HydrogenBondAnalysis(base.AnalysisBase): - """Perform a hydrogen bond analysis - - The analysis of the trajectory is performed with the - :meth:`HydrogenBondAnalysis.run` method. The result is stored in - :attr:`HydrogenBondAnalysis.timeseries`. See - :meth:`~HydrogenBondAnalysis.run` for the format. - - The default atom names are taken from the CHARMM 27 force field files, which - will also work for e.g. OPLS/AA in Gromacs, and GLYCAM06. - - *Donors* (associated hydrogens are deduced from topology) - *CHARMM 27* - N of the main chain, water OH2/OW, ARG NE/NH1/NH2, ASN ND2, HIS ND1/NE2, - SER OG, TYR OH, CYS SG, THR OG1, GLN NE2, LYS NZ, TRP NE1 - *GLYCAM06* - N,NT,N3,OH,OW - - *Acceptors* - *CHARMM 27* - O of the main chain, water OH2/OW, ASN OD1, ASP OD1/OD2, CYH SG, GLN OE1, - GLU OE1/OE2, HIS ND1/NE2, MET SD, SER OG, THR OG1, TYR OH - *GLYCAM06* - N,NT,O,O2,OH,OS,OW,OY,P,S,SM - *amber99sb-ildn(Gromacs)* - OC1, OC2 of the main chain - - See Also - -------- - :ref:`Default atom names for hydrogen bonding analysis` - - - .. versionchanged:: 0.7.6 - DEFAULT_DONORS/ACCEPTORS is now embedded in a dict to switch between - default values for different force fields. - - .. versionchanged:: 1.0.0 - Added autocorrelation (MDAnalysis.lib.correlations.py) for calculating hydrogen bond lifetime - ``save_table()`` method has been removed. You can use ``np.save()`` or - ``cPickle.dump()`` on :attr:`HydrogenBondAnalysis.table` instead. - """ - - # use tuple(set()) here so that one can just copy&paste names from the - # table; set() takes care for removing duplicates. At the end the - # DEFAULT_DONORS and DEFAULT_ACCEPTORS should simply be tuples. - - #: default heavy atom names whose hydrogens are treated as *donors* - #: (see :ref:`Default atom names for hydrogen bonding analysis`); - #: use the keyword `donors` to add a list of additional donor names. - DEFAULT_DONORS = { - 'CHARMM27': tuple(set([ - 'N', 'OH2', 'OW', 'NE', 'NH1', 'NH2', 'ND2', 'SG', 'NE2', 'ND1', 'NZ', 'OG', 'OG1', 'NE1', 'OH'])), - 'GLYCAM06': tuple(set(['N', 'NT', 'N3', 'OH', 'OW'])), - 'other': tuple(set([]))} - - #: default atom names that are treated as hydrogen *acceptors* - #: (see :ref:`Default atom names for hydrogen bonding analysis`); - #: use the keyword `acceptors` to add a list of additional acceptor names. - DEFAULT_ACCEPTORS = { - 'CHARMM27': tuple(set([ - 'O', 'OC1', 'OC2', 'OH2', 'OW', 'OD1', 'OD2', 'SG', 'OE1', 'OE1', 'OE2', 'ND1', 'NE2', 'SD', 'OG', 'OG1', 'OH'])), - 'GLYCAM06': tuple(set(['N', 'NT', 'O', 'O2', 'OH', 'OS', 'OW', 'OY', 'SM'])), - 'other': tuple(set([]))} - - #: A :class:`collections.defaultdict` of covalent radii of common donors - #: (used in :meth`_get_bonded_hydrogens_list` to check if a hydrogen is - #: sufficiently close to its donor heavy atom). Values are stored for - #: N, O, P, and S. Any other heavy atoms are assumed to have hydrogens - #: covalently bound at a maximum distance of 1.5 Å. - r_cov = defaultdict(lambda: 1.5, # default value - N=1.31, O=1.31, P=1.58, S=1.55) - - def __init__(self, universe, selection1='protein', selection2='all', selection1_type='both', - update_selection1=True, update_selection2=True, filter_first=True, distance_type='hydrogen', - distance=3.0, angle=120.0, - forcefield='CHARMM27', donors=None, acceptors=None, - debug=False, detect_hydrogens='distance', verbose=False, pbc=False, **kwargs): - """Set up calculation of hydrogen bonds between two selections in a universe. - - The timeseries is accessible as the attribute :attr:`HydrogenBondAnalysis.timeseries`. - - Some initial checks are performed. If there are no atoms selected by - `selection1` or `selection2` or if no donor hydrogens or acceptor atoms - are found then a :exc:`SelectionError` is raised for any selection that - does *not* update (`update_selection1` and `update_selection2` - keywords). For selections that are set to update, only a warning is - logged because it is assumed that the selection might contain atoms at - a later frame (e.g. for distance based selections). - - If no hydrogen bonds are detected or if the initial check fails, look - at the log output (enable with :func:`MDAnalysis.start_logging` and set - `verbose` ``=True``). It is likely that the default names for donors - and acceptors are not suitable (especially for non-standard - ligands). In this case, either change the `forcefield` or use - customized `donors` and/or `acceptors`. - - Parameters - ---------- - universe : Universe - Universe object - selection1 : str (optional) - Selection string for first selection ['protein'] - selection2 : str (optional) - Selection string for second selection ['all'] - selection1_type : {"donor", "acceptor", "both"} (optional) - Selection 1 can be 'donor', 'acceptor' or 'both'. Note that the - value for `selection1_type` automatically determines how - `selection2` handles donors and acceptors: If `selection1` contains - 'both' then `selection2` will also contain 'both'. If `selection1` - is set to 'donor' then `selection2` is 'acceptor' (and vice versa). - ['both']. - update_selection1 : bool (optional) - Update selection 1 at each frame? Setting to ``False`` is recommended - for any static selection to increase performance. [``True``] - update_selection2 : bool (optional) - Update selection 2 at each frame? Setting to ``False`` is recommended - for any static selection to increase performance. [``True``] - filter_first : bool (optional) - Filter selection 2 first to only atoms 3 * `distance` away [``True``] - distance : float (optional) - Distance cutoff for hydrogen bonds; only interactions with a H-A distance - <= `distance` (and the appropriate D-H-A angle, see `angle`) are - recorded. (Note: `distance_type` can change this to the D-A distance.) [3.0] - angle : float (optional) - Angle cutoff for hydrogen bonds; an ideal H-bond has an angle of - 180º. A hydrogen bond is only recorded if the D-H-A angle is - >= `angle`. The default of 120º also finds fairly non-specific - hydrogen interactions and a possibly better value is 150º. [120.0] - forcefield : {"CHARMM27", "GLYCAM06", "other"} (optional) - Name of the forcefield used. Switches between different - :attr:`~HydrogenBondAnalysis.DEFAULT_DONORS` and - :attr:`~HydrogenBondAnalysis.DEFAULT_ACCEPTORS` values. - ["CHARMM27"] - donors : sequence (optional) - Extra H donor atom types (in addition to those in - :attr:`~HydrogenBondAnalysis.DEFAULT_DONORS`), must be a sequence. - acceptors : sequence (optional) - Extra H acceptor atom types (in addition to those in - :attr:`~HydrogenBondAnalysis.DEFAULT_ACCEPTORS`), must be a sequence. - detect_hydrogens : {"distance", "heuristic"} (optional) - Determine the algorithm to find hydrogens connected to donor - atoms. Can be "distance" (default; finds all hydrogens in the - donor's residue within a cutoff of the donor) or "heuristic" - (looks for the next few atoms in the atom list). "distance" should - always give the correct answer but "heuristic" is faster, - especially when the donor list is updated each - for each frame. ["distance"] - distance_type : {"hydrogen", "heavy"} (optional) - Measure hydrogen bond lengths between donor and acceptor heavy - attoms ("heavy") or between donor hydrogen and acceptor heavy - atom ("hydrogen"). If using "heavy" then one should set the *distance* - cutoff to a higher value such as 3.5 Å. ["hydrogen"] - debug : bool (optional) - If set to ``True`` enables per-frame debug logging. This is disabled - by default because it generates a very large amount of output in - the log file. (Note that a logger must have been started to see - the output, e.g. using :func:`MDAnalysis.start_logging`.) [``False``] - verbose : bool (optional) - Toggle progress output. (Can also be given as keyword argument to - :meth:`run`.) - pbc : bool (optional) - Whether to consider periodic boundaries in calculations [``False``] - - Raises - ------ - :exc:`SelectionError` - is raised for each static selection without the required - donors and/or acceptors. - - Notes - ----- - In order to speed up processing, atoms are filtered by a coarse - distance criterion before a detailed hydrogen bonding analysis is - performed (`filter_first` = ``True``). If one of your selections is - e.g. the solvent then `update_selection1` (or `update_selection2`) must - also be ``True`` so that the list of candidate atoms is updated at each - step: this is now the default. - - If your selections will essentially remain the same for all time steps - (i.e. residues are not moving farther than 3 x `distance`), for - instance, if neither water nor large conformational changes are - involved or if the optimization is disabled (`filter_first` = - ``False``) then you can improve performance by setting the - `update_selection1` and/or `update_selection2` keywords to ``False``. - - - .. versionchanged:: 0.7.6 - New `verbose` keyword (and per-frame debug logging disabled by - default). - - New `detect_hydrogens` keyword to switch between two different - algorithms to detect hydrogens bonded to donor. "distance" is a new, - rigorous distance search within the residue of the donor atom, - "heuristic" is the previous list scan (improved with an additional - distance check). - - New `forcefield` keyword to switch between different values of - DEFAULT_DONORS/ACCEPTORS to accomodate different force fields. - Also has an option "other" for no default values. - - .. versionchanged:: 0.8 - The new default for `update_selection1` and `update_selection2` is now - ``True`` (see `Issue 138`_). Set to ``False`` if your selections only - need to be determined once (will increase performance). - - .. versionchanged:: 0.9.0 - New keyword `distance_type` to select between calculation between - heavy atoms or hydrogen-acceptor. It defaults to the previous - behavior (i.e. "hydrogen"). - - .. versionchanged:: 0.11.0 - Initial checks for selections that potentially raise :exc:`SelectionError`. - - .. versionchanged:: 0.17.0 - use 0-based indexing - - .. deprecated:: 0.16 - The previous `verbose` keyword argument was replaced by - `debug`. Note that the `verbose` keyword argument is now - consistently used to toggle progress meters throughout the library. - - .. _`Issue 138`: https://github.com/MDAnalysis/mdanalysis/issues/138 - - """ - super(HydrogenBondAnalysis, self).__init__(universe.trajectory, **kwargs) - - warnings.warn( - "This class is deprecated as of MDAnalysis version 1.0 and will " - "be removed in version 2.0." - "Please use MDAnalysis.analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalysis instead.", - category=DeprecationWarning - ) - - # per-frame debugging output? - self.debug = debug - - self._get_bonded_hydrogens_algorithms = { - "distance": self._get_bonded_hydrogens_dist, # 0.7.6 default - "heuristic": self._get_bonded_hydrogens_list, # pre 0.7.6 - } - if not detect_hydrogens in self._get_bonded_hydrogens_algorithms: - raise ValueError("detect_hydrogens must be one of {0!r}".format( - self._get_bonded_hydrogens_algorithms.keys())) - self.detect_hydrogens = detect_hydrogens - - self.u = universe - self.selection1 = selection1 - self.selection2 = selection2 - self.selection1_type = selection1_type - self.update_selection1 = update_selection1 - self.update_selection2 = update_selection2 - self.filter_first = filter_first - self.distance = distance - self.distance_type = distance_type # note: everything except 'heavy' will give the default behavior - self.angle = angle - self.pbc = pbc and all(self.u.dimensions[:3]) - - # set up the donors/acceptors lists - if donors is None: - donors = [] - if acceptors is None: - acceptors = [] - self.forcefield = forcefield - self.donors = tuple(set(self.DEFAULT_DONORS[forcefield]).union(donors)) - self.acceptors = tuple(set(self.DEFAULT_ACCEPTORS[forcefield]).union(acceptors)) - - if not (self.selection1 and self.selection2): - raise ValueError('HydrogenBondAnalysis: invalid selections') - elif self.selection1_type not in ('both', 'donor', 'acceptor'): - raise ValueError('HydrogenBondAnalysis: Invalid selection type {0!s}'.format(self.selection1_type)) - - self._timeseries = None # final result accessed as self.timeseries - self.timesteps = None # time for each frame - - self.table = None # placeholder for output table - - self._update_selection_1() - self._update_selection_2() - - self._log_parameters() - - if self.selection1_type == 'donor': - self._sanity_check(1, 'donors') - self._sanity_check(2, 'acceptors') - elif self.selection1_type == 'acceptor': - self._sanity_check(1, 'acceptors') - self._sanity_check(2, 'donors') - else: # both - self._sanity_check(1, 'donors') - self._sanity_check(1, 'acceptors') - self._sanity_check(2, 'acceptors') - self._sanity_check(2, 'donors') - logger.info("HBond analysis: initial checks passed.") - - def _sanity_check(self, selection, htype): - """sanity check the selections 1 and 2 - - *selection* is 1 or 2, *htype* is "donors" or "acceptors" - - If selections do not update and the required donor and acceptor - selections are empty then a :exc:`SelectionError` is immediately - raised. - - If selections update dynamically then it is possible that the selection - will yield donors/acceptors at a later step and we only issue a - warning. - - .. versionadded:: 0.11.0 - """ - assert selection in (1, 2) - assert htype in ("donors", "acceptors") - # horrible data organization: _s1_donors, _s2_acceptors, etc, update_selection1, ... - atoms = getattr(self, "_s{0}_{1}".format(selection, htype)) - update = getattr(self, "update_selection{0}".format(selection)) - if not atoms: - errmsg = "No {1} found in selection {0}. " \ - "You might have to specify a custom '{1}' keyword.".format( - selection, htype) - if not update: - logger.error(errmsg) - raise SelectionError(errmsg) - else: - errmsg += " Selection will update so continuing with fingers crossed." - warnings.warn(errmsg, category=SelectionWarning) - logger.warning(errmsg) - - def _log_parameters(self): - """Log important parameters to the logfile.""" - logger.info("HBond analysis: selection1 = %r (update: %r)", self.selection1, self.update_selection1) - logger.info("HBond analysis: selection2 = %r (update: %r)", self.selection2, self.update_selection2) - logger.info("HBond analysis: criterion: donor %s atom and acceptor atom distance <= %.3f A", self.distance_type, - self.distance) - logger.info("HBond analysis: criterion: angle D-H-A >= %.3f degrees", self.angle) - logger.info("HBond analysis: force field %s to guess donor and acceptor names", self.forcefield) - logger.info("HBond analysis: bonded hydrogen detection algorithm: %r", self.detect_hydrogens) - - def _get_bonded_hydrogens(self, atom, **kwargs): - """Find hydrogens bonded to `atom`. - - This method is typically not called by a user but it is documented to - facilitate understanding of the internals of - :class:`HydrogenBondAnalysis`. - - Parameters - ---------- - atom : groups.Atom - heavy atom - **kwargs - passed through to the calculation method that was selected with - the `detect_hydrogens` kwarg of :class:`HydrogenBondAnalysis`. - - Returns - ------- - hydrogen_atoms : AtomGroup or [] - list of hydrogens (can be a :class:`~MDAnalysis.core.groups.AtomGroup`) - or empty list ``[]`` if none were found. - - See Also - -------- - :meth:`_get_bonded_hydrogens_dist` - :meth:`_get_bonded_hydrogens_list` - - - .. versionchanged:: 0.7.6 - Can switch algorithm by using the `detect_hydrogens` keyword to the - constructor. *kwargs* can be used to supply arguments for algorithm. - - """ - return self._get_bonded_hydrogens_algorithms[self.detect_hydrogens](atom, **kwargs) - - def _get_bonded_hydrogens_dist(self, atom): - """Find hydrogens bonded within cutoff to `atom`. - - Hydrogens are detected by either name ("H*", "[123]H*") or type ("H"); - this is not fool-proof as the atom type is not always a character but - the name pattern should catch most typical occurrences. - - The distance from `atom` is calculated for all hydrogens in the residue - and only those within a cutoff are kept. The cutoff depends on the - heavy atom (more precisely, on its element, which is taken as the first - letter of its name ``atom.name[0]``) and is parameterized in - :attr:`HydrogenBondAnalysis.r_cov`. If no match is found then the - default of 1.5 Å is used. - - - Parameters - ---------- - atom : groups.Atom - heavy atom - - Returns - ------- - hydrogen_atoms : AtomGroup or [] - list of hydrogens (can be a :class:`~MDAnalysis.core.groups.AtomGroup`) - or empty list ``[]`` if none were found. - - Notes - ----- - The performance of this implementation could be improved once the - topology always contains bonded information; it currently uses the - selection parser with an "around" selection. - - - .. versionadded:: 0.7.6 - - """ - try: - return atom.residue.atoms.select_atoms( - "(name H* 1H* 2H* 3H* or type H) and around {0:f} name {1!s}" - "".format(self.r_cov[atom.name[0]], atom.name)) - except NoDataError: - return [] - - def _get_bonded_hydrogens_list(self, atom, **kwargs): - """Find "bonded" hydrogens to the donor *atom*. - - At the moment this relies on the **assumption** that the - hydrogens are listed directly after the heavy atom in the - topology. If this is not the case then this function will - fail. - - Hydrogens are detected by name ``H*``, ``[123]H*`` and they have to be - within a maximum distance from the heavy atom. The cutoff distance - depends on the heavy atom and is parameterized in - :attr:`HydrogenBondAnalysis.r_cov`. - - Parameters - ---------- - atom : groups.Atom - heavy atom - **kwargs - ignored - - Returns - ------- - hydrogen_atoms : AtomGroup or [] - list of hydrogens (can be a :class:`~MDAnalysis.core.groups.AtomGroup`) - or empty list ``[]`` if none were found. - - - .. versionchanged:: 0.7.6 - - Added detection of ``[123]H`` and additional check that a - selected hydrogen is bonded to the donor atom (i.e. its - distance to the donor is less than the covalent radius - stored in :attr:`HydrogenBondAnalysis.r_cov` or the default - 1.5 Å). - - Changed name to - :meth:`~HydrogenBondAnalysis._get_bonded_hydrogens_list` - and added *kwargs* so that it can be used instead of - :meth:`~HydrogenBondAnalysis._get_bonded_hydrogens_dist`. - - """ - warnings.warn("_get_bonded_hydrogens_list() (heuristic detection) does " - "not always find " - "all hydrogens; Using detect_hydrogens='distance', when " - "constructing the HydrogenBondAnalysis class is safer. " - "Removal of this feature is targeted for 1.0", - category=DeprecationWarning) - box = self.u.dimensions if self.pbc else None - try: - hydrogens = [ - a for a in self.u.atoms[atom.index + 1:atom.index + 4] - if (a.name.startswith(('H', '1H', '2H', '3H')) and - distances.calc_bonds(atom.position, a.position, box=box) < self.r_cov[atom.name[0]]) - ] - except IndexError: - hydrogens = [] # weird corner case that atom is the last one in universe - return hydrogens - - def _update_selection_1(self): - self._s1 = self.u.select_atoms(self.selection1) - self.logger_debug("Size of selection 1: {0} atoms".format(len(self._s1))) - if not self._s1: - logger.warning("Selection 1 '{0}' did not select any atoms." - .format(str(self.selection1)[:80])) - self._s1_donors = {} - self._s1_donors_h = {} - self._s1_acceptors = {} - if self.selection1_type in ('donor', 'both'): - self._s1_donors = self._s1.select_atoms( - 'name {0}'.format(' '.join(self.donors))) - self._s1_donors_h = {} - for i, d in enumerate(self._s1_donors): - tmp = self._get_bonded_hydrogens(d) - if tmp: - self._s1_donors_h[i] = tmp - self.logger_debug("Selection 1 donors: {0}".format(len(self._s1_donors))) - self.logger_debug("Selection 1 donor hydrogens: {0}".format(len(self._s1_donors_h))) - if self.selection1_type in ('acceptor', 'both'): - self._s1_acceptors = self._s1.select_atoms( - 'name {0}'.format(' '.join(self.acceptors))) - self.logger_debug("Selection 1 acceptors: {0}".format(len(self._s1_acceptors))) - - def _update_selection_2(self): - box = self.u.dimensions if self.pbc else None - self._s2 = self.u.select_atoms(self.selection2) - if self.filter_first and self._s2: - self.logger_debug('Size of selection 2 before filtering:' - ' {} atoms'.format(len(self._s2))) - ns_selection_2 = AtomNeighborSearch(self._s2, box) - self._s2 = ns_selection_2.search(self._s1, 3. * self.distance) - self.logger_debug('Size of selection 2: {0} atoms'.format(len(self._s2))) - if not self._s2: - logger.warning('Selection 2 "{0}" did not select any atoms.' - .format(str(self.selection2)[:80])) - self._s2_donors = {} - self._s2_donors_h = {} - self._s2_acceptors = {} - if not self._s2: - return None - if self.selection1_type in ('donor', 'both'): - self._s2_acceptors = self._s2.select_atoms( - 'name {0}'.format(' '.join(self.acceptors))) - self.logger_debug("Selection 2 acceptors: {0:d}".format(len(self._s2_acceptors))) - if self.selection1_type in ('acceptor', 'both'): - self._s2_donors = self._s2.select_atoms( - 'name {0}'.format(' '.join(self.donors))) - self._s2_donors_h = {} - for i, d in enumerate(self._s2_donors): - tmp = self._get_bonded_hydrogens(d) - if tmp: - self._s2_donors_h[i] = tmp - self.logger_debug("Selection 2 donors: {0:d}".format(len(self._s2_donors))) - self.logger_debug("Selection 2 donor hydrogens: {0:d}".format(len(self._s2_donors_h))) - - def logger_debug(self, *args): - if self.debug: - logger.debug(*args) - - def run(self, start=None, stop=None, step=None, verbose=None, **kwargs): - """Analyze trajectory and produce timeseries. - - Stores the hydrogen bond data per frame as - :attr:`HydrogenBondAnalysis.timeseries` (see there for output - format). - - Parameters - ---------- - start : int (optional) - starting frame-index for analysis, ``None`` is the first one, 0. - `start` and `stop` are 0-based frame indices and are used to slice - the trajectory (if supported) [``None``] - stop : int (optional) - last trajectory frame for analysis, ``None`` is the last one [``None``] - step : int (optional) - read every `step` between `start` (included) and `stop` (excluded), - ``None`` selects 1. [``None``] - verbose : bool (optional) - toggle progress meter output :class:`~MDAnalysis.lib.log.ProgressMeter` - [``True``] - debug : bool (optional) - enable detailed logging of debugging information; this can create - *very big* log files so it is disabled (``False``) by default; setting - `debug` toggles the debug status for :class:`HydrogenBondAnalysis`, - namely the value of :attr:`HydrogenBondAnalysis.debug`. - - Other Parameters - ---------------- - remove_duplicates : bool (optional) - duplicate hydrogen bonds are removed from output if set to the - default value ``True``; normally, this should not be changed. - - See Also - -------- - :meth:`HydrogenBondAnalysis.generate_table` : - processing the data into a different format. - - - .. versionchanged:: 0.7.6 - Results are not returned, only stored in - :attr:`~HydrogenBondAnalysis.timeseries` and duplicate hydrogen bonds - are removed from output (can be suppressed with `remove_duplicates` = - ``False``) - - .. versionchanged:: 0.11.0 - Accept `quiet` keyword. Analysis will now proceed through frames even if - no donors or acceptors were found in a particular frame. - - .. deprecated:: 0.16 - The `quiet` keyword argument is deprecated in favor of the `verbose` - one. Previous use of `verbose` now corresponds to the new keyword - argument `debug`. - - """ - # sets self.start/stop/step and _pm - self._setup_frames(self._trajectory, start, stop, step) - - logger.info("HBond analysis: starting") - logger.debug("HBond analysis: donors %r", self.donors) - logger.debug("HBond analysis: acceptors %r", self.acceptors) - - remove_duplicates = kwargs.pop('remove_duplicates', True) # False: old behaviour - if not remove_duplicates: - logger.warning("Hidden feature remove_duplicates=False activated: you will probably get duplicate H-bonds.") - - debug = kwargs.pop('debug', None) - if debug is not None and debug != self.debug: - self.debug = debug - logger.debug("Toggling debug to %r", self.debug) - if not self.debug: - logger.debug("HBond analysis: For full step-by-step debugging output use debug=True") - - self._timeseries = [] - self.timesteps = [] - - try: - self.u.trajectory.time - def _get_timestep(): - return self.u.trajectory.time - logger.debug("HBond analysis is recording time step") - except NotImplementedError: - # chained reader or xyz(?) cannot do time yet - def _get_timestep(): - return self.u.trajectory.frame - logger.warning("HBond analysis is recording frame number instead of time step") - - logger.info("Starting analysis (frame index start=%d stop=%d, step=%d)", - self.start, self.stop, self.step) - - for ts in ProgressBar(self.u.trajectory[self.start:self.stop:self.step], - desc="HBond analysis", - verbose=kwargs.get('verbose', False)): - # all bonds for this timestep - frame_results = [] - # dict of tuples (atom.index, atom.index) for quick check if - # we already have the bond (to avoid duplicates) - already_found = {} - - frame = ts.frame - timestep = _get_timestep() - self.timesteps.append(timestep) - - self.logger_debug("Analyzing frame %(frame)d, timestep %(timestep)f ps", vars()) - if self.update_selection1: - self._update_selection_1() - if self.update_selection2: - self._update_selection_2() - - box = self.u.dimensions if self.pbc else None - if self.selection1_type in ('donor', 'both') and self._s2_acceptors: - self.logger_debug("Selection 1 Donors <-> Acceptors") - ns_acceptors = AtomNeighborSearch(self._s2_acceptors, box) - for i, donor_h_set in self._s1_donors_h.items(): - d = self._s1_donors[i] - for h in donor_h_set: - res = ns_acceptors.search(h, self.distance) - for a in res: - angle = distances.calc_angles(d.position, - h.position, - a.position, box=box) - angle = np.rad2deg(angle) - donor_atom = h if self.distance_type != 'heavy' else d - dist = distances.calc_bonds(donor_atom.position, a.position, box=box) - if angle >= self.angle and dist <= self.distance: - self.logger_debug( - "S1-D: {0!s} <-> S2-A: {1!s} {2:f} A, {3:f} DEG".format(h.index, a.index, dist, angle)) - frame_results.append( - [h.index, a.index, - (h.resname, h.resid, h.name), - (a.resname, a.resid, a.name), - dist, angle]) - - already_found[(h.index, a.index)] = True - if self.selection1_type in ('acceptor', 'both') and self._s1_acceptors: - self.logger_debug("Selection 1 Acceptors <-> Donors") - ns_acceptors = AtomNeighborSearch(self._s1_acceptors, box) - for i, donor_h_set in self._s2_donors_h.items(): - d = self._s2_donors[i] - for h in donor_h_set: - res = ns_acceptors.search(h, self.distance) - for a in res: - if remove_duplicates and ( - (h.index, a.index) in already_found or - (a.index, h.index) in already_found): - continue - angle = distances.calc_angles(d.position, - h.position, - a.position, box=box) - angle = np.rad2deg(angle) - donor_atom = h if self.distance_type != 'heavy' else d - dist = distances.calc_bonds(donor_atom.position, a.position, box=box) - if angle >= self.angle and dist <= self.distance: - self.logger_debug( - "S1-A: {0!s} <-> S2-D: {1!s} {2:f} A, {3:f} DEG".format(a.index, h.index, dist, angle)) - frame_results.append( - [h.index, a.index, - (h.resname, h.resid, h.name), - (a.resname, a.resid, a.name), - dist, angle]) - - self._timeseries.append(frame_results) - logger.info("HBond analysis: complete; timeseries %s.timeseries", - self.__class__.__name__) - - @property - def timeseries(self): - """Time series of hydrogen bonds. - - The results of the hydrogen bond analysis can be accessed as a `list` of `list` of `list`: - - 1. `timeseries[i]`: data for the i-th trajectory frame (at time - `timesteps[i]`, see :attr:`timesteps`) - 2. `timeseries[i][j]`: j-th hydrogen bond that was detected at the i-th - frame. - 3. ``donor_index, acceptor_index, - donor_name_str, acceptor_name_str, distance, angle = - timeseries[i][j]``: structure of one hydrogen bond data item - - - In the following description, ``#`` indicates comments that are not - part of the output:: - - results = [ - [ # frame 1 - [ # hbond 1 - , , - , , , - ], - [ # hbond 2 - , , - , , , - ], - .... - ], - [ # frame 2 - [ ... ], [ ... ], ... - ], - ... - ] - - The time of each step is not stored with each hydrogen bond frame but in - :attr:`~HydrogenBondAnalysis.timesteps`. - - - Note - ---- - For instance, to find an acceptor atom in :attr:`Universe.atoms` by - *index* one would use ``u.atoms[acceptor_index]``. - - The :attr:`timeseries` is a managed attribute and it is generated - from the underlying data in :attr:`_timeseries` every time the - attribute is accessed. It is therefore costly to call and if - :attr:`timeseries` is needed repeatedly it is recommended that you - assign to a variable:: - - h = HydrogenBondAnalysis(u) - h.run() - timeseries = h.timeseries - - - See Also - -------- - :attr:`table` : structured array of the data - - - .. versionchanged:: 0.16.1 - :attr:`timeseries` has become a managed attribute and is generated from the stored - :attr:`_timeseries` when needed. :attr:`_timeseries` contains the donor atom and - acceptor atom specifiers as tuples `(resname, resid, atomid)` instead of strings. - - .. versionchanged:: 0.17.0 - The 1-based indices "donor_idx" and "acceptor_idx" are being - removed in favor of the 0-based indices "donor_index" and - "acceptor_index". - - """ - return [[self._reformat_hb(hb) for hb in hframe] for hframe in self._timeseries] - - @staticmethod - def _reformat_hb(hb, atomformat="{0[0]!s}{0[1]!s}:{0[2]!s}"): - """Convert 0.16.1 _timeseries hbond item to 0.16.0 hbond item. - - In 0.16.1, donor and acceptor are stored as a tuple(resname, - resid, atomid). In 0.16.0 and earlier they were stored as a string. - - """ - return (hb[:2] - + [atomformat.format(hb[2]), atomformat.format(hb[3])] - + hb[4:]) - - def generate_table(self): - """Generate a normalised table of the results. - - The table is stored as a :class:`numpy.recarray` in the - attribute :attr:`~HydrogenBondAnalysis.table`. - - See Also - -------- - HydrogenBondAnalysis.table - - """ - if self._timeseries is None: - msg = "No timeseries computed, do run() first." - warnings.warn(msg, category=MissingDataWarning) - logger.warning(msg) - return - - num_records = np.sum([len(hframe) for hframe in self._timeseries]) - # build empty output table - dtype = [ - ("time", float), - ("donor_index", int), ("acceptor_index", int), - ("donor_resnm", "|U4"), ("donor_resid", int), ("donor_atom", "|U4"), - ("acceptor_resnm", "|U4"), ("acceptor_resid", int), ("acceptor_atom", "|U4"), - ("distance", float), ("angle", float)] - # according to Lukas' notes below, using a recarray at this stage is ineffective - # and speedups of ~x10 can be achieved by filling a standard array, like this: - out = np.empty((num_records,), dtype=dtype) - cursor = 0 # current row - for t, hframe in zip(self.timesteps, self._timeseries): - for (donor_index, acceptor_index, donor, - acceptor, distance, angle) in hframe: - # donor|acceptor = (resname, resid, atomid) - out[cursor] = (t, donor_index, acceptor_index) + \ - donor + acceptor + (distance, angle) - cursor += 1 - assert cursor == num_records, "Internal Error: Not all HB records stored" - self.table = out.view(np.recarray) - logger.debug("HBond: Stored results as table with %(num_records)d entries.", vars()) - - def _has_timeseries(self): - has_timeseries = self._timeseries is not None - if not has_timeseries: - msg = "No timeseries computed, do run() first." - warnings.warn(msg, category=MissingDataWarning) - logger.warning(msg) - return has_timeseries - - def count_by_time(self): - """Counts the number of hydrogen bonds per timestep. - - Processes :attr:`HydrogenBondAnalysis._timeseries` into the time series - ``N(t)`` where ``N`` is the total number of observed hydrogen bonds at - time ``t``. - - Returns - ------- - counts : numpy.recarray - The resulting array can be thought of as rows ``(time, N)`` where - ``time`` is the time (in ps) of the time step and ``N`` is the - total number of hydrogen bonds. - - """ - if not self._has_timeseries(): - return - - out = np.empty((len(self.timesteps),), dtype=[('time', float), ('count', int)]) - for cursor, time_count in enumerate(zip(self.timesteps, - (len(series) for series in self._timeseries))): - out[cursor] = time_count - return out.view(np.recarray) - - def count_by_type(self): - """Counts the frequency of hydrogen bonds of a specific type. - - Processes :attr:`HydrogenBondAnalysis._timeseries` and returns a - :class:`numpy.recarray` containing atom indices, residue names, residue - numbers (for donors and acceptors) and the fraction of the total time - during which the hydrogen bond was detected. - - Returns - ------- - counts : numpy.recarray - Each row of the array contains data to define a unique hydrogen - bond together with the frequency (fraction of the total time) that - it has been observed. - - - .. versionchanged:: 0.17.0 - The 1-based indices "donor_idx" and "acceptor_idx" are being - deprecated in favor of zero-based indices. - """ - if not self._has_timeseries(): - return - - hbonds = defaultdict(int) - for hframe in self._timeseries: - for (donor_index, acceptor_index, donor, - acceptor, distance, angle) in hframe: - donor_resnm, donor_resid, donor_atom = donor - acceptor_resnm, acceptor_resid, acceptor_atom = acceptor - # generate unambigous key for current hbond \ - # (the donor_heavy_atom placeholder '?' is added later) - # idx_zero is redundant for an unambigous key, but included for - # consistency. - hb_key = ( - donor_index, acceptor_index, - donor_resnm, donor_resid, "?", donor_atom, - acceptor_resnm, acceptor_resid, acceptor_atom) - - hbonds[hb_key] += 1 - - # build empty output table - dtype = [ - ("donor_index", int), ("acceptor_index", int), ('donor_resnm', 'U4'), - ('donor_resid', int), ('donor_heavy_atom', 'U4'), ('donor_atom', 'U4'), - ('acceptor_resnm', 'U4'), ('acceptor_resid', int), ('acceptor_atom', 'U4'), - ('frequency', float) - ] - out = np.empty((len(hbonds),), dtype=dtype) - - # float because of division later - tsteps = float(len(self.timesteps)) - for cursor, (key, count) in enumerate(six.iteritems(hbonds)): - out[cursor] = key + (count / tsteps,) - - # return array as recarray - # The recarray has not been used within the function, because accessing the - # the elements of a recarray (3.65 us) is much slower then accessing those - # of a ndarray (287 ns). - r = out.view(np.recarray) - - # patch in donor heavy atom names (replaces '?' in the key) - h2donor = self._donor_lookup_table_byindex() - r.donor_heavy_atom[:] = [h2donor[idx] for idx in r.donor_index] - - return r - - def timesteps_by_type(self): - """Frames during which each hydrogen bond existed, sorted by hydrogen bond. - - Processes :attr:`HydrogenBondAnalysis._timeseries` and returns a - :class:`numpy.recarray` containing atom indices, residue names, residue - numbers (for donors and acceptors) and each timestep at which the - hydrogen bond was detected. - - In principle, this is the same as :attr:`~HydrogenBondAnalysis.table` - but sorted by hydrogen bond and with additional data for the - *donor_heavy_atom* and angle and distance omitted. - - - Returns - ------- - data : numpy.recarray - - - .. versionchanged:: 0.17.0 - The 1-based indices "donor_idx" and "acceptor_idx" are being - replaced in favor of zero-based indices. - - """ - if not self._has_timeseries(): - return - - hbonds = defaultdict(list) - for (t, hframe) in zip(self.timesteps, self._timeseries): - for (donor_index, acceptor_index, donor, - acceptor, distance, angle) in hframe: - donor_resnm, donor_resid, donor_atom = donor - acceptor_resnm, acceptor_resid, acceptor_atom = acceptor - # generate unambigous key for current hbond - # (the donor_heavy_atom placeholder '?' is added later) - # idx_zero is redundant for key but added for consistency - hb_key = ( - donor_index, acceptor_index, - donor_resnm, donor_resid, "?", donor_atom, - acceptor_resnm, acceptor_resid, acceptor_atom) - hbonds[hb_key].append(t) - - out_nrows = 0 - # count number of timesteps per key to get length of output table - for ts_list in six.itervalues(hbonds): - out_nrows += len(ts_list) - - # build empty output table - dtype = [ - ('donor_index', int), - ('acceptor_index', int), ('donor_resnm', 'U4'), ('donor_resid', int), - ('donor_heavy_atom', 'U4'), ('donor_atom', 'U4'),('acceptor_resnm', 'U4'), - ('acceptor_resid', int), ('acceptor_atom', 'U4'), ('time', float)] - out = np.empty((out_nrows,), dtype=dtype) - - out_row = 0 - for (key, times) in six.iteritems(hbonds): - for tstep in times: - out[out_row] = key + (tstep,) - out_row += 1 - - # return array as recarray - # The recarray has not been used within the function, because accessing the - # the elements of a recarray (3.65 us) is much slower then accessing those - # of a ndarray (287 ns). - r = out.view(np.recarray) - - # patch in donor heavy atom names (replaces '?' in the key) - h2donor = self._donor_lookup_table_byindex() - r.donor_heavy_atom[:] = [h2donor[idx] for idx in r.donor_index] - - return r - - def _donor_lookup_table_byres(self): - """Look-up table to identify the donor heavy atom from resid and hydrogen name. - - Assumptions: - * resids are unique - * hydrogen atom names are unique within a residue - * selections have not changed (because we are simply looking at the last content - of the donors and donor hydrogen lists) - - Donors from `selection1` and `selection2` are merged. - - Output dictionary ``h2donor`` can be used as:: - - heavy_atom_name = h2donor[resid][hydrogen_name] - - """ - s1d = self._s1_donors # list of donor Atom instances - s1h = self._s1_donors_h # dict indexed by donor position in donor list, containg AtomGroups of H - s2d = self._s2_donors - s2h = self._s2_donors_h - - def _make_dict(donors, hydrogens): - # two steps so that entry for one residue can be UPDATED for multiple donors - d = dict((donors[k].resid, {}) for k in range(len(donors)) if k in hydrogens) - for k in range(len(donors)): - if k in hydrogens: - d[donors[k].resid].update(dict((atom.name, donors[k].name) for atom in hydrogens[k])) - return d - - h2donor = _make_dict(s2d, s2h) # 2 is typically the larger group - # merge (in principle h2donor.update(_make_dict(s1d, s1h) should be sufficient - # with our assumptions but the following should be really safe) - for resid, names in _make_dict(s1d, s1h).items(): - if resid in h2donor: - h2donor[resid].update(names) - else: - h2donor[resid] = names - - return h2donor - - def _donor_lookup_table_byindex(self): - """Look-up table to identify the donor heavy atom from hydrogen atom index. - - Assumptions: - * selections have not changed (because we are simply looking at the last content - of the donors and donor hydrogen lists) - - Donors from `selection1` and `selection2` are merged. - - Output dictionary ``h2donor`` can be used as:: - - heavy_atom_name = h2donor[index] - - """ - s1d = self._s1_donors # list of donor Atom instances - s1h = self._s1_donors_h # dict indexed by donor position in donor list, containg AtomGroups of H - s2d = self._s2_donors - s2h = self._s2_donors_h - - def _make_dict(donors, hydrogens): - #return dict(flatten_1([(atom.id, donors[k].name) for atom in hydrogens[k]] for k in range(len(donors)) - # if k in hydrogens)) - x = [] - for k in range(len(donors)): - if k in hydrogens: - x.extend([(atom.index, donors[k].name) for atom in hydrogens[k]]) - return dict(x) - - h2donor = _make_dict(s2d, s2h) # 2 is typically the larger group - h2donor.update(_make_dict(s1d, s1h)) - - return h2donor - - - def autocorrelation(self, tau_max=20, window_step=1, intermittency=0): - """ - Computes and returns the autocorrelation (HydrogenBondLifetimes ) on the computed hydrogen bonds. - - Parameters - ---------- - window_step : int, optional - Jump every `step`-th frame. This is compatible but independant of the taus used. - Note that `step` and `tau_max` work consistently with intermittency. - tau_max : int, optional - Survival probability is calculated for the range 1 <= `tau` <= `tau_max` - intermittency : int, optional - The maximum number of consecutive frames for which a bond can disappear but be counted as present if it - returns at the next frame. An intermittency of `0` is equivalent to a continuous autocorrelation, which does - not allow for the hydrogen bond disappearance. For example, for `intermittency=2`, any given hbond may - disappear for up to two consecutive frames yet be treated as being present at all frames. - The default is continuous (0). - - Returns - ------- - tau_timeseries : list - tau from 1 to `tau_max`. Saved in the field tau_timeseries. - timeseries : list - autcorrelation value for each value of `tau`. Saved in the field timeseries. - timeseries_data: list - raw datapoints from which the average is taken (timeseries). - Time dependancy and distribution can be extracted. - """ - - if self._timeseries is None: - logging.error("Autocorrelation analysis of hydrogen bonds cannot be done before the hydrogen bonds are found") - logging.error("Autocorrelation: Please use the .run() before calling this function") - return - - if self.step != 1: - logging.warning("Autocorrelation function should be carried out on consecutive frames. ") - logging.warning("Autocorrelation: if you would like to allow bonds to break and reform, please use 'intermittency'") - - # Extract the hydrogen bonds IDs only in the format [set(superset(x1,x2), superset(x3,x4)), ..] - hydrogen_bonds = self.timeseries - found_hydrogen_bonds = [{frozenset(bond[0:2]) for bond in frame} for frame in hydrogen_bonds] - - intermittent_hbonds = correct_intermittency(found_hydrogen_bonds, intermittency=intermittency) - tau_timeseries, timeseries, timeseries_data = autocorrelation(intermittent_hbonds, tau_max, - window_step=window_step) - self.acf_tau_timeseries = tau_timeseries - self.acf_timeseries = timeseries - # user can investigate the distribution and sample size - self.acf_timeseries_data = timeseries_data - - return self diff --git a/package/MDAnalysis/analysis/hydrogenbonds/__init__.py b/package/MDAnalysis/analysis/hydrogenbonds/__init__.py index e77893e5d2e..b63f86db917 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/__init__.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/__init__.py @@ -20,10 +20,13 @@ # 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 __all__ = [ - 'HydrogenBondAnalysis' + 'HydrogenBondAnalysis', + 'HydrogenBondAutoCorrel', + 'WaterBridgeAnalysis', ] from .hbond_analysis import HydrogenBondAnalysis +from .hbond_autocorrel import HydrogenBondAutoCorrel, find_hydrogen_donors +from .wbridge_analysis import WaterBridgeAnalysis diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index bab6066b3ac..d662bfd1eca 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -166,8 +166,6 @@ .. autoclass:: HydrogenBondAnalysis :members: """ -from __future__ import absolute_import, division - import warnings import numpy as np diff --git a/package/MDAnalysis/analysis/hbonds/hbond_autocorrel.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_autocorrel.py similarity index 97% rename from package/MDAnalysis/analysis/hbonds/hbond_autocorrel.py rename to package/MDAnalysis/analysis/hydrogenbonds/hbond_autocorrel.py index 1d064e5c34e..409751737a8 100644 --- a/package/MDAnalysis/analysis/hbonds/hbond_autocorrel.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_autocorrel.py @@ -21,8 +21,8 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # """ -Hydrogen bond autocorrelation --- :mod:`MDAnalysis.analysis.hbonds.hbond_autocorrel` -==================================================================================== +Hydrogen bond autocorrelation --- :mod:`MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel` +=========================================================================================== :Author: Richard J. Gowers :Year: 2014 @@ -204,10 +204,6 @@ """ -from __future__ import division, absolute_import -from six.moves import zip -from six import raise_from - import numpy as np import scipy.optimize @@ -220,7 +216,7 @@ from MDAnalysis.due import due, Doi due.cite(Doi("10.1063/1.4922445"), description="Hydrogen bonding autocorrelation time", - path='MDAnalysis.analysis.hbonds.hbond_autocorrel', + path='MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel', ) del Doi @@ -302,16 +298,12 @@ def __init__(self, universe, nsamples=50, # number of different points to sample in a run pbc=True): - #warnings.warn("This class is deprecated, use analysis.hbonds.HydrogenBondAnalysis " - # "which has .autocorrelation function", - # category=DeprecationWarning) - self.u = universe # check that slicing is possible try: self.u.trajectory[0] except Exception: - raise_from(ValueError("Trajectory must support slicing"), None) + raise ValueError("Trajectory must support slicing") from None self.h = hydrogens self.a = acceptors diff --git a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py similarity index 99% rename from package/MDAnalysis/analysis/hbonds/wbridge_analysis.py rename to package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py index 3e7453f0edf..6a43fb8f85c 100644 --- a/package/MDAnalysis/analysis/hbonds/wbridge_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/wbridge_analysis.py @@ -559,12 +559,10 @@ def analysis(current, output, u, **kwargs): of a water bridge existence. """ -from __future__ import print_function, absolute_import, division from collections import defaultdict import logging import warnings -import six import numpy as np from ..base import AnalysisBase @@ -1499,7 +1497,7 @@ def timesteps_by_type(self, analysis_func=None, **kwargs): link_func=self._full_link, time=time, **kwargs) result_list = [] - for key, time_list in six.iteritems(result): + for key, time_list in result.items(): for time in time_list: if output == 'combined': key = list(key) diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py index b7f0e206677..b448ba64ac9 100644 --- a/package/MDAnalysis/analysis/waterdynamics.py +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -449,14 +449,13 @@ :inherited-members: """ -from __future__ import print_function, division, absolute_import +import MDAnalysis.analysis.hydrogenbonds from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency -import MDAnalysis.analysis.hbonds -from six.moves import range, zip_longest import logging import warnings import numpy as np +from itertools import zip_longest logger = logging.getLogger('MDAnalysis.analysis.waterdynamics') @@ -640,11 +639,11 @@ def _getGraphics(self, HBP, t0, tf, maxdt): def run(self, **kwargs): """Analyze trajectory and produce timeseries""" - h_list = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(self.universe, - self.selection1, - self.selection2, - distance=3.5, - angle=120.0) + h_list = MDAnalysis.analysis.hydrogenbonds.HydrogenBondAnalysis(self.universe, + self.selection1, + self.selection2, + d_h_cutoff=3.5, + d_h_a_angle_cutoff=120.0) h_list.run(**kwargs) self.timeseries = self._getGraphics(h_list.timeseries, self.t0, self.tf, self.dtmax) diff --git a/package/MDAnalysis/lib/correlations.py b/package/MDAnalysis/lib/correlations.py index 494e5482a11..84adf24cfeb 100644 --- a/package/MDAnalysis/lib/correlations.py +++ b/package/MDAnalysis/lib/correlations.py @@ -55,10 +55,6 @@ Calculates the continuous or intermittent survival probability of an atom group in a region of interest. - * :class:`MDAnalysis.analysis.hbonds.hbond_analysis` - Calculates the continuous or intermittent hydrogen bond - lifetime. - .. rubric:: References .. [Gu2019] Gu, R.-X.; Baoukina, S.; Tieleman, D. P. (2019) @@ -68,7 +64,6 @@ """ -from __future__ import absolute_import, division import numpy as np from copy import deepcopy diff --git a/package/doc/sphinx/source/documentation_pages/analysis/hbond_analysis.rst b/package/doc/sphinx/source/documentation_pages/analysis/hbond_analysis.rst deleted file mode 100644 index c44a4a8e424..00000000000 --- a/package/doc/sphinx/source/documentation_pages/analysis/hbond_analysis.rst +++ /dev/null @@ -1 +0,0 @@ -.. automodule:: MDAnalysis.analysis.hbonds.hbond_analysis diff --git a/package/doc/sphinx/source/documentation_pages/analysis/wbridge_analysis.rst b/package/doc/sphinx/source/documentation_pages/analysis/wbridge_analysis.rst index b022b17c4a4..ea9ab08ee18 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis/wbridge_analysis.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis/wbridge_analysis.rst @@ -1 +1 @@ -.. automodule:: MDAnalysis.analysis.hbonds.wbridge_analysis +.. automodule:: MDAnalysis.analysis.hydrogenbonds.wbridge_analysis diff --git a/testsuite/MDAnalysisTests/analysis/test_hbonds.py b/testsuite/MDAnalysisTests/analysis/test_hbonds.py deleted file mode 100644 index 93979aa44ea..00000000000 --- a/testsuite/MDAnalysisTests/analysis/test_hbonds.py +++ /dev/null @@ -1,556 +0,0 @@ -# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- -# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 -# -# MDAnalysis --- https://www.mdanalysis.org -# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors -# (see the file AUTHORS for the full list of names) -# -# Released under the GNU Public Licence, v2 or any higher version -# -# Please cite your use of MDAnalysis in published work: -# -# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, -# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. -# MDAnalysis: A Python package for the rapid analysis of molecular dynamics -# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th -# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. -# doi: 10.25080/majora-629e541a-00e -# -# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. -# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. -# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 -# -# Note: to be removed with MDAnalysis.analysis.hbonds.hbond_analysis in 2.0 - -from __future__ import print_function, absolute_import - - -import MDAnalysis -import MDAnalysis.analysis.hbonds -import itertools -import pytest -from MDAnalysis import SelectionError - -from numpy.testing import ( - assert_equal, assert_array_equal, assert_almost_equal, - assert_array_almost_equal, assert_allclose, -) -import numpy as np - - -from six import StringIO - -from MDAnalysisTests.datafiles import PDB_helix, GRO, XTC, waterPSF, waterDCD - -# For type guessing: -from MDAnalysis.topology.core import guess_atom_type -from MDAnalysis.core.topologyattrs import Atomtypes - - -def guess_types(names): - """GRO doesn't supply types, this returns an Attr""" - return Atomtypes(np.array([guess_atom_type(name) for name in names], dtype=object)) - - -class TestHydrogenBondAnalysis(object): - @staticmethod - @pytest.fixture(scope='class') - def universe(): - return MDAnalysis.Universe(PDB_helix) - - @staticmethod - @pytest.fixture(scope='class') - def values(universe): - return { - 'num_bb_hbonds': universe.atoms.n_residues - universe.select_atoms('resname PRO').n_residues - 4, - 'donor_resid': np.array([5, 6, 8, 9, 10, 11, 12, 13]), - 'acceptor_resnm': np.array(['ALA', 'ALA', 'ALA', 'ALA', 'ALA', 'PRO', 'ALA', 'ALA'], dtype='U4'), - } - - kwargs = { - 'selection1': 'protein', - 'selection2': 'protein', - 'detect_hydrogens': "distance", - 'distance': 3.0, - 'angle': 150.0, - } - - @pytest.fixture(scope='class') - def h(self, universe): - kw = self.kwargs.copy() - # kw.update(kwargs) - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(universe, **kw) - if kw['detect_hydrogens'] == 'heuristic': - with pytest.warns(DeprecationWarning): - h.run(verbose=False) - else: - h.run(verbose=False) - return h - - def test_helix_backbone(self, values, h): - - assert len(h.timeseries[0]) == values['num_bb_hbonds'], "wrong number of backbone hydrogen bonds" - assert h.timesteps, [0.0] - - def test_generate_table(self, values, h): - - h.generate_table() - assert len(h.table) == values['num_bb_hbonds'], "wrong number of backbone hydrogen bonds in table" - assert isinstance(h.table, np.core.records.recarray) - assert_array_equal(h.table.donor_resid, values['donor_resid']) - assert_array_equal(h.table.acceptor_resnm, values['acceptor_resnm']) - - def test_atoms_too_far(self): - pdb = ''' -ATOM 1 N LEU 1 32.310 13.778 14.372 1.00 0.00 SYST N 0 -ATOM 2 OW SOL 2 3.024 4.456 4.147 1.00 0.00 SYST H 0''' - - u = MDAnalysis.Universe(StringIO(pdb), format="pdb") - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u, 'resname SOL', 'protein') - h.run(verbose=False) - assert h.timeseries == [[]] - - def test_acceptor_OC1_OC2(self): - gro = '''test -3 - 1ALA OC1 1 0.000 0.000 0.000 - 2ALA N 2 0.300 0.000 0.000 - 2ALA H1 3 0.200 0.000 0.000 -7.29748 7.66094 9.82962''' - - u = MDAnalysis.Universe(StringIO(gro), format="gro") - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u) - h.run(verbose=False) - assert h.timeseries[0][0][2] == 'ALA2:H1' - - def test_true_traj(self): - u = MDAnalysis.Universe(GRO, XTC) - u.add_TopologyAttr(guess_types(u.atoms.names)) - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(u,'protein','resname ASP', distance=3.0, angle=120.0) - h.run() - assert len(h.timeseries) == 10 - - def test_count_by_time(self, values, h): - - c = h.count_by_time() - assert c.tolist(), [(0.0, values['num_bb_hbonds'])] - - def test_count_by_type(self, values, h): - - c = h.count_by_type() - assert_equal(c.frequency, values['num_bb_hbonds'] * [1.0]) - - # TODO: Rename - def test_count_by_type2(self, values, h): - - t = h.timesteps_by_type() - assert_equal(t.time, values['num_bb_hbonds'] * [0.0]) - -class TestHydrogenBondAnalysisPBC(TestHydrogenBondAnalysis): - # This system is identical to above class - # But has a box introduced, and atoms moved into neighbouring images - # The results however should remain identical if PBC is used - # If pbc:True in kwargs is changed, these tests should fail - @staticmethod - @pytest.fixture(scope='class') - def universe(): - u = MDAnalysis.Universe(PDB_helix) - # transfer to memory to changes to coordinates are reset - u.transfer_to_memory() - # place in huge oversized box - # real system is about 30A wide at most - boxsize = 150. - box = np.array([boxsize, boxsize, boxsize, 90., 90., 90.]) - u.dimensions = box - - # then scatter the atoms throughout various images of this box - u.atoms[::4].translate([boxsize * 2, 0, 0]) - u.atoms[1::4].translate([0, boxsize *4, 0]) - u.atoms[2::4].translate([-boxsize * 5, 0, -boxsize * 2]) - - return u - - kwargs = { - 'selection1': 'protein', - 'selection2': 'protein', - 'detect_hydrogens': "distance", - 'distance': 3.0, - 'angle': 150.0, - 'pbc': True, - } - - - -class TestHydrogenBondAnalysisHeuristic(TestHydrogenBondAnalysis): - kwargs = { - 'selection1': 'protein', - 'selection2': 'protein', - 'detect_hydrogens': "heuristic", - 'distance': 3.0, - 'angle': 150.0, - } - - -class TestHydrogenBondAnalysisHeavy(TestHydrogenBondAnalysis): - - kwargs = { - 'selection1': 'protein', - 'selection2': 'protein', - 'detect_hydrogens': "heuristic", - 'distance': 3.5, - 'angle': 150.0, - 'distance_type': 'heavy', - } - - -class TestHydrogenBondAnalysisHeavyFail(TestHydrogenBondAnalysisHeavy): - kwargs = { - 'selection1': 'protein', - 'selection2': 'protein', - 'detect_hydrogens': "heuristic", - 'distance': 3.0, - 'angle': 150.0, - 'distance_type': 'heavy', - } - - @staticmethod - @pytest.fixture() - def values(universe): - return { - 'num_bb_hbonds': 0, # no H-bonds with a D-A distance < 3.0 A (they start at 3.05 A) - 'donor_resid': np.array([]), - 'acceptor_resnm': np.array([], dtype="= framePlus1 for frame, framePlus1 in zip(hbonds.acf_timeseries, hbonds.acf_timeseries[1:])]) - - - -class TestHydrogenBondAnalysisTIP3PHeavyPBC(object): - @staticmethod - @pytest.fixture(scope='class') - def universe(): - return MDAnalysis.Universe(waterPSF, waterDCD) - - kwargs = { - 'selection1': 'all', - 'selection2': 'all', - 'detect_hydrogens': "distance", - 'distance': 3.0, - 'angle': 120.0, - 'distance_type': 'heavy', - 'pbc': True - } - - @pytest.fixture(scope='class') - def h(self, universe): - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(universe, **self.kwargs) - h.run(verbose=False) - h.generate_table() - return h - - @pytest.fixture(scope='class') - def normalized_timeseries(self, h): - # timeseries in normalized form: (t, d_indx1, a_indx1, d_index0, a_index0, donor, acceptor, dist, angle) - # array index: 0 1 2 3 4 5 6 7 8 - timeseries = [[t] + item - for t, hframe in zip(h.timesteps, h.timeseries) - for item in hframe] - return timeseries - - def test_timeseries(self, h, normalized_timeseries): - - assert len(h.timeseries) == 10 - assert len(normalized_timeseries) == 32 - - reference = { - 'distance': {'mean': 2.7627309, 'std': 0.0905052}, - 'angle': {'mean': 158.9038039, 'std': 12.0362826}, - } - - assert_allclose(np.mean([hbond[5] for hbond in normalized_timeseries]), reference['distance']['mean']) - assert_allclose(np.std([hbond[5] for hbond in normalized_timeseries]), reference['distance']['std']) - assert_allclose(np.mean([hbond[6] for hbond in normalized_timeseries]), reference['angle']['mean']) - assert_allclose(np.std([hbond[6] for hbond in normalized_timeseries]), reference['angle']['std']) - - def test_count_by_time(self, h): - - ref_times = np.arange(0.02, 0.21, 0.02) - ref_counts = np.array([3, 2, 4, 4, 4, 4, 3, 2, 3, 3]) - - times, counts = np.array([[count[0], count[1]] for count in h.count_by_time()]).T - - assert_array_almost_equal(times, ref_times) - assert_array_equal(counts, ref_counts) - - def test_count_by_type(self, h): - - ref_values = [1.0, 1.0, 0.5, 0.4, 0.2, 0.1] - - c = h.count_by_type() - assert_equal(np.sort(c.frequency)[::-1], ref_values) - - -class TestHydrogenBondAnalysisTIP3PHeavyPBCStartStep(object): - @staticmethod - @pytest.fixture(scope='class') - def universe(): - return MDAnalysis.Universe(waterPSF, waterDCD) - - kwargs = { - 'selection1': 'all', - 'selection2': 'all', - 'detect_hydrogens': "distance", - 'distance': 3.0, - 'angle': 120.0, - 'distance_type': 'heavy', - 'pbc': True - } - - @pytest.fixture(scope='class') - def h(self, universe): - h = MDAnalysis.analysis.hbonds.HydrogenBondAnalysis(universe, **self.kwargs) - h.run(start=1, step=2, verbose=False) - h.generate_table() - return h - - @pytest.fixture(scope='class') - def normalized_timeseries(self, h): - # timeseries in normalized form: (t, d_indx1, a_indx1, d_index0, a_index0, donor, acceptor, dist, angle) - # array index: 0 1 2 3 4 5 6 7 8 - timeseries = [[t] + item - for t, hframe in zip(h.timesteps, h.timeseries) - for item in hframe] - return timeseries - - def test_timeseries(self, h, normalized_timeseries): - - assert len(h.timeseries) == 5 - assert len(normalized_timeseries) == 15 - - reference = { - 'distance': {'mean': 2.73942464, 'std': 0.05867924}, - 'angle': {'mean': 157.07768079, 'std': 9.72636682}, - } - - assert_allclose(np.mean([hbond[5] for hbond in normalized_timeseries]), reference['distance']['mean']) - assert_allclose(np.std([hbond[5] for hbond in normalized_timeseries]), reference['distance']['std']) - assert_allclose(np.mean([hbond[6] for hbond in normalized_timeseries]), reference['angle']['mean']) - assert_allclose(np.std([hbond[6] for hbond in normalized_timeseries]), reference['angle']['std']) - - def test_count_by_time(self, h): - - ref_times = np.array([0.04, 0.08, 0.12, 0.16, 0.20]) - ref_counts = np.array([2, 4, 4, 2, 3]) - - times, counts = np.array([[count[0], count[1]] for count in h.count_by_time()]).T - - assert_array_almost_equal(times, ref_times) - assert_array_equal(counts, ref_counts) diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbondautocorrel.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbondautocorrel.py index acb0331c244..e5e72c3fe0d 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbondautocorrel.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbondautocorrel.py @@ -20,11 +20,8 @@ # 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 division, absolute_import import pytest -import six -from six.moves import range from MDAnalysisTests.datafiles import ( TRZ, TRZ_psf, @@ -277,7 +274,7 @@ def test_repr(self, u, hydrogens, oxygens, nitrogens): bond_type='continuous', sample_time=0.06, ) - assert isinstance(repr(hbond), six.string_types) + assert isinstance(repr(hbond), str) def test_find_donors(): u = mda.Universe(waterPSF, waterDCD) diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index 56bddc47784..8a55578da23 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -21,9 +21,6 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # - -from __future__ import absolute_import, division - import numpy as np import MDAnalysis from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis diff --git a/testsuite/MDAnalysisTests/analysis/test_wbridge.py b/testsuite/MDAnalysisTests/analysis/test_wbridge.py index ec2904ac720..8639bf656d7 100644 --- a/testsuite/MDAnalysisTests/analysis/test_wbridge.py +++ b/testsuite/MDAnalysisTests/analysis/test_wbridge.py @@ -1,5 +1,4 @@ -from __future__ import print_function, absolute_import -from six import StringIO +from io import StringIO from collections import defaultdict from numpy.testing import ( @@ -7,16 +6,17 @@ import pytest import MDAnalysis -import MDAnalysis.analysis.hbonds -from MDAnalysis.analysis.hbonds.wbridge_analysis import WaterBridgeAnalysis +from MDAnalysis.analysis.hydrogenbonds.wbridge_analysis import WaterBridgeAnalysis + def test_import_from_hbonds(): try: - from MDAnalysis.analysis.hbonds import WaterBridgeAnalysis + from MDAnalysis.analysis.hydrogenbonds import WaterBridgeAnalysis except ImportError: raise AssertionError("Issue #2064 not fixed: " "importing WaterBridgeAnalysis from " - "MDAnalysis.analysis.hbonds failed.'") + "MDAnalysis.analysis.hydrogenbonds failed.'") + class TestWaterBridgeAnalysis(object): @staticmethod diff --git a/testsuite/MDAnalysisTests/utils/test_duecredit.py b/testsuite/MDAnalysisTests/utils/test_duecredit.py index cec267f5cab..7cde28ef8bf 100644 --- a/testsuite/MDAnalysisTests/utils/test_duecredit.py +++ b/testsuite/MDAnalysisTests/utils/test_duecredit.py @@ -58,8 +58,8 @@ def test_duecredit_collector_primary(self, module, path, citekey): ("MDAnalysis.analysis.psa", "MDAnalysis.analysis.psa", "10.1371/journal.pcbi.1004568"), - ("MDAnalysis.analysis.hbonds.hbond_autocorrel", - "MDAnalysis.analysis.hbonds.hbond_autocorrel", + ("MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel", + "MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel", "10.1063/1.4922445"), ("MDAnalysis.analysis.leaflet", "MDAnalysis.analysis.leaflet",