Skip to content

Commit

Permalink
Amended hbond test for full coverage, added _get_dh_pairs ability to …
Browse files Browse the repository at this point in the history
…find hydrogen atoms.
  • Loading branch information
jaclark5 committed Sep 27, 2022
1 parent f7f8bf5 commit 299b898
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 1 deletion.
7 changes: 6 additions & 1 deletion package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down
34 changes: 34 additions & 0 deletions testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -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 =\
Expand Down

0 comments on commit 299b898

Please sign in to comment.