Hbond update

.gitignore

@@ -21,7 +21,8 @@ examples/output.txt
 # ignore Vagrant virtual machines
 .vagrant
 # ignore coverage files
-.coverage
+.coverage*
+!.coveragerc
 .noseids
 htmlcov
 # ignore trajectory offset caches

package/CHANGELOG

@@ -19,6 +19,7 @@ The rules for this file:
  * 2.4.0
 
 Fixes
+  * Update hbond analysis doc string to use exclusive bound language (#3847)
   * LAMMPSDump Reader translates the box to the origin (#3741)
   * hbond analysis: access hbonds results through new results member in count_by_ids() and count_by_type()
   * Added isolayer selection method (Issue #3845)
@@ -31,6 +32,10 @@ Fixes
     (e.g. bonds, angles) (PR #3779).
 
 Enhancements
+  * Added ability for hbond analysis to use types when resnames are not 
+    present (#3847)
+  * Added explanatory warnings when hbond analysis doesn't find any hbonds
+    (#3847)
   * LAMMPSDump Reader optionally unwraps trajectories with image flags upon 
     loading (#3843
   * LAMMPSDump Reader now imports velocities and forces (#3843)

package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py

@@ -174,6 +174,16 @@
 protein-water and protein-protein hydrogen bonds will be found, but
 no water-water hydrogen bonds.
 
+One can also define hydrogen bonds with atom types::
+
+  from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA
+  hbonds = HBA(
+               universe=u,
+               donors_sel='type 2',
+               hydrogens_sel='type 1',
+               acceptors_sel='type 2',
+              )
+
 In order to compute the hydrogen bond lifetime, after finding hydrogen bonds
 one can use the :attr:`lifetime` function::
 
@@ -236,6 +246,7 @@
 from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency
 from MDAnalysis.exceptions import NoDataError
 from MDAnalysis.core.groups import AtomGroup
+from MDAnalysis.analysis.hydrogenbonds.hbond_autocorrel import find_hydrogen_donors
 
 from ...due import due, Doi
 
@@ -260,6 +271,10 @@ def __init__(self, universe,
         """Set up atom selections and geometric criteria for finding hydrogen
         bonds in a Universe.
 
+        Hydrogen bond selections with `donors_sel` , `hydrogens_sel`, and
+        `acceptors_sel` may be achieved with either a *resname*, atom *name* 
+        combination, or when those are absent, with atom *type* selections.
+
         Parameters
         ----------
         universe : Universe
@@ -310,6 +325,8 @@ def __init__(self, universe,
         information is used, as this is the only way that guarantees the
         correct identification of donor-hydrogen pairs.
 
+        .. versionchanged:: 2.4.0
+            Added use of atom types
         .. versionadded:: 2.0.0
             Added `between` keyword
         """
@@ -370,11 +387,14 @@ def guess_hydrogens(self,
         Parameters
         ----------
         select: str (optional)
-            Selection string for atom group from which hydrogens will be identified.
+            :doc:`Selection string </documentation_pages/selections>` for atom group from which hydrogens 
+            will be identified. (e.g., ``(resname X and name H1)`` or ``type 2``)
         max_mass: float (optional)
-            Maximum allowed mass of a hydrogen atom.
+            The mass of a hydrogen atom must be less than this value.
+        min_mass: float (optional)
+            The mass of a hydrogen atom must be greater than this value.
         min_charge: float (optional)
-            Minimum allowed charge of a hydrogen atom.
+            The charge of a hydrogen atom must be greater than this value.
 
         Returns
         -------
@@ -384,6 +404,9 @@ def guess_hydrogens(self,
 
         Notes
         -----
+        Hydrogen selections may be achieved with either a resname, atom 
+        name combination, or when those are absent, atom types.
+
         This function makes use of atomic masses and atomic charges to identify which atoms are hydrogen atoms that are
         capable of participating in hydrogen bonding. If an atom has a mass less than :attr:`max_mass` and an atomic
         charge greater than :attr:`min_charge` then it is considered capable of participating in hydrogen bonds.
@@ -393,6 +416,10 @@ def guess_hydrogens(self,
         Alternatively, this function may be used to quickly generate a :class:`str` of potential hydrogen atoms involved
         in hydrogen bonding. This str may then be modified before being used to set the attribute
         :attr:`hydrogens_sel`.
+
+        .. versionchanged:: 2.4.0
+            Added ability to use atom types
+
         """
 
         if min_mass > max_mass:
@@ -407,13 +434,7 @@ def guess_hydrogens(self,
             ))
         ]
 
-        hydrogens_list = np.unique(
-            [
-                '(resname {} and name {})'.format(r, p) for r, p in zip(hydrogens_ag.resnames, hydrogens_ag.names)
-            ]
-        )
-
-        return " or ".join(hydrogens_list)
+        return self.__group_categories(hydrogens_ag)
 
     def guess_donors(self, select='all', max_charge=-0.5):
         """Guesses which atoms could be considered donors in the analysis. Only use if the universe topology does not
@@ -422,9 +443,10 @@ def guess_donors(self, select='all', max_charge=-0.5):
         Parameters
         ----------
         select: str (optional)
-            Selection string for atom group from which donors will be identified.
+            :doc:`Selection string </documentation_pages/selections>` for atom group from which donors 
+            will be identified. (e.g., ``(resname X and name O1)`` or ``type 2``)
         max_charge: float (optional)
-            Maximum allowed charge of a donor atom.
+            The charge of a donor atom must be less than this value.
 
         Returns
         -------
@@ -434,6 +456,9 @@ def guess_donors(self, select='all', max_charge=-0.5):
 
         Notes
         -----
+        Donor selections may be achieved with either a resname, atom 
+        name combination, or when those are absent, atom types.
+
         This function makes use of and atomic charges to identify which atoms could be considered donor atoms in the
         hydrogen bond analysis. If an atom has an atomic charge less than :attr:`max_charge`, and it is within
         :attr:`d_h_cutoff` of a hydrogen atom, then it is considered capable of participating in hydrogen bonds.
@@ -445,6 +470,9 @@ def guess_donors(self, select='all', max_charge=-0.5):
         in hydrogen bonding. This :class:`str` may then be modified before being used to set the attribute
         :attr:`donors_sel`.
 
+        .. versionchanged:: 2.4.0
+            Added ability to use atom types
+
         """
 
         # We need to know `hydrogens_sel` before we can find donors
@@ -455,31 +483,38 @@ def guess_donors(self, select='all', max_charge=-0.5):
             hydrogens_sel = self.hydrogens_sel
         hydrogens_ag = self.u.select_atoms(hydrogens_sel)
 
-        ag = hydrogens_ag.residues.atoms.select_atoms(
-            "({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format(
-                donors_sel=select,
-                d_h_cutoff=self.d_h_cutoff,
-                hydrogens_sel=hydrogens_sel
+        # 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 (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0):
+            donors_ag = find_hydrogen_donors(hydrogens_ag)
+            donors_ag = donors_ag.intersection(self.u.select_atoms(select))
+        else:
+            donors_ag = hydrogens_ag.residues.atoms.select_atoms(
+                "({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format(
+                    donors_sel=select,
+                    d_h_cutoff=self.d_h_cutoff,
+                    hydrogens_sel=hydrogens_sel
+                )
             )
-        )
-        donors_ag = ag[ag.charges < max_charge]
-        donors_list = np.unique(
-            [
-                '(resname {} and name {})'.format(r, p) for r, p in zip(donors_ag.resnames, donors_ag.names)
-            ]
-        )
 
-        return " or ".join(donors_list)
+        donors_ag = donors_ag[donors_ag.charges < max_charge]
+
+        return self.__group_categories(donors_ag)
 
     def guess_acceptors(self, select='all', max_charge=-0.5):
         """Guesses which atoms could be considered acceptors in the analysis.
 
+        Acceptor selections may be achieved with either a resname, atom 
+        name combination, or when those are absent, atom types.
+
         Parameters
         ----------
         select: str (optional)
-            Selection string for atom group from which acceptors will be identified.
+            :doc:`Selection string </documentation_pages/selections>` for atom group from which acceptors will
+            be identified. (e.g., ``(resname X and name O1)`` or ``type 2``)
         max_charge: float (optional)
-            Maximum allowed charge of an acceptor atom.
+            The charge of an acceptor atom must be less than this value.
 
         Returns
         -------
@@ -489,6 +524,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5):
 
         Notes
         -----
+        Acceptor selections may be achieved with either a resname, atom 
+        name combination, or when those are absent, atom types.
+
         This function makes use of and atomic charges to identify which atoms could be considered acceptor atoms in the
         hydrogen bond analysis. If an atom has an atomic charge less than :attr:`max_charge` then it is considered
         capable of participating in hydrogen bonds.
@@ -499,17 +537,50 @@ def guess_acceptors(self, select='all', max_charge=-0.5):
         in hydrogen bonding. This :class:`str` may then be modified before being used to set the attribute
         :attr:`acceptors_sel`.
 
+        .. versionchanged:: 2.4.0
+            Added ability to use atom types
+
         """
 
         ag = self.u.select_atoms(select)
         acceptors_ag = ag[ag.charges < max_charge]
-        acceptors_list = np.unique(
-            [
-                '(resname {} and name {})'.format(r, p) for r, p in zip(acceptors_ag.resnames, acceptors_ag.names)
-            ]
-        )
 
-        return " or ".join(acceptors_list)
+        return self.__group_categories(acceptors_ag)
+
+    @staticmethod
+    def __group_categories(group):
+        """ Find categories according to universe constraints
+
+        Parameters
+        ----------
+        group : AtomGroup
+            AtomGroups corresponding to either hydrogen bond acceptors, 
+            donors, or hydrogen atoms that meet their respective charge
+            and mass constraints. 
+
+        Returns
+        -------
+        select : str
+            String for each hydrogen bond acceptor/donor/hydrogen atom category.
+
+        .. versionadded:: 2.4.0
+
+        """
+
+        if hasattr(group, "resnames") and hasattr(group, "names"):
+            group_list = np.unique(
+                [
+                    '(resname {} and name {})'.format(r, p) for r, p in zip(group.resnames, group.names)
+                ]
+            )
+        else:
+            group_list = np.unique(
+                [
+                    'type {}'.format(tp) for tp in group.types
+                ]
+            )
+
+        return " or ".join(group_list)
 
     def _get_dh_pairs(self):
         """Finds donor-hydrogen pairs.
@@ -622,6 +693,13 @@ def _single_frame(self):
             return_distances=True,
         )
 
+        if np.size(d_a_indices) == 0:
+            warnings.warn(
+                f"No hydrogen bonds were found given d-a cutoff of "
+                "{self.d_a_cutoff} between Donor, {self.donors_sel}, and "
+                "Acceptor, {self.acceptors_sel}."
+            )
+
         # Remove D-A pairs more than d_a_cutoff away from one another
         tmp_donors = self._donors[d_a_indices.T[0]]
         tmp_hydrogens = self._hydrogens[d_a_indices.T[0]]
@@ -644,6 +722,13 @@ def _single_frame(self):
         )
         hbond_indices = np.where(d_h_a_angles > self.d_h_a_angle)[0]
 
+        if np.size(hbond_indices) == 0:
+            warnings.warn(
+                "No hydrogen bonds were found given angle of "
+                "{self.d_h_a_angle} between Donor, {self.donors_sel}, and "
+                "Acceptor, {self.acceptors_sel}."
+            )
+
         # Retrieve atoms, distances and angles of hydrogen bonds
         hbond_donors = tmp_donors[hbond_indices]
         hbond_hydrogens = tmp_hydrogens[hbond_indices]
@@ -800,8 +885,14 @@ def count_by_type(self):
         d = self.u.atoms[self.results.hbonds[:, 1].astype(np.intp)]
         a = self.u.atoms[self.results.hbonds[:, 3].astype(np.intp)]
 
-        tmp_hbonds = np.array([d.resnames, d.types, a.resnames, a.types],
-                              dtype=str).T
+        if hasattr(d, "resnames"):
+            d_res = d.resnames
+            a_res = a.resnames
+        else:
+            d_res = [None for x in range(len(d.types))]
+            a_res = [None for x in range(len(a.types))]
+
+        tmp_hbonds = np.array([d_res, d.types, a_res, a.types], dtype=str).T
         hbond_type, type_counts = np.unique(
             tmp_hbonds, axis=0, return_counts=True)
         hbond_type_list = []