diff --git a/package/AUTHORS b/package/AUTHORS index 3d87a0651c..ed19d3ac17 100644 --- a/package/AUTHORS +++ b/package/AUTHORS @@ -286,6 +286,7 @@ Chronological list of authors - Aryaman Chaudhri - Francesco Siciliani - Akshit Boora + - Sumit Saraswat External code diff --git a/package/CHANGELOG b/package/CHANGELOG index 72b880c821..e70fe474f9 100644 --- a/package/CHANGELOG +++ b/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). diff --git a/package/MDAnalysis/analysis/rdf.py b/package/MDAnalysis/analysis/rdf.py index 6a33abc840..902dbb7100 100644 --- a/package/MDAnalysis/analysis/rdf.py +++ b/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.12.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.12.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.12.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` + and :ref:`radial densities`. + This attribute only exists if `norm` is ``"rdf"`` or ``"density"``. + + .. versionadded:: 2.12.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.12.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.12.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 = [] diff --git a/testsuite/MDAnalysisTests/analysis/test_rdf.py b/testsuite/MDAnalysisTests/analysis/test_rdf.py index 36da9be327..dfc84b1377 100644 --- a/testsuite/MDAnalysisTests/analysis/test_rdf.py +++ b/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 diff --git a/testsuite/MDAnalysisTests/analysis/test_rdf_s.py b/testsuite/MDAnalysisTests/analysis/test_rdf_s.py index 8e23410404..441fe4b87e 100644 --- a/testsuite/MDAnalysisTests/analysis/test_rdf_s.py +++ b/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]