From dbd39c2900d8b84184d61ff1adfd5250616994e8 Mon Sep 17 00:00:00 2001
From: Oliver Beckstein <orbeckst@gmail.com>
Date: Fri, 10 Jul 2020 14:28:03 -0700
Subject: [PATCH] changed InterRDF_s density keyword semantics

- 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
---
 package/CHANGELOG                             |  6 ++--
 package/MDAnalysis/analysis/rdf.py            | 29 +++++++++++++------
 .../MDAnalysisTests/analysis/test_rdf_s.py    | 11 ++++---
 3 files changed, 31 insertions(+), 15 deletions(-)

diff --git a/package/CHANGELOG b/package/CHANGELOG
index d0d568b052f..4e7ac1d120d 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -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,
diff --git a/package/MDAnalysis/analysis/rdf.py b/package/MDAnalysis/analysis/rdf.py
index 41f7d1b5534..b78326a29f8 100644
--- a/package/MDAnalysis/analysis/rdf.py
+++ b/package/MDAnalysis/analysis/rdf.py
@@ -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}
@@ -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)
@@ -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
@@ -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
@@ -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
diff --git a/testsuite/MDAnalysisTests/analysis/test_rdf_s.py b/testsuite/MDAnalysisTests/analysis/test_rdf_s.py
index 874f414a698..2f0babab54b 100644
--- a/testsuite/MDAnalysisTests/analysis/test_rdf_s.py
+++ b/testsuite/MDAnalysisTests/analysis/test_rdf_s.py
@@ -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])
+