diff --git a/package/CHANGELOG b/package/CHANGELOG index 678f0fcf293..d4106cdeb2e 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -39,6 +39,8 @@ Enhancements * Added computation of Mean Squared Displacements (#2438, PR #2619) * Improved performances when parsing TPR files (PR #2804) * Added converter between Cartesian and Bond-Angle-Torsion coordinates (PR #2668) + * Added Hydrogen Bond Lifetime via existing autocorrelation features (PR #2791) + * Added Hydrogen Bond Lifetime keyword "between" (PR #2791) Changes * Changes development status from Beta to Mature (Issue #2773) diff --git a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py index d662bfd1eca..0ce9062ded2 100644 --- a/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py +++ b/package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py @@ -141,13 +141,57 @@ hbonds.acceptors_sel = f"({protein_acceptors_sel}) or ({water_acceptors_sel} and around 10 not resname TIP3})" hbonds.run() -It is highly recommended that a topology with bonding information is used to generate the universe, e.g `PSF`, `TPR`, or -`PRMTOP` files. This is the only method by which it can be guaranteed that donor-hydrogen pairs are correctly identified. -However, if, for example, a `PDB` file is used instead, a :attr:`donors_sel` may be provided along with a -:attr:`hydrogens_sel` and the donor-hydrogen pairs will be identified via a distance cutoff, :attr:`d_h_cutoff`:: +To calculate the hydrogen bonds between different groups, for example a +protein and water, one can use the :attr:`between` keyword. The +following will find protein-water hydrogen bonds but not protein-protein +or water-water hydrogen bonds:: import MDAnalysis - from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis as HBA + from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import ( + HydrogenBondAnalysis as HBA) + + u = MDAnalysis.Universe(psf, trajectory) + + hbonds = HBA( + universe=u, + between=['resname TIP3', 'protein'] + ) + + protein_hydrogens_sel = hbonds.guess_hydrogens("protein") + protein_acceptors_sel = hbonds.guess_acceptors("protein") + + water_hydrogens_sel = "resname TIP3 and name H1 H2" + water_acceptors_sel = "resname TIP3 and name OH2" + + hbonds.hydrogens_sel = f"({protein_hydrogens_sel}) or ({water_hydrogens_sel}" + hbonds.acceptors_sel = f"({protein_acceptors_sel}) or ({water_acceptors_sel}" + + hbonds.run() + +It is further possible to compute hydrogen bonds between several groups with +with use of :attr:`between`. If in the above example, +`between=[['resname TIP3', 'protein'], ['protein', 'protein']]`, all +protein-water and protein-protein hydrogen bonds will be found, but +no water-water hydrogen bonds. + +In order to compute the hydrogen bond lifetime, after finding hydrogen bonds +one can use the :attr:`lifetime` function: + + ... + hbonds.run() + tau_timeseries, timeseries = hbonds.lifetime() + +It is **highly recommended** that a topology with bond information is used to +generate the universe, e.g `PSF`, `TPR`, or `PRMTOP` files. This is the only +method by which it can be guaranteed that donor-hydrogen pairs are correctly +identified. However, if, for example, a `PDB` file is used instead, a +:attr:`donors_sel` may be provided along with a :attr:`hydrogens_sel` and the +donor-hydrogen pairs will be identified via a distance cutoff, +:attr:`d_h_cutoff`:: + + import MDAnalysis + from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import ( + HydrogenBondAnalysis as HBA) u = MDAnalysis.Universe(pdb, trajectory) @@ -166,11 +210,14 @@ .. autoclass:: HydrogenBondAnalysis :members: """ -import warnings +import logging +from collections.abc import Iterable + import numpy as np from .. import base from MDAnalysis.lib.distances import capped_distance, calc_angles +from MDAnalysis.lib.correlations import autocorrelation, correct_intermittency from MDAnalysis.exceptions import NoDataError from ...due import due, Doi @@ -188,24 +235,45 @@ class HydrogenBondAnalysis(base.AnalysisBase): Perform an analysis of hydrogen bonds in a Universe. """ - def __init__(self, universe, donors_sel=None, hydrogens_sel=None, acceptors_sel=None, - d_h_cutoff=1.2, d_a_cutoff=3.0, d_h_a_angle_cutoff=150, update_selections=True): - """Set up atom selections and geometric criteria for finding hydrogen bonds in a Universe. + def __init__(self, universe, + donors_sel=None, hydrogens_sel=None, acceptors_sel=None, + between=None, d_h_cutoff=1.2, + d_a_cutoff=3.0, d_h_a_angle_cutoff=150, + update_selections=True): + """Set up atom selections and geometric criteria for finding hydrogen + bonds in a Universe. Parameters ---------- universe : Universe MDAnalysis Universe object donors_sel : str - Selection string for the hydrogen bond donor atoms. If the universe topology contains bonding information, - leave :attr:`donors_sel` as `None` so that donor-hydrogen pairs can be correctly identified. + Selection string for the hydrogen bond donor atoms. If the + universe topology contains bonding information, leave + :attr:`donors_sel` as `None` so that donor-hydrogen pairs can be + correctly identified. hydrogens_sel : str - Selection string for the hydrogen bond hydrogen atoms. Leave as `None` to guess which hydrogens to use in - the analysis using :attr:`guess_hydrogens`. If :attr:`hydrogens_sel` is left as `None`, also leave - :attr:`donors_sel` as None so that donor-hydrogen pairs can be correctly identified. + Selection string for the hydrogen bond hydrogen atoms. Leave as + `None` to guess which hydrogens to use in the analysis using + :attr:`guess_hydrogens`. If :attr:`hydrogens_sel` is left as + `None`, also leave :attr:`donors_sel` as None so that + donor-hydrogen pairs can be correctly identified. acceptors_sel : str - Selection string for the hydrogen bond acceptor atoms. Leave as `None` to guess which atoms to use in the - analysis using :attr:`guess_acceptors` + Selection string for the hydrogen bond acceptor atoms. Leave as + `None` to guess which atoms to use in the analysis using + :attr:`guess_acceptors` + between : List (optional), + Specify two selection strings for non-updating atom groups between + which hydrogen bonds will be calculated. For example, if the donor + and acceptor selections include both protein and water, it is + possible to find only protein-water hydrogen bonds - and not + protein-protein or water-water - by specifying + between=["protein", "SOL"]`. If a two-dimensional list is + passed, hydrogen bonds between each pair will be found. For + example, between=[["protein", "SOL"], ["protein", "protein"]]` + will calculate all protein-water and protein-protein hydrogen + bonds but not water-water hydrogen bonds. If `None`, hydrogen + bonds between all donors and acceptors will be calculated. d_h_cutoff : float (optional) Distance cutoff used for finding donor-hydrogen pairs [1.2]. Only used to find donor-hydrogen pairs if the universe topology does not contain bonding information @@ -219,8 +287,12 @@ def __init__(self, universe, donors_sel=None, hydrogens_sel=None, acceptors_sel= Note ---- - It is highly recommended that a universe topology with bonding information is used, as this is the only way - that guarantees the correct identification of donor-hydrogen pairs. + It is highly recommended that a universe topology with bond + information is used, as this is the only way that guarantees the + correct identification of donor-hydrogen pairs. + + .. versionadded:: 2.0.0 + Added `between` keyword """ self.u = universe @@ -228,10 +300,34 @@ def __init__(self, universe, donors_sel=None, hydrogens_sel=None, acceptors_sel= self.donors_sel = donors_sel self.hydrogens_sel = hydrogens_sel self.acceptors_sel = acceptors_sel + + # If hydrogen bonding groups are selected, then generate + # corresponding atom groups + if between is not None: + if not isinstance(between, Iterable) or len(between) == 0: + raise ValueError("between must be a non-empty list/iterable") + if isinstance(between[0], str): + between = [between] + + between_ags = [] + for group1, group2 in between: + between_ags.append( + [ + self.u.select_atoms(group1, updating=False), + self.u.select_atoms(group2, updating=False) + ] + ) + + self.between_ags = between_ags + else: + self.between_ags = None + + self.d_h_cutoff = d_h_cutoff self.d_a_cutoff = d_a_cutoff self.d_h_a_angle = d_h_a_angle_cutoff self.update_selections = update_selections + self.hbonds = None def guess_hydrogens(self, select='all', @@ -428,6 +524,40 @@ def _get_dh_pairs(self): return donors, hydrogens + def _filter_atoms(self, donors, hydrogens, acceptors): + """Filter donor, hydrogen and acceptor atoms to consider only hydrogen + bonds between two or more specified groups. + + Groups are specified with the `between` keyword when creating the + HydrogenBondAnalysis object. + + Returns + ------- + donors, hydrogens, acceptors: Filtered AtomGroups + """ + + mask = np.full(donors.n_atoms, fill_value=False) + for group1, group2 in self.between_ags: + + # Find donors in G1 and acceptors in G2 + mask[ + np.logical_and( + np.in1d(donors.indices, group1.indices), + np.in1d(acceptors.indices, group2.indices) + ) + ] = True + + # Find acceptors in G1 and donors in G2 + mask[ + np.logical_and( + np.in1d(acceptors.indices, group1.indices), + np.in1d(donors.indices, group2.indices) + ) + ] = True + + return donors[mask], hydrogens[mask], acceptors[mask] + + def _prepare(self): self.hbonds = [[], [], [], [], [], []] @@ -466,6 +596,12 @@ def _single_frame(self): tmp_hydrogens = self._hydrogens[d_a_indices.T[0]] tmp_acceptors = self._acceptors[d_a_indices.T[1]] + # Remove donor-acceptor pairs between pairs of AtomGroups we are not + # interested in + if self.between_ags is not None: + tmp_donors, tmp_hydrogens, tmp_acceptors = \ + self._filter_atoms(tmp_donors, tmp_hydrogens, tmp_acceptors) + # Find D-H-A angles greater than d_h_a_angle_cutoff d_h_a_angles = np.rad2deg( calc_angles( @@ -496,6 +632,95 @@ def _conclude(self): self.hbonds = np.asarray(self.hbonds).T + def lifetime(self, tau_max=20, window_step=1, intermittency=0): + """ + Computes and returns the time-autocorrelation (HydrogenBondLifetimes) + of hydrogen bonds. + + Before calling this method, the hydrogen bonds must first be computed + with the `run()` function. The same `start`, `stop` and `step` + parameters used in finding hydrogen bonds will be used here for + calculating hydrogen bond lifetimes. That is, the same frames will be + used in the analysis. + + Unique hydrogen bonds are identified using hydrogen-acceptor pairs. + This means an acceptor switching to a different hydrogen atom - with + the same donor - from one frame to the next is considered a different + hydrogen bond. + + Please see :func:`MDAnalysis.lib.correlations.autocorrelation` and + :func:`MDAnalysis.lib.correlations.intermittency` functions for more + details. + + + Parameters + ---------- + window_step : int, optional + The number of frames between each t(0). + tau_max : int, optional + Hydrogen bond lifetime is calculated for frames in the range + 1 <= `tau` <= `tau_max` + intermittency : int, optional + The maximum number of consecutive frames for which a bond can + disappear but be counted as present if it returns at the next + frame. An intermittency of `0` is equivalent to a continuous + autocorrelation, which does not allow for hydrogen bond + disappearance. For example, for `intermittency=2`, any given + hydrogen bond may disappear for up to two consecutive frames yet + be treated as being present at all frames. The default is + continuous (intermittency=0). + + Returns + ------- + tau_timeseries : np.array + tau from 1 to `tau_max` + timeseries : np.array + autcorrelation value for each value of `tau` + """ + + if self.hbonds is None: + logging.error( + "Autocorrelation analysis of hydrogen bonds cannot be done" + "before the hydrogen bonds are found" + ) + logging.error( + "Autocorrelation: Please use the .run() before calling this" + "function" + ) + raise NoDataError(".hbonds attribute is None: use .run() first") + + if self.step != 1: + logging.warning( + "Autocorrelation: Hydrogen bonds were computed with step > 1." + ) + logging.warning( + "Autocorrelation: We recommend recomputing hydrogen bonds with" + " step = 1." + ) + logging.warning( + "Autocorrelation: if you would like to allow bonds to break" + " and reform, please use 'intermittency'" + ) + + # Extract the hydrogen bonds IDs only in the format + # [set(superset(x1,x2), superset(x3,x4)), ..] + found_hydrogen_bonds = [set() for _ in self.frames] + for frame_index, frame in enumerate(self.frames): + for hbond in self.hbonds[self.hbonds[:, 0] == frame]: + found_hydrogen_bonds[frame_index].add(frozenset(hbond[2:4])) + + intermittent_hbonds = correct_intermittency( + found_hydrogen_bonds, + intermittency=intermittency + ) + tau_timeseries, timeseries, timeseries_data = autocorrelation( + intermittent_hbonds, + tau_max, + window_step=window_step + ) + + return np.vstack([tau_timeseries, timeseries]) + def count_by_time(self): """Counts the number of hydrogen bonds per timestep. diff --git a/package/MDAnalysis/analysis/waterdynamics.py b/package/MDAnalysis/analysis/waterdynamics.py index c5d32264b9b..fd26da5f021 100644 --- a/package/MDAnalysis/analysis/waterdynamics.py +++ b/package/MDAnalysis/analysis/waterdynamics.py @@ -508,11 +508,23 @@ class HydrogenBondLifetimes(object): .. versionchanged:: 1.0.0 The ``nproc`` keyword was removed as it linked to a portion of code that may have failed in some cases. + .. deprecated:: 1.0.1 + ``waterdynamics.HydrogenBondLifetimes`` is deprecated and will be + removed in 2.0.0. Instead, please use (available in 2.0.0) + MDAnalysis.analysis.hydrogenbonds.HydrogenBondAnalysis.lifetime + """ def __init__(self, universe, selection1, selection2, t0, tf, dtmax, nproc=1): + warnings.warn( + "This class is deprecated. " + "Instrad, please use" + "MDAnalysis.analysis.hydrogenbonds.HydrogenBondAnalysis.lifetime", + category=DeprecationWarning + ) + self.universe = universe self.selection1 = selection1 self.selection2 = selection2 diff --git a/package/MDAnalysis/lib/correlations.py b/package/MDAnalysis/lib/correlations.py index bf52dd61707..9b12823b576 100644 --- a/package/MDAnalysis/lib/correlations.py +++ b/package/MDAnalysis/lib/correlations.py @@ -223,36 +223,40 @@ def correct_intermittency(list_of_sets, intermittency): list_of_sets = deepcopy(list_of_sets) - for i, ids in enumerate(list_of_sets): + # an element (a superset) takes the form of: + # - an atom pair when computing hydrogen bonds lifetime, + # - atom ID in the case of water survival probability, + for i, elements in enumerate(list_of_sets): # initially update each frame as seen 0 ago (now) - seen_frames_ago = {i: 0 for i in ids} + seen_frames_ago = {i: 0 for i in elements} for j in range(1, intermittency + 2): - for atomid in seen_frames_ago.keys(): + for element in seen_frames_ago.keys(): # no more frames if i + j >= len(list_of_sets): continue - # if the atom is absent now - if not atomid in list_of_sets[i + j]: + # if the element is absent now + if element not in list_of_sets[i + j]: # increase its absence counter - seen_frames_ago[atomid] += 1 + seen_frames_ago[element] += 1 continue - # the atom is found - if seen_frames_ago[atomid] == 0: - # the atom was present in the last frame + # the element is present + if seen_frames_ago[element] == 0: + # the element was present in the last frame continue - # it was absent more times than allowed - if seen_frames_ago[atomid] > intermittency: + # element was absent more times than allowed + if seen_frames_ago[element] > intermittency: continue - # the atom was absent but returned (within <= intermittency_value) - # add it to the frames where it was absent. + # The element was absent but returned, i.e. + # within <= intermittency_value. + # Add it to the frames where it was absent. # Introduce the corrections. - for k in range(seen_frames_ago[atomid], 0, -1): - list_of_sets[i + j - k].add(atomid) + for k in range(seen_frames_ago[element], 0, -1): + list_of_sets[i + j - k].add(element) - seen_frames_ago[atomid] = 0 + seen_frames_ago[element] = 0 return list_of_sets diff --git a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py index 80356bcd7a6..d107f9932fe 100644 --- a/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py +++ b/testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py @@ -21,14 +21,17 @@ # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 # import numpy as np -import MDAnalysis -from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import HydrogenBondAnalysis -from MDAnalysis.exceptions import NoDataError - +import logging import pytest + from numpy.testing import (assert_allclose, assert_equal, assert_array_almost_equal, assert_array_equal, assert_almost_equal) + +import MDAnalysis +from MDAnalysis.analysis.hydrogenbonds.hbond_analysis import ( + HydrogenBondAnalysis) +from MDAnalysis.exceptions import NoDataError from MDAnalysisTests.datafiles import waterPSF, waterDCD @@ -101,7 +104,7 @@ def test_count_by_ids(self, h, universe): assert_allclose(counts, ref_counts) -class TestHydrogenBondAnalysisMock(object): +class TestHydrogenBondAnalysisIdeal(object): kwargs = { 'donors_sel': 'name O', @@ -166,6 +169,17 @@ def universe(): return u + @staticmethod + @pytest.fixture(scope='class') + def hydrogen_bonds(universe): + h = HydrogenBondAnalysis( + universe, + **TestHydrogenBondAnalysisIdeal.kwargs + ) + h.run() + return h + + def test_no_bond_info_exception(self, universe): kwargs = { @@ -181,32 +195,181 @@ def test_no_bond_info_exception(self, universe): 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 + def test_first(self, hydrogen_bonds): + assert len(hydrogen_bonds.hbonds) == 2 - frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, dha_angle = h.hbonds[0] + frame_no, donor_index, hydrogen_index, acceptor_index, da_dst, angle =\ + hydrogen_bonds.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() + 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 = h.count_by_time() - assert_array_almost_equal(h.times, ref_times) + 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): + + kwargs = { + 'donors_sel': 'name O P', + 'hydrogens_sel': 'name H1 H2 PH', + 'acceptors_sel': 'name O P', + '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 and two "protein" molecules + # P1-PH1 are the two atoms that comprise the toy protein PROT1 + """ + H4 + \ + O1-H2 .... O2-H3 ... P1-PH1 ... P2-PH2 + / + H1 + """ + n_residues = 4 + n_sol_residues = 2 + n_prot_residues = 2 + u = MDAnalysis.Universe.empty( + n_atoms=n_sol_residues*3 + n_prot_residues*2, + n_residues=n_residues, + atom_resindex=[0, 0, 0, 1, 1, 1, 2, 2, 3, 3], + residue_segindex=[0, 0, 1, 1], + trajectory=True, # necessary for adding coordinates + ) + + u.add_TopologyAttr( + 'name', + ['O', 'H1', 'H2'] * n_sol_residues + ['P', 'PH'] * n_prot_residues + ) + u.add_TopologyAttr( + 'type', + ['O', 'H', 'H'] * n_sol_residues + ['P', 'PH'] * n_prot_residues + ) + u.add_TopologyAttr( + 'resname', + ['SOL'] * n_sol_residues + ['PROT'] * n_prot_residues + ) + u.add_TopologyAttr( + 'resid', + list(range(1, n_residues + 1)) + ) + u.add_TopologyAttr( + 'id', + list(range(1, (n_sol_residues * 3 + n_prot_residues * 2) + 1)) + ) + + # Atomic coordinates with hydrogen bonds between: + # O1-H2---O2 + # O2-H3---P1 + # P1-PH1---P2 + pos = 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 + [5.5, 0, 0], # P1 + [6.5, 0, 0], # PH1 + [8.5, 0, 0], # P2 + [9.5, 0, 0], # PH2 + ]) + + coordinates = np.empty((1, # number of frames + u.atoms.n_atoms, + 3)) + coordinates[0] = pos + u.load_new(coordinates, order='fac') + + return u + + def test_between_all(self, universe): + # don't specify groups between which to find hydrogen bonds + hbonds = HydrogenBondAnalysis(universe, between=None, **self.kwargs) + hbonds.run() + + # indices of [donor, hydrogen, acceptor] for each hydrogen bond + expected_hbond_indices = [ + [0, 2, 3], # water-water + [3, 4, 6], # protein-water + [6, 7, 8] # protein-protein + ] + assert_array_equal(hbonds.hbonds[:, 1:4], expected_hbond_indices) + + def test_between_PW(self, universe): + # Find only protein-water hydrogen bonds + hbonds = HydrogenBondAnalysis( + universe, + between=["resname PROT", "resname SOL"], + **self.kwargs + ) + hbonds.run() + + # indices of [donor, hydrogen, acceptor] for each hydrogen bond + expected_hbond_indices = [ + [3, 4, 6] # protein-water + ] + assert_array_equal(hbonds.hbonds[:, 1:4], expected_hbond_indices) + + def test_between_PW_PP(self, universe): + # Find protein-water and protein-protein hydrogen bonds (not + # water-water) + hbonds = HydrogenBondAnalysis( + universe, + between=[ + ["resname PROT", "resname SOL"], + ["resname PROT", "resname PROT"] + ], + **self.kwargs + ) + hbonds.run() + + # indices of [donor, hydrogen, acceptor] for each hydrogen bond + expected_hbond_indices = [ + [3, 4, 6], # protein-water + [6, 7, 8] # protein-protein + ] + assert_array_equal(hbonds.hbonds[:, 1:4], expected_hbond_indices) + class TestHydrogenBondAnalysisTIP3P_GuessAcceptors_GuessHydrogens_UseTopology_(TestHydrogenBondAnalysisTIP3P): """Uses the same distance and cutoff hydrogen bond criteria as :class:`TestHydrogenBondAnalysisTIP3P`, so the @@ -358,4 +521,4 @@ def test_count_by_type(self, h): ref_count = 15 counts = h.count_by_type() - assert int(counts[0, 2]) == ref_count \ No newline at end of file + assert int(counts[0, 2]) == ref_count