Merge pull request #572 from ISISNeutronMuon/chi/van-hove-normalization

Update the normalization of the van hove functions

MDANSE/Extensions/van_hove.pyx

@@ -116,6 +116,7 @@ def van_hove_distinct(
 def van_hove_self(
     double[:,:] xyz,
     double[:,:] histograms,
+    double[:] cell_vols,
     double rmin,
     double dr,
     int n_configs,
@@ -131,6 +132,8 @@ def van_hove_self(
         The trajectory of an atom.
     histograms : np.ndarray
         The histograms to be updated.
+    cell_vols : np.ndarray
+        The cell volumes.
     rmin : float
         The minimum distance of the histogram.
     dr : float
@@ -159,4 +162,4 @@ def van_hove_self(
             if ((bin < 0) or (bin >= nbins)):
                 continue
 
-            histograms[bin, j] += 1.0
+            histograms[bin, j] += cell_vols[i+j]

MDANSE/Src/MDANSE/Framework/Jobs/VanHoveFunctionSelf.py

@@ -105,6 +105,13 @@ def initialize(self):
             axis="r|time",
             units="au",
         )
+        self._outputData.add(
+            "4_pi_r2_g(r,t)_total",
+            "SurfaceOutputVariable",
+            (self.n_mid_points, self.n_frames),
+            axis="r|time",
+            units="au",
+        )
         for element in self.selectedElements:
             self._outputData.add(
                 "g(r,t)_%s" % element,
@@ -113,6 +120,29 @@ def initialize(self):
                 axis="r|time",
                 units="au",
             )
+            self._outputData.add(
+                "4_pi_r2_g(r,t)_%s" % element,
+                "SurfaceOutputVariable",
+                (self.n_mid_points, self.n_frames),
+                axis="r|time",
+                units="au",
+            )
+
+        # usually the normalization is 4 * pi * r^2 * dr which is
+        # correct for small values of dr or large values of r.
+        # unlike the PDF, g(r, t) may not be zero around r=0 we will use
+        # the actual shell volume instead.
+        self.shell_volumes = []
+        for i in range(self.n_mid_points):
+            self.shell_volumes.append(
+                (
+                    self.configuration["r_values"]["value"][i]
+                    + self.configuration["r_values"]["step"]
+                )
+                ** 3
+                - self.configuration["r_values"]["value"][i] ** 3
+            )
+        self.shell_volumes = (4 / 3) * np.pi * np.array(self.shell_volumes)
 
     def run_step(self, atm_index: int) -> tuple[int, tuple[np.ndarray, np.ndarray]]:
         """Calculates a distance histograms of an atoms displacement.
@@ -131,18 +161,30 @@ def run_step(self, atm_index: int) -> tuple[int, tuple[np.ndarray, np.ndarray]]:
             A tuple containing the atom index and distance histograms.
         """
         histograms = np.zeros((self.n_mid_points, self.n_frames))
+        first = self.configuration["frames"]["first"]
+        last = self.configuration["frames"]["last"] + 1
+        step = self.configuration["frames"]["step"]
 
         idx = self.configuration["atom_selection"]["indexes"][atm_index][0]
         series = self.configuration["trajectory"]["instance"].read_atomic_trajectory(
             idx,
-            first=self.configuration["frames"]["first"],
-            last=self.configuration["frames"]["last"] + 1,
-            step=self.configuration["frames"]["step"],
+            first=first,
+            last=last,
+            step=step,
+        )
+        cell_vols = np.array(
+            [
+                self.configuration["trajectory"]["instance"]
+                .configuration(i)
+                .unit_cell.volume
+                for i in range(first, last, step)
+            ]
         )
 
         van_hove.van_hove_self(
             series,
             histograms,
+            cell_vols,
             self.configuration["r_values"]["first"],
             self.configuration["r_values"]["step"],
             self.n_configs,
@@ -165,6 +207,7 @@ def combine(self, atm_index: int, histogram: np.ndarray):
         """
         element = self.configuration["atom_selection"]["names"][atm_index]
         self._outputData["g(r,t)_{}".format(element)][:] += histogram
+        self._outputData["4_pi_r2_g(r,t)_{}".format(element)][:] += histogram
 
     def finalize(self):
         """Using the distance histograms calculate, normalize and save the
@@ -173,7 +216,12 @@ def finalize(self):
 
         nAtomsPerElement = self.configuration["atom_selection"].get_natoms()
         for element, number in list(nAtomsPerElement.items()):
-            self._outputData["g(r,t)_%s" % element][:] /= number * self.n_configs
+            self._outputData["g(r,t)_%s" % element][:] /= (
+                self.shell_volumes[:, np.newaxis] * number**2 * self.n_configs
+            )
+            self._outputData["4_pi_r2_g(r,t)_%s" % element][:] /= (
+                number**2 * self.n_configs * self.configuration["r_values"]["step"]
+            )
 
         weights = self.configuration["weights"].get_weights()
         self._outputData["g(r,t)_total"][:] = weight(
@@ -184,6 +232,14 @@ def finalize(self):
             "g(r,t)_%s",
             update_partials=True,
         )
+        self._outputData["4_pi_r2_g(r,t)_total"][:] = weight(
+            weights,
+            self._outputData,
+            nAtomsPerElement,
+            1,
+            "4_pi_r2_g(r,t)_%s",
+            update_partials=True,
+        )
 
         self._outputData.write(
             self.configuration["output_files"]["root"],