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.

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

add ignore_same_residue kwarg to InterRDF #4161

Merged
merged 15 commits into from
Jul 19, 2023
Merged
Show file tree
Hide file tree
Changes from 9 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: 2 additions & 1 deletion package/CHANGELOG
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ The rules for this file:
* release numbers follow "Semantic Versioning" http://semver.org

------------------------------------------------------------------------------
??/??/?? IAlibay, jaclark5, MohitKumar020291
??/??/?? IAlibay, jaclark5, MohitKumar020291, orionarcher

* 2.6.0

Expand All @@ -23,6 +23,7 @@ Fixes
(PR #4163, Issue #4159)

Enhancements
* Added a `exclude_same` kwarg to InterRDF (PR #4161)
* LAMMPSDump Reader use of `continuous` option to stitch trajectories
(Issue #3546)

Expand Down
23 changes: 23 additions & 0 deletions package/MDAnalysis/analysis/rdf.py
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,10 @@ class InterRDF(AnalysisBase):

exclusion_block : tuple
A tuple representing the tile to exclude from the distance array.
exclude_same : str
Will exclude pairs of atoms that share the same "residue", "segment", or "chain".
Those are the only valid values. This is intended to remove atoms that are
spatially correlated due to direct bonded connections.
verbose : bool
Show detailed progress of the calculation if set to `True`

Expand Down Expand Up @@ -220,6 +224,7 @@ def __init__(self,
range=(0.0, 15.0),
norm="rdf",
exclusion_block=None,
exclude_same=None,
**kwargs):
super(InterRDF, self).__init__(g1.universe.trajectory, **kwargs)
self.g1 = g1
Expand All @@ -229,6 +234,17 @@ def __init__(self,
self.rdf_settings = {'bins': nbins,
'range': range}
self._exclusion_block = exclusion_block
if exclude_same is not None and exclude_same not in ['residue', 'segment', 'chain']:
orionarcher marked this conversation as resolved.
Show resolved Hide resolved
raise ValueError(
"The exclude_same argument to InterRDF must be None, 'residue', 'segment' "
"or 'chain'."
)
if exclude_same is not None and exclusion_block is not None:
raise ValueError(
"The exclude_same argument to InterRDF cannot be used with exclusion_block."
)
name_to_attr = {'residue': 'resindices', 'segment': 'segindices', 'chain': 'chainIDs'}
self.exclude_same = name_to_attr.get(exclude_same)

if self.norm not in ['rdf', 'density', 'none']:
raise ValueError(f"'{self.norm}' is an invalid norm. "
Expand Down Expand Up @@ -261,6 +277,13 @@ def _single_frame(self):
mask = np.where(idxA != idxB)[0]
dist = dist[mask]

if self.exclude_same is not None:
# Ignore distances between atoms in the same attribute
attr_ix_a = getattr(self.g1, self.exclude_same)[pairs[:, 0]]
attr_ix_b = getattr(self.g2, self.exclude_same)[pairs[:, 1]]
mask = np.where(attr_ix_a != attr_ix_b)[0]
dist = dist[mask]

count, _ = np.histogram(dist, **self.rdf_settings)
self.results.count += count

Expand Down
15 changes: 15 additions & 0 deletions testsuite/MDAnalysisTests/analysis/test_rdf.py
Original file line number Diff line number Diff line change
Expand Up @@ -94,6 +94,21 @@ def test_exclusion(sels):
assert rdf.results.count.sum() == 4


@pytest.mark.parametrize("attr", ("residue"))
orionarcher marked this conversation as resolved.
Show resolved Hide resolved
def test_ignore_same_residues(sels, attr):
# should see two distances with 4 counts each
s1, s2 = sels
rdf = InterRDF(s2, s2, exclude_same=attr).run()
assert rdf.rdf[0] == 0
assert rdf.results.count.sum() == 8


def test_ignore_same_residues_fails(sels):
# should see two distances with 4 counts each
orionarcher marked this conversation as resolved.
Show resolved Hide resolved
orionarcher marked this conversation as resolved.
Show resolved Hide resolved
s1, s2 = sels
with pytest.raises(ValueError):
orionarcher marked this conversation as resolved.
Show resolved Hide resolved
rdf = InterRDF(s2, s2, exclude_same="unsupported").run()

orionarcher marked this conversation as resolved.
Show resolved Hide resolved
@pytest.mark.parametrize("attr", ("rdf", "bins", "edges", "count"))
def test_rdf_attr_warning(sels, attr):
s1, s2 = sels
Expand Down