From c47513d5b17a371038a7ad1832cde8b6332c19c8 Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 21 Sep 2022 17:46:03 -0400 Subject: [PATCH 01/17] Added warning explain hbond nan --- .../analysis/hydrogenbonds/hbond_analysis.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 088d5999c66..a6670888bc0 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -622,6 +622,14 @@ def _single_frame(self): return_distances=True, ) + if not d_a_indices: + warnings.warn( + "No hydrogen bonds were found given d-a cutoff of {} between "\ + "Donor, {}, and Acceptor, {}.".format(self.d_a_cutoff, + self.donors_sel, + 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 +652,14 @@ def _single_frame(self): ) hbond_indices = np.where(d_h_a_angles > self.d_h_a_angle)[0] + if not hbond_indices: + warnings.warn( + "No hydrogen bonds were found given angle of {} between "\ + "Donor, {}, and Acceptor, {}.".format(self.d_h_a_cutoff, + self.donors_sel, + self.acceptors_sel) + ) + # Retrieve atoms, distances and angles of hydrogen bonds hbond_donors = tmp_donors[hbond_indices] hbond_hydrogens = tmp_hydrogens[hbond_indices] From 014000512fd09c3069052b87f800ea5b0c47858a Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 21 Sep 2022 18:05:40 -0400 Subject: [PATCH 02/17] Allow atom types for hbond, inclusive to LAMMPS --- .../analysis/hydrogenbonds/hbond_analysis.py | 66 ++++++++++++++----- 1 file changed, 49 insertions(+), 17 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index a6670888bc0..f700a47632e 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -407,11 +407,19 @@ 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) - ] - ) + if hasattr(hydrogens_ag,"resnames") and hasattr(hydrogens_ag,"names"): + hydrogens_list = np.unique( + [ + '(resname {} and name {})'.format(r, p) for r, p in zip(hydrogens_ag.resnames, hydrogens_ag.names) + ] + ) + else: + hydrogens_list = np.unique( + [ + 'type {}'.format(tp) for tp in hydrogens_ag.types + ] + ) + return " or ".join(hydrogens_list) @@ -463,11 +471,19 @@ def guess_donors(self, select='all', max_charge=-0.5): ) ) 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) - ] - ) + + if hasattr(donors_ag,"resnames") and hasattr(donors_ag,"names"): + donors_list = np.unique( + [ + '(resname {} and name {})'.format(r, p) for r, p in zip(donors_ag.resnames, donors_ag.names) + ] + ) + else: + donors_list = np.unique( + [ + 'type {}'.format(tp) for tp in donors_ag.types if tp not in hydrogens_ag.types + ] + ) return " or ".join(donors_list) @@ -503,11 +519,18 @@ def guess_acceptors(self, select='all', max_charge=-0.5): 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) - ] - ) + if hasattr(acceptors_ag,"resnames") and hasattr(acceptors_ag,"names"): + acceptors_list = np.unique( + [ + '(resname {} and name {})'.format(r, p) for r, p in zip(acceptors_ag.resnames, acceptors_ag.names) + ] + ) + else: + acceptors_list = np.unique( + [ + 'type {}'.format(tp) for tp in acceptors_ag.types + ] + ) return " or ".join(acceptors_list) @@ -816,8 +839,17 @@ def count_by_type(self): d = self.u.atoms[self.hbonds[:, 1].astype(np.intp)] a = self.u.atoms[self.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 + else: + d_res = [None for x in range(len(d.types))] + + if hasattr(a,"resnames"): + a_res = a.resnames + else: + 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 = [] From 340cd3fedf7c7bd202de53f7446d30e420d159ba Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 21 Sep 2022 18:06:31 -0400 Subject: [PATCH 03/17] Use bond info if available for guess_donor --- .../analysis/hydrogenbonds/hbond_analysis.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index f700a47632e..6e64bc9b441 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -463,14 +463,17 @@ 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 + if hasattr(hydrogens_ag[0],"bonded_atoms") and hydrogens_ag[0].bonded_atoms: + donors_ag = self.u.atoms[[x.bonded_atoms[0].ix for x in hydrogens_ag]] + else: + 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_ag = ag[ag.charges < max_charge] if hasattr(donors_ag,"resnames") and hasattr(donors_ag,"names"): donors_list = np.unique( From d1a7182b22e36aa254671558ebcf236ffcd5a824 Mon Sep 17 00:00:00 2001 From: jac16 Date: Thu, 22 Sep 2022 09:14:30 -0400 Subject: [PATCH 04/17] Added hbond PR tests and CHANGELOG with bug fixes --- package/CHANGELOG | 6 +- .../analysis/hydrogenbonds/hbond_analysis.py | 10 +- .../source/documentation_pages/references.rst | 4 + .../analysis/test_hydrogenbonds_analysis.py | 131 ++++++++++++++++++ 4 files changed, 146 insertions(+), 5 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 0cce35f9bd9..bff03629f96 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -13,7 +13,8 @@ The rules for this file: * release numbers follow "Semantic Versioning" http://semver.org ------------------------------------------------------------------------------ -??/??/?? IAlibay, Luthaf, hmacdope, rafaelpap, jbarnoud, BFedder, aya9aladdin +??/??/?? IAlibay, Luthaf, hmacdope, rafaelpap, jbarnoud, BFedder, aya9aladdin, + jaclark5 * 2.4.0 @@ -25,6 +26,9 @@ Fixes * fixes guessed_attributes and read_attributes methods of the Topology class, as those methods were giving unexpected results for some topology attributes (e.g. bonds, angles) (PR #3779). + * Make cutoffs in hbond analysis inclusive to match language (#3847) + * Added ability for hbond anaylsis to use types when resnames aren't there + * Added explanatory warnings when hbond analysis doesn't find any hbonds diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 6e64bc9b441..c453e793a25 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -236,6 +236,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 @@ -464,7 +465,8 @@ def guess_donors(self, select='all', max_charge=-0.5): hydrogens_ag = self.u.select_atoms(hydrogens_sel) if hasattr(hydrogens_ag[0],"bonded_atoms") and hydrogens_ag[0].bonded_atoms: - donors_ag = self.u.atoms[[x.bonded_atoms[0].ix for x in hydrogens_ag]] + donors_ag = find_hydrogen_donors(hydrogens_ag) +# donors_ag = self.u.atoms[[x.bonded_atoms[0].ix for x in hydrogens_ag]] else: ag = hydrogens_ag.residues.atoms.select_atoms( "({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format( @@ -648,7 +650,7 @@ def _single_frame(self): return_distances=True, ) - if not d_a_indices: + if np.size(d_a_indices) == 0: warnings.warn( "No hydrogen bonds were found given d-a cutoff of {} between "\ "Donor, {}, and Acceptor, {}.".format(self.d_a_cutoff, @@ -678,10 +680,10 @@ def _single_frame(self): ) hbond_indices = np.where(d_h_a_angles > self.d_h_a_angle)[0] - if not hbond_indices: + if np.size(hbond_indices) == 0: warnings.warn( "No hydrogen bonds were found given angle of {} between "\ - "Donor, {}, and Acceptor, {}.".format(self.d_h_a_cutoff, + "Donor, {}, and Acceptor, {}.".format(self.d_h_a_angle, self.donors_sel, self.acceptors_sel) ) diff --git a/package/doc/sphinx/source/documentation_pages/references.rst b/package/doc/sphinx/source/documentation_pages/references.rst index 166ddeb27ca..9f6dd797b2b 100644 --- a/package/doc/sphinx/source/documentation_pages/references.rst +++ b/package/doc/sphinx/source/documentation_pages/references.rst @@ -180,6 +180,10 @@ If you use the Mean Squared Displacement analysis code in .. _`10.21105/joss.00877`: https://doi.org/10.21105/joss.00877 +For additional information on non-gaussian parameter, consider [Fuchs1998]_. + +.. [Fuchs1998] Fuchs M., Gotze W., and Mayr M. R. Asymptotic laws for tagged-particle motion in glassy systems. *Phys. Rev. E*, 58(3), 3384-3399. doi:`10.1103/PhysRevE.58.3384 `_ + If you calculate shape parameters using :meth:`~MDAnalysis.core.group.AtomGroup.shape_parameter`, :meth:`~MDAnalysis.core.group.ResidueGroup.shape_parameter`, diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index 2fd191073c3..be434af6940 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -250,6 +250,137 @@ def test_logging_step_not_1(self, universe, caplog): "Autocorrelation: Hydrogen bonds were computed with step > 1." assert any(warning in rec.getMessage() for rec in caplog.records) +class TestHydrogenBondAnalysisNoRes(object): + + kwargs = { + 'donors_sel': 'type O', + 'hydrogens_sel': 'type H H', + 'acceptors_sel': 'type 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, + trajectory=True, # necessary for adding coordinates + ) + + u.add_TopologyAttr('type', ['O', 'H', 'H'] * n_residues) + 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 + + @staticmethod + @pytest.fixture(scope='class') + def hydrogen_bonds(universe): + h = HydrogenBondAnalysis( + universe, + **TestHydrogenBondAnalysisNoRes.kwargs + ) + h.run() + return h + + + 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_hbond(self, hydrogen_bonds): + assert len(hydrogen_bonds.results.hbonds) == 2 + frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, angle =\ + hydrogen_bonds.results.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(angle, 180) + + def test_count_by_time(self, hydrogen_bonds): + ref_times = np.array([0, 1, 2]) # u.trajectory.dt is 1 + ref_counts = np.array([1, 0, 1]) + + counts = hydrogen_bonds.count_by_time() + assert_array_almost_equal(hydrogen_bonds.times, ref_times) + assert_array_equal(counts, ref_counts) + + def test_hydrogen_bond_lifetime(self, hydrogen_bonds): + tau_timeseries, timeseries = hydrogen_bonds.lifetime(tau_max=2) + assert_array_equal(timeseries, [1, 0, 0]) + + def test_hydrogen_bond_lifetime_intermittency(self, hydrogen_bonds): + tau_timeseries, timeseries = hydrogen_bonds.lifetime( + tau_max=2, intermittency=1 + ) + assert_array_equal(timeseries, 1) + + def test_no_attr_hbonds(self, universe): + hbonds = HydrogenBondAnalysis(universe, **self.kwargs) + # hydrogen bonds are not computed + with pytest.raises(NoDataError, match=".hbonds attribute is None"): + hbonds.lifetime(tau_max=2, intermittency=1) + + def test_logging_step_not_1(self, universe, caplog): + hbonds = HydrogenBondAnalysis(universe, **self.kwargs) + # using step 2 + hbonds.run(step=2) + + caplog.set_level(logging.WARNING) + hbonds.lifetime(tau_max=2, intermittency=1) + + warning = \ + "Autocorrelation: Hydrogen bonds were computed with step > 1." + assert any(warning in rec.getMessage() for rec in caplog.records) + class TestHydrogenBondAnalysisBetween(object): From f7f8bf5eef9c71b9994da62296191b3aacc4df37 Mon Sep 17 00:00:00 2001 From: jac16 Date: Tue, 27 Sep 2022 10:57:52 -0400 Subject: [PATCH 05/17] Hydrogen Bonding: Update docs --- .../analysis/hydrogenbonds/hbond_analysis.py | 26 ++++++++++++++++++- .../source/documentation_pages/references.rst | 4 --- 2 files changed, 25 insertions(+), 5 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index c453e793a25..94ec0ab80d9 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -174,6 +174,17 @@ protein-water and protein-protein hydrogen bonds will be found, but no water-water hydrogen bonds. +One can now also define hydrogen bonds with atom types:: + + import MDAnalysis + 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:: @@ -261,6 +272,9 @@ def __init__(self, universe, """Set up atom selections and geometric criteria for finding hydrogen bonds in a Universe. + Hydrogen bond selections may be achieved with either a resname, atom + name combination, or when those are absent, atom types. + Parameters ---------- universe : Universe @@ -311,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 """ @@ -368,6 +384,9 @@ def guess_hydrogens(self, ): """Guesses which hydrogen atoms should be used in the analysis. + Hydrogen selections may be achieved with either a resname, atom + name combination, or when those are absent, atom types. + Parameters ---------- select: str (optional) @@ -428,6 +447,9 @@ 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 contain bonding information, otherwise donor-hydrogen pairs may be incorrectly assigned. + Donor selections may be achieved with either a resname, atom + name combination, or when those are absent, atom types. + Parameters ---------- select: str (optional) @@ -466,7 +488,6 @@ def guess_donors(self, select='all', max_charge=-0.5): if hasattr(hydrogens_ag[0],"bonded_atoms") and hydrogens_ag[0].bonded_atoms: donors_ag = find_hydrogen_donors(hydrogens_ag) -# donors_ag = self.u.atoms[[x.bonded_atoms[0].ix for x in hydrogens_ag]] else: ag = hydrogens_ag.residues.atoms.select_atoms( "({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format( @@ -495,6 +516,9 @@ def guess_donors(self, select='all', max_charge=-0.5): 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) diff --git a/package/doc/sphinx/source/documentation_pages/references.rst b/package/doc/sphinx/source/documentation_pages/references.rst index 9f6dd797b2b..166ddeb27ca 100644 --- a/package/doc/sphinx/source/documentation_pages/references.rst +++ b/package/doc/sphinx/source/documentation_pages/references.rst @@ -180,10 +180,6 @@ If you use the Mean Squared Displacement analysis code in .. _`10.21105/joss.00877`: https://doi.org/10.21105/joss.00877 -For additional information on non-gaussian parameter, consider [Fuchs1998]_. - -.. [Fuchs1998] Fuchs M., Gotze W., and Mayr M. R. Asymptotic laws for tagged-particle motion in glassy systems. *Phys. Rev. E*, 58(3), 3384-3399. doi:`10.1103/PhysRevE.58.3384 `_ - If you calculate shape parameters using :meth:`~MDAnalysis.core.group.AtomGroup.shape_parameter`, :meth:`~MDAnalysis.core.group.ResidueGroup.shape_parameter`, From 299b89889661a179b24893791c9dbb2dcfff7852 Mon Sep 17 00:00:00 2001 From: jac16 Date: Tue, 27 Sep 2022 11:32:16 -0400 Subject: [PATCH 06/17] Amended hbond test for full coverage, added _get_dh_pairs ability to find hydrogen atoms. --- .../analysis/hydrogenbonds/hbond_analysis.py | 7 +++- .../analysis/test_hydrogenbonds_analysis.py | 34 +++++++++++++++++++ 2 files changed, 40 insertions(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 94ec0ab80d9..f99b98b80bf 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -592,7 +592,12 @@ def _get_dh_pairs(self): # Otherwise, use d_h_cutoff as a cutoff distance else: - hydrogens = self.u.select_atoms(self.hydrogens_sel) + if self.hydrogens_sel is None: + hydrogens_sel = self.guess_hydrogens() + else: + hydrogens_sel = self.hydrogens_sel + + hydrogens = self.u.select_atoms(hydrogens_sel) donors = self.u.select_atoms(self.donors_sel) donors_indices, hydrogen_indices = capped_distance( donors.positions, diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index be434af6940..eb1a87c1af1 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -335,6 +335,40 @@ def test_no_bond_info_exception(self, universe): h = HydrogenBondAnalysis(universe, **kwargs) h._get_dh_pairs() + def test_no_bond_no_mass_donor_sel(self, universe): + + kwargs = { + 'donors_sel': "type O", + '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="not contain mass information"): + h = HydrogenBondAnalysis(universe, **kwargs) + h._get_dh_pairs() + + def test_no_bond_donor_sel(self, universe): + + kwargs = { + 'donors_sel': "type O", + 'hydrogens_sel': None, + 'acceptors_sel': None, + 'd_h_cutoff': 1.2, + 'd_a_cutoff': 3.0, + 'd_h_a_angle_cutoff': 120.0 + } + u = universe.copy() + n_residues = 2 + u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) + u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) + h = HydrogenBondAnalysis(u, **kwargs) + pairs = h._get_dh_pairs() + + assert len(pairs) == 2 + def test_first_hbond(self, hydrogen_bonds): assert len(hydrogen_bonds.results.hbonds) == 2 frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, angle =\ From bd9590012bfae8523755c1d7b9ff6ee12dd3cf4d Mon Sep 17 00:00:00 2001 From: jac16 Date: Wed, 28 Sep 2022 20:55:16 -0400 Subject: [PATCH 07/17] hbond_analysis complete codecov --- .gitignore | 2 +- .../analysis/hydrogenbonds/hbond_analysis.py | 5 ++- .../analysis/test_hydrogenbonds_analysis.py | 31 +++++++++++++++++-- 3 files changed, 34 insertions(+), 4 deletions(-) diff --git a/.gitignore b/.gitignore index 0760f8129e4..ab758cf45e1 100644 --- a/.gitignore +++ b/.gitignore @@ -21,7 +21,7 @@ examples/output.txt # ignore Vagrant virtual machines .vagrant # ignore coverage files -.coverage +.coverage* .noseids htmlcov # ignore trajectory offset caches diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index f99b98b80bf..9201f04fe7a 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -486,7 +486,10 @@ def guess_donors(self, select='all', max_charge=-0.5): hydrogens_sel = self.hydrogens_sel hydrogens_ag = self.u.select_atoms(hydrogens_sel) - if hasattr(hydrogens_ag[0],"bonded_atoms") and hydrogens_ag[0].bonded_atoms: + # 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) else: ag = hydrogens_ag.residues.atoms.select_atoms( diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index eb1a87c1af1..bf44af95de0 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -204,6 +204,23 @@ def test_no_bond_info_exception(self, universe): h = HydrogenBondAnalysis(universe, **kwargs) h._get_dh_pairs() + def test_no_bond_donor_sel(self, universe): + + kwargs = { + 'donors_sel': "type O", + 'hydrogens_sel': None, + 'acceptors_sel': None, + 'd_h_cutoff': 1.2, + 'd_a_cutoff': 3.0, + 'd_h_a_angle_cutoff': 120.0 + } + u = universe.copy() + n_residues = 2 + u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) + u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) + h = HydrogenBondAnalysis(u, **kwargs) + donors = u.select_atoms(h.guess_donors()) + def test_first_hbond(self, hydrogen_bonds): assert len(hydrogen_bonds.results.hbonds) == 2 frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, angle =\ @@ -319,6 +336,16 @@ def hydrogen_bonds(universe): h.run() return h + def test_count_by_type(self, universe): + + h = HydrogenBondAnalysis( + universe, + **TestHydrogenBondAnalysisNoRes.kwargs + ) + h.run() + counts = h.count_by_type() + ref_count = 2 + assert int(counts[0, 2]) == ref_count def test_no_bond_info_exception(self, universe): @@ -365,9 +392,9 @@ def test_no_bond_donor_sel(self, universe): u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) h = HydrogenBondAnalysis(u, **kwargs) - pairs = h._get_dh_pairs() + donors = u.select_atoms(h.guess_donors()) - assert len(pairs) == 2 + assert len(donors) == 2 def test_first_hbond(self, hydrogen_bonds): assert len(hydrogen_bonds.results.hbonds) == 2 From 5d27895503ed29c299e472513e59c48e87f9a2ff Mon Sep 17 00:00:00 2001 From: jac16 Date: Thu, 29 Sep 2022 08:37:43 -0400 Subject: [PATCH 08/17] Adjusted hbond doc strings for max_*/min_* bounds --- .../analysis/hydrogenbonds/hbond_analysis.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 9201f04fe7a..d586f22b870 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -392,9 +392,11 @@ def guess_hydrogens(self, select: str (optional) Selection string for atom group from which hydrogens will be identified. 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 ------- @@ -455,7 +457,7 @@ def guess_donors(self, select='all', max_charge=-0.5): select: str (optional) Selection string for atom group from which donors will be identified. max_charge: float (optional) - Maximum allowed charge of a donor atom. + The charge of a donor atom must be less than this value. Returns ------- @@ -527,7 +529,7 @@ def guess_acceptors(self, select='all', max_charge=-0.5): select: str (optional) Selection string for atom group from which acceptors will be identified. max_charge: float (optional) - Maximum allowed charge of an acceptor atom. + The charge of an acceptor atom must be less than this value. Returns ------- From 7af74d84a1a48413d3228b897b9ddc4363808073 Mon Sep 17 00:00:00 2001 From: jac16 Date: Tue, 18 Oct 2022 10:47:58 -0400 Subject: [PATCH 09/17] Incorporate reviewer comments --- .gitignore | 1 + package/CHANGELOG | 2 +- .../analysis/hydrogenbonds/hbond_analysis.py | 82 ++++++++----------- .../analysis/test_hydrogenbonds_analysis.py | 9 +- 4 files changed, 42 insertions(+), 52 deletions(-) diff --git a/.gitignore b/.gitignore index ab758cf45e1..e9f3571ca8e 100644 --- a/.gitignore +++ b/.gitignore @@ -22,6 +22,7 @@ examples/output.txt .vagrant # ignore coverage files .coverage* +!.coveragerc .noseids htmlcov # ignore trajectory offset caches diff --git a/package/CHANGELOG b/package/CHANGELOG index 4b783e120ec..9417f481484 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -19,7 +19,6 @@ The rules for this file: * 2.4.0 Fixes - * Update hbond analysis doc string to use exclusive bound language (#3847) * Added isolayer selection method (Issue #3845) * Auxiliary; determination of representative frames: Removed undefined behaviour for cutoff values < -1 (PR # 3749) @@ -32,6 +31,7 @@ Fixes Enhancements + * Update hbond analysis doc string to use exclusive bound language (#3847) * Added ability for hbond anaylsis to use types when resnames aren't there * Added explanatory warnings when hbond analysis doesn't find any hbonds * AuxReaders now have a memory_limit parameter to control when memory usage diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index d586f22b870..e527d99946b 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -429,21 +429,7 @@ def guess_hydrogens(self, )) ] - if hasattr(hydrogens_ag,"resnames") and hasattr(hydrogens_ag,"names"): - hydrogens_list = np.unique( - [ - '(resname {} and name {})'.format(r, p) for r, p in zip(hydrogens_ag.resnames, hydrogens_ag.names) - ] - ) - else: - hydrogens_list = np.unique( - [ - 'type {}'.format(tp) for tp in hydrogens_ag.types - ] - ) - - - 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 @@ -503,20 +489,7 @@ def guess_donors(self, select='all', max_charge=-0.5): ) donors_ag = ag[ag.charges < max_charge] - if hasattr(donors_ag,"resnames") and hasattr(donors_ag,"names"): - donors_list = np.unique( - [ - '(resname {} and name {})'.format(r, p) for r, p in zip(donors_ag.resnames, donors_ag.names) - ] - ) - else: - donors_list = np.unique( - [ - 'type {}'.format(tp) for tp in donors_ag.types if tp not in hydrogens_ag.types - ] - ) - - return " or ".join(donors_list) + 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. @@ -553,20 +526,40 @@ def guess_acceptors(self, select='all', max_charge=-0.5): ag = self.u.select_atoms(select) acceptors_ag = ag[ag.charges < max_charge] - if hasattr(acceptors_ag,"resnames") and hasattr(acceptors_ag,"names"): - acceptors_list = np.unique( + + 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 meant their respective charge + and mass constraints. + + Returns + ------- + select : str + String for each hydrogen bond acceptor/donor/hydrogen atom category. + """ + + if hasattr(group,"resnames") and hasattr(group,"names"): + group_list = np.unique( [ - '(resname {} and name {})'.format(r, p) for r, p in zip(acceptors_ag.resnames, acceptors_ag.names) + '(resname {} and name {})'.format(r, p) for r, p in zip(group.resnames, group.names) ] ) else: - acceptors_list = np.unique( + group_list = np.unique( [ - 'type {}'.format(tp) for tp in acceptors_ag.types + 'type {}'.format(tp) for tp in group.types ] ) - return " or ".join(acceptors_list) + return " or ".join(group_list) def _get_dh_pairs(self): """Finds donor-hydrogen pairs. @@ -686,10 +679,9 @@ def _single_frame(self): if np.size(d_a_indices) == 0: warnings.warn( - "No hydrogen bonds were found given d-a cutoff of {} between "\ - "Donor, {}, and Acceptor, {}.".format(self.d_a_cutoff, - self.donors_sel, - self.acceptors_sel) + 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 @@ -716,10 +708,9 @@ def _single_frame(self): if np.size(hbond_indices) == 0: warnings.warn( - "No hydrogen bonds were found given angle of {} between "\ - "Donor, {}, and Acceptor, {}.".format(self.d_h_a_angle, - self.donors_sel, - self.acceptors_sel) + "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 @@ -880,12 +871,9 @@ def count_by_type(self): if hasattr(d,"resnames"): d_res = d.resnames - else: - d_res = [None for x in range(len(d.types))] - - if hasattr(a,"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 diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index bf44af95de0..e8b0e8006ad 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -263,10 +263,11 @@ def test_logging_step_not_1(self, universe, caplog): caplog.set_level(logging.WARNING) hbonds.lifetime(tau_max=2, intermittency=1) - warning = \ - "Autocorrelation: Hydrogen bonds were computed with step > 1." + warning = ("Autocorrelation: Hydrogen bonds were computed with " + "step > 1.") assert any(warning in rec.getMessage() for rec in caplog.records) + class TestHydrogenBondAnalysisNoRes(object): kwargs = { @@ -438,8 +439,8 @@ def test_logging_step_not_1(self, universe, caplog): caplog.set_level(logging.WARNING) hbonds.lifetime(tau_max=2, intermittency=1) - warning = \ - "Autocorrelation: Hydrogen bonds were computed with step > 1." + warning = ("Autocorrelation: Hydrogen bonds were computed with " + "step > 1.") assert any(warning in rec.getMessage() for rec in caplog.records) From 77af51e7cc87d926fbf0c6bb0541d23af749b1b3 Mon Sep 17 00:00:00 2001 From: jac16 Date: Fri, 4 Nov 2022 10:44:45 -0400 Subject: [PATCH 10/17] Review comment integration --- .../analysis/hydrogenbonds/hbond_analysis.py | 8 +- .../analysis/test_hydrogenbonds_analysis.py | 113 +----------------- 2 files changed, 10 insertions(+), 111 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index c84acfade9d..f8dab7b996d 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -479,15 +479,17 @@ def guess_donors(self, select='all', max_charge=-0.5): # 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(u.select_atoms(select)) else: - ag = hydrogens_ag.residues.atoms.select_atoms( + 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_ag = donors_ag[donors_ag.charges < max_charge] return self.__group_categories(donors_ag) @@ -537,7 +539,7 @@ def __group_categories(group): ---------- group : AtomGroup AtomGroups corresponding to either hydrogen bond acceptors, - donors, or hydrogen atoms that meant their respective charge + donors, or hydrogen atoms that meet their respective charge and mass constraints. Returns diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index e8b0e8006ad..e450adac70d 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -268,7 +268,7 @@ def test_logging_step_not_1(self, universe, caplog): assert any(warning in rec.getMessage() for rec in caplog.records) -class TestHydrogenBondAnalysisNoRes(object): +class TestHydrogenBondAnalysisNoRes(TestHydrogenBondAnalysisIdeal): kwargs = { 'donors_sel': 'type O', @@ -280,7 +280,7 @@ class TestHydrogenBondAnalysisNoRes(object): } @staticmethod - @pytest.fixture(scope='class') + @pytest.fixture(scope='class', autouse=True) def universe(): # create two water molecules """ @@ -293,6 +293,9 @@ def universe(): 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 ) @@ -337,112 +340,6 @@ def hydrogen_bonds(universe): h.run() return h - def test_count_by_type(self, universe): - - h = HydrogenBondAnalysis( - universe, - **TestHydrogenBondAnalysisNoRes.kwargs - ) - h.run() - counts = h.count_by_type() - ref_count = 2 - assert int(counts[0, 2]) == ref_count - - 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_no_bond_no_mass_donor_sel(self, universe): - - kwargs = { - 'donors_sel': "type O", - '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="not contain mass information"): - h = HydrogenBondAnalysis(universe, **kwargs) - h._get_dh_pairs() - - def test_no_bond_donor_sel(self, universe): - - kwargs = { - 'donors_sel': "type O", - 'hydrogens_sel': None, - 'acceptors_sel': None, - 'd_h_cutoff': 1.2, - 'd_a_cutoff': 3.0, - 'd_h_a_angle_cutoff': 120.0 - } - u = universe.copy() - n_residues = 2 - u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) - u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) - h = HydrogenBondAnalysis(u, **kwargs) - donors = u.select_atoms(h.guess_donors()) - - assert len(donors) == 2 - - def test_first_hbond(self, hydrogen_bonds): - assert len(hydrogen_bonds.results.hbonds) == 2 - frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, angle =\ - hydrogen_bonds.results.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(angle, 180) - - def test_count_by_time(self, hydrogen_bonds): - ref_times = np.array([0, 1, 2]) # u.trajectory.dt is 1 - ref_counts = np.array([1, 0, 1]) - - counts = hydrogen_bonds.count_by_time() - assert_array_almost_equal(hydrogen_bonds.times, ref_times) - assert_array_equal(counts, ref_counts) - - def test_hydrogen_bond_lifetime(self, hydrogen_bonds): - tau_timeseries, timeseries = hydrogen_bonds.lifetime(tau_max=2) - assert_array_equal(timeseries, [1, 0, 0]) - - def test_hydrogen_bond_lifetime_intermittency(self, hydrogen_bonds): - tau_timeseries, timeseries = hydrogen_bonds.lifetime( - tau_max=2, intermittency=1 - ) - assert_array_equal(timeseries, 1) - - def test_no_attr_hbonds(self, universe): - hbonds = HydrogenBondAnalysis(universe, **self.kwargs) - # hydrogen bonds are not computed - with pytest.raises(NoDataError, match=".hbonds attribute is None"): - hbonds.lifetime(tau_max=2, intermittency=1) - - def test_logging_step_not_1(self, universe, caplog): - hbonds = HydrogenBondAnalysis(universe, **self.kwargs) - # using step 2 - hbonds.run(step=2) - - caplog.set_level(logging.WARNING) - hbonds.lifetime(tau_max=2, intermittency=1) - - warning = ("Autocorrelation: Hydrogen bonds were computed with " - "step > 1.") - assert any(warning in rec.getMessage() for rec in caplog.records) - class TestHydrogenBondAnalysisBetween(object): From 4ae3b7bcea5fa7d0891c0bc8237c0927360db9e0 Mon Sep 17 00:00:00 2001 From: jac16 Date: Fri, 4 Nov 2022 15:02:27 -0400 Subject: [PATCH 11/17] bug fix --- package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index f8dab7b996d..7b3e584fa13 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -479,7 +479,7 @@ def guess_donors(self, select='all', max_charge=-0.5): # 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(u.select_atoms(select)) + 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( From 333ef6bd0a97f318e0f8ff0ff8935548da68b47c Mon Sep 17 00:00:00 2001 From: jac16 Date: Fri, 4 Nov 2022 17:50:03 -0400 Subject: [PATCH 12/17] Updated tests --- .../MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index e450adac70d..2484171cc9a 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -272,7 +272,7 @@ class TestHydrogenBondAnalysisNoRes(TestHydrogenBondAnalysisIdeal): kwargs = { 'donors_sel': 'type O', - 'hydrogens_sel': 'type H H', +# 'hydrogens_sel': 'type H H', 'acceptors_sel': 'type O', 'd_h_cutoff': 1.2, 'd_a_cutoff': 3.0, @@ -301,6 +301,8 @@ def universe(): u.add_TopologyAttr('type', ['O', 'H', 'H'] * n_residues) u.add_TopologyAttr('id', list(range(1, (n_residues * 3) + 1))) + u.add_TopologyAttr('mass', [15.999, 1.008, 1.008] * n_residues) + u.add_TopologyAttr('charge', [-1.04, 0.52, 0.52] * n_residues) # Atomic coordinates with a single hydrogen bond between O1-H2---O2 pos1 = np.array([[0, 0, 0], # O1 From 049c0eef9fdc4d2810deddbbe5b930f35123e739 Mon Sep 17 00:00:00 2001 From: jac16 Date: Fri, 4 Nov 2022 18:58:18 -0400 Subject: [PATCH 13/17] Remove unused line --- .../MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 7b3e584fa13..3caea9b5fa4 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -592,12 +592,7 @@ def _get_dh_pairs(self): # Otherwise, use d_h_cutoff as a cutoff distance else: - if self.hydrogens_sel is None: - hydrogens_sel = self.guess_hydrogens() - else: - hydrogens_sel = self.hydrogens_sel - - hydrogens = self.u.select_atoms(hydrogens_sel) + hydrogens = self.u.select_atoms(self.hydrogens_sel) donors = self.u.select_atoms(self.donors_sel) donors_indices, hydrogen_indices = capped_distance( donors.positions, From 278893562cb6c91d5679da9a14030b089ebc1b33 Mon Sep 17 00:00:00 2001 From: jac16 Date: Sat, 5 Nov 2022 11:28:22 -0400 Subject: [PATCH 14/17] Integrate reviewer comments --- package/CHANGELOG | 6 ++- .../analysis/hydrogenbonds/hbond_analysis.py | 49 +++++++++++++------ 2 files changed, 38 insertions(+), 17 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index 75c2eea38de..1ee3097f8b4 100644 --- a/package/CHANGELOG +++ b/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,9 +32,10 @@ Fixes (e.g. bonds, angles) (PR #3779). Enhancements - * Update hbond analysis doc string to use exclusive bound language (#3847) - * Added ability for hbond anaylsis to use types when resnames aren't there + * 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) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 3caea9b5fa4..6a2bf6bb592 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -174,9 +174,8 @@ protein-water and protein-protein hydrogen bonds will be found, but no water-water hydrogen bonds. -One can now also define hydrogen bonds with atom types:: +One can also define hydrogen bonds with atom types:: - import MDAnalysis from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA hbonds = HBA( universe=u, @@ -272,8 +271,9 @@ def __init__(self, universe, """Set up atom selections and geometric criteria for finding hydrogen bonds in a Universe. - Hydrogen bond selections may be achieved with either a resname, atom - name combination, or when those are absent, atom types. + 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 ---------- @@ -384,13 +384,11 @@ def guess_hydrogens(self, ): """Guesses which hydrogen atoms should be used in the analysis. - Hydrogen 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 hydrogens will be identified. + :doc:`Selection string ` for atom group from which hydrogens + will be identified. (e.g., ``(resname X and name H1)`` or ``type 2``) max_mass: float (optional) The mass of a hydrogen atom must be less than this value. min_mass: float (optional) @@ -406,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. @@ -415,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: @@ -435,13 +440,11 @@ 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 contain bonding information, otherwise donor-hydrogen pairs may be incorrectly assigned. - Donor 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 donors will be identified. + :doc:`Selection string ` for atom group from which donors + will be identified. (e.g., ``(resname X and name O1)`` or ``type 2``) max_charge: float (optional) The charge of a donor atom must be less than this value. @@ -453,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. @@ -464,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 @@ -502,7 +511,8 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Parameters ---------- select: str (optional) - Selection string for atom group from which acceptors will be identified. + :doc:`Selection string ` for atom group from which acceptors will + be identified. (e.g., ``(resname X and name O1)`` or ``type 2``) max_charge: float (optional) The charge of an acceptor atom must be less than this value. @@ -514,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. @@ -524,6 +537,9 @@ 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) @@ -546,9 +562,12 @@ def __group_categories(group): ------- select : str String for each hydrogen bond acceptor/donor/hydrogen atom category. + + .. versionadded:: 2.4.0 + """ - if hasattr(group,"resnames") and hasattr(group,"names"): + 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) @@ -866,7 +885,7 @@ 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)] - if hasattr(d,"resnames"): + if hasattr(d, "resnames"): d_res = d.resnames a_res = a.resnames else: From a60760519a229ee6706c7a5268765c776a6abae7 Mon Sep 17 00:00:00 2001 From: jac16 Date: Sat, 5 Nov 2022 11:44:45 -0400 Subject: [PATCH 15/17] Add test for full coverage --- .../analysis/test_hydrogenbonds_analysis.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index 2484171cc9a..22a1ee0ddc6 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -188,6 +188,13 @@ def hydrogen_bonds(universe): h.run() return h + def test_count_by_type(self, hydrogen_bonds): + + # Only one type of hydrogen bond in this system + ref_count = 2 + + counts = hydrogen_bonds.count_by_type() + assert int(counts[0, 2]) == ref_count def test_no_bond_info_exception(self, universe): From 50dd8326f22e86a26a9738bb536eebb3bd0eb593 Mon Sep 17 00:00:00 2001 From: jac16 Date: Sun, 6 Nov 2022 07:35:55 -0500 Subject: [PATCH 16/17] Integrate reviewer comments --- .../analysis/hydrogenbonds/hbond_analysis.py | 141 +++++++++++------- 1 file changed, 83 insertions(+), 58 deletions(-) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index 6a2bf6bb592..ca09e71447a 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -325,10 +325,12 @@ 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 + .. versionchanged:: 2.4.0 + Added use of atom types + """ self.u = universe @@ -387,8 +389,9 @@ def guess_hydrogens(self, Parameters ---------- select: str (optional) - :doc:`Selection string ` for atom group from which hydrogens - will be identified. (e.g., ``(resname X and name H1)`` or ``type 2``) + :ref:`Selection string ` for atom group + from which hydrogens will be identified. (e.g., ``(resname X and + name H1)`` or ``type 2``) max_mass: float (optional) The mass of a hydrogen atom must be less than this value. min_mass: float (optional) @@ -399,24 +402,29 @@ def guess_hydrogens(self, Returns ------- potential_hydrogens: str - String containing the :attr:`resname` and :attr:`name` of all hydrogen atoms potentially capable of forming - hydrogen bonds. + String containing the :attr:`resname` and :attr:`name` of all + hydrogen atoms potentially capable of forming hydrogen bonds. 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. + 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. - If :attr:`hydrogens_sel` is `None`, this function is called to guess the selection. + If :attr:`hydrogens_sel` is `None`, this function is called to guess + the selection. - 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 + 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 @@ -434,41 +442,48 @@ def guess_hydrogens(self, )) ] - return self.__group_categories(hydrogens_ag) + 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 - contain bonding information, otherwise donor-hydrogen pairs may be incorrectly assigned. + """Guesses which atoms could be considered donors in the analysis. Only + use if the universe topology does not contain bonding information, + otherwise donor-hydrogen pairs may be incorrectly assigned. Parameters ---------- select: str (optional) - :doc:`Selection string ` for atom group from which donors - will be identified. (e.g., ``(resname X and name O1)`` or ``type 2``) + :ref:`Selection string ` for atom group + from which donors will be identified. (e.g., ``(resname X and name + O1)`` or ``type 2``) max_charge: float (optional) The charge of a donor atom must be less than this value. Returns ------- potential_donors: str - String containing the :attr:`resname` and :attr:`name` of all atoms that potentially capable of forming - hydrogen bonds. + String containing the :attr:`resname` and :attr:`name` of all atoms + that potentially capable of forming hydrogen bonds. 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. + 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. - If :attr:`donors_sel` is `None`, and the universe topology does not have bonding information, this function is - called to guess the selection. + If :attr:`donors_sel` is `None`, and the universe topology does not + have bonding information, this function is called to guess the + selection. + + Alternatively, this function may be used to quickly generate a + :class:`str` of potential donor atoms involved in hydrogen bonding. + This :class:`str` may then be modified before being used to set the + attribute :attr:`donors_sel`. - Alternatively, this function may be used to quickly generate a :class:`str` of potential donor atoms involved - 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 @@ -476,17 +491,20 @@ def guess_donors(self, select='all', max_charge=-0.5): """ # We need to know `hydrogens_sel` before we can find donors - # Use a new variable `hydrogens_sel` so that we do not set `self.hydrogens_sel` if it is currently `None` + # Use a new variable `hydrogens_sel` so that we do not set + # `self.hydrogens_sel` if it is currently `None` if self.hydrogens_sel is None: hydrogens_sel = self.guess_hydrogens() else: hydrogens_sel = self.hydrogens_sel hydrogens_ag = self.u.select_atoms(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): + # 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: @@ -500,7 +518,7 @@ def guess_donors(self, select='all', max_charge=-0.5): donors_ag = donors_ag[donors_ag.charges < max_charge] - return self.__group_categories(donors_ag) + 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. @@ -511,31 +529,36 @@ def guess_acceptors(self, select='all', max_charge=-0.5): Parameters ---------- select: str (optional) - :doc:`Selection string ` for atom group from which acceptors will - be identified. (e.g., ``(resname X and name O1)`` or ``type 2``) + :ref:`Selection string ` for atom group + from which acceptors will be identified. (e.g., ``(resname X and + name O1)`` or ``type 2``) max_charge: float (optional) The charge of an acceptor atom must be less than this value. Returns ------- potential_acceptors: str - String containing the :attr:`resname` and :attr:`name` of all atoms that potentially capable of forming - hydrogen bonds. + String containing the :attr:`resname` and :attr:`name` of all atoms + that potentially capable of forming hydrogen bonds. 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. + 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. - If :attr:`acceptors_sel` is `None`, this function is called to guess the selection. + If :attr:`acceptors_sel` is `None`, this function is called to guess + the selection. + + Alternatively, this function may be used to quickly generate a + :class:`str` of potential acceptor atoms involved in hydrogen bonding. + This :class:`str` may then be modified before being used to set the + attribute :attr:`acceptors_sel`. - Alternatively, this function may be used to quickly generate a :class:`str` of potential acceptor atoms involved - 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 @@ -545,10 +568,10 @@ def guess_acceptors(self, select='all', max_charge=-0.5): ag = self.u.select_atoms(select) acceptors_ag = ag[ag.charges < max_charge] - return self.__group_categories(acceptors_ag) + return self._group_categories(acceptors_ag) @staticmethod - def __group_categories(group): + def _group_categories(group): """ Find categories according to universe constraints Parameters @@ -563,16 +586,16 @@ def __group_categories(group): 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) - ] - ) + group_list = np.unique([ + '(resname {} and name {})'.format(r, + p) for r, p in zip(group.resnames, group.names) + ]) else: group_list = np.unique( [ @@ -588,7 +611,8 @@ def _get_dh_pairs(self): Returns ------- donors, hydrogens: AtomGroup, AtomGroup - AtomGroups corresponding to all donors and all hydrogens. AtomGroups are ordered such that, if zipped, will + AtomGroups corresponding to all donors and all hydrogens. + AtomGroups are ordered such that, if zipped, will produce a list of donor-hydrogen pairs. """ @@ -873,13 +897,14 @@ def count_by_type(self): Returns ------- counts : numpy.ndarray - Each row of the array contains the donor resname, donor atom type, acceptor resname, acceptor atom type and - the total number of times the hydrogen bond was found. + Each row of the array contains the donor resname, donor atom type, + acceptor resname, acceptor atom type and the total number of times + the hydrogen bond was found. Note ---- - Unique hydrogen bonds are determined through a consideration of the resname and atom type of the donor and - acceptor atoms in a hydrogen bond. + Unique hydrogen bonds are determined through a consideration of the + resname and atom type of the donor and acceptor atoms in a hydrogen bond. """ d = self.u.atoms[self.results.hbonds[:, 1].astype(np.intp)] @@ -889,8 +914,8 @@ def count_by_type(self): 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))] + d_res = len(d.types) * ["None"] + a_res = len(a.types) * ["None"] tmp_hbonds = np.array([d_res, d.types, a_res, a.types], dtype=str).T hbond_type, type_counts = np.unique( From d5fbe48785b37822ff5b85ffdf35f9a2c7d79dc1 Mon Sep 17 00:00:00 2001 From: jac16 Date: Fri, 11 Nov 2022 21:13:02 -0500 Subject: [PATCH 17/17] Update test and docs --- package/CHANGELOG | 11 ++++++----- .../analysis/hydrogenbonds/hbond_analysis.py | 15 ++++++++------- .../analysis/test_hydrogenbonds_analysis.py | 12 ++++++++++++ 3 files changed, 26 insertions(+), 12 deletions(-) diff --git a/package/CHANGELOG b/package/CHANGELOG index c91c192e6ba..9b087ce637b 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -19,7 +19,8 @@ The rules for this file: * 2.4.0 Fixes - * Update hbond analysis doc string to use exclusive bound language (#3847) + * Update hbond analysis doc string to use exclusive bound language + (Issue #3847) * NoDataError raised on empty atomgroup provided to `DCDReader.timeseries` was changed to ValueError (see #3901). A matching ValueError was added to `MemoryReader.timeseries`(PR #3890). @@ -37,15 +38,15 @@ Fixes Enhancements * Added ability for hbond analysis to use types when resnames are not - present (#3847) + present (Issue #3847) * Added explanatory warnings when hbond analysis doesn't find any hbonds - (#3847) + (Issue #3847) * The timeseries method for exporting coordinates from multiple frames to a NumPy array was added to ProtoReader (PR #3890) * MDAnalysis now officially supports py3.11 (Issue #3878) * LAMMPSDump Reader optionally unwraps trajectories with image flags upon - loading (#3843 - * LAMMPSDump Reader now imports velocities and forces (#3843) + loading (Issue #3843) + * LAMMPSDump Reader now imports velocities and forces (Issue #3843) * Minimum NumPy version for py3.11 is now set to 1.23.2. * Added decorator to check for an empty atom group, applied in topological attributes(Issue #3837) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index ca09e71447a..21b08a3d457 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -329,7 +329,8 @@ def __init__(self, universe, .. versionadded:: 2.0.0 Added `between` keyword .. versionchanged:: 2.4.0 - Added use of atom types + Added use of atom types in selection strings for hydrogen atoms, + bond donors, or bond acceptors """ @@ -462,7 +463,7 @@ def guess_donors(self, select='all', max_charge=-0.5): ------- potential_donors: str String containing the :attr:`resname` and :attr:`name` of all atoms - that potentially capable of forming hydrogen bonds. + that are potentially capable of forming hydrogen bonds. Notes ----- @@ -719,9 +720,9 @@ def _single_frame(self): 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}." + "No hydrogen bonds were found given d-a cutoff of " + f"{self.d_a_cutoff} between Donor, {self.donors_sel}, and " + f"Acceptor, {self.acceptors_sel}." ) # Remove D-A pairs more than d_a_cutoff away from one another @@ -749,8 +750,8 @@ def _single_frame(self): 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}." + f"{self.d_h_a_angle} between Donor, {self.donors_sel}, and " + f"Acceptor, {self.acceptors_sel}." ) # Retrieve atoms, distances and angles of hydrogen bonds diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index 22a1ee0ddc6..c7295631862 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -23,6 +23,7 @@ import numpy as np import logging import pytest +import copy from numpy.testing import (assert_allclose, assert_equal, assert_array_almost_equal, assert_array_equal, @@ -349,6 +350,17 @@ def hydrogen_bonds(universe): h.run() return h + def test_no_hydrogen_bonds(self, universe): + tmp_kwargs = copy.deepcopy(self.kwargs) + tmp_kwargs["d_h_a_angle_cutoff"] = 50 + hbonds = HydrogenBondAnalysis(universe, **tmp_kwargs) + + with pytest.warns(UserWarning, + match=("No hydrogen bonds were found given angle " + "of 50 between Donor, type O, and Acceptor," + " type O.")): + hbonds.run(step=1) + class TestHydrogenBondAnalysisBetween(object):