Static structure factor S(q)

.circleci/ci-oldest-reqs.txt

@@ -1,3 +1,4 @@
+ase==3.16.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,108 @@
+// 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_min < 0)
+    {
+        throw std::invalid_argument("StaticStructureFactorDebye requires k_min to be non-negative.");
+    }
+    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
+    const auto axes = S_kHistogram::Axes {std::make_shared<util::RegularAxis>(bins, k_min, k_max)};
+    m_histogram = S_kHistogram(axes);
+    m_local_histograms = S_kHistogram::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)
+{
+    const auto& box = neighbor_query->getBox();
+
+    // The r_max should be just less than half of the smallest side length of the box
+    const auto box_L = box.getL();
+    const auto 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));
+    const auto r_max = std::nextafter(float(0.5) * min_box_length, float(0));
+    const auto* const points = neighbor_query->getPoints();
+    const auto 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);
+
+    std::vector<float> distances(n_points * n_query_points);
+    box.computeAllDistances(points, n_points, query_points, n_query_points, distances.data());
+
+    const auto k_bin_centers = m_histogram.getBinCenters()[0];
+
+    util::forLoopWrapper(0, m_histogram.getAxisSizes()[0], [&](size_t begin, size_t end) {
+        for (size_t k_index = begin; k_index < end; ++k_index)
+        {
+            const auto k = k_bin_centers[k_index];
+            double S_k = 0.0;
+            for (const auto& distance : distances)
+            {
+                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(), [&](size_t begin, size_t end) {
+            for (size_t k_index = begin; k_index < end; ++k_index)
+            {
+                m_structure_factor[k_index] /= static_cast<float>(m_frame_counter);
+            }
+        });
+    }
+}
+
+}; }; // namespace freud::diffraction

cpp/diffraction/StaticStructureFactorDebye.h

@@ -0,0 +1,91 @@
+// 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 S_kHistogram = util::Histogram<float>;
+
+public:
+    //! Constructor
+    StaticStructureFactorDebye(unsigned int bins, float k_max, float k_min = 0);
+
+    //! 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:
+    S_kHistogram m_histogram; //!< Histogram to hold computed structure factor
+    S_kHistogram::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

@@ -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;

cpp/util/utils.h

@@ -3,6 +3,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <stdexcept>
 #include <tbb/blocked_range.h>
 #include <tbb/blocked_range2d.h>
 #include <tbb/parallel_for.h>
@@ -21,11 +22,64 @@ inline float clamp(float v, float lo, float hi)
     \returns The remainder of a/b, between min(0, b) and max(0, b).
     \note This is the same behavior of the modulus operator % in Python (but not C++).
 */
-template<class Scalar> inline Scalar modulusPositive(Scalar a, Scalar b)
+template<typename Scalar> inline Scalar modulusPositive(Scalar a, Scalar b)
 {
     return std::fmod(std::fmod(a, b) + b, b);
 }
 
+//! Sinc function, sin(x)/x
+/*! \param x Argument.
+    \returns Sin of x divided by x.
+    \note There is no factor of pi in this definition (some conventions include pi).
+*/
+template<typename Scalar> inline Scalar sinc(Scalar x)
+{
+    if (x == 0)
+    {
+        return 1;
+    }
+    return std::sin(x) / x;
+}
+
+//! Simpson's Rule numerical integration
+/*! \param integrand Callable that returns the integrand value for a provided bin index.
+    \param num_bins Number of bins to integrate over. Must be odd.
+    \param dx Step size between bins.
+    \note The integration summation is performed in double-precision regardless of Scalar type.
+*/
+template<typename Integrand> inline double simpson_integrate(Integrand& integrand, size_t num_bins, double dx)
+{
+    if (num_bins % 2 != 1)
+    {
+        // This only implements the easiest case of Simpson's rule.
+        // Even numbers of bins require additional logic.
+        throw std::invalid_argument("The number of integration bins must be odd.");
+    }
+
+    double integral = 0.0;
+
+    // Simpson's rule uses prefactors 1, 4, 2, 4, 2, ..., 4, 1
+    auto simpson_prefactor = [=](size_t bin) {
+        if (bin == 0 || bin == num_bins - 1)
+        {
+            return 1.0;
+        }
+        if (bin % 2 == 0)
+        {
+            return 2.0;
+        }
+        return 4.0;
+    };
+
+    for (size_t bin_index = 0; bin_index < num_bins; bin_index++)
+    {
+        integral += simpson_prefactor(bin_index) * integrand(bin_index);
+    }
+
+    integral *= dx / 3.0;
+    return integral;
+}
+
 //! Wrapper for for-loop to allow the execution in parallel or not.
 /*! \param begin Beginning index.
  *  \param end Ending index.

doc/source/modules/diffraction.rst

@@ -8,6 +8,7 @@ Diffraction Module
     :nosignatures:
 
     freud.diffraction.DiffractionPattern
+    freud.diffraction.StaticStructureFactorDebye
 
 .. rubric:: Details
 

doc/source/reference/freud.bib

@@ -1,3 +1,15 @@
+@article{Farrow2009,
+author = "Farrow, Christopher L. and Billinge, Simon J. L.",
+title = "{Relationship between the atomic pair distribution function and small-angle scattering: implications for modeling of nanoparticles}",
+journal = "Acta Crystallographica Section A",
+year = "2009",
+volume = "65",
+number = "3",
+pages = "232--239",
+doi = {10.1107/S0108767309009714},
+url = {https://doi.org/10.1107/S0108767309009714}
+}
+
 @article{Filion_2010,
   author        = {Filion, L. and Hermes, M. and Ni, R. and Dijkstra, M.},
   date-added    = {2019-10-03 00:25:57 -0400},
@@ -205,3 +217,16 @@ @article{Mickel2013
   url     = {https://doi.org/10.1063/1.4774084},
   eprint  = {arXiv:1209.6180}
 }
+
+@article{Liu2016,
+  author    = {Hongjun Liu and Stephen J. Paddison},
+  title     = {Direct calculation of the X-ray structure factor of ionic liquids},
+  journal   = {Phys. Chem. Chem. Phys.},
+  year      = {2016},
+  volume    = {18},
+  issue     = {16},
+  pages     = {11000-11007},
+  publisher = {The Royal Society of Chemistry},
+  doi       = {10.1039/C5CP06199G},
+  url       = {https://dx.doi.org/10.1039/C5CP06199G}
+}