Skip to content

Commit

Permalink
updated-hbond-analysis-with-atomgroups
Browse files Browse the repository at this point in the history
  • Loading branch information
@lunamorrow
lunamorrow committed Apr 16, 2024
1 parent 670781e commit 24b5371
Show file tree
Hide file tree
Showing 2 changed files with 77 additions and 68 deletions.
135 changes: 71 additions & 64 deletions package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -406,9 +406,9 @@ def guess_hydrogens(self,
Returns
-------
potential_hydrogens: AtomGroup
AtomGroup corresponding to all hydrogen atoms potentially capable
of forming hydrogen bonds.
potential_hydrogens: str
String containing the :attr:`resname` and :attr:`name` of all
hydrogen atoms potentially capable of forming hydrogen bonds.
Notes
-----
Expand All @@ -425,9 +425,9 @@ def guess_hydrogens(self,
the selection.
Alternatively, this function may be used to quickly generate a
:class:`AtomGroup` of potential hydrogen atoms involved in hydrogen
bonding. This :class:`AtomGroup` may then be modified before being used
to set the attribute :attr:`hydrogens_sel`.
: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
Expand All @@ -439,13 +439,17 @@ def guess_hydrogens(self,
raise ValueError("min_mass is higher than (or equal to) max_mass")

ag = self.u.select_atoms(select)
hydrogens_ag = ag[
np.logical_and.reduce((
ag.masses < max_mass,
ag.charges > min_charge,
ag.masses > min_mass,
))
]
# hydrogens_ag = ag[
# np.logical_and.reduce((
# ag.masses < max_mass,
# ag.charges > min_charge,
# ag.masses > min_mass,
# ))
# ]
hydrogens_ag = ag.select_atoms(f"prop charge > {min_charge} and prop \
mass > {min_mass} and prop mass < \
{max_mass}",
updating=self.update_selections)

return hydrogens_ag

Expand All @@ -465,9 +469,9 @@ def guess_donors(self, select='all', max_charge=-0.5):
Returns
-------
potential_donors: AtomGroup
AtomGroup corresponding to all atoms potentially capable of forming
hydrogen bonds.
potential_donors: str
String containing the :attr:`resname` and :attr:`name` of all atoms
that are potentially capable of forming hydrogen bonds.
Notes
-----
Expand All @@ -485,9 +489,9 @@ def guess_donors(self, select='all', max_charge=-0.5):
selection.
Alternatively, this function may be used to quickly generate a
:class:`AtomGroup` of potential donor atoms involved in hydrogen
bonding. This :class:`AtomGroup` may then be modified before being used
to set the attribute :attr:`donors_sel`.
: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
Expand All @@ -498,11 +502,12 @@ 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`
hydrogens_ag = self.guess_hydrogens()
if self.hydrogens_sel is None:
hydrogens_sel = self._group_categories(hydrogens_ag)
hydrogens_sel = self._group_categories(self.guess_hydrogens())
else:
hydrogens_sel = self.hydrogens_sel
hydrogens_ag = self.u.select_atoms(hydrogens_sel,
updating=self.update_selections)

# 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
Expand All @@ -521,7 +526,12 @@ def guess_donors(self, select='all', max_charge=-0.5):
)
)

return donors_ag[donors_ag.charges < max_charge]
#donors_ag = donors_ag[donors_ag.charges < max_charge]

donors_ag = donors_ag.select_atoms(f"prop charge < {max_charge}",
updating=self.update_selections)

return donors_ag

def guess_acceptors(self, select='all', max_charge=-0.5):
"""Guesses which atoms could be considered acceptors in the analysis.
Expand All @@ -540,9 +550,9 @@ def guess_acceptors(self, select='all', max_charge=-0.5):
Returns
-------
potential_acceptors: AtomGroup
AtomGroup corresponding to all atoms potentially capable of forming
hydrogen bonds.
potential_acceptors: str
String containing the :attr:`resname` and :attr:`name` of all atoms
that potentially capable of forming hydrogen bonds.
Notes
-----
Expand All @@ -558,19 +568,22 @@ def guess_acceptors(self, select='all', max_charge=-0.5):
the selection.
Alternatively, this function may be used to quickly generate a
:class:`AtomGroup` of potential acceptor atoms involved in hydrogen
bonding. This :class:`AtomGroup` may then be modified before being used
to set the attribute :attr:`acceptors_sel`.
: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
"""

ag = self.u.select_atoms(select, updating=self.update_selections)
ag = self.u.select_atoms(select)
#acceptors_ag = ag[ag.charges < max_charge]
acceptors_ag = ag.select_atoms(f"prop charge < {max_charge}",
updating=self.update_selections)

return ag[ag.charges < max_charge]
return acceptors_ag

@staticmethod
def _group_categories(group):
Expand Down Expand Up @@ -629,18 +642,24 @@ def _get_dh_pairs(self):
'Please either: load a topology file with bond information; use the guess_bonds() '
'topology guesser; or set HydrogenBondAnalysis.donors_sel so that a distance cutoff '
'can be used.')

hydrogens = self.hydrogens_ag
if self.hydrogens_sel is None:
hydrogens = self.guess_hydrogens()
else:
hydrogens = self.u.select_atoms(self.hydrogens_sel,
updating=self.update_selections)
donors = sum(h.bonded_atoms[0] for h in hydrogens) if hydrogens \
else AtomGroup([], self.u)

# Otherwise, use d_h_cutoff as a cutoff distance
else:

if self.hydrogens_sel is None:
hydrogens = self.guess_hydrogens()
else:
hydrogens = self.u.select_atoms(self.hydrogens_sel,
updating=self.update_selections)
donors = self.u.select_atoms(self.donors_sel)
# hydrogens = self.guess_hydrogens()
# donors = self.guess_donors()
hydrogens = self.hydrogens_ag #NOTE: do I need this?
donors = self.donors_ag #NOTE: do I need this?
donors_indices, hydrogen_indices = capped_distance(
donors.positions,
hydrogens.positions,
Expand All @@ -652,19 +671,6 @@ def _get_dh_pairs(self):
donors = donors[donors_indices]
hydrogens = hydrogens[hydrogen_indices]

# hydrogens = self.hydrogens_ag #NOTE: do I need this?
# donors = self.donors_ag #NOTE: do I need this?
# donors_indices, hydrogen_indices = capped_distance(
# donors.positions,
# hydrogens.positions,
# max_cutoff=self.d_h_cutoff,
# box=self.u.dimensions,
# return_distances=False
# ).T

# donors = donors[donors_indices]
# hydrogens = hydrogens[hydrogen_indices]

return donors, hydrogens

def _filter_atoms(self, donors, acceptors):
Expand Down Expand Up @@ -711,25 +717,18 @@ def _filter_atoms(self, donors, acceptors):
def _prepare(self):
self.results.hbonds = [[], [], [], [], [], []]

# Set unpaired acceptor and hydrogen AtomGroups
# self.acceptors_ag = self.guess_acceptors() #NOTE: do I need this?
# self.hydrogens_ag = self.guess_hydrogens()
# self.donors_ag = self.guess_donors()

# Set atom selections if they have not been provided
if self.acceptors_sel is None:
self.acceptors_ag = self.guess_acceptors()
self.acceptors_sel = self._group_categories(self.acceptors_ag) #NOTE: do I need this?
if self.hydrogens_sel is None:
self.hydrogens_ag = self.guess_hydrogens()
self.hydrogens_sel = self._group_categories(self.hydrogens_ag)
#NOTE: add in self.donors_sel? Is this not done for a reason?
if self.donors_sel is None:
self.donors_ag = self.guess_donors()
self.donors_sel = self._group_categories(self.donors_ag)
# if self.acceptors_sel is None:
# self.acceptors_sel = self.guess_acceptors()
# if self.hydrogens_sel is None:
# self.hydrogens_sel = self.guess_hydrogens()

# Select atom groups
self._acceptors = self.guess_acceptors()
if self.acceptors_sel == None:
self._acceptors = self.guess_acceptors()
else:
self._acceptors = self.u.select_atoms(self.acceptors_sel,
updating=self.update_selections)
self._donors, self._hydrogens = self._get_dh_pairs()

def _single_frame(self):
Expand All @@ -742,6 +741,8 @@ def _single_frame(self):

# find D and A within cutoff distance of one another
# min_cutoff = 1.0 as an atom cannot form a hydrogen bond with itself
print(f"D POS: {self._donors.positions}")
print(f"A POS: {self._acceptors.positions}")
d_a_indices, d_a_distances = capped_distance(
self._donors.positions,
self._acceptors.positions,
Expand All @@ -760,6 +761,9 @@ def _single_frame(self):

# Remove D-A pairs more than d_a_cutoff away from one another
tmp_donors = self._donors[d_a_indices.T[0]]
print(f"TMP DONORS INITIAL: {tmp_donors}")
print(f"DONORS: {self._donors}")
print(f"DA INDICIES: {d_a_indices.T[0]}")
tmp_hydrogens = self._hydrogens[d_a_indices.T[0]]
tmp_acceptors = self._acceptors[d_a_indices.T[1]]

Expand Down Expand Up @@ -792,6 +796,9 @@ def _single_frame(self):

# Retrieve atoms, distances and angles of hydrogen bonds
hbond_donors = tmp_donors[hbond_indices]
print(f"HBOND DONORS: {hbond_donors}")
print(f"TMP DONORS: {tmp_donors}")
print(f"HBOND INDICIES: {hbond_indices}")
hbond_hydrogens = tmp_hydrogens[hbond_indices]
hbond_acceptors = tmp_acceptors[hbond_indices]
hbond_distances = d_a_distances[hbond_indices]
Expand Down
10 changes: 6 additions & 4 deletions testsuite/MDAnalysisTests/analysis/test_hydrogenbonds_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -227,7 +227,7 @@ 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)
donors = u.select_atoms(h.guess_donors())
donors = h.guess_donors()

def test_first_hbond(self, hydrogen_bonds):
assert len(hydrogen_bonds.results.hbonds) == 2
Expand Down Expand Up @@ -554,7 +554,8 @@ def test_guess_donors(self, h):

ref_donors = "(resname TIP3 and name OH2)"
donors = h.guess_donors(select='all', max_charge=-0.5)
assert donors == ref_donors
donors_sel = h._group_categories(donors)
assert donors_sel == ref_donors


class TestHydrogenBondAnalysisTIP3P_GuessHydrogens_NoTopology(object):
Expand Down Expand Up @@ -586,7 +587,8 @@ def test_guess_hydrogens(self, h):

ref_hydrogens = "(resname TIP3 and name H1) or (resname TIP3 and name H2)"
hydrogens = h.guess_hydrogens(select='all')
assert hydrogens == ref_hydrogens
hydrogens_sel = h._group_categories(hydrogens)
assert hydrogens_sel == ref_hydrogens

pytest.mark.parametrize(
"min_mass, max_mass, min_charge",
Expand All @@ -599,7 +601,7 @@ def test_guess_hydrogens(self, h):
def test_guess_hydrogens_empty_selection(self, h):

hydrogens = h.guess_hydrogens(select='all', min_charge=1.0)
assert hydrogens == ""
assert len(hydrogens) == 0

def test_guess_hydrogens_min_max_mass(self, h):

Expand Down

0 comments on commit 24b5371

Please sign in to comment.