From cfda4ed3d096d381af67ff8bee8529a37f9920cf Mon Sep 17 00:00:00 2001 From: Mikael Lund Date: Fri, 30 Jan 2026 06:55:54 +0100 Subject: [PATCH] Add atom-specific scattering form factors MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add scattering_f0 property to AtomData (defaults to 1.0) - Create FormFactorAtomicConstant class using atom-specific form factors - Normalize scattering intensity by sum of squared form factors (Σfᵢ²) - Update ScatteringFunction to use Scatterer struct with position and atom id - Add tests for new form factor functionality - Update documentation in analysis.md and topology.md Backwards compatible: when scattering_f0=1.0 (default), normalization by Σfᵢ² equals normalization by N, giving identical results. --- docs/_docs/analysis.md | 13 ++-- docs/_docs/topology.md | 1 + src/analysis.cpp | 12 ++-- src/analysis.h | 8 +-- src/atomdata.cpp | 7 +- src/atomdata.h | 1 + src/scatter.cpp | 37 ++++++++-- src/scatter.h | 158 +++++++++++++++++++++++++++-------------- 8 files changed, 165 insertions(+), 72 deletions(-) diff --git a/docs/_docs/analysis.md b/docs/_docs/analysis.md index 5dfa88319..cee0118c6 100644 --- a/docs/_docs/analysis.md +++ b/docs/_docs/analysis.md @@ -159,17 +159,20 @@ and as a function of separation, _r_. In addition, the radial distribution funct ### Structure Factor -The isotropically averaged static structure factor between $N$ point scatterers is calculated using +The isotropically averaged static structure factor between $N$ scatterers is calculated using the [Debye formula](http://doi.org/dmb9wm), $$ - S(q) = 1 + \frac{2}{N} \left \langle - \sum_{i=1}^{N-1}\sum_{j=i+1}^N \frac{\sin(qr_{ij})}{qr_{ij}} + I(q) = \frac{1}{\sum_i f_i^2} \left \langle + \sum_{i}^{N}\sum_{j}^N f_i f_j \frac{\sin(qr_{ij})}{qr_{ij}} \right \rangle $$ -The selected `molecules` can be treated either as single point scatterers (`com=true`) or as a group of individual -point scatterers of equal intensity, i.e., with a form factor of unity. +where $f_i$ is the (q-independent) atomic form factor for particle $i$, defined by `scattering_f0` in the +[atom properties](topology/#atoms). If not specified, `scattering_f0` defaults to 1. + +The selected `molecules` can be treated either as single point scatterers (`com=true`) or as a group of +individual scatterers with atom-specific form factors. The computation of the structure factor is rather computationally intensive task, scaling quadratically with the number of particles and linearly with the number of scattering vector mesh points. If OpenMP is available, multiple threads diff --git a/docs/_docs/topology.md b/docs/_docs/topology.md index 2b27b93e5..47a0ba7fa 100644 --- a/docs/_docs/topology.md +++ b/docs/_docs/topology.md @@ -72,6 +72,7 @@ Atoms are the smallest possible particle entities with properties defined below. `sigma=0` | `2r` [Å] (overrides radius) `tension=0` | Surface tension [kJ/mol/Å$^2$] `tfe=0` | Transfer free energy [kJ/mol/Å$^2$/M] +`scattering_f0=1` | Atomic scattering form factor (q-independent) `psc` | Patchy sphero-cylinders properties (object) A filename (`.json`) may be given instead of an atom definition to load diff --git a/src/analysis.cpp b/src/analysis.cpp index f316bd848..a5d17de2c 100644 --- a/src/analysis.cpp +++ b/src/analysis.cpp @@ -2676,11 +2676,12 @@ void ScatteringFunction::_sample() std::ranges::for_each(groups, [&](auto& group) { if (mass_center_scattering && group.isMolecular()) { - scatter_positions.push_back(group.mass_center); + scatter_positions.push_back({group.mass_center, -1}); // id=-1 for unity form factor } else { - auto positions = group.positions(); - std::copy(positions.begin(), positions.end(), std::back_inserter(scatter_positions)); + for (const auto& particle : group) { + scatter_positions.push_back({particle.pos, particle.id}); + } } }); @@ -2759,11 +2760,12 @@ try : Analysis(spc, "scatter") { const int pmax = j.value("pmax", 15); if (ipbc) { scheme = Schemes::EXPLICIT_IPBC; - explicit_average_ipbc = std::make_unique>(pmax); + explicit_average_ipbc = + std::make_unique>(pmax); } else { scheme = Schemes::EXPLICIT_PBC; - explicit_average_pbc = std::make_unique>(pmax); + explicit_average_pbc = std::make_unique>(pmax); } } else { diff --git a/src/analysis.h b/src/analysis.h index 695177915..4835fc876 100644 --- a/src/analysis.h +++ b/src/analysis.h @@ -841,14 +841,14 @@ class ScatteringFunction : public Analysis bool mass_center_scattering; //!< scatter from mass center, only? bool save_after_sample = false; //!< if true, save average S(q) after each sample point std::string filename; //!< output file name - std::vector scatter_positions; //!< vector of scattering points + std::vector scatter_positions; //!< vector of scatterers with position and atom id std::vector molecule_ids; //!< Molecule ids std::vector molecule_names; //!< Molecule names corresponding to `molecule_ids` - using Tformfactor = Scatter::FormFactorUnity; + using Tformfactor = Scatter::FormFactorAtomicConstant; std::unique_ptr> debye; - std::unique_ptr> explicit_average_pbc; - std::unique_ptr> explicit_average_ipbc; + std::unique_ptr> explicit_average_pbc; + std::unique_ptr> explicit_average_ipbc; void _sample() override; void _to_disk() override; void _to_json(json& j) const override; diff --git a/src/atomdata.cpp b/src/atomdata.cpp index aaaae798b..bf11c99a4 100644 --- a/src/atomdata.cpp +++ b/src/atomdata.cpp @@ -99,6 +99,7 @@ void to_json(json& j, const AtomData& a) {"mu", a.mu}, {"mulen", a.mulen}, {"psc", a.sphero_cylinder}, + {"scattering_f0", a.scattering_f0}, {"id", a.id()}}; if (a.dp.has_value()) { _j["dp"] = a.dp.value() / 1.0_angstrom; @@ -155,6 +156,7 @@ void from_json(const json& j, AtomData& a) a.tfe = val.value("tfe", a.tfe) * 1.0_kJmol / (1.0_angstrom * 1.0_angstrom * 1.0_molar); a.hydrophobic = val.value("hydrophobic", false); a.implicit = val.value("implicit", false); + a.scattering_f0 = val.value("scattering_f0", 1.0); set_dp_and_dprot(val, a); if (val.contains("activity")) { a.activity = val.at("activity").get() * 1.0_molar; @@ -247,7 +249,7 @@ TEST_CASE("[Faunus] AtomData") json j = R"({ "atomlist" : [ { "A": { "sigma": 2.5, "pactivity":2, "eps_custom": 0.1 } }, - { "B": { "r":1.1, "activity":0.2, "eps":0.05, "dp":9.8, "dprot":3.14, "mw":1.1, "tfe":0.98, "tension":0.023 } } + { "B": { "r":1.1, "activity":0.2, "eps":0.05, "dp":9.8, "dprot":3.14, "mw":1.1, "tfe":0.98, "tension":0.023, "scattering_f0": 7.5 } } ]})"_json; pc::temperature = 298.15_K; @@ -265,6 +267,8 @@ TEST_CASE("[Faunus] AtomData") // "unknown atom property"); CHECK_EQ(v.front().sigma, Approx(2.5e-10_m)); CHECK_EQ(v.front().activity, Approx(0.01_molar)); + CHECK_EQ(v.front().scattering_f0, Approx(1.0)); // default value + CHECK_EQ(v.back().scattering_f0, Approx(7.5)); // explicit value CHECK_EQ(v.back().tfe, Approx(0.98_kJmol / (1.0_angstrom * 1.0_angstrom * 1.0_molar))); AtomData a = json(v.back()); // AtomData -> JSON -> AtomData @@ -277,6 +281,7 @@ TEST_CASE("[Faunus] AtomData") CHECK_EQ(a.dp, Approx(9.8)); CHECK_EQ(a.dprot, Approx(3.14)); CHECK_EQ(a.mw, Approx(1.1)); + CHECK_EQ(a.scattering_f0, Approx(7.5)); // check JSON round-trip CHECK_EQ(a.tfe, Approx(0.98_kJmol / 1.0_angstrom / 1.0_angstrom / 1.0_molar)); CHECK_EQ(a.tension, Approx(0.023_kJmol / 1.0_angstrom / 1.0_angstrom)); diff --git a/src/atomdata.h b/src/atomdata.h index 10094ee9a..51042cb76 100644 --- a/src/atomdata.h +++ b/src/atomdata.h @@ -93,6 +93,7 @@ class AtomData double mulen = 0; //!< Dipole moment length bool hydrophobic = false; //!< Is the particle hydrophobic? bool implicit = false; //!< Is the particle implicit (e.g. proton)? + double scattering_f0 = 1.0; //!< Atomic scattering form factor (q-independent) InteractionData interaction; //!< Arbitrary interaction parameters, e.g., epsilons in various potentials SpheroCylinderData sphero_cylinder; //!< Data for patchy sphero cylinders (PSCs) diff --git a/src/scatter.cpp b/src/scatter.cpp index 1557da306..7db6e760c 100644 --- a/src/scatter.cpp +++ b/src/scatter.cpp @@ -9,8 +9,10 @@ namespace Faunus::Scatter { using doctest::Approx; -TEST_CASE_TEMPLATE("[Faunus] StructureFactorPBC", T, StructureFactorPBC, - StructureFactorPBC, StructureFactorPBC) +TEST_CASE_TEMPLATE("[Faunus] StructureFactorPBC", T, + StructureFactorPBC, float, SIMD>, + StructureFactorPBC, float, EIGEN>, + StructureFactorPBC, float, GENERIC>) { size_t cnt = 0; Point box = {80.0, 80.0, 80.0}; @@ -39,9 +41,9 @@ TEST_CASE("Benchmark") } ankerl::nanobench::Bench bench; bench.minEpochIterations(100); - bench.run("SIMD", [&] { StructureFactorPBC(10).sample(positions, box); }); - bench.run("EIGEN", [&] { StructureFactorPBC(10).sample(positions, box); }); - bench.run("GENERIC", [&] { StructureFactorPBC(10).sample(positions, box); }); + bench.run("SIMD", [&] { StructureFactorPBC, double, SIMD>(10).sample(positions, box); }); + bench.run("EIGEN", [&] { StructureFactorPBC, double, EIGEN>(10).sample(positions, box); }); + bench.run("GENERIC", [&] { StructureFactorPBC, double, GENERIC>(10).sample(positions, box); }); } #endif @@ -63,4 +65,29 @@ TEST_CASE("[Faunus] StructureFactorIPBC") CHECK_EQ(cnt, result.size()); } +TEST_CASE("[Faunus] FormFactorAtomicConstant") +{ + // Setup atom types with different scattering_f0 values + Faunus::atoms = R"([ + { "A": { "sigma": 2.0, "scattering_f0": 6.0 } }, + { "B": { "sigma": 3.0, "scattering_f0": 7.5 } } + ])"_json.get(); + + Scatterer s1{{0, 0, 0}, 0}; // type A + Scatterer s2{{1, 1, 1}, 1}; // type B + Scatterer s_unity{{2, 2, 2}, -1}; // special id for unity form factor + + FormFactorAtomicConstant ff; + CHECK_EQ(ff(0.1, s1), Approx(6.0)); // q-independent + CHECK_EQ(ff(0.5, s1), Approx(6.0)); // same for different q + CHECK_EQ(ff(0.1, s2), Approx(7.5)); // different atom type + CHECK_EQ(ff(0.5, s2), Approx(7.5)); // q-independent for type B + CHECK_EQ(ff(0.1, s_unity), Approx(1.0)); // id=-1 returns unity + + // Verify FormFactorUnity still returns 1 + FormFactorUnity ff_unity; + CHECK_EQ(ff_unity(0.1, s1), Approx(1.0)); + CHECK_EQ(ff_unity(0.1, s2), Approx(1.0)); +} + } // namespace Faunus::Scatter diff --git a/src/scatter.h b/src/scatter.h index add1b8b3f..354889e64 100644 --- a/src/scatter.h +++ b/src/scatter.h @@ -3,6 +3,7 @@ #include #include #include +#include "atomdata.h" namespace Faunus { @@ -49,7 +50,57 @@ template class FormFactorSphere */ template struct FormFactorUnity { - template T operator()(T, const Tparticle&) const { return 1; } + template [[nodiscard]] constexpr T operator()(T, const Tparticle&) const + { + return T{1}; + } +}; + +/** + * @brief Lightweight scatterer holding position and atom type id. + * + * Used for scattering calculations where both position and atom-specific + * form factors are needed. + */ +struct Scatterer +{ + Point pos; + int id = 0; //!< Atom type id for looking up form factor +}; + +/** + * @brief Helper to extract position from a scatterer or point. + * + * For Scatterer types, returns the `.pos` member. + * For Point types (Eigen vectors), returns the point directly. + */ +template +constexpr const auto& getPosition(const T& scatterer) +{ + if constexpr (requires { scatterer.pos; }) { + return scatterer.pos; + } + else { + return scatterer; + } +} + +/** + * @brief Atom-specific constant form factor (q independent). + * + * Returns the scattering_f0 property from AtomData for each scatterer. + * This allows atom-type specific scattering weights. + * For scatterers with id < 0, returns 1.0 (unity form factor). + */ +template struct FormFactorAtomicConstant +{ + template [[nodiscard]] T operator()(T, const Tscatterer& scatterer) const + { + if (scatterer.id < 0) { + return T{1}; + } + return static_cast(atoms.at(scatterer.id).scattering_f0); + } }; /** @@ -160,8 +211,8 @@ template class DebyeFormula #pragma omp for schedule(dynamic) for (int i = 0; i < N - 1; ++i) { for (int j = i + 1; j < N; ++j) { - T r = - T(Faunus::Geometry::Sphere::sqdist(p[i], p[j])); // the square root follows + T r = T(Faunus::Geometry::Sphere::sqdist( + getPosition(p[i]), getPosition(p[j]))); // the square root follows if (r < r_cutoff * r_cutoff) { r = std::sqrt(r); // Black magic: The q_mesh function must be inlineable otherwise the loop @@ -205,7 +256,8 @@ template class DebyeFormula } sampling[m] += weight; intensity[m] += - ((2 * intensity_sum[m] + intensity_self_sum) / N + intensity_corr) * weight; + ((2 * intensity_sum[m] + intensity_self_sum) / intensity_self_sum + intensity_corr) * + weight; } } @@ -275,19 +327,19 @@ template class SamplingPolicy }; /** - * @brief Calculate structure factor using explicit q averaging. + * @brief Calculate scattering intensity using explicit q averaging. * * This averages over the thirteen permutations of the Miller index [100], [110], [101] using: * - * @f[ S(\mathbf{q}) = \frac{1}{N} \left < - * \left ( \sum_i^N \sin(\mathbf{qr}_i) \right )^2 + - * \left ( \sum_j^N \cos(\mathbf{qr}_j) \right )^2 + * @f[ I(\mathbf{q}) = \frac{1}{\sum_i f_i^2} \left < + * \left ( \sum_i^N f_i \sin(\mathbf{qr}_i) \right )^2 + + * \left ( \sum_j^N f_j \cos(\mathbf{qr}_j) \right )^2 * \right > * @f] * * For more information, see @see http://doi.org/d8zgw5 and @see http://doi.org/10.1063/1.449987. */ -template , typename T = double, Algorithm method = SIMD, typename TSamplingPolicy = SamplingPolicy> class StructureFactorPBC : private TSamplingPolicy { @@ -299,6 +351,7 @@ class StructureFactorPBC : private TSamplingPolicy }; const int p_max; //!< multiples of q to be sampled + Tformfactor form_factor; //!< form factor functor using TSamplingPolicy::addSampling; public: @@ -310,57 +363,52 @@ class StructureFactorPBC : private TSamplingPolicy /** * https://gcc.gnu.org/gcc-9/porting_to.html#ompdatasharing * #pragma omp parallel for collapse(2) default(none) shared(directions, p_max, boxlength) - * shared(positions) + * shared(scatterers) */ - template - void sample(const Tpositions& positions, const Point& boxlength) + template + void sample(const Tscatterers& scatterers, const Point& boxlength) { #pragma omp parallel for collapse(2) default(shared) for (size_t i = 0; i < directions.size(); ++i) { // openmp req. tradional loop for (int p = 1; p <= p_max; ++p) { // loop over multiples of q const Point q = 2.0 * pc::pi * p * directions[i].cwiseQuotient(boxlength); // scattering vector - const auto s_of_q = calculateStructureFactor(positions, q); + const auto intensity = calculateIntensity(scatterers, q); #pragma omp critical // avoid race conditions when updating the map - addSampling(q.norm(), s_of_q, 1.0); + addSampling(q.norm(), intensity, 1.0); } } } - template - T calculateStructureFactor(const Tpositions& positions, const Point& q) const + template + T calculateIntensity(const Tscatterers& scatterers, const Point& q) const { T sum_cos = 0.0; T sum_sin = 0.0; - if constexpr (method == SIMD) { - // When sine and cosine is computed in separate loops, sine and cosine SIMD - // instructions may be used to get at least 4 times performance boost. - // Note January 2020: only GCC exploits this using libmvec library if --ffast-math is - // enabled. - auto dot_product = [q](const auto& pos) { return static_cast(q.dot(pos)); }; - auto qdotr = positions | std::views::transform(dot_product) | ranges::to; - std::for_each(qdotr.begin(), qdotr.end(), [&](auto qr) { sum_cos += cos(qr); }); - std::for_each(qdotr.begin(), qdotr.end(), [&](auto qr) { sum_sin += sin(qr); }); - } - else if constexpr (method == EIGEN) { - // Map is a Nx3 matrix facade into original positions (std::vector) - using namespace Eigen; - static_assert(std::is_same_v>); - auto qdotr = - (Map>((double*)positions.data(), positions.size(), 3) * q) - .array() - .eval(); - sum_cos = qdotr.cast().cos().sum(); - sum_sin = qdotr.cast().sin().sum(); + T sum_f_squared = 0.0; + if constexpr (method == SIMD || method == EIGEN) { + // For SIMD/EIGEN with form factors, fall back to GENERIC approach + // as vectorization with varying weights is more complex + for (const auto& scatterer : scatterers) { + const auto& pos = getPosition(scatterer); + const auto f = form_factor(q.norm(), scatterer); + const auto qr = static_cast(q.dot(pos)); + sum_cos += f * cos(qr); + sum_sin += f * sin(qr); + sum_f_squared += f * f; + } } else if constexpr (method == GENERIC) { - for (const auto& r : positions) { - const auto qr = static_cast(q.dot(r)); - sum_cos += cos(qr); // sine and cosine in same loop obstructs - sum_sin += sin(qr); // vectorization on most compilers... + for (const auto& scatterer : scatterers) { + const auto& pos = getPosition(scatterer); + const auto f = form_factor(q.norm(), scatterer); + const auto qr = static_cast(q.dot(pos)); + sum_cos += f * cos(qr); + sum_sin += f * sin(qr); + sum_f_squared += f * f; } - }; - return std::norm(std::complex(sum_cos, sum_sin)) / static_cast(positions.size()); + } + return std::norm(std::complex(sum_cos, sum_sin)) / sum_f_squared; } int getQMultiplier() { return p_max; } @@ -369,13 +417,14 @@ class StructureFactorPBC : private TSamplingPolicy }; /** - * @brief Calculate structure factor using explicit q averaging in isotropic periodic boundary + * @brief Calculate scattering intensity using explicit q averaging in isotropic periodic boundary * conditions (IPBC). * * The sample directions reduce to 3 compared to 13 in regular periodic boundary conditions. Overall * simplification shall yield roughly 10 times faster computation. */ -template > +template , std::floating_point T = float, + typename TSamplingPolicy = SamplingPolicy> class StructureFactorIPBC : private TSamplingPolicy { //! Sample directions (h,k,l). @@ -383,6 +432,7 @@ class StructureFactorIPBC : private TSamplingPolicy std::vector directions = {{1, 0, 0}, {1, 1, 0}, {1, 1, 1}}; int p_max; //!< multiples of q to be sampled + Tformfactor form_factor; //!< form factor functor using TSamplingPolicy::addSampling; public: @@ -391,19 +441,22 @@ class StructureFactorIPBC : private TSamplingPolicy { } - template - void sample(const Tpositions& positions, const Point& boxlength) + template + void sample(const Tscatterers& scatterers, const Point& boxlength) { // https://gcc.gnu.org/gcc-9/porting_to.html#ompdatasharing -// #pragma omp parallel for collapse(2) default(none) shared(directions, p_max, positions, +// #pragma omp parallel for collapse(2) default(none) shared(directions, p_max, scatterers, // boxlength) #pragma omp parallel for collapse(2) default(shared) for (size_t i = 0; i < directions.size(); ++i) { for (int p = 1; p <= p_max; ++p) { // loop over multiples of q const Point q = 2.0 * pc::pi * p * directions[i].cwiseQuotient(boxlength); // scattering vector - T sum_cos = 0; - for (auto& r : positions) { // loop over positions + T sum_f_cos = 0; + T sum_f_squared = 0; + for (const auto& scatterer : scatterers) { + const auto& r = getPosition(scatterer); + const auto f = form_factor(q.norm(), scatterer); // if q[i] == 0 then its cosine == 1 hence we can avoid cosine computation for // performance reasons T product = std::cos(T(q[0] * r[0])); @@ -411,15 +464,16 @@ class StructureFactorIPBC : private TSamplingPolicy product *= std::cos(T(q[1] * r[1])); if (q[2] != 0) product *= std::cos(T(q[2] * r[2])); - sum_cos += product; + sum_f_cos += f * product; + sum_f_squared += f * f; } // collect average, `norm()` gives the scattering vector length const T ipbc_factor = std::pow(2, directions[i].count()); // 2 ^ number of non-zero elements - const T sf = (sum_cos * sum_cos) / (float)(positions.size()) * ipbc_factor; + const T intensity = (sum_f_cos * sum_f_cos) / sum_f_squared * ipbc_factor; #pragma omp critical // avoid race conditions when updating the map - addSampling(q.norm(), sf, 1.0); + addSampling(q.norm(), intensity, 1.0); } } }