Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
3 changes: 3 additions & 0 deletions package/CHANGELOG
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,9 @@ Fixes
array warnings (PR #3139, backports PRs #2845 and #2834)
* Fixed several issues with NSGrid and triclinic boxes not finding some pairs.
(Issues #2229 #2345 #2919, PR #2937).
* Removed deprecation warning from numpy in hbond_autocorrel
(Issue #2987 PR #3242)
* Fixed bug in exclusion matrix of hbond_autocorrel (PR #3242)

Changes
* Maximum pinned versions in setup.py removed for python 3.6+ (PR #3139)
Expand Down
20 changes: 11 additions & 9 deletions package/MDAnalysis/analysis/hbonds/hbond_autocorrel.py
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@


The keyword **exclusions** allows a tuple of array addresses to be provided,
(Hidx, Aidx),these pairs of hydrogen-acceptor are then not permitted to be
(Hidx, Aidx), these pairs of hydrogen-acceptor are then not permitted to be
counted as part of the analysis. This could be used to exclude the
consideration of hydrogen bonds within the same functional group, or to perform
analysis on strictly intermolecular hydrogen bonding.
Expand Down Expand Up @@ -215,6 +215,7 @@

from MDAnalysis.lib.log import ProgressBar
from MDAnalysis.lib.distances import capped_distance, calc_angles, calc_bonds
from MDAnalysis.lib._cutil import _in2d
from MDAnalysis.core.groups import requires

from MDAnalysis.due import due, Doi
Expand Down Expand Up @@ -320,11 +321,13 @@ def __init__(self, universe,
raise ValueError("Donors and Hydrogen groups must be identical "
"length. Try using `find_hydrogen_donors`.")

self.exclusions = exclusions
if self.exclusions:
if not len(self.exclusions[0]) == len(self.exclusions[1]):
if exclusions is not None:
if len(exclusions[0]) != len(exclusions[1]):
raise ValueError(
"'exclusion' must be two arrays of identical length")
self.exclusions = np.column_stack((exclusions[0], exclusions[1])).astype(np.intp)
else:
self.exclusions = None

self.bond_type = bond_type
if self.bond_type not in ['continuous', 'intermittent']:
Expand Down Expand Up @@ -427,11 +430,10 @@ def _single_run(self, start, stop):
box = self.u.dimensions if self.pbc else None

# 2d array of all distances
pair, d = capped_distance(self.h.positions, self.a.positions, max_cutoff=self.d_crit, box=box)
if self.exclusions:
# set to above dist crit to exclude
exclude = np.column_stack((self.exclusions[0], self.exclusions[1]))
pair = np.delete(pair, np.where(pair==exclude), 0)
pair = capped_distance(self.h.positions, self.a.positions, max_cutoff=self.d_crit, box=box,
return_distances=False)
if self.exclusions is not None:
pair = pair[~ _in2d(pair, self.exclusions)]
Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

I'm surprised there was no suitable numpy operation for this :/

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

The below passes the analysis suite locally without using a custom Cython function. Maybe double check with Richard if his mixture of Cython/C++ is more efficient -- probably any advantage is related to the integer only nature of the data structure?

--- a/package/MDAnalysis/analysis/hbonds/hbond_autocorrel.py
+++ b/package/MDAnalysis/analysis/hbonds/hbond_autocorrel.py
@@ -210,6 +210,7 @@ from six import raise_from
 
 import numpy as np
 import scipy.optimize
+from scipy.spatial.distance import cdist
 
 import warnings
 
@@ -433,7 +434,7 @@ class HydrogenBondAutoCorrel(object):
         pair = capped_distance(self.h.positions, self.a.positions, max_cutoff=self.d_crit, box=box,
                                return_distances=False)
         if not self.exclusions is None:
-            pair = pair[~ _in2d(pair, self.exclusions)]
+            pair = pair[~ (cdist(pair, self.exclusions)==0).any(axis=1)]
 
         hidx, aidx = np.transpose(pair)

Copy link
Copy Markdown
Member Author

Choose a reason for hiding this comment

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

You could use cdist but it's odd to use an O(n2) solution when you can just use a set

Copy link
Copy Markdown
Member

@tylerjereddy tylerjereddy Apr 26, 2021

Choose a reason for hiding this comment

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

I don't have a strong view on it--if this is genuinely a performance-critical part of the code maybe the custom maintenance burden is worth it.


hidx, aidx = np.transpose(pair)

Expand Down
Loading