Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

NEW feature: intermediate scattering function #1042

Open
wants to merge 24 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
24 commits
Select commit Hold shift + click to select a range
af1da72
boilerplate for intermediate scattering
Nov 30, 2022
565a0be
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Nov 30, 2022
acc9b40
complete the framework of isf, need to debug and test
Dec 1, 2022
bec9619
fix conflict
Dec 1, 2022
fcf006c
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Dec 1, 2022
32c340d
fix lots of bugs, compile pass
Dec 3, 2022
9138a86
new problem: We cant set new attribute in diffraction.pyx, so we cant…
Dec 3, 2022
a495664
fix conflicts
Dec 3, 2022
18ac4fb
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Dec 3, 2022
96f59d0
restructure previous code, comment out intermediate scattering function
tommy-waltmann Dec 7, 2022
41b57de
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Dec 7, 2022
c09064c
linting
tommy-waltmann Dec 7, 2022
3323202
Merge branch 'dev-isf' of github.com:roy-kid/freud into dev-isf
tommy-waltmann Dec 7, 2022
9575ee6
Merge branch 'master' into dev-isf
tommy-waltmann Dec 7, 2022
aabf02b
put boilerplate in place for implementing intermediate scattering fun…
tommy-waltmann Dec 7, 2022
9671b5f
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Dec 7, 2022
0402e0f
more code restructuring
tommy-waltmann Dec 13, 2022
99d4b6e
finish code restructure, implementation is wrong
tommy-waltmann Dec 14, 2022
82cbbcc
linting
tommy-waltmann Dec 14, 2022
d2a74c1
get fatal error at IntermediateScattering.cc line 121 and 122: Attem…
Dec 25, 2022
53642c4
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Dec 25, 2022
fd08212
Merge branch 'master' into dev-isf
tommy-waltmann Jan 3, 2023
78e2bb6
fix oob crashes
tommy-waltmann Jan 3, 2023
86374b0
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Jan 3, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -22,3 +22,5 @@ freud/*.cpp
freud/*.cxx
freud/*.dylib
freud/*.so

.vscode/*
14 changes: 11 additions & 3 deletions cpp/diffraction/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,8 +1,16 @@
add_library(
_diffraction OBJECT
StaticStructureFactor.h StaticStructureFactor.cc StaticStructureFactorDebye.h
StaticStructureFactorDebye.cc StaticStructureFactorDirect.h
StaticStructureFactorDirect.cc)
StructureFactor.h
StaticStructureFactor.h
StaticStructureFactor.cc
StaticStructureFactorDebye.h
StaticStructureFactorDebye.cc
StructureFactorDirect.h
StructureFactorDirect.cc
StaticStructureFactorDirect.h
StaticStructureFactorDirect.cc
IntermediateScattering.h
IntermediateScattering.cc)

target_link_libraries(_diffraction PUBLIC TBB::tbb)

Expand Down
167 changes: 167 additions & 0 deletions cpp/diffraction/IntermediateScattering.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,167 @@
// 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 "IntermediateScattering.h"
#include "Box.h"
#include "ManagedArray.h"
#include "NeighborQuery.h"
#include "utils.h"

/*! \file IntermediateScattering.cc
\brief Routines for computing intermediate scattering function.
*/

namespace freud { namespace diffraction {

IntermediateScattering::IntermediateScattering(const box::Box& box, unsigned int bins, float k_max,
float k_min, unsigned int num_sampled_k_points)
: StructureFactorDirect(bins, k_max, k_min, num_sampled_k_points), StructureFactor(bins, k_max, k_min),
m_box(box)
{
if (bins == 0)
{
throw std::invalid_argument("IntermediateScattering requires a nonzero number of bins.");
}
if (k_max <= 0)
{
throw std::invalid_argument("IntermediateScattering requires k_max to be positive.");
}
if (k_min < 0)
{
throw std::invalid_argument("IntermediateScattering requires k_min to be non-negative.");
}
if (k_max <= k_min)
{
throw std::invalid_argument("IntermediateScattering requires that k_max must be greater than k_min.");
}

m_k_points = IntermediateScattering::reciprocal_isotropic(box, k_max, k_min, m_num_sampled_k_points);
}

void IntermediateScattering::compute(const vec3<float>* points, unsigned int num_points,
const vec3<float>* query_points, unsigned int num_query_points,
unsigned int num_frames, unsigned int n_total)
{
// initialize the histograms, now that the size of both axes are known
m_self_function = StructureFactorHistogram({
std::make_shared<util::RegularAxis>(num_frames, -0.5, num_frames - 0.5),
std::make_shared<util::RegularAxis>(m_nbins, m_k_min, m_k_max),
});
m_local_self_function = StructureFactorHistogram::ThreadLocalHistogram(m_self_function);
m_distinct_function = StructureFactorHistogram({
std::make_shared<util::RegularAxis>(num_frames, -0.5, num_frames - 0.5),
std::make_shared<util::RegularAxis>(m_nbins, m_k_min, m_k_max),
});
m_local_distinct_function = StructureFactorHistogram::ThreadLocalHistogram(m_distinct_function);
// Compute k vectors by sampling reciprocal space.
if (m_box.is2D())
{
throw std::invalid_argument("2D boxes are not currently supported.");
}

// The minimum valid k value is 2 * pi / L, where L is the smallest side length.
const auto box_L = m_box.getL();
const auto min_box_length = std::min(box_L.x, std::min(box_L.y, box_L.z));
m_min_valid_k = freud::constants::TWO_PI / min_box_length;

// record the point at t=0
const vec3<float>* query_r0 = &query_points[0];
const vec3<float>* r0 = &points[0];

util::forLoopWrapper(0, num_frames, [&](size_t begin, size_t end) {
for (size_t t = begin; t < end; t++)
{
size_t offset = t * num_points;
// Compute self-part
const auto self_part
= IntermediateScattering::compute_self(&points[offset], r0, num_points, n_total, m_k_points);

// Compute distinct-part
const auto distinct_part = IntermediateScattering::compute_distinct(
&points[offset], query_r0, num_query_points, n_total, m_k_points);

std::vector<float> S_k_self_part = StaticStructureFactorDirect::compute_S_k(self_part, self_part);
std::vector<float> S_k_distinct_part
= StaticStructureFactorDirect::compute_S_k(distinct_part, distinct_part);

// Bin the S_k values and track the number of k values in each bin.
util::forLoopWrapper(0, m_k_points.size(), [&](size_t begin, size_t end) {
for (size_t k_index = begin; k_index < end; ++k_index)
{
const auto& k_vec = m_k_points[k_index];
const auto k_magnitude = std::sqrt(dot(k_vec, k_vec));
const auto k_bin = m_k_histogram.bin({k_magnitude});
m_local_self_function.increment(k_bin, S_k_self_part[k_index]);
m_local_distinct_function.increment(k_bin, S_k_distinct_part[k_index]);
m_local_k_histograms.increment(k_bin);
}
});
}
});

m_reduce = true;
}

void IntermediateScattering::reduce()
{
const auto k_axis_size = m_k_histogram.getAxisSizes();
m_k_histogram.prepare(k_axis_size);
const auto self_axis_size = m_self_function.getAxisSizes();
m_self_function.prepare(self_axis_size);
const auto distinct_axis_size = m_distinct_function.getAxisSizes();
m_distinct_function.prepare(distinct_axis_size);

// Reduce the bin counts over all threads, then use them to normalize the
// structure factor when computing. This computes a binned mean over all k
// points.
m_k_histogram.reduceOverThreads(m_local_k_histograms);
m_self_function.reduceOverThreadsPerBin(m_local_self_function, [&](size_t i) {
m_self_function[i] /= m_k_histogram[util::modulusPositive(i, m_nbins)];
});
m_distinct_function.reduceOverThreadsPerBin(m_local_distinct_function, [&](size_t i) {
m_distinct_function[i] /= m_k_histogram[util::modulusPositive(i, m_nbins)];
});
}

std::vector<std::complex<float>>
IntermediateScattering::compute_self(const vec3<float>* rt, const vec3<float>* r0, unsigned int n_points,
unsigned int n_total, const std::vector<vec3<float>>& k_points)
{
//
std::vector<vec3<float>> r_i_t0(n_points); // rt - rt0 element-wisely
util::forLoopWrapper(0, n_points, [&](size_t begin, size_t end) {
for (size_t i = begin; i < end; ++i)
{
r_i_t0[i] = rt[i] - rt[0];
}
});

return StructureFactorDirect::compute_F_k(r_i_t0.data(), n_points, n_total, m_k_points);
}

std::vector<std::complex<float>>
IntermediateScattering::compute_distinct(const vec3<float>* rt, const vec3<float>* r0, unsigned int n_points,
unsigned int n_total, const std::vector<vec3<float>>& k_points)
{
const auto n_rij = n_points * (n_points - 1);
std::vector<vec3<float>> r_ij(n_rij);
size_t i = 0;

util::forLoopWrapper(0, n_rij, [&](size_t begin, size_t end) {
for (size_t rt_index = begin; rt_index < end; ++rt_index)
{
for (size_t r0_index = 0; r0_index < n_rij; ++r0_index)
{
if (rt_index != r0_index)
{
r_ij[i] = rt[rt_index] - r0[r0_index];
++i;
}
}
};
});

return StructureFactorDirect::compute_F_k(r_ij.data(), n_rij, n_total, m_k_points);
}

}} // namespace freud::diffraction
101 changes: 101 additions & 0 deletions cpp/diffraction/IntermediateScattering.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,101 @@
// 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 INTERMEDIATE_SCATTERING_H
#define INTERMEDIATE_SCATTERING_H

#include <complex>
#include <limits>
#include <vector>

#include "Box.h"
#include "Histogram.h"
#include "ManagedArray.h"
#include "NeighborQuery.h"
#include "StaticStructureFactorDirect.h"

/*! \file IntermediateScattering.h
\brief Routines for computing intermediate scattering function.

This computes the intermediate scattering function F(k, t) between sets of
time-dependent points and query points.

The k-vectors are used to compute complex scattering amplitudes F(k, t) by
summing over the scattering contribution of all particle positions (atomic
form factors are assumed to be 1) at every k-vector. The complex conjugate
of the scattering amplitudes of the points at each k-vector are multiplied
by the scattering amplitudes of the query points at each k-vector. Finally,
the results are binned according to their k-vector's magnitude and normalized
by the number of samples in each radial bin in k-space.

Note that k-vectors are in the physics convention, and q-vectors are in the
crystallographic convention. These conventions differ by a factor of 2\pi.

See also:
https://en.wikipedia.org/wiki/Reciprocal_lattice#Arbitrary_collection_of_atoms
*/

namespace freud { namespace diffraction {

class IntermediateScattering : public StructureFactorDirect
{
public:
//! Constructor
IntermediateScattering(const box::Box& box, unsigned int bins, float k_max, float k_min = 0,
unsigned int num_sampled_k_points = 0);

void compute(const vec3<float>* points, unsigned int num_points, const vec3<float>* query_points,
unsigned int num_query_points, unsigned int num_frames, unsigned int n_total);

std::vector<float> getBinEdges() const
{
// the second axis of the histogram is the k-axis
return m_self_function.getBinEdges()[1];
}

std::vector<float> getBinCenters() const
{
// the second axis of the histogram is the k-axis
return m_self_function.getBinCenters()[1];
}

// more intuitively named alias for getStructureFactor
const util::ManagedArray<float>& getSelfFunction()
{
return reduceAndReturn(m_self_function.getBinCounts());
}

const util::ManagedArray<float>& getDistinctFunction()
{
return reduceAndReturn(m_distinct_function.getBinCounts());
}

private:
void reduce() override;

//!< box for the calculation, we assume the box is constant over the time interval
box::Box m_box;

//!< Histogram to hold self part of the ISF at each frame
StructureFactorHistogram m_self_function {};
//!< Thread local histograms for TBB parallelism
StructureFactorHistogram::ThreadLocalHistogram m_local_self_function {};

//!< Histogram to hold distinct part of the ISF at each frame
StructureFactorHistogram m_distinct_function {};
//!< Thread local histograms for TBB parallelism
StructureFactorHistogram::ThreadLocalHistogram m_local_distinct_function {};

//!< Helpers to compute self and distinct parts
std::vector<std::complex<float>> compute_self(const vec3<float>* rt, const vec3<float>* r0,
unsigned int n_points, unsigned int n_total,
const std::vector<vec3<float>>& k_points);

std::vector<std::complex<float>> compute_distinct(const vec3<float>* rt, const vec3<float>* r0,
unsigned int n_points, unsigned int n_total,
const std::vector<vec3<float>>& k_points);
};

}; }; // namespace freud::diffraction

#endif // INTERMEDIATE_SCATTERING_H
3 changes: 2 additions & 1 deletion cpp/diffraction/StaticStructureFactor.cc
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,8 @@
namespace freud { namespace diffraction {

StaticStructureFactor::StaticStructureFactor(unsigned int bins, float k_max, float k_min)
: m_structure_factor({std::make_shared<util::RegularAxis>(bins, k_min, k_max)}),
: StructureFactor(bins, k_max, k_min),
m_structure_factor({std::make_shared<util::RegularAxis>(bins, k_min, k_max)}),
m_local_structure_factor(m_structure_factor)
{
// Validation logic is not shared in the parent StaticStructureFactor
Expand Down
52 changes: 23 additions & 29 deletions cpp/diffraction/StaticStructureFactor.h
Original file line number Diff line number Diff line change
Expand Up @@ -10,28 +10,43 @@
#include "Histogram.h"
#include "ManagedArray.h"
#include "NeighborQuery.h"
#include "StructureFactor.h"

/*! \file StaticStructureFactor.h
\brief Base class for structure factor classes.
\brief Base class for static tructure factor classes.
*/

namespace freud { namespace diffraction {

class StaticStructureFactor
/* Abstract base class for all static structure factors.
*
* A static structure factor is a structure factor which can be computed using
* only the data from one frame of a simulation. A typical use case is to compute
* the static structure factor for each frame over many frames of a simulation,
* and get an average for better statistics/curve smoothness. To support this use
* case, all static structure factors must have logic for either continuing to
* accumulate histogram data or resetting the data on each successive call to
* accumulate().
*
* This class does not have any members, it simply defines the interface for
* all static structure factors.
*
* */
class StaticStructureFactor : virtual public StructureFactor
{
protected:
using StructureFactorHistogram = util::Histogram<float>;

StaticStructureFactor(unsigned int bins, float k_max, float k_min = 0);

public:
virtual ~StaticStructureFactor() = default;

//!< do the histogram calculation, and accumulate if necessary
virtual void accumulate(const freud::locality::NeighborQuery* neighbor_query,
const vec3<float>* query_points, unsigned int n_query_points,
unsigned int n_total)
= 0;

//<! reset the structure factor if the user doesn't wish to accumulate
virtual void reset() = 0;

//! Get the structure factor
Expand All @@ -52,32 +67,11 @@ class StaticStructureFactor
return m_structure_factor.getBinCenters()[0];
}

//! Get the minimum valid k value
float getMinValidK() const
{
return m_min_valid_k;
}

protected:
virtual void reduce() = 0;

//! 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;
}

StructureFactorHistogram m_structure_factor; //!< Histogram to hold computed structure factor
StructureFactorHistogram::ThreadLocalHistogram
m_local_structure_factor; //!< Thread local histograms for TBB parallelism

bool m_reduce {true}; //! Whether to reduce local histograms
float m_min_valid_k {std::numeric_limits<float>::infinity()}; //! Minimum valid k-vector magnitude
//!< Histogram to hold computed structure factor
StructureFactorHistogram m_structure_factor;
//!< Thread local histograms for TBB parallelism
StructureFactorHistogram::ThreadLocalHistogram m_local_structure_factor;
};

}; }; // namespace freud::diffraction
Expand Down
Loading