Skip to content

Replace the distance search of finding hydrogens with a dictionary lookup #2519

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

Merged
richardjgowers merged 3 commits into from
Feb 15, 2020
Merged

Replace the distance search of finding hydrogens with a dictionary lookup #2519

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2b2d497
Replace the distance search of finding hydrogens with a dictionary lo…
Feb 9, 2020
3cc2620
bug fix
Feb 9, 2020
2791289
small fix
Feb 10, 2020
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
50 changes: 35 additions & 15 deletions package/MDAnalysis/analysis/hbonds/wbridge_analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -820,6 +820,29 @@ def _log_parameters(self):
logger.info("WaterBridgeAnalysis: force field %s to guess donor and \
acceptor names", self.forcefield)

def _build_residue_dict(self, selection):
# Build the residue_dict where the key is the residue name
# The content is a dictionary where hydrogen bond donor heavy atom names is the key
# The content is the hydrogen bond donor hydrogen atom names
atom_group = self.u.select_atoms(selection)
for residue in atom_group.residues:
if not residue.resname in self._residue_dict:
self._residue_dict[residue.resname] = defaultdict(set)
for atom in residue.atoms:
if atom.name in self.donors:
self._residue_dict[residue.resname][atom.name].update(self._get_bonded_hydrogens(atom).names)

def _update_donor_h(self, atom_ix, h_donors, donors_h):
Copy link
Member

@richardjgowers richardjgowers Feb 12, 2020

Choose a reason for hiding this comment

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

this needs a docstring, it's not clear who/what/where atom_ix is

atom = self.u.atoms[atom_ix]
residue = atom.residue
hydrogen_names = self._residue_dict[residue.resname][atom.name]
if hydrogen_names:
hydrogens = residue.atoms.select_atoms('name {0}'.format(
' '.join(hydrogen_names))).ix
for atom in hydrogens:
h_donors[atom] = atom_ix
donors_h[atom_ix].append(atom)

def _update_selection(self):
self._s1_donors = []
self._s1_h_donors = {}
Expand Down Expand Up @@ -861,11 +884,8 @@ def _update_selection(self):
if self.selection1_type in ('donor', 'both'):
self._s1_donors = self.u.atoms[self._s1].select_atoms(
'name {0}'.format(' '.join(self.donors))).ix
for i, atom_ix in enumerate(self._s1_donors):
tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix
for atom in tmp:
self._s1_h_donors[atom] = atom_ix
self._s1_donors_h[atom_ix].append(atom)
for atom_ix in self._s1_donors:
self._update_donor_h(atom_ix, self._s1_h_donors, self._s1_donors_h)
self.logger_debug("Selection 1 donors: {0}".format(len(self._s1_donors)))
self.logger_debug("Selection 1 donor hydrogens: {0}".format(len(self._s1_h_donors)))
if self.selection1_type in ('acceptor', 'both'):
Expand All @@ -882,11 +902,8 @@ def _update_selection(self):
if self.selection1_type in ('acceptor', 'both'):
self._s2_donors = self.u.atoms[self._s2].select_atoms(
'name {0}'.format(' '.join(self.donors))).ix
for i, atom_ix in enumerate(self._s2_donors):
tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix
for atom in tmp:
self._s2_h_donors[atom] = atom_ix
self._s2_donors_h[atom_ix].append(atom)
for atom_ix in self._s2_donors:
self._update_donor_h(atom_ix, self._s2_h_donors, self._s2_donors_h)
self.logger_debug("Selection 2 donors: {0:d}".format(len(self._s2_donors)))
self.logger_debug("Selection 2 donor hydrogens: {0:d}".format(len(self._s2_h_donors)))

Expand Down Expand Up @@ -914,11 +931,8 @@ def _update_water_selection(self):
else:
self._water_donors = self.u.atoms[self._water].select_atoms(
'name {0}'.format(' '.join(self.donors))).ix
for i, atom_ix in enumerate(self._water_donors):
tmp = self._get_bonded_hydrogens(self.u.atoms[atom_ix]).ix
for atom in tmp:
self._water_h_donors[atom] = atom_ix
self._water_donors_h[atom_ix].append(atom)
for atom_ix in self._water_donors:
self._update_donor_h(atom_ix, self._water_h_donors, self._water_donors_h)
self.logger_debug("Water donors: {0}".format(len(self._water_donors)))
self.logger_debug("Water donor hydrogens: {0}".format(len(self._water_h_donors)))
self._water_acceptors = self.u.atoms[self._water].select_atoms(
Expand Down Expand Up @@ -969,11 +983,17 @@ def _prepare(self):
self.selection_distance = (2 * 2 + self.order * (2 + self.distance))

self.box = self.u.dimensions if self.pbc else None
self._residue_dict = {}
self._build_residue_dict(self.selection1)
self._build_residue_dict(self.selection2)
self._build_residue_dict(self.water_selection)

self._update_selection()

self.timesteps = []
if len(self._s1) and len(self._s2):
self._update_water_selection()
else:
logger.info("WaterBridgeAnalysis: no atoms found in the selection.")

logger.info("WaterBridgeAnalysis: initial checks passed.")
Expand Down