From 73d33f4df6cf0d962e935130f90bcc0a36b56c4d Mon Sep 17 00:00:00 2001 From: Vyas Ramasubramani Date: Thu, 17 Oct 2019 16:59:04 -0400 Subject: [PATCH 1/5] Fix repr_png for RDF --- freud/density.pyx | 2 +- tests/test_density_RDF.py | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/freud/density.pyx b/freud/density.pyx index b548e3849..8cba5727d 100644 --- a/freud/density.pyx +++ b/freud/density.pyx @@ -545,7 +545,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)$", diff --git a/tests/test_density_RDF.py b/tests/test_density_RDF.py index 4dec29b6e..19363f97b 100644 --- a/tests/test_density_RDF.py +++ b/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): From bd46c732495af4a50a95e7920d17a0771c2d446f Mon Sep 17 00:00:00 2001 From: Vyas Ramasubramani Date: Thu, 17 Oct 2019 17:12:58 -0400 Subject: [PATCH 2/5] Add optional constructor flag for normalization --- cpp/density/RDF.cc | 6 +++--- cpp/density/RDF.h | 3 ++- freud/_density.pxd | 2 +- freud/density.pyx | 32 ++++++++++++++++++-------------- 4 files changed, 24 insertions(+), 19 deletions(-) diff --git a/cpp/density/RDF.cc b/cpp/density/RDF.cc index 9f29cf79d..ad6046cd8 100644 --- a/cpp/density/RDF.cc +++ b/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,9 @@ 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 = (m_normalize ? float(m_n_points) : float(m_n_query_points)) / m_box.getVolume(); float np = static_cast(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 vol_array = m_box.is2D() ? m_vol_array2D : m_vol_array3D; m_histogram.reduceOverThreadsPerBin(m_local_histograms, diff --git a/cpp/density/RDF.h b/cpp/density/RDF.h index 3ebc5d7fd..64a8b7b52 100644 --- a/cpp/density/RDF.h +++ b/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 m_pcf; //!< The computed pair correlation function. util::ManagedArray m_N_r; //!< Cumulative bin sum N(r) (the average number of points in a ball of radius r). util::ManagedArray m_vol_array2D; //!< Areas of concentric rings corresponding to the histogram bins in 2D. diff --git a/freud/_density.pxd b/freud/_density.pxd index 8222ace3d..57a6230d5 100644 --- a/freud/_density.pxd +++ b/freud/_density.pxd @@ -50,7 +50,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]*, diff --git a/freud/density.pyx b/freud/density.pyx index 8cba5727d..024096c88 100644 --- a/freud/density.pyx +++ b/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:: + + g(r) = V\frac{N-1}{N} \langle \delta(r) \rangle - 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. + 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. @@ -451,10 +454,11 @@ cdef class RDF(_SpatialHistogram1D): """ 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, bool + 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 From 47be969975a0f90d4eae63ef72f9c8ce50a41946 Mon Sep 17 00:00:00 2001 From: Vyas Ramasubramani Date: Thu, 17 Oct 2019 17:51:22 -0400 Subject: [PATCH 3/5] Fix normalization logic --- cpp/density/RDF.cc | 6 +++++- freud/density.pyx | 14 ++++++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/cpp/density/RDF.cc b/cpp/density/RDF.cc index ad6046cd8..02cc676ab 100644 --- a/cpp/density/RDF.cc +++ b/cpp/density/RDF.cc @@ -48,7 +48,11 @@ void RDF::reduce() m_N_r.prepare(getAxisSizes()[0]); // Define prefactors with appropriate types to simplify and speed later code. - float number_density = (m_normalize ? float(m_n_points) : 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(m_n_query_points-1)/(m_n_query_points); + } float np = static_cast(m_n_points); float prefactor = float(1.0)/(np*number_density*m_frame_counter); diff --git a/freud/density.pyx b/freud/density.pyx index 024096c88..ad48b969d 100644 --- a/freud/density.pyx +++ b/freud/density.pyx @@ -15,6 +15,7 @@ from cython.operator cimport dereference from freud.util cimport _Compute from freud.locality cimport _PairCompute, _SpatialHistogram1D from freud.util cimport vec3 +from libcpp cimport bool from collections.abc import Sequence @@ -451,6 +452,19 @@ 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 From 9c2c6c8b1fc1a079357ab3b443b91c1a360c5077 Mon Sep 17 00:00:00 2001 From: Vyas Ramasubramani Date: Thu, 17 Oct 2019 18:15:40 -0400 Subject: [PATCH 4/5] Fix boolean passing to C++ --- freud/_density.pxd | 1 + freud/density.pyx | 3 +-- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/freud/_density.pxd b/freud/_density.pxd index 57a6230d5..56d7fcf3a 100644 --- a/freud/_density.pxd +++ b/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 diff --git a/freud/density.pyx b/freud/density.pyx index ad48b969d..0950a33e7 100644 --- a/freud/density.pyx +++ b/freud/density.pyx @@ -15,7 +15,6 @@ from cython.operator cimport dereference from freud.util cimport _Compute from freud.locality cimport _PairCompute, _SpatialHistogram1D from freud.util cimport vec3 -from libcpp cimport bool from collections.abc import Sequence @@ -468,7 +467,7 @@ cdef class RDF(_SpatialHistogram1D): """ cdef freud._density.RDF * thisptr - def __cinit__(self, unsigned int bins, float r_max, float r_min=0, bool + 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( From a7c25f43b2f234f3adc4aa7e762dd6de2de1c7c0 Mon Sep 17 00:00:00 2001 From: Vyas Ramasubramani Date: Thu, 17 Oct 2019 18:20:53 -0400 Subject: [PATCH 5/5] Update credits and changelog --- ChangeLog.md | 1 + doc/source/credits.rst | 1 + 2 files changed, 2 insertions(+) diff --git a/ChangeLog.md b/ChangeLog.md index 331807f97..ea4ff3cab 100644 --- a/ChangeLog.md +++ b/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. diff --git a/doc/source/credits.rst b/doc/source/credits.rst index 5b3a4061b..2cc55b710 100644 --- a/doc/source/credits.rst +++ b/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**