Static structure factor S(q)

.circleci/ci-oldest-reqs.txt

@@ -1,3 +1,4 @@
+ase==3.22.0
 cmake==3.6.3.post1
 cython==0.29.14
 garnett==0.7.1

cpp/CMakeLists.txt

@@ -8,6 +8,7 @@ endif()
 
 add_subdirectory(cluster)
 add_subdirectory(density)
+add_subdirectory(diffraction)
 add_subdirectory(environment)
 add_subdirectory(locality)
 add_subdirectory(order)
@@ -19,6 +20,7 @@ add_library(
   libfreud SHARED
   $<TARGET_OBJECTS:_cluster>
   $<TARGET_OBJECTS:_density>
+  $<TARGET_OBJECTS:_diffraction>
   $<TARGET_OBJECTS:_environment>
   $<TARGET_OBJECTS:_locality>
   $<TARGET_OBJECTS:_order>

cpp/diffraction/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_library(_diffraction OBJECT StaticStructureFactorDebye.h
+                                StaticStructureFactorDebye.cc)
+
+target_include_directories(_diffraction
+                           PUBLIC ${PROJECT_SOURCE_DIR}/cpp/density)

cpp/diffraction/StaticStructureFactorDebye.cc

@@ -0,0 +1,104 @@
+// Copyright (c) 2010-2020 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+
+#include "Box.h"
+#include "ManagedArray.h"
+#include "NeighborComputeFunctional.h"
+#include "NeighborQuery.h"
+#include "RDF.h"
+#include "StaticStructureFactorDebye.h"
+#include "utils.h"
+
+/*! \file StaticStructureFactorDebye.cc
+    \brief Routines for computing static structure factors.
+*/
+
+namespace freud { namespace diffraction {
+
+StaticStructureFactorDebye::StaticStructureFactorDebye(unsigned int bins, float k_max, float k_min)
+{
+    if (bins == 0)
+    {
+        throw std::invalid_argument("StaticStructureFactorDebye requires a nonzero number of bins.");
+    }
+    if (k_max <= 0)
+    {
+        throw std::invalid_argument("StaticStructureFactorDebye requires k_max to be positive.");
+    }
+    if (k_max <= k_min)
+    {
+        throw std::invalid_argument("StaticStructureFactorDebye requires that k_max must be greater than k_min.");
+    }
+
+    // Construct the Histogram object that will be used to track the structure factor
+    auto axes
+        = StaticStructureFactorDebyeHistogram::Axes {std::make_shared<util::RegularAxis>(bins, k_min, k_max)};
+    m_histogram = StaticStructureFactorDebyeHistogram(axes);
+    m_local_histograms = StaticStructureFactorDebyeHistogram::ThreadLocalHistogram(m_histogram);
+    m_min_valid_k = std::numeric_limits<float>::infinity();
+    m_structure_factor.prepare(bins);
+}
+
+void StaticStructureFactorDebye::accumulate(const freud::locality::NeighborQuery* neighbor_query,
+                                             const vec3<float>* query_points, unsigned int n_query_points, unsigned int n_total)
+{
+    auto const& box = neighbor_query->getBox();
+
+    // The r_max should be just less than half of the smallest side length of the box
+    auto const box_L = box.getL();
+    auto const min_box_length
+        = box.is2D() ? std::min(box_L.x, box_L.y) : std::min(box_L.x, std::min(box_L.y, box_L.z));
+    auto const r_max = std::nextafter(float(0.5) * min_box_length, float(0));
+    auto const points= neighbor_query->getPoints();
+    auto const n_points= neighbor_query->getNPoints();
+
+    // The minimum k value of validity is 4 * pi / L, where L is the smallest side length.
+    // This is equal to 2 * pi / r_max.
+    m_min_valid_k = std::min(m_min_valid_k, freud::constants::TWO_PI / r_max);
+
+    auto distances = std::vector<float>(n_points * n_query_points);
+    box.computeAllDistances(points, n_points, query_points, n_query_points, distances.data());
+
+    auto const k_bin_centers = m_histogram.getBinCenters()[0];
+
+    util::forLoopWrapper(0, m_histogram.getAxisSizes()[0], [&](size_t begin_k_index, size_t end_k_index) {
+        for (size_t k_index = begin_k_index; k_index < end_k_index; ++k_index)
+        {
+            auto const k = k_bin_centers[k_index];
+            double S_k = 0.0;
+            for (size_t distance_index = 0; distance_index < distances.size(); ++distance_index)
+            {
+                auto const distance = distances[distance_index];
+                S_k += util::sinc(k * distance);
+            }
+            S_k /= static_cast<double>(n_total);
+            m_local_histograms.increment(k_index, S_k);
+        };
+    });
+    m_frame_counter++;
+    m_reduce = true;
+}
+
+void StaticStructureFactorDebye::reduce()
+{
+    m_structure_factor.prepare(m_histogram.shape());
+    m_local_histograms.reduceInto(m_structure_factor);
+
+    // Normalize by the frame count if necessary
+    if (m_frame_counter > 1)
+    {
+        util::forLoopWrapper(0, m_structure_factor.size(), [this](size_t begin, size_t end) {
+            for (size_t i = begin; i < end; ++i)
+            {
+                m_structure_factor[i] /= static_cast<float>(m_frame_counter);
+            }
+        });
+    }
+}
+
+}; }; // namespace freud::diffraction

cpp/diffraction/StaticStructureFactorDebye.h

@@ -0,0 +1,94 @@
+// Copyright (c) 2010-2020 The Regents of the University of Michigan
+// This file is from the freud project, released under the BSD 3-Clause License.
+
+#ifndef STATIC_STRUCTURE_FACTOR_DEBYE_H
+#define STATIC_STRUCTURE_FACTOR_DEBYE_H
+
+#include "Box.h"
+#include "Histogram.h"
+#include "NeighborQuery.h"
+
+/*! \file StaticStructureFactorDebye.h
+    \brief Routines for computing static structure factors.
+
+    Computes structure factors from the Fourier transform of a radial distribution function (RDF).
+    This method is not capable of resolving k-vectors smaller than the magnitude 4 * pi / L, where L is the
+    smallest side length of the system's periodic box.
+*/
+
+namespace freud { namespace diffraction {
+
+class StaticStructureFactorDebye
+{
+    using StaticStructureFactorDebyeHistogram = util::Histogram<float>;
+
+public:
+    //! Constructor
+    StaticStructureFactorDebye(unsigned int bins, float k_max, float k_min = 0);
+
+    //! Destructor
+    virtual ~StaticStructureFactorDebye() = default;
+
+    //! Compute the structure factor S(k) using the Debye formula
+    void accumulate(const freud::locality::NeighborQuery* neighbor_query, const vec3<float>* query_points,
+                    unsigned int n_query_points, unsigned int n_total);
+
+    //! Reduce thread-local arrays onto the primary data arrays.
+    void reduce();
+
+    //! Return thing_to_return after reducing if necessary.
+    template<typename U> U& reduceAndReturn(U& thing_to_return)
+    {
+        if (m_reduce)
+        {
+            reduce();
+        }
+        m_reduce = false;
+        return thing_to_return;
+    }
+
+    //! Reset the histogram to all zeros
+    void reset()
+    {
+        m_local_histograms.reset();
+        m_frame_counter = 0;
+        m_reduce = true;
+    }
+
+    //! Get the structure factor
+    const util::ManagedArray<float>& getStructureFactor()
+    {
+        return reduceAndReturn(m_structure_factor);
+    }
+
+    //! Get the k bin edges
+    std::vector<float> getBinEdges() const
+    {
+        return m_histogram.getBinEdges()[0];
+    }
+
+    //! Get the k bin centers
+    std::vector<float> getBinCenters() const
+    {
+        return m_histogram.getBinCenters()[0];
+    }
+
+    //! Get the minimum valid k value
+    float getMinValidK() const
+    {
+        return m_min_valid_k;
+    }
+
+private:
+    StaticStructureFactorDebyeHistogram m_histogram; //!< Histogram to hold computed structure factor
+    StaticStructureFactorDebyeHistogram::ThreadLocalHistogram
+        m_local_histograms;                       //!< Thread local histograms for TBB parallelism
+    util::ManagedArray<float> m_structure_factor; //!< The computed structure factor
+    unsigned int m_frame_counter {0};                 //!< Number of frames calculated
+    float m_min_valid_k {0};                          //!< The minimum valid k-value based on the computed box
+    bool m_reduce {true};                                //!< Whether to reduce
+};
+
+}; }; // namespace freud::diffraction
+
+#endif // STATIC_STRUCTURE_FACTOR_DEBYE_H

cpp/locality/NeighborQuery.h

@@ -51,6 +51,20 @@ struct QueryArgs
 {
     QueryArgs() = default;
 
+   //! Constructor
+    QueryArgs(QueryType mode, unsigned int num_neighbors, float r_max, float r_min, float r_guess,
+              float scale, bool exclude_ii)
+        : mode(mode), num_neighbors(num_neighbors), r_max(r_max), r_min(r_min), r_guess(r_guess),
+          scale(scale), exclude_ii(exclude_ii)
+    {}
+
+    //! Ball query factory
+    static QueryArgs make_ball(float r_max, float r_min = 0, bool exclude_ii = DEFAULT_EXCLUDE_II)
+    {
+        return QueryArgs(QueryType::ball, DEFAULT_NUM_NEIGHBORS, r_max, r_min, DEFAULT_R_GUESS, DEFAULT_SCALE,
+                         exclude_ii);
+    }
+
     QueryType mode {DEFAULT_MODE}; //! Whether to perform a ball or k-nearest neighbor query.
     unsigned int num_neighbors {DEFAULT_NUM_NEIGHBORS}; //! The number of nearest neighbors to find.
     float r_max {DEFAULT_R_MAX};          //! The cutoff distance within which to find neighbors.

cpp/util/Histogram.h

@@ -72,13 +72,13 @@ class Axis
     virtual size_t bin(const float& value) const = 0;
 
     //! Return the boundaries of bins.
-    std::vector<float> getBinEdges() const
+    const std::vector<float> getBinEdges() const
     {
         return m_bin_edges;
     }
 
     //! Return the centers of bins.
-    std::vector<float> getBinCenters() const
+    const std::vector<float> getBinCenters() const
     {
         std::vector<float> bin_centers(m_nbins);
         for (size_t i = 0; i < m_nbins; i++)
@@ -369,7 +369,7 @@ template<typename T> class Histogram
     }
 
     //! Get the shape of the computed histogram.
-    std::vector<size_t> shape() const
+    const std::vector<size_t> shape() const
     {
         return m_bin_counts.shape();
     }
@@ -383,7 +383,7 @@ template<typename T> class Histogram
     //! Return the edges of bins.
     /*! This vector will be of size axis.size()+1 for each axis.
      */
-    std::vector<std::vector<float>> getBinEdges() const
+    const std::vector<std::vector<float>> getBinEdges() const
     {
         std::vector<std::vector<float>> bins(m_axes.size());
         for (unsigned int i = 0; i < m_axes.size(); ++i)
@@ -396,7 +396,7 @@ template<typename T> class Histogram
     //! Return the bin centers.
     /*! This vector will be of size axis.size() for each axis.
      */
-    std::vector<std::vector<float>> getBinCenters() const
+    const std::vector<std::vector<float>> getBinCenters() const
     {
         std::vector<std::vector<float>> bins(m_axes.size());
         for (unsigned int i = 0; i < m_axes.size(); ++i)
@@ -407,7 +407,7 @@ template<typename T> class Histogram
     }
 
     //! Return a vector of tuples (min, max) indicating the bounds of each axis.
-    std::vector<std::pair<float, float>> getBounds() const
+    const std::vector<std::pair<float, float>> getBounds() const
     {
         std::vector<std::pair<float, float>> bounds(m_axes.size());
         for (unsigned int i = 0; i < m_axes.size(); ++i)
@@ -418,7 +418,7 @@ template<typename T> class Histogram
     }
 
     //! Return a vector indicating the number of bins in each axis.
-    std::vector<size_t> getAxisSizes() const
+    const std::vector<size_t> getAxisSizes() const
     {
         std::vector<size_t> sizes(m_axes.size());
         for (unsigned int i = 0; i < m_axes.size(); ++i)
@@ -489,7 +489,7 @@ template<typename T> class Histogram
      * variadic templating to accept an arbitrary set of float values and
      * construct a vector out of them.
      */
-    std::pair<std::vector<float>, Weight<T>> getValueVector(float value) const
+    const std::pair<std::vector<float>, Weight<T>> getValueVector(float value) const
     {
         return {{value}, Weight<T>()};
     }
@@ -499,14 +499,14 @@ template<typename T> class Histogram
      * variadic templating to accept an arbitrary set of float values and
      * construct a vector out of them.
      */
-    std::pair<std::vector<float>, Weight<T>> getValueVector(Weight<T> weight) const
+    const std::pair<std::vector<float>, Weight<T>> getValueVector(Weight<T> weight) const
     {
         return {{}, weight};
     }
 
     //! The recursive case for constructing a vector of values (see base-case function docs).
     template<typename... FloatsOrWeight>
-    std::pair<std::vector<float>, Weight<T>> getValueVector(float value, FloatsOrWeight... values) const
+    const std::pair<std::vector<float>, Weight<T>> getValueVector(float value, FloatsOrWeight... values) const
     {
         std::pair<std::vector<float>, Weight<T>> tmp = getValueVector(values...);
         tmp.first.insert(tmp.first.begin(), value);
@@ -515,7 +515,8 @@ template<typename T> class Histogram
 
     //! The recursive case for constructing a vector of values (see base-case function docs).
     template<typename... FloatsOrWeight>
-    std::pair<std::vector<float>, Weight<T>> getValueVector(Weight<T> weight, FloatsOrWeight... values) const
+    const std::pair<std::vector<float>, Weight<T>> getValueVector(Weight<T> weight,
+                                                                  FloatsOrWeight... values) const
     {
         std::pair<std::vector<float>, Weight<T>> tmp = getValueVector(values...);
         tmp.second = weight;