Fixed the run() method of analysis.hydrogenbonds.hbond_analysis method to use atom indices instead of ids

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

package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py

@@ -165,14 +165,14 @@
 
 .. autoclass:: HydrogenBondAnalysis
    :members:
-
 """
 from __future__ import absolute_import, division
 
 import numpy as np
 
 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)
 

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.
     """