Skip to content

Commit

Permalink
changed InterRDF_s density keyword semantics
Browse files Browse the repository at this point in the history
- fix (because it was not documented before)
- density=True now means "calculate the single site density";
  default is still "calculate RDF" (now meaning density=False)
- update docs
- update tests (add test for default behavior)
- update CHANGELOG
- backported issue #2811 fix (directly rebased onto master
  instead of added to PR #2798
  • Loading branch information
orbeckst committed Jul 12, 2020
1 parent a2a9a5c commit dbd39c2
Show file tree
Hide file tree
Showing 3 changed files with 31 additions and 15 deletions.
6 changes: 4 additions & 2 deletions package/CHANGELOG
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,10 @@ The rules for this file:
Fixes
* Development status changed from beta to mature (Issue #2773)
* pip installation only requests Python 2.7-compatible packages (#2736)
* fixed InterRDF_s function now gives correct results if density=True.
(Issue #2811, PR #2812)
* rdf.InterRDF_s density keyword documented and now gives correct results for
density=True; the keyword was available since 0.19.0 but with incorrect
semantics and not documented and did not produce correct results (Issue
#2811, PR #2812)


06/09/20 richardjgowers, kain88-de, lilyminium, p-j-smith, bdice, joaomcteixeira,
Expand Down
29 changes: 20 additions & 9 deletions package/MDAnalysis/analysis/rdf.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,8 @@
`pair distribution functions`_ (`radial distribution functions`_ or "RDF").
The RDF :math:`g_{ab}(r)` between types of particles :math:`a` and :math:`b` is
.. _equation-gab:
.. math::
g_{ab}(r) = (N_{a} N_{b})^{-1} \sum_{i=1}^{N_a} \sum_{j=1}^{N_b}
Expand Down Expand Up @@ -58,9 +60,11 @@
the first minimum in :math:`g(r)`.
In :class:`InterRDF_s`, we provide an option `density`. When `density` is
`True`, it will return the RDF :math:`g_{ab}(r)`; when `density` is False`,
it will return the density of particle :math:`b` in a shell at distance
:math:`r` around a :math:`a` particle, which is
``False``, it will return the RDF :math:`g_{ab}(r)`; when `density` is
``True``, it will return the density of particle :math:`b` in a shell at
distance :math:`r` around a :math:`a` particle, which is
.. _equation-nab:
.. math::
n_{ab}(r) = \rho g_{ab}(r)
Expand Down Expand Up @@ -317,9 +321,16 @@ class InterRDF_s(AnalysisBase):
range : tuple or list (optional)
The size of the RDF [0.0, 15.0]
density : bool (optional)
Calculate :math:`g_{ab}(r)` if set to ``True``; or calculate
:math:`n_{ab}(r)` if set to ``false``; the default is
``True``.
``False``: calculate :math:`g_{ab}(r)`; ``True``: calculate
the true :ref:`single particle density<equation-nab>`
:math:`n_{ab}(r)`. The default is ``False``.
.. versionadded:: 1.0.1
This keyword was available since 0.19.0 but was not
documented. Furthermore, it had the opposite
meaning. Since 1.0.1 it is officially supported as
documented.
Example
Expand Down Expand Up @@ -370,7 +381,7 @@ class InterRDF_s(AnalysisBase):
"""
def __init__(self, u, ags,
nbins=75, range=(0.0, 15.0), density=True, **kwargs):
nbins=75, range=(0.0, 15.0), density=False, **kwargs):
super(InterRDF_s, self).__init__(u.universe.trajectory, **kwargs)

# List of pairs of AtomGroups
Expand Down Expand Up @@ -428,9 +439,9 @@ def _conclude(self):
density = 1 / box_vol

if self._density:
rdf.append(self.count[i] / (density * vol * self.n_frames))
else:
rdf.append(self.count[i] / (vol * self.n_frames))
else:
rdf.append(self.count[i] / (density * vol * self.n_frames))

self.rdf = rdf
self.indices = indices
Expand Down
11 changes: 7 additions & 4 deletions testsuite/MDAnalysisTests/analysis/test_rdf_s.py
Original file line number Diff line number Diff line change
Expand Up @@ -102,14 +102,17 @@ def test_cdf(rdf):


@pytest.mark.parametrize("density, value", [
(True, 26551.55088100731),
(False, 0.021915460340071267)])
(None, 26551.55088100731), # default, like False (no kwarg, see below)
(False, 26551.55088100731),
(True, 0.021915460340071267)])
def test_density(u, sels, density, value):
rdf = InterRDF_s(u, sels, density=density).run()
kwargs = {'density': density} if density is not None else {}
rdf = InterRDF_s(u, sels, **kwargs).run()
assert_almost_equal(max(rdf.rdf[0][0][0]), value)
if density:
if not density:
s1 = u.select_atoms('name ZND and resid 289')
s2 = u.select_atoms(
'name OD1 and resid 51 and sphzone 5.0 (resid 289)')
rdf_ref = InterRDF(s1, s2).run()
assert_almost_equal(rdf_ref.rdf, rdf.rdf[0][0][0])

0 comments on commit dbd39c2

Please sign in to comment.