Merge pull request #1075 from glotzerlab/feature/rad

RAD Neighborlist Filter
glotzerlab · Apr 21, 2023 · bde386f · bde386f
2 parents fcdd0bc + f54897b
commit bde386f
Show file tree

Hide file tree

Showing 14 changed files with 528 additions and 111 deletions.
diff --git a/ChangeLog.md b/ChangeLog.md
@@ -7,7 +7,7 @@ and this project adheres to
 ## v2.13.0 -- YYYY-MM-DD
 
 ### Added
-* Filter neighborlists with `freud.locality.FilterSANN`.
+* Filter neighborlists with `freud.locality.FilterSANN` and `freud.locality.FilterRAD`.
 
 ## v2.12.1 -- 2022-12-05
 

diff --git a/cpp/locality/CMakeLists.txt b/cpp/locality/CMakeLists.txt
@@ -11,6 +11,8 @@ add_library(
   Filter.h
   FilterSANN.cc
   FilterSANN.h
+  FilterRAD.cc
+  FilterRAD.h
   LinkCell.cc
   LinkCell.h
   NeighborBond.h

diff --git a/cpp/locality/Filter.h b/cpp/locality/Filter.h
@@ -6,6 +6,7 @@
 
 #include "NeighborList.h"
 #include "NeighborQuery.h"
+#include <iostream>
 
 namespace freud { namespace locality {
 
@@ -22,9 +23,9 @@ namespace freud { namespace locality {
 class Filter
 {
 public:
-    Filter()
+    Filter(bool allow_incomplete_shell)
         : m_unfiltered_nlist(std::make_shared<NeighborList>()),
-          m_filtered_nlist(std::make_shared<NeighborList>())
+          m_filtered_nlist(std::make_shared<NeighborList>()), m_allow_incomplete_shell(allow_incomplete_shell)
     {}
 
     virtual ~Filter() = default;
@@ -46,8 +47,51 @@ class Filter
 protected:
     //!< The unfiltered neighborlist
     std::shared_ptr<NeighborList> m_unfiltered_nlist;
+
     //!< The filtered neighborlist
     std::shared_ptr<NeighborList> m_filtered_nlist;
+
+    //<! whether a warning (true) or error (false) should be raised if the filter
+    //<! algorithm implementation cannot guarantee that all neighbors have completely filled shells
+    bool m_allow_incomplete_shell;
+
+    /*! Output the appropriate warning/error message for particles with unfilled shells
+     *
+     * In general, the filter concept cannot guarantee that each query point will have
+     * a completely filled shell according to the implemented algorithm. This happens
+     * when the initial unfiltered neighbor list doesn't have enough neighbors to
+     * guarantee this criterion.
+     *
+     * \param unfilled_qps Vector of query points which may have unfilled neighbor shells.
+     *                     The vector should have the value
+     *                     ``unfilled_qps[i] = std::numeric_limits<unsigned int>::max()``
+     *                     for all query point indices ``i`` which have filled shells, and
+     *                     should have ``unfilled_qps[i] = i`` for all query point indices
+     *                     ``i`` which may not have a filled neighbor shell.
+     * */
+    void warnAboutUnfilledNeighborShells(const std::vector<unsigned int>& unfilled_qps) const
+    {
+        std::string indices;
+        for (const auto& idx : unfilled_qps)
+        {
+            if (idx != std::numeric_limits<unsigned int>::max())
+            {
+                indices += std::to_string(idx);
+                indices += ", ";
+            }
+        }
+        indices = indices.substr(0, indices.size() - 2);
+        if (!indices.empty())
+        {
+            std::ostringstream error_str;
+            error_str << "Query point indices " << indices << " do not have full neighbor shells.";
+            if (!m_allow_incomplete_shell)
+            {
+                throw std::runtime_error(error_str.str());
+            }
+            std::cout << "WARNING: " << error_str.str() << std::endl;
+        }
+    }
 };
 
 }; };  // namespace freud::locality

diff --git a/cpp/locality/FilterRAD.cc b/cpp/locality/FilterRAD.cc
@@ -0,0 +1,110 @@
+// Copyright (c) 2010-2023 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
+#include "FilterRAD.h"
+#include "NeighborBond.h"
+#include "NeighborComputeFunctional.h"
+#include "utils.h"
+#include <tbb/enumerable_thread_specific.h>
+#include <vector>
+
+namespace freud { namespace locality {
+
+void FilterRAD::compute(const NeighborQuery* nq, const vec3<float>* query_points,
+                        unsigned int num_query_points, const NeighborList* nlist, const QueryArgs& qargs)
+{
+    // make the unfiltered neighborlist from the arguments
+    m_unfiltered_nlist = std::make_shared<NeighborList>(
+        std::move(makeDefaultNlist(nq, nlist, query_points, num_query_points, qargs)));
+
+    // work with nlist sorted by distance
+    NeighborList sorted_nlist(*m_unfiltered_nlist);
+    sorted_nlist.sort(true);
+
+    // hold index of query point for a thread if its RAD shell isn't filled
+    std::vector<unsigned int> unfilled_qps(sorted_nlist.getNumQueryPoints(),
+                                           std::numeric_limits<unsigned int>::max());
+
+    const auto& points = nq->getPoints();
+    const auto& box = nq->getBox();
+    const auto& sorted_neighbors = sorted_nlist.getNeighbors();
+    const auto& sorted_dist = sorted_nlist.getDistances();
+    const auto& sorted_weights = sorted_nlist.getWeights();
+    const auto& sorted_counts = sorted_nlist.getCounts();
+
+    using BondVector = tbb::enumerable_thread_specific<std::vector<NeighborBond>>;
+    BondVector filtered_bonds;
+
+    // parallelize over query_point_index
+    util::forLoopWrapper(0, sorted_nlist.getNumQueryPoints(), [&](size_t begin, size_t end) {
+        // grab thread-local vector
+        BondVector::reference local_bonds(filtered_bonds.local());
+        for (auto i = begin; i < end; i++)
+        {
+            const auto num_unfiltered_neighbors = sorted_counts(i);
+            const auto first_idx = sorted_nlist.find_first_index(i);
+            bool good_neighbor = true;
+
+            // loop over each potential neighbor particle j
+            for (unsigned int j = 0; j < num_unfiltered_neighbors; j++)
+            {
+                const auto first_neighbor_idx = sorted_neighbors(first_idx + j, 1);
+                const auto v1 = box.wrap(query_points[i] - points[first_neighbor_idx]);
+                good_neighbor = true;
+
+                // loop over particles which may be blocking the neighbor j
+                for (unsigned int k = 0; k < j; k++)
+                {
+                    const auto second_neighbor_idx = sorted_neighbors(first_idx + k, 1);
+                    const auto v2 = box.wrap(query_points[i] - points[second_neighbor_idx]);
+
+                    // check if k blocks j
+                    if ((dot(v2, v2) * sorted_dist(first_idx + j) * sorted_dist(first_idx + k))
+                        < (dot(v1, v2) * dot(v1, v1)))
+                    {
+                        good_neighbor = false;
+                        break;
+                    }
+                }
+
+                // if no k blocks j, add a bond from i to j
+                if (good_neighbor)
+                {
+                    local_bonds.emplace_back(i, first_neighbor_idx, sorted_dist(first_idx + j));
+                }
+                else if (m_terminate_after_blocked)
+                {
+                    // if a particle blocks j and "RAD-closed" is requested
+                    // stop looking for more neighbors
+                    break;
+                }
+            }
+
+            // if we have searched over all potential neighbors j and still have
+            // not found one that is blocked, the neighbor shell may be incomplete.
+            // This only applies to the RAD-closed case, because RAD-open will
+            // never terminate prematurely.
+            if (good_neighbor && m_terminate_after_blocked)
+            {
+                // in principle, the incomplete shell exception can be raised here,
+                // but the error is more informative if the exception raised
+                // includes each query point with an unfilled neighbor shell
+                unfilled_qps[i] = i;
+            }
+        }
+    });
+
+    // print warning/exception about query point indices with unfilled neighbor shells
+    Filter::warnAboutUnfilledNeighborShells(unfilled_qps);
+
+    // combine thread-local arrays
+    tbb::flattened2d<BondVector> flat_filtered_bonds = tbb::flatten2d(filtered_bonds);
+    std::vector<NeighborBond> rad_bonds(flat_filtered_bonds.begin(), flat_filtered_bonds.end());
+
+    // sort final bonds array by distance
+    tbb::parallel_sort(rad_bonds.begin(), rad_bonds.end(), compareNeighborDistance);
+
+    m_filtered_nlist = std::make_shared<NeighborList>(rad_bonds);
+};
+
+}; }; // namespace freud::locality
diff --git a/cpp/locality/FilterRAD.h b/cpp/locality/FilterRAD.h
@@ -0,0 +1,34 @@
+// Copyright (c) 2010-2023 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
+#ifndef __FILTERRAD_H__
+#define __FILTERRAD_H__
+
+#include "Filter.h"
+
+#include "NeighborList.h"
+#include "NeighborQuery.h"
+
+namespace freud { namespace locality {
+
+/* Class for the RAD method of filtering a neighborlist.
+ * */
+class FilterRAD : public Filter
+{
+public:
+    //<! Construct with an empty NeighborList, fill it upon calling compute
+    FilterRAD(bool allow_incomplete_shell, bool terminate_after_blocked)
+        : Filter(allow_incomplete_shell), m_terminate_after_blocked(terminate_after_blocked)
+    {}
+
+    void compute(const NeighborQuery* nq, const vec3<float>* query_points, unsigned int num_query_points,
+                 const NeighborList* nlist, const QueryArgs& qargs) override;
+
+private:
+    //<! variable that determines if RAD-open (False) or RAD-closed (True) is computed
+    bool m_terminate_after_blocked;
+};
+
+}; };  // namespace freud::locality
+
+#endif // __FILTERRAD_H__
diff --git a/cpp/locality/FilterSANN.cc b/cpp/locality/FilterSANN.cc
@@ -5,7 +5,6 @@
 #include "NeighborBond.h"
 #include "NeighborComputeFunctional.h"
 #include "utils.h"
-#include <iostream>
 #include <tbb/enumerable_thread_specific.h>
 #include <vector>
 
@@ -78,7 +77,7 @@ void FilterSANN::compute(const NeighborQuery* nq, const vec3<float>* query_point
     });
 
     // print warning/exception about query point indices with unfilled neighbor shells
-    warnAboutUnfilledNeighborShells(unfilled_qps);
+    Filter::warnAboutUnfilledNeighborShells(unfilled_qps);
 
     // combine thread-local NeighborBond vectors into a single vector
     tbb::flattened2d<BondVector> flat_filtered_bonds = tbb::flatten2d(filtered_bonds);
@@ -90,28 +89,4 @@ void FilterSANN::compute(const NeighborQuery* nq, const vec3<float>* query_point
     m_filtered_nlist = std::make_shared<NeighborList>(sann_bonds);
 };
 
-void FilterSANN::warnAboutUnfilledNeighborShells(const std::vector<unsigned int>& unfilled_qps) const
-{
-    std::string indices;
-    for (const auto& idx : unfilled_qps)
-    {
-        if (idx != std::numeric_limits<unsigned int>::max())
-        {
-            indices += std::to_string(idx);
-            indices += ", ";
-        }
-    }
-    indices = indices.substr(0, indices.size() - 2);
-    if (!indices.empty())
-    {
-        std::ostringstream error_str;
-        error_str << "Query point indices " << indices << " do not have full neighbor shells.";
-        if (!m_allow_incomplete_shell)
-        {
-            throw std::runtime_error(error_str.str());
-        }
-        std::cout << "WARNING: " << error_str.str() << std::endl;
-    }
-}
-
 }; }; // namespace freud::locality
diff --git a/cpp/locality/FilterSANN.h b/cpp/locality/FilterSANN.h
@@ -21,19 +21,14 @@ class FilterSANN : public Filter
      * \param allow_incomplete_shell whether incomplete neighbor shells should
      *                               result in a warning or error
      * */
-    explicit FilterSANN(bool allow_incomplete_shell)
-        : Filter(), m_allow_incomplete_shell(allow_incomplete_shell)
-    {}
+    explicit FilterSANN(bool allow_incomplete_shell) : Filter(allow_incomplete_shell) {}
 
     void compute(const NeighborQuery* nq, const vec3<float>* query_points, unsigned int num_query_points,
                  const NeighborList* nlist, const QueryArgs& qargs) override;
 
 private:
     //! warn/raise exception about unfilled shells
     void warnAboutUnfilledNeighborShells(const std::vector<unsigned int>& unfilled_qps) const;
-
-    //! whether incomplete neighbor shell should result in a warning or error
-    bool m_allow_incomplete_shell;
 };
 
 }; };  // namespace freud::locality

diff --git a/doc/source/modules/locality.rst b/doc/source/modules/locality.rst
@@ -9,6 +9,7 @@ Locality Module
 
     freud.locality.AABBQuery
     freud.locality.Filter
+    freud.locality.FilterRAD
     freud.locality.FilterSANN
     freud.locality.LinkCell
     freud.locality.NeighborList

diff --git a/doc/source/reference/credits.rst b/doc/source/reference/credits.rst
@@ -334,6 +334,7 @@ Tommy Waltmann
 * Add support for compilation with the C++17 standard.
 * Normalize n(r) in ``RDF`` class by number of query points.
 * Contributed code, design, documentation, and testing for ``freud.locality.FilterSANN`` class.
+* Contributed code, design, documentation, and testing for ``freud.locality.FilterRAD`` class.
 
 Maya Martirossyan
 
@@ -346,12 +347,14 @@ Pavel Buslaev
 Charlotte Zhao
 
 * Worked with Vyas Ramasubramani and Bradley Dice to add the ``out`` option for ``box.Box.wrap``.
+* Contributed code, design, documentation, and testing for ``freud.locality.FilterRAD`` class.
 
 Domagoj Fijan
 
 * Contributed code, design, documentation, and testing for ``StaticStructureFactorDebye`` class.
 * Contributed code, design, documentation, and testing for ``StaticStructureFactorDirect`` class.
 * Contributed code, design, documentation, and testing for ``freud.locality.FilterSANN`` class.
+* Contributed code, design, documentation, and testing for ``freud.locality.FilterRAD`` class.
 
 Andrew Kerr
 
@@ -373,6 +376,10 @@ Alain Kadar
 
 * Introduced mass dependence for ``ClusterProperties`` class, inertia tensor and radius of gyration.
 
+Philipp Schönhöfer
+
+* Contributed code, design, documentation, and testing for ``freud.locality.FilterRAD`` class.
+
 Source code
 -----------
 

diff --git a/doc/source/reference/freud.bib b/doc/source/reference/freud.bib
@@ -295,3 +295,15 @@ @article{vanMeel2012
 doi = {10.1063/1.4729313},
 url = {https://doi.org/10.1063/1.4729313},
 }
+
+@article{Higham2016,
+author = {Higham,Jonathan and Henchman,Richard H. },
+title = {Locally adaptive method to define coordination shell},
+journal = {The Journal of Chemical Physics},
+volume = {145},
+number = {8},
+pages = {084108},
+year = {2016},
+doi = {10.1063/1.4961439},
+url = {https://doi.org/10.1063/1.4961439},
+}
diff --git a/freud/_locality.pxd b/freud/_locality.pxd
@@ -161,3 +161,7 @@ cdef extern from "Filter.h" namespace "freud::locality":
 cdef extern from "FilterSANN.h" namespace "freud::locality":
     cdef cppclass FilterSANN(Filter):
         FilterSANN(bool)
+
+cdef extern from "FilterRAD.h" namespace "freud::locality":
+    cdef cppclass FilterRAD(Filter):
+        FilterRAD(bool, bool)
diff --git a/freud/locality.pxd b/freud/locality.pxd
@@ -79,3 +79,6 @@ cdef class Filter(_PairCompute):
 
 cdef class FilterSANN(Filter):
     cdef freud._locality.FilterSANN *_thisptr
+
+cdef class FilterRAD(Filter):
+    cdef freud._locality.FilterRAD *_thisptr