Feature/rdf normalization

ChangeLog.md

@@ -22,6 +22,7 @@ and this project adheres to
 * Extensive new documentation including tutorial for new users and reference sections on crucial topics.
 * ClusterProperties computes radius of gyration from the gyration tensor for each cluster.
 * `freud.data` module.
+* Optional normalization for RDF, useful for small systems.
 
 ### Changed
 * All compute objects that perform neighbor computations now use NeighborQuery internally.

cpp/density/RDF.cc

@@ -11,7 +11,7 @@
 
 namespace freud { namespace density {
 
-RDF::RDF(unsigned int bins, float r_max, float r_min) : BondHistogramCompute()
+RDF::RDF(unsigned int bins, float r_max, float r_min, bool normalize) : BondHistogramCompute(), m_normalize(normalize)
 {
     if (bins == 0)
         throw std::invalid_argument("RDF requires a nonzero number of bins.");
@@ -48,9 +48,13 @@ void RDF::reduce()
     m_N_r.prepare(getAxisSizes()[0]);
 
     // Define prefactors with appropriate types to simplify and speed later code.
-    float ndens = float(m_n_query_points) / m_box.getVolume();
+    float number_density = float(m_n_query_points) / m_box.getVolume();
+    if (m_normalize)
+    {
+        number_density *= static_cast<float>(m_n_query_points-1)/(m_n_query_points);
+    }
     float np = static_cast<float>(m_n_points);
-    float prefactor = float(1.0)/(np*ndens*m_frame_counter);
+    float prefactor = float(1.0)/(np*number_density*m_frame_counter);
 
     util::ManagedArray<float> vol_array = m_box.is2D() ? m_vol_array2D : m_vol_array3D;
     m_histogram.reduceOverThreadsPerBin(m_local_histograms,

cpp/density/RDF.h

@@ -17,7 +17,7 @@ class RDF : public locality::BondHistogramCompute
 {
 public:
     //! Constructor
-    RDF(unsigned int bins, float r_max, float r_min = 0);
+    RDF(unsigned int bins, float r_max, float r_min = 0, bool normalize = false);
 
     //! Destructor
     virtual ~RDF() {};
@@ -51,6 +51,7 @@ class RDF : public locality::BondHistogramCompute
     }
 
 private:
+    bool m_normalize;                          //!< Whether to enforce that the RDF should tend to 1 (instead of num_query_points/num_points).
     util::ManagedArray<float> m_pcf;         //!< The computed pair correlation function.
     util::ManagedArray<float> m_N_r;         //!< Cumulative bin sum N(r) (the average number of points in a ball of radius r).
     util::ManagedArray<float> m_vol_array2D; //!< Areas of concentric rings corresponding to the histogram bins in 2D.

doc/source/credits.rst

@@ -45,6 +45,7 @@ Vyas Ramasubramani - **Lead developer**
 * Moved all citations into Bibtex format.
 * Created data module.
 * Standardized PMFT normalization.
+* Enabled optional normalization of RDF.
 
 Bradley Dice - **Lead developer**
 

freud/_density.pxd

@@ -6,6 +6,7 @@ from libcpp.memory cimport shared_ptr
 from libcpp.complex cimport complex
 from libcpp.vector cimport vector
 from freud._locality cimport BondHistogramCompute
+from libcpp cimport bool
 
 cimport freud._box
 cimport freud._locality
@@ -50,7 +51,7 @@ cdef extern from "LocalDensity.h" namespace "freud::density":
 
 cdef extern from "RDF.h" namespace "freud::density":
     cdef cppclass RDF(BondHistogramCompute):
-        RDF(float, float, float) except +
+        RDF(float, float, float, bool) except +
         const freud._box.Box & getBox() const
         void accumulate(const freud._locality.NeighborQuery*,
                         const vec3[float]*,

freud/density.pyx

@@ -421,20 +421,23 @@ cdef class LocalDensity(_PairCompute):
 
 
 cdef class RDF(_SpatialHistogram1D):
-    R"""Computes RDF for supplied data.
+    R"""Computes the RDF :math:`g \left( r \right)` for supplied data.
 
-    The RDF (:math:`g \left( r \right)`) is computed and averaged for a given
-    set of query points in a sea of data points. Providing the same points
-    calculates them against themselves. Computing the RDF results in an RDF
-    array listing the value of the RDF at each given :math:`r`, listed in the
-    :code:`R` array.
+    Note that the RDF is defined strictly according to the pair correlation
+    function, i.e.
+
+     .. math::
 
-    The values of :math:`r` to compute the RDF are set by the values of
-    :code:`r_min`, :code:`r_max`, :code:`bins` in the constructor.
-    :code:`r_max` sets the maximum distance at which to calculate the :math:`g
-    \left( r \right)`, :code:`r_min` sets the minimum distance at which to
-    calculate the :math:`g \left( r \right)`, and :code:`bins` determines the
-    number of bins.
+         g(r) = V\frac{N-1}{N} \langle \delta(r) \rangle
+
+    In the thermodynamic limit, the fraction tends to unity and the limiting
+    behavior of :math:`\lim_{r->\infty} g(r)=1` is recovered. However, for very
+    small systems the long range behavior of the radial distribution will
+    instead tend to :math:`\frac{N-1}{N}`. If you are analyzing a very small
+    system but wish to recover the more familiar behavior, you may use the
+    `normalize` flag to enforce this requirement upon construction of this
+    object. Note that this will have little to no effect on larger systems (for
+    example, for systems of 100 particles the RDF will differ by 1%).
 
     .. note::
         **2D:** :class:`freud.density.RDF` properly handles 2D boxes.
@@ -448,13 +451,27 @@ cdef class RDF(_SpatialHistogram1D):
         r_min (float, optional):
             Minimum interparticle distance to include in the calculation
             (Default value = :code:`0`).
+        normalize (bool, optional):
+            Scale the RDF values by
+            :math:`\frac{N_{query\_points}+1}{N_{query\_points}+1}`. This
+            argument primarily exists to deal with standard RDF calculations
+            where no special ``query_points`` or ``neighbors`` are provided,
+            but where the number of ``query_points`` is small enough that the
+            long-ranged limit of :math:`g(r)` deviates significantly from
+            :math:`1`. It should not be used if :code:`query_points` is
+            provided as a different set of points, or if unusual query
+            arguments are provided to :meth:`~.compute`, specifically if
+            :code`exclude_ii` is set to :code:`False`. This normalization is
+            not meaningful in such cases and will simply convolute the data.
+
     """
     cdef freud._density.RDF * thisptr
 
-    def __cinit__(self, unsigned int bins, float r_max, float r_min=0):
+    def __cinit__(self, unsigned int bins, float r_max, float r_min=0,
+                  normalize=False):
         if type(self) == RDF:
             self.thisptr = self.histptr = new freud._density.RDF(
-                bins, r_max, r_min)
+                bins, r_max, r_min, normalize)
 
             # r_max is left as an attribute rather than a property for now
             # since that change needs to happen at the _SpatialHistogram level
@@ -545,7 +562,7 @@ cdef class RDF(_SpatialHistogram1D):
             (:class:`matplotlib.axes.Axes`): Axis with the plot.
         """
         import freud.plot
-        return freud.plot.line_plot(self.R, self.RDF,
+        return freud.plot.line_plot(self.bin_centers, self.rdf,
                                     title="RDF",
                                     xlabel=r"$r$",
                                     ylabel=r"$g(r)$",

tests/test_density_RDF.py

@@ -114,6 +114,7 @@ def test_repr_png(self):
         self.assertEqual(rdf._repr_png_(), None)
 
         rdf.compute((box, points), reset=False)
+        rdf.plot()
         rdf._repr_png_()
 
     def test_points_ne_query_points(self):