Expose computed density in InterRDF results

package/AUTHORS

@@ -286,6 +286,7 @@ Chronological list of authors
   - Aryaman Chaudhri
   - Francesco Siciliani
   - Akshit Boora
+  - Sumit Saraswat
 
 
 External code

package/CHANGELOG

@@ -18,7 +18,8 @@ The rules for this file:
          spyke7, talagayev, tanii1125, BradyAJohnston, hejamu, jeremyleung521,
          harshitgajjela-droid, kunjsinha, aygarwal, jauy123, Dreamstick9,
          ollyfutur, Amarendra22, charity-g, ParthUppal523, apoorva-01, RMeli,
-         raulloiscuns, Aryaman-Chaudhri, sici17, AkshitBoora
+         raulloiscuns, Aryaman-Chaudhri, sici17, AkshitBoora,
+         sumitsaraswat362
 
  * 2.11.0
 
@@ -109,6 +110,10 @@ Fixes
    DSSP by porting upstream PyDSSP 0.9.1 fix (Issue #4913)
 
 Enhancements
+ * InterRDF and InterRDF_s now expose computed average box volume and densities
+   as results.box_volume, results.normalization_density, and (for InterRDF)
+   results.density_1 and results.density_2 when norm is "rdf" or "density"
+   (Issue #5362, PR #5482)
  * Add reading single-frame files from mmCIF, also known as PDBx format, using
    crystallographic library `gemmi` as a backend (Issue #2367 and extension of
    #4303, also solves #5089, PR #4712).

package/MDAnalysis/analysis/rdf.py

@@ -217,6 +217,39 @@ class InterRDF(AnalysisBase):
            This attribute will be removed in 3.0.0.
            Use :attr:`results.count` instead.
 
+    results.normalization_density : float
+       The pair normalization density :math:`N_{\text{pairs}} / \langle V \rangle`
+       used to normalize the RDF (or ``None`` if box volume is not available),
+       where :math:`\langle V \rangle` is the average box volume and
+       :math:`N_{\text{pairs}} = N_1 N_2 - N_{\text{excl}}`.
+       This attribute only exists if `norm` is ``"rdf"`` or ``"density"``.
+
+       .. versionadded:: 2.11.0
+
+    results.box_volume : float
+       The average box volume :math:`\langle V \rangle` across analyzed frames
+       (or ``None`` if box volume is not available).
+       This attribute only exists if `norm` is ``"rdf"`` or ``"density"``.
+
+       .. versionadded:: 2.11.0
+
+    results.density_1 : float
+       The bulk number density :math:`\rho_1 = N_1 / \langle V \rangle` of the
+       first selection `g1` (or ``None`` if box volume is not available).
+       This attribute only exists if `norm` is ``"rdf"`` or ``"density"``.
+
+       .. versionadded:: 2.11.0
+
+    results.density_2 : float
+       The bulk number density :math:`\rho_2 = N_2 / \langle V \rangle` of the
+       second selection `g2` (or ``None`` if box volume is not available).
+       This corresponds to the bulk density :math:`\rho` appearing in the
+       equations for :ref:`cumulative coordination counts<equation-countab>`
+       and :ref:`radial densities<equation-nab>`.
+       This attribute only exists if `norm` is ``"rdf"`` or ``"density"``.
+
+       .. versionadded:: 2.11.0
+
     Example
     -------
     First create the :class:`InterRDF` object, by supplying two
@@ -318,7 +351,7 @@ def _prepare(self):
         self.results.edges = edges
         self.results.bins = 0.5 * (edges[:-1] + edges[1:])
 
-        if self.norm == "rdf":
+        if self.norm in ["rdf", "density"]:
             # Cumulative volume for rdf normalization
             self.results.volume_cum = 0
         # Set the max range to filter the search radius
@@ -349,7 +382,7 @@ def _single_frame(self):
         count, _ = np.histogram(dist, **self.rdf_settings)
         self.results.count += count
 
-        if self.norm == "rdf":
+        if self.norm in ["rdf", "density"]:
             self.results.volume_cum += self._ts.volume
 
     def _get_aggregator(self):
@@ -369,7 +402,7 @@ def _conclude(self):
             vols = np.power(self.results.edges, 3)
             norm *= 4 / 3 * np.pi * np.diff(vols)
 
-        if self.norm == "rdf":
+        if self.norm in ["rdf", "density"]:
             # Number of each selection
             nA = len(self.g1)
             nB = len(self.g2)
@@ -381,10 +414,23 @@ def _conclude(self):
                 nblocks = nA / xA
                 N -= xA * xB * nblocks
 
-            # Average number density
+            # Average box volume and densities
             self.volume_cum = self.results.volume_cum
             box_vol = self.volume_cum / self.n_frames
-            norm *= N / box_vol
+
+            if box_vol > 0:
+                self.results.box_volume = box_vol
+                self.results.normalization_density = N / box_vol
+                self.results.density_1 = nA / box_vol
+                self.results.density_2 = nB / box_vol
+            else:
+                self.results.box_volume = None
+                self.results.normalization_density = None
+                self.results.density_1 = None
+                self.results.density_2 = None
+
+            if self.norm == "rdf":
+                norm *= self.results.normalization_density
 
         self.results.rdf = self.results.count / norm
 
@@ -569,6 +615,21 @@ class InterRDF_s(AnalysisBase):
             This attribute will be removed in 3.0.0.
             Use :attr:`results.cdf` instead.
 
+    results.normalization_density : float
+        The normalization density :math:`1 / \langle V \rangle` used for
+        normalizing site-specific pair RDFs (or ``None`` if box volume is not
+        available), where :math:`\langle V \rangle` is the average box volume.
+        This attribute only exists if `norm` is ``"rdf"`` or ``"density"``.
+
+        .. versionadded:: 2.11.0
+
+    results.box_volume : float
+        The average box volume :math:`\langle V \rangle` across analyzed frames
+        (or ``None`` if box volume is not available).
+        This attribute only exists if `norm` is ``"rdf"`` or ``"density"``.
+
+        .. versionadded:: 2.11.0
+
     Example
     -------
     First create the :class:`InterRDF_s` object, by supplying one Universe and
@@ -704,7 +765,7 @@ def _prepare(self):
         self.results.edges = edges
         self.results.bins = 0.5 * (edges[:-1] + edges[1:])
 
-        if self.norm == "rdf":
+        if self.norm in ["rdf", "density"]:
             # Cumulative volume for rdf normalization
             self.results.volume_cum = 0
         self._maxrange = self.rdf_settings["range"][1]
@@ -730,7 +791,7 @@ def _single_frame(self):
             idx2s = pairs[:, 1]
             self.results.count[i][idx1s, idx2s, bin_indices] += counts
 
-        if self.norm == "rdf":
+        if self.norm in ["rdf", "density"]:
             self.results.volume_cum += self._ts.volume
 
     def _conclude(self):
@@ -740,10 +801,19 @@ def _conclude(self):
             vols = np.power(self.results.edges, 3)
             norm *= 4 / 3 * np.pi * np.diff(vols)
 
-        if self.norm == "rdf":
-            # Average number density
+        if self.norm in ["rdf", "density"]:
+            # Average box volume and normalization density
             self.volume_cum = self.results.volume_cum
-            norm *= 1 / (self.volume_cum / self.n_frames)
+            if self.volume_cum > 0:
+                box_vol = self.volume_cum / self.n_frames
+                self.results.box_volume = box_vol
+                self.results.normalization_density = 1.0 / box_vol
+            else:
+                self.results.box_volume = None
+                self.results.normalization_density = None
+
+            if self.norm == "rdf":
+                norm *= self.results.normalization_density
 
         # Empty lists to restore indices, RDF
         self.results.indices = []

testsuite/MDAnalysisTests/analysis/test_rdf.py

@@ -171,3 +171,49 @@ def test_backend(sels, backend):
     s1, s2 = sels
     rdf = InterRDF(s1, s2, norm="none", backend=backend).run()
     assert_allclose(max(rdf.results.rdf), 4)
+
+
+@pytest.mark.parametrize("norm", ("rdf", "density"))
+def test_exposed_density(sels, norm, client_InterRDF):
+    s1, s2 = sels
+    expected_volume = 1_000_000.0
+    rdf = InterRDF(s1, s2, norm=norm).run(**client_InterRDF)
+
+    assert hasattr(rdf.results, "box_volume")
+    assert hasattr(rdf.results, "normalization_density")
+    assert hasattr(rdf.results, "density_1")
+    assert hasattr(rdf.results, "density_2")
+
+    assert_allclose(rdf.results.box_volume, expected_volume)
+    assert_allclose(rdf.results.density_1, len(s1) / expected_volume)
+    assert_allclose(rdf.results.density_2, len(s2) / expected_volume)
+    assert_allclose(
+        rdf.results.normalization_density,
+        (len(s1) * len(s2)) / expected_volume,
+    )
+
+
+def test_exposed_density_norm_none(sels, client_InterRDF):
+    s1, s2 = sels
+    rdf = InterRDF(s1, s2, norm="none").run(**client_InterRDF)
+    assert not hasattr(rdf.results, "box_volume")
+    assert not hasattr(rdf.results, "normalization_density")
+    assert not hasattr(rdf.results, "density_1")
+    assert not hasattr(rdf.results, "density_2")
+
+
+def test_density_zero_box(u, client_InterRDF):
+    trajectory_data = np.tile(u.atoms.positions, (NFRAMES, 1, 1))
+    dimensions_data = np.zeros((NFRAMES, 6))
+    u.load_new(trajectory_data, format="MEMORY", dimensions=dimensions_data)
+    s1 = u.select_atoms("name OW")
+    s2 = u.select_atoms("name HW1 HW2")
+    rdf = InterRDF(s1, s2, norm="density").run(**client_InterRDF)
+    assert hasattr(rdf.results, "box_volume")
+    assert hasattr(rdf.results, "normalization_density")
+    assert hasattr(rdf.results, "density_1")
+    assert hasattr(rdf.results, "density_2")
+    assert rdf.results.box_volume is None
+    assert rdf.results.normalization_density is None
+    assert rdf.results.density_1 is None
+    assert rdf.results.density_2 is None

testsuite/MDAnalysisTests/analysis/test_rdf_s.py

@@ -116,6 +116,13 @@ def test_density(u, sels, density, value, client_InterRDF_s):
     kwargs = {"density": density} if density is not None else {}
     rdf = InterRDF_s(u, sels, **kwargs).run(**client_InterRDF_s)
     assert_almost_equal(max(rdf.results.rdf[0][0][0]), value)
+    assert hasattr(rdf.results, "box_volume")
+    assert hasattr(rdf.results, "normalization_density")
+    expected_box_vol = sum(ts.volume for ts in u.trajectory) / len(
+        u.trajectory
+    )
+    assert_allclose(rdf.results.box_volume, expected_box_vol)
+    assert_allclose(rdf.results.normalization_density, 1.0 / expected_box_vol)
     if not density:
         s1 = u.select_atoms("name ZND and resid 289")
         s2 = u.select_atoms(
@@ -196,6 +203,21 @@ def test_nested_array_sum():
     assert np.array_equal(result[1], arr_2 + arr_4)
 
 
+def test_density_zero_box(client_InterRDF_s):
+    u = mda.Universe(GRO_MEMPROT, in_memory=True)
+    trajectory_data = np.tile(u.atoms.positions, (2, 1, 1))
+    dimensions_data = np.zeros((2, 6))
+    u.load_new(trajectory_data, format="MEMORY", dimensions=dimensions_data)
+    s1 = u.select_atoms("name ZND and resid 289")
+    s2 = u.select_atoms("name OD1 and resid 51")
+    ags = [[s1, s2]]
+    rdf = InterRDF_s(u, ags, norm="density").run(**client_InterRDF_s)
+    assert hasattr(rdf.results, "box_volume")
+    assert hasattr(rdf.results, "normalization_density")
+    assert rdf.results.box_volume is None
+    assert rdf.results.normalization_density is None
+
+
 @pytest.mark.parametrize("n_pairs", [1, 2, 3])
 def test_n_pairs(u, sels, n_pairs, client_InterRDF_s):
     ags = sels[:n_pairs]