Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Hbond analysis update #4107

Closed
wants to merge 5 commits into from
Closed

Hbond analysis update #4107

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
f894064
changed comparisons of scalars to use pytest.approx function and comp…
yickysan Mar 24, 2023
d455e4e
reverting previous commit
yickysan Apr 1, 2023
1c49ead
Revert "reverting previous commit"
yickysan Apr 1, 2023
7c98e9f
Revert "changed comparisons of scalars to use pytest.approx function …
yickysan Apr 1, 2023
bb25dcd
hbond analysis update
yickysan Apr 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
74 changes: 51 additions & 23 deletions package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -380,7 +380,7 @@ def __init__(self, universe,
self.results.hbonds = None

def guess_hydrogens(self,
select='all',
select=None,
max_mass=1.1,
min_charge=0.3,
min_mass=0.9
Expand Down Expand Up @@ -434,18 +434,29 @@ def guess_hydrogens(self,
if min_mass > max_mass:
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,
))
]
if select:

ag = self.u.select_atoms(select)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ratther than doing a completely new code branch for each option, why not just do this selection solely in an if/else and leaving everything else the same?

hydrogens_ag = ag[
np.logical_and.reduce((
ag.masses < max_mass,
ag.charges > min_charge,
ag.masses > min_mass,
))
]
else:
ag = self.u.atoms
hydrogens_ag = ag[
np.logical_and.reduce((
ag.masses < max_mass,
ag.charges > min_charge,
ag.masses > min_mass,
))
]

return self._group_categories(hydrogens_ag)

def guess_donors(self, select='all', max_charge=-0.5):
def guess_donors(self, select=None, 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.
Expand Down Expand Up @@ -504,24 +515,37 @@ def guess_donors(self, select='all', max_charge=-0.5):
# 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))
if select:
if (hasattr(self.u._topology, 'bonds')
and len(self.u._topology.bonds.values) != 0):
donors_ag = find_hydrogen_donors(hydrogens_ag)
donors_ag = donors_ag.intersection(self.u.select_atoms(select))
else:
donors_ag = hydrogens_ag.residues.atoms.select_atoms(
"({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format(
donors_sel=select,
d_h_cutoff=self.d_h_cutoff,
hydrogens_sel=hydrogens_sel
)
)
else:
donors_ag = hydrogens_ag.residues.atoms.select_atoms(
"({donors_sel}) and around {d_h_cutoff} {hydrogens_sel}".format(
donors_sel=select,
d_h_cutoff=self.d_h_cutoff,
hydrogens_sel=hydrogens_sel
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.atoms)
else:
donors_ag = hydrogens_ag.residues.atoms.select_atoms(
"(all) and around {d_h_cutoff} {hydrogens_sel}".format(
d_h_cutoff=self.d_h_cutoff,
hydrogens_sel=hydrogens_sel
)
)
)

donors_ag = donors_ag[donors_ag.charges < max_charge]

return self._group_categories(donors_ag)

def guess_acceptors(self, select='all', max_charge=-0.5):
def guess_acceptors(self, select=None, max_charge=-0.5):
"""Guesses which atoms could be considered acceptors in the analysis.

Acceptor selections may be achieved with either a resname, atom
Expand Down Expand Up @@ -566,8 +590,12 @@ def guess_acceptors(self, select='all', max_charge=-0.5):

"""

ag = self.u.select_atoms(select)
acceptors_ag = ag[ag.charges < max_charge]
if select:
ag = self.u.select_atoms(select)
acceptors_ag = ag[ag.charges < max_charge]
else:
ag = self.u.atoms
acceptors_ag = ag[ag.charges < max_charge]

return self._group_categories(acceptors_ag)

Expand Down