diff --git a/package/CHANGELOG b/package/CHANGELOG index b14f58ebe32..cf733f35cfe 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -68,6 +68,8 @@ Fixes match TPRParser. (PR #2408) * Added parmed to setup.py * Fixed example docs for polymer persistence length (#2582) + * Fixed HydrogenBondAnalysis to return atom indices rather than atom ids (PR #2572). + Fixed the check for bond information in the _get_dh_pairs method (Issue #2396). * Added missing selection module to leaflet.py (Issue #2612) Enhancements diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 443277ee2bd..ef5c746cb88 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -165,7 +165,6 @@ .. autoclass:: HydrogenBondAnalysis :members: - """ from __future__ import absolute_import, division @@ -173,6 +172,7 @@ from .. import base from MDAnalysis.lib.distances import capped_distance, calc_angles +from MDAnalysis.exceptions import NoDataError from ...due import due, Doi @@ -409,9 +409,12 @@ def _get_dh_pairs(self): # If donors_sel is not provided, use topology to find d-h pairs if not self.donors_sel: - if len(self.u.bonds) == 0: - raise Exception('Cannot assign donor-hydrogen pairs via topology as no bonded information is present. ' - 'Please either: load a topology file with bonded information; use the guess_bonds() ' + # We're using u._topology.bonds rather than u.bonds as it is a million times faster to access. + # This is because u.bonds also calculates properties of each bond (e.g bond length). + # See https://github.com/MDAnalysis/mdanalysis/issues/2396#issuecomment-596251787 + if not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0): + raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. ' + 'Please either: load a topology file with bond information; use the guess_bonds() ' 'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff ' 'can be used.') @@ -496,9 +499,9 @@ def _single_frame(self): # 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.ids) - self.hbonds[2].extend(hbond_hydrogens.ids) - self.hbonds[3].extend(hbond_acceptors.ids) + 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) diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index 934ee4e8bd6..f1fd34cf0c4 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -27,9 +27,11 @@ import numpy as np import MDAnalysis from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis +from MDAnalysis.exceptions import NoDataError import pytest -from numpy.testing import assert_allclose, assert_array_almost_equal, assert_array_equal +from numpy.testing import assert_allclose, assert_equal, assert_array_almost_equal, assert_array_equal, \ + assert_almost_equal from MDAnalysisTests.datafiles import waterPSF, waterDCD @@ -87,15 +89,125 @@ def test_count_by_type(self, h): counts = h.count_by_type() assert int(counts[0, 2]) == ref_count - def test_count_by_ids(self, h): + def test_count_by_ids(self, h, universe): ref_counts = [1.0, 1.0, 0.5, 0.4, 0.2, 0.1] unique_hbonds = h.count_by_ids() + most_common_hbond_ids = [12, 14, 9] + assert_equal(unique_hbonds[0,:3], most_common_hbond_ids) + # count_by_ids() returns raw counts # convert to fraction of time that bond was observed counts = unique_hbonds[:, 3] / len(h.timesteps) + assert_allclose(counts, ref_counts) + + +class TestHydrogenBondAnalysisMock(object): + + kwargs = { + 'donors_sel': 'name O', + 'hydrogens_sel': 'name H1 H2', + 'acceptors_sel': 'name O', + 'd_h_cutoff': 1.2, + 'd_a_cutoff': 3.0, + 'd_h_a_angle_cutoff': 120.0 + } + + @staticmethod + @pytest.fixture(scope='class') + def universe(): + # create two water molecules + """ + H4 + \ + O1-H2 .... O2-H3 + / + H1 + """ + n_residues = 2 + u = MDAnalysis.Universe.empty( + n_atoms=n_residues*3, + n_residues=n_residues, + atom_resindex=np.repeat(range(n_residues), 3), + residue_segindex=[0] * n_residues, + trajectory=True, # necessary for adding coordinates + ) + + u.add_TopologyAttr('name', ['O', 'H1', 'H2'] * n_residues) + u.add_TopologyAttr('type', ['O', 'H', 'H'] * n_residues) + u.add_TopologyAttr('resname', ['SOL'] * n_residues) + u.add_TopologyAttr('resid', list(range(1, n_residues + 1))) + u.add_TopologyAttr('id', list(range(1, (n_residues * 3) + 1))) + + # Atomic coordinates with a single hydrogen bond between O1-H2---O2 + pos1 = np.array([[0, 0, 0], # O1 + [-0.249, -0.968, 0], # H1 + [1, 0, 0], # H2 + [2.5, 0, 0], # O2 + [3., 0, 0], # H3 + [2.250, 0.968, 0] # H4 + ]) + + # Atomic coordinates with no hydrogen bonds + pos2 = np.array([[0, 0, 0], # O1 + [-0.249, -0.968, 0], # H1 + [1, 0, 0], # H2 + [4.5, 0, 0], # O2 + [5., 0, 0], # H3 + [4.250, 0.968, 0] # H4 + ]) + + coordinates = np.empty((3, # number of frames + u.atoms.n_atoms, + 3)) + coordinates[0] = pos1 + coordinates[1] = pos2 + coordinates[2] = pos1 + u.load_new(coordinates, order='fac') + + return u + + def test_no_bond_info_exception(self, universe): + + kwargs = { + 'donors_sel': None, + 'hydrogens_sel': None, + 'acceptors_sel': None, + 'd_h_cutoff': 1.2, + 'd_a_cutoff': 3.0, + 'd_h_a_angle_cutoff': 120.0 + } + + with pytest.raises(NoDataError, match="no bond information"): + h = HydrogenBondAnalysis(universe, **kwargs) + h._get_dh_pairs() + + def test_first(self, universe): + + h = HydrogenBondAnalysis(universe, **self.kwargs) + h.run() + + assert len(h.hbonds) == 2 + + frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, dha_angle = h.hbonds[0] + assert_equal(donor_index, 0) + assert_equal(hydrogen_index, 2) + assert_equal(acceptor_index, 3) + assert_almost_equal(da_dst, 2.5) + assert_almost_equal(dha_angle, 180) + + def test_count_by_time(self, universe): + + h = HydrogenBondAnalysis(universe, **self.kwargs) + h.run() + + ref_times = np.array([0, 1, 2]) # u.trajectory.dt is 1 + ref_counts = np.array([1, 0, 1]) + + counts = h.count_by_time() + assert_array_almost_equal(h.timesteps, ref_times) assert_array_equal(counts, ref_counts) @@ -112,7 +224,6 @@ class TestHydrogenBondAnalysisTIP3P_GuessAcceptors_GuessHydrogens_UseTopology_(T 'd_h_a_angle_cutoff': 120.0 } - class TestHydrogenBondAnalysisTIP3P_GuessDonors_NoTopology(object): """Guess the donor atoms involved in hydrogen bonds using the partial charges of the atoms. """