Add atom-specific scattering form factors

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

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

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<Scatter::StructureFactorIPBC<>>(pmax);
+            explicit_average_ipbc =
+                std::make_unique<Scatter::StructureFactorIPBC<Tformfactor, float>>(pmax);
         }
         else {
             scheme = Schemes::EXPLICIT_PBC;
-            explicit_average_pbc = std::make_unique<Scatter::StructureFactorPBC<>>(pmax);
+            explicit_average_pbc = std::make_unique<Scatter::StructureFactorPBC<Tformfactor>>(pmax);
         }
     }
     else {

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<Point> scatter_positions; //!< vector of scattering points
+    std::vector<Scatter::Scatterer> scatter_positions; //!< vector of scatterers with position and atom id
     std::vector<MoleculeData::index_type> molecule_ids; //!< Molecule ids
     std::vector<std::string> molecule_names; //!< Molecule names corresponding to `molecule_ids`
-    using Tformfactor = Scatter::FormFactorUnity<double>;
+    using Tformfactor = Scatter::FormFactorAtomicConstant<double>;
 
     std::unique_ptr<Scatter::DebyeFormula<Tformfactor>> debye;
-    std::unique_ptr<Scatter::StructureFactorPBC<>> explicit_average_pbc;
-    std::unique_ptr<Scatter::StructureFactorIPBC<>> explicit_average_ipbc;
+    std::unique_ptr<Scatter::StructureFactorPBC<Tformfactor>> explicit_average_pbc;
+    std::unique_ptr<Scatter::StructureFactorIPBC<Tformfactor, float>> explicit_average_ipbc;
     void _sample() override;
     void _to_disk() override;
     void _to_json(json& j) const override;

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<double>() * 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));
 

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)

src/scatter.cpp

@@ -9,8 +9,10 @@ namespace Faunus::Scatter {
 
 using doctest::Approx;
 
-TEST_CASE_TEMPLATE("[Faunus] StructureFactorPBC", T, StructureFactorPBC<float, SIMD>,
-                   StructureFactorPBC<float, EIGEN>, StructureFactorPBC<float, GENERIC>)
+TEST_CASE_TEMPLATE("[Faunus] StructureFactorPBC", T,
+                   StructureFactorPBC<FormFactorUnity<float>, float, SIMD>,
+                   StructureFactorPBC<FormFactorUnity<float>, float, EIGEN>,
+                   StructureFactorPBC<FormFactorUnity<float>, 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<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); });
+    bench.run("SIMD", [&] { StructureFactorPBC<FormFactorUnity<double>, double, SIMD>(10).sample(positions, box); });
+    bench.run("EIGEN", [&] { StructureFactorPBC<FormFactorUnity<double>, double, EIGEN>(10).sample(positions, box); });
+    bench.run("GENERIC", [&] { StructureFactorPBC<FormFactorUnity<double>, 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<decltype(Faunus::atoms)>();
+
+    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<double> 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<double> 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