Skip to content

Commit

Permalink
gsoc-1
Browse files Browse the repository at this point in the history
  • Loading branch information
lunamorrow committed Mar 27, 2024
1 parent 4daea21 commit 670781e
Show file tree
Hide file tree
Showing 2 changed files with 53 additions and 47 deletions.
1 change: 1 addition & 0 deletions package/AUTHORS
Original file line number Diff line number Diff line change
Expand Up @@ -236,6 +236,7 @@ Chronological list of authors
2024
- Aditya Keshari
- Philipp Stärk
- Luna Morrow

External code
-------------
Expand Down
99 changes: 52 additions & 47 deletions package/MDAnalysis/analysis/hydrogenbonds/hbond_analysis.py
Original file line number Diff line number Diff line change
Expand Up @@ -618,50 +618,52 @@ def _get_dh_pairs(self):
produce a list of donor-hydrogen pairs.
"""

# # If donors_sel is not provided, use topology to find d-h pairs
# if self.donors_sel is None:

# # 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 not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0):
# raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. '
# '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.u.select_atoms(self.hydrogens_sel)
# 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:

# hydrogens = self.guess_hydrogens()
# donors = self.guess_donors()
# 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]

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]
# If donors_sel is not provided, use topology to find d-h pairs
if self.donors_sel is None:

# 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 not (hasattr(self.u._topology, 'bonds') and len(self.u._topology.bonds.values) != 0):
raise NoDataError('Cannot assign donor-hydrogen pairs via topology as no bond information is present. '
'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
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:

# 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,
max_cutoff=self.d_h_cutoff,
box=self.u.dimensions,
return_distances=False
).T

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

Expand Down Expand Up @@ -710,17 +712,20 @@ 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()
# 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)

# Select atom groups
Expand Down

0 comments on commit 670781e

Please sign in to comment.