Skip to content

Switch hbonds to results.hbonds for HydrogenBondAnalysis #3283

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
orbeckst merged 1 commit into from
May 7, 2021
Merged

Switch hbonds to results.hbonds for HydrogenBondAnalysis #3283

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 5 additions & 0 deletions package/CHANGELOG
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,8 @@ Enhancements
checking if it can be used in parallel analysis. (Issue #2996, PR #2950)

Changes
* `analysis.hydrogenbonds.hbond_analysis.HydrogenBondAnalyis` now stores
`hbonds` data using the `analysis.base.Results class (Issues #3271, #3261)
* `analysis.rms.RMSD` and `analysis.rms.RMSF` classes now store `rmsd` and
`rmsf` data using the `analysis.base.Results` class (Issues #3274 #3261)
* `analysis.dihedrals` classes now store angle data using the
Expand Down Expand Up @@ -235,6 +237,9 @@ Changes
* Added OpenMM coordinate and topology converters (Issue #2863, PR #2917)

Deprecations
* The `hbonds` attribute of
`hydrogenbonds.hbond_analysis.HydrogenBondAnalysis.hbonds` is now
deprecated in favour of `results.hbonds` (Issues #3271, #3261)
* The `analysis.rms.RMSD.rmsd` and `analysis.rms.RMSF.rmsf` attributes are
now deprecated in favour of `analysis.rms.RMSD.results.rmsd` and
`analysis.rms.RMSF.results.rmsf` respectively. They will be removed in
Expand Down
66 changes: 46 additions & 20 deletions package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@
- *angle* (degrees): angle of the hydrogen bond

Hydrogen bond data are returned in a :class:`numpy.ndarray` on a "one line, one observation" basis
and can be accessed via :attr:`HydrogenBondAnalysis.hbonds`::
and can be accessed via :attr:`HydrogenBondAnalysis.results.hbonds`::

results = [
[
Expand Down Expand Up @@ -209,13 +209,29 @@

.. autoclass:: HydrogenBondAnalysis
:members:

.. attribute:: results.hbonds

A :class:`numpy.ndarray` which contains a list of all observed hydrogen
bond interactions. See `Output`_ for more information.

.. versionadded:: 2.0.0

.. attribute:: hbonds

Alias to the :attr:`results.hbonds` attribute.

.. deprecated:: 2.0.0
Will be removed in MDAnalysis 3.0.0. Please use
:attr:`results.hbonds` instead.
"""
import logging
import warnings
from collections.abc import Iterable

import numpy as np

from .. import base
from ..base import AnalysisBase, Results
from MDAnalysis.lib.distances import capped_distance, calc_angles
from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency
from MDAnalysis.exceptions import NoDataError
Expand All @@ -231,7 +247,7 @@
del Doi


class HydrogenBondAnalysis(base.AnalysisBase):
class HydrogenBondAnalysis(AnalysisBase):
"""
Perform an analysis of hydrogen bonds in a Universe.
"""
Expand Down Expand Up @@ -330,7 +346,8 @@ def __init__(self, universe,
self.d_a_cutoff = d_a_cutoff
self.d_h_a_angle = d_h_a_angle_cutoff
self.update_selections = update_selections
self.hbonds = None
self.results = Results()
Copy link
Member

@orbeckst orbeckst May 7, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, that's necessary when not calling super().__init__()

self.results.hbonds = None

def guess_hydrogens(self,
select='all',
Expand Down Expand Up @@ -563,7 +580,7 @@ def _filter_atoms(self, donors, hydrogens, acceptors):


def _prepare(self):
self.hbonds = [[], [], [], [], [], []]
self.results.hbonds = [[], [], [], [], [], []]

# Set atom selections if they have not been provided
if not self.acceptors_sel:
Expand Down Expand Up @@ -625,16 +642,25 @@ def _single_frame(self):
hbond_angles = d_h_a_angles[hbond_indices]

# Store data on hydrogen bonds found at this frame
self.hbonds[0].extend(np.full_like(hbond_donors, self._ts.frame))
self.hbonds[1].extend(hbond_donors.indices)
self.hbonds[2].extend(hbond_hydrogens.indices)
self.hbonds[3].extend(hbond_acceptors.indices)
self.hbonds[4].extend(hbond_distances)
self.hbonds[5].extend(hbond_angles)
self.results.hbonds[0].extend(np.full_like(hbond_donors,
self._ts.frame))
self.results.hbonds[1].extend(hbond_donors.indices)
self.results.hbonds[2].extend(hbond_hydrogens.indices)
self.results.hbonds[3].extend(hbond_acceptors.indices)
self.results.hbonds[4].extend(hbond_distances)
self.results.hbonds[5].extend(hbond_angles)

def _conclude(self):

self.hbonds = np.asarray(self.hbonds).T
self.results.hbonds = np.asarray(self.results.hbonds).T

@property
def hbonds(self):
wmsg = ("The `hbonds` attribute was deprecated in MDAnalysis 2.0.0 "
"and will be removed in MDAnalysis 3.0.0. Please use "
"`results.hbonds` instead.")
warnings.warn(wmsg, DeprecationWarning)
return self.results.hbonds

def lifetime(self, tau_max=20, window_step=1, intermittency=0):
"""Computes and returns the time-autocorrelation
Expand Down Expand Up @@ -681,7 +707,7 @@ def lifetime(self, tau_max=20, window_step=1, intermittency=0):
autcorrelation value for each value of `tau`
"""

if self.hbonds is None:
if self.results.hbonds is None:
logging.error(
"Autocorrelation analysis of hydrogen bonds cannot be done"
"before the hydrogen bonds are found"
Expand Down Expand Up @@ -709,7 +735,7 @@ def lifetime(self, tau_max=20, window_step=1, intermittency=0):
# [set(superset(x1,x2), superset(x3,x4)), ..]
found_hydrogen_bonds = [set() for _ in self.frames]
for frame_index, frame in enumerate(self.frames):
for hbond in self.hbonds[self.hbonds[:, 0] == frame]:
for hbond in self.results.hbonds[self.results.hbonds[:, 0] == frame]:
found_hydrogen_bonds[frame_index].add(frozenset(hbond[2:4]))

intermittent_hbonds = correct_intermittency(
Expand All @@ -735,7 +761,7 @@ def count_by_time(self):
the number of hydrogen bonds over time.
"""

indices, tmp_counts = np.unique(self.hbonds[:, 0], axis=0,
indices, tmp_counts = np.unique(self.results.hbonds[:, 0], axis=0,
return_counts=True)

indices -= self.start
Expand All @@ -760,8 +786,8 @@ def count_by_type(self):
acceptor atoms in a hydrogen bond.
"""

d = self.u.atoms[self.hbonds[:, 1].astype(int)]
a = self.u.atoms[self.hbonds[:, 3].astype(int)]
d = self.u.atoms[self.results.hbonds[:, 1].astype(int)]
a = self.u.atoms[self.results.hbonds[:, 3].astype(int)]

tmp_hbonds = np.array([d.resnames, d.types, a.resnames, a.types],
dtype=str).T
Expand Down Expand Up @@ -789,9 +815,9 @@ def count_by_ids(self):
in a hydrogen bond.
"""

d = self.u.atoms[self.hbonds[:, 1].astype(int)]
h = self.u.atoms[self.hbonds[:, 2].astype(int)]
a = self.u.atoms[self.hbonds[:, 3].astype(int)]
d = self.u.atoms[self.results.hbonds[:, 1].astype(int)]
h = self.u.atoms[self.results.hbonds[:, 2].astype(int)]
a = self.u.atoms[self.results.hbonds[:, 3].astype(int)]

tmp_hbonds = np.array([d.ids, h.ids, a.ids]).T
hbond_ids, ids_counts = np.unique(tmp_hbonds, axis=0,
Expand Down
52 changes: 34 additions & 18 deletions testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,18 +59,22 @@ def h(self, universe):

def test_hbond_analysis(self, h):

assert len(np.unique(h.hbonds[:, 0])) == 10
assert len(h.hbonds) == 32
assert len(np.unique(h.results.hbonds[:, 0])) == 10
assert len(h.results.hbonds) == 32

reference = {
'distance': {'mean': 2.7627309, 'std': 0.0905052},
'angle': {'mean': 158.9038039, 'std': 12.0362826},
}

assert_allclose(np.mean(h.hbonds[:, 4]), reference['distance']['mean'])
assert_allclose(np.std(h.hbonds[:, 4]), reference['distance']['std'])
assert_allclose(np.mean(h.hbonds[:, 5]), reference['angle']['mean'])
assert_allclose(np.std(h.hbonds[:, 5]), reference['angle']['std'])
assert_allclose(np.mean(h.results.hbonds[:, 4]),
reference['distance']['mean'])
assert_allclose(np.std(h.results.hbonds[:, 4]),
reference['distance']['std'])
assert_allclose(np.mean(h.results.hbonds[:, 5]),
reference['angle']['mean'])
assert_allclose(np.std(h.results.hbonds[:, 5]),
reference['angle']['std'])

def test_count_by_time(self, h):

Expand Down Expand Up @@ -103,6 +107,11 @@ def test_count_by_ids(self, h, universe):

assert_allclose(counts, ref_counts)

def test_hbonds_deprecated_attr(self, h):
wmsg = "The `hbonds` attribute was deprecated in MDAnalysis 2.0.0"
with pytest.warns(DeprecationWarning, match=wmsg):
assert_equal(h.hbonds, h.results.hbonds)


class TestHydrogenBondAnalysisIdeal(object):

Expand Down Expand Up @@ -196,10 +205,10 @@ def test_no_bond_info_exception(self, universe):
h._get_dh_pairs()

def test_first_hbond(self, hydrogen_bonds):
assert len(hydrogen_bonds.hbonds) == 2
assert len(hydrogen_bonds.results.hbonds) == 2

frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, angle =\
hydrogen_bonds.hbonds[0]
hydrogen_bonds.results.hbonds[0]
assert_equal(donor_index, 0)
assert_equal(hydrogen_index, 2)
assert_equal(acceptor_index, 3)
Expand Down Expand Up @@ -333,7 +342,8 @@ def test_between_all(self, universe):
[3, 4, 6], # protein-water
[6, 7, 8] # protein-protein
]
assert_array_equal(hbonds.hbonds[:, 1:4], expected_hbond_indices)
assert_array_equal(hbonds.results.hbonds[:, 1:4],
expected_hbond_indices)

def test_between_PW(self, universe):
# Find only protein-water hydrogen bonds
Expand All @@ -348,7 +358,8 @@ def test_between_PW(self, universe):
expected_hbond_indices = [
[3, 4, 6] # protein-water
]
assert_array_equal(hbonds.hbonds[:, 1:4], expected_hbond_indices)
assert_array_equal(hbonds.results.hbonds[:, 1:4],
expected_hbond_indices)

def test_between_PW_PP(self, universe):
# Find protein-water and protein-protein hydrogen bonds (not
Expand All @@ -368,7 +379,8 @@ def test_between_PW_PP(self, universe):
[3, 4, 6], # protein-water
[6, 7, 8] # protein-protein
]
assert_array_equal(hbonds.hbonds[:, 1:4], expected_hbond_indices)
assert_array_equal(hbonds.results.hbonds[:, 1:4],
expected_hbond_indices)


class TestHydrogenBondAnalysisTIP3P_GuessAcceptors_GuessHydrogens_UseTopology_(TestHydrogenBondAnalysisTIP3P):
Expand Down Expand Up @@ -396,7 +408,7 @@ def test_no_hydrogens(self, universe):

assert h._hydrogens.n_atoms == 0
assert h._donors.n_atoms == 0
assert h.hbonds.size == 0
assert h.results.hbonds.size == 0

class TestHydrogenBondAnalysisTIP3P_GuessDonors_NoTopology(object):
"""Guess the donor atoms involved in hydrogen bonds using the partial charges of the atoms.
Expand Down Expand Up @@ -507,18 +519,22 @@ def h(self, universe):

def test_hbond_analysis(self, h):

assert len(np.unique(h.hbonds[:, 0])) == 5
assert len(h.hbonds) == 15
assert len(np.unique(h.results.hbonds[:, 0])) == 5
assert len(h.results.hbonds) == 15

reference = {
'distance': {'mean': 2.73942464, 'std': 0.05867924},
'angle': {'mean': 157.07768079, 'std': 9.72636682},
}

assert_allclose(np.mean(h.hbonds[:, 4]), reference['distance']['mean'])
assert_allclose(np.std(h.hbonds[:, 4]), reference['distance']['std'])
assert_allclose(np.mean(h.hbonds[:, 5]), reference['angle']['mean'])
assert_allclose(np.std(h.hbonds[:, 5]), reference['angle']['std'])
assert_allclose(np.mean(h.results.hbonds[:, 4]),
reference['distance']['mean'])
assert_allclose(np.std(h.results.hbonds[:, 4]),
reference['distance']['std'])
assert_allclose(np.mean(h.results.hbonds[:, 5]),
reference['angle']['mean'])
assert_allclose(np.std(h.results.hbonds[:, 5]),
reference['angle']['std'])

def test_count_by_time(self, h):

Expand Down