Refactor and test Structure Factor

src/core/statistics.cpp

@@ -39,8 +39,10 @@
 #include <utils/contains.hpp>
 #include <utils/math/sqr.hpp>
 
+#include <cassert>
 #include <cstdlib>
 #include <limits>
+#include <stdexcept>
 
 /****************************************************************************************
  *                                 basic observables calculation
@@ -263,84 +265,64 @@ void calc_part_distribution(PartCfg &partCfg, std::vector<int> const &p1_types,
     dist[i] /= (double)cnt;
 }
 
-std::vector<double> calc_structurefactor(PartCfg &partCfg,
-                                         std::vector<int> const &p_types,
-                                         int order) {
-  auto const order2 = order * order;
-  std::vector<double> ff;
-  ff.resize(2 * order2);
-  ff[2 * order2] = 0;
+void calc_structurefactor(PartCfg &partCfg, std::vector<int> const &p_types,
+                          int order, std::vector<double> &wavevectors,
+                          std::vector<double> &intensities) {
+
+  if (order < 1)
+    throw std::domain_error("order has to be a strictly positive number");
+
+  auto const order_sq = order * order;
+  std::vector<double> ff(2 * order_sq + 1);
   auto const twoPI_L = 2 * Utils::pi() / box_geo.length()[0];
 
-  if (order < 1) {
-    fprintf(stderr,
-            "WARNING: parameter \"order\" has to be a whole positive number");
-    fflush(nullptr);
-    errexit();
-  } else {
-    for (int qi = 0; qi < 2 * order2; qi++) {
-      ff[qi] = 0.0;
-    }
-    for (int i = 0; i <= order; i++) {
-      for (int j = -order; j <= order; j++) {
-        for (int k = -order; k <= order; k++) {
-          auto const n = i * i + j * j + k * k;
-          if ((n <= order2) && (n >= 1)) {
-            double C_sum = 0.0, S_sum = 0.0;
-            for (auto const &p : partCfg) {
-              for (int t : p_types) {
-                if (p.p.type == t) {
-                  auto const qr =
-                      twoPI_L * (Utils::Vector3i{{i, j, k}} * p.r.p);
-                  C_sum += cos(qr);
-                  S_sum += sin(qr);
-                }
+  for (int i = 0; i <= order; i++) {
+    for (int j = -order; j <= order; j++) {
+      for (int k = -order; k <= order; k++) {
+        auto const n = i * i + j * j + k * k;
+        if ((n <= order_sq) && (n >= 1)) {
+          double C_sum = 0.0, S_sum = 0.0;
+          for (auto const &p : partCfg) {
+            for (int t : p_types) {
+              if (p.p.type == t) {
+                auto const qr = twoPI_L * (Utils::Vector3i{{i, j, k}} * p.r.p);
+                C_sum += cos(qr);
+                S_sum += sin(qr);
               }
             }
-            ff[2 * n - 2] += C_sum * C_sum + S_sum * S_sum;
-            ff[2 * n - 1]++;
           }
+          ff[2 * n - 2] += C_sum * C_sum + S_sum * S_sum;
+          ff[2 * n - 1]++;
         }
       }
     }
-    int n = 0;
-    for (auto const &p : partCfg) {
-      for (int t : p_types) {
-        if (p.p.type == t)
-          n++;
-      }
+  }
+
+  long n_particles = 0l;
+  for (auto const &p : partCfg) {
+    for (int t : p_types) {
+      if (p.p.type == t)
+        n_particles++;
     }
-    for (int qi = 0; qi < order2; qi++)
-      if (ff[2 * qi + 1] != 0)
-        ff[2 * qi] /= n * ff[2 * qi + 1];
   }
-  return ff;
-}
 
-std::vector<std::vector<double>> modify_stucturefactor(int order,
-                                                       double const *sf) {
   int length = 0;
-
-  for (int i = 0; i < order * order; i++) {
-    if (sf[2 * i + 1] > 0) {
+  for (int qi = 0; qi < order_sq; qi++) {
+    if (ff[2 * qi + 1] != 0) {
+      ff[2 * qi] /= static_cast<double>(n_particles) * ff[2 * qi + 1];
       length++;
     }
   }
 
-  auto const qfak = 2.0 * Utils::pi() / box_geo.length()[0];
-  std::vector<double> intern;
-  intern.assign(2, 0.0);
-  std::vector<std::vector<double>> structure_factor;
-  structure_factor.assign(length, intern);
+  wavevectors.resize(length);
+  intensities.resize(length);
 
   int cnt = 0;
-  for (int i = 0; i < order * order; i++) {
-    if (sf[2 * i + 1] > 0) {
-      structure_factor[cnt][0] = qfak * sqrt(i + 1);
-      structure_factor[cnt][1] = sf[2 * i];
+  for (int i = 0; i < order_sq; i++) {
+    if (ff[2 * i + 1] != 0) {
+      wavevectors[cnt] = twoPI_L * sqrt(i + 1);
+      intensities[cnt] = ff[2 * i];
       cnt++;
     }
   }
-
-  return structure_factor;
 }

src/core/statistics.hpp

@@ -98,16 +98,15 @@ void calc_part_distribution(PartCfg &partCfg, std::vector<int> const &p1_types,
  *  and sf[1]=1. For q=7, there are no possible wave vectors, so
  *  sf[2*(7-1)]=sf[2*(7-1)+1]=0.
  *
- *  @param partCfg   particle collection
- *  @param p_types   list with types of particles to be analyzed
- *  @param order     the maximum wave vector length in 2PI/L
+ *  @param[in]  partCfg   particle collection
+ *  @param[in]  p_types   list with types of particles to be analyzed
+ *  @param[in]  order     the maximum wave vector length in units of 2PI/L
+ *  @param[out] wavevectors  the scattering vectors q
+ *  @param[out] intensities  the structure factor S(q)
  */
-std::vector<double> calc_structurefactor(PartCfg &partCfg,
-                                         std::vector<int> const &p_types,
-                                         int order);
-
-std::vector<std::vector<double>> modify_stucturefactor(int order,
-                                                       double const *sf);
+void calc_structurefactor(PartCfg &partCfg, std::vector<int> const &p_types,
+                          int order, std::vector<double> &wavevectors,
+                          std::vector<double> &intensities);
 
 /** Calculate the center of mass of a special type of the current configuration.
  *  \param part_type  type of the particle

src/python/espressomd/analyze.pxd

@@ -61,8 +61,7 @@ cdef extern from "Observable_stat.hpp":
         size_t chunk_size()
 
 cdef extern from "statistics.hpp":
-    cdef vector[double] calc_structurefactor(PartCfg & , const vector[int] & p_types, int order)
-    cdef vector[vector[double]] modify_stucturefactor(int order, double * sf)
+    cdef void calc_structurefactor(PartCfg & , const vector[int] & p_types, int order, vector[double] & wavevectors, vector[double] & intensities) except +
     cdef double mindist(PartCfg & , const vector[int] & set1, const vector[int] & set2)
     cdef vector[int] nbhood(PartCfg & , const Vector3d & pos, double r_catch, const Vector3i & planedims)
     cdef vector[double] calc_linear_momentum(int include_particles, int include_lbfluid)

src/python/espressomd/analyze.pyx

@@ -629,15 +629,17 @@ class Analysis:
 
         """
 
-        if (sf_types is None) or (not hasattr(sf_types, '__iter__')):
+        if sf_types is None or not hasattr(sf_types, '__iter__'):
             raise ValueError("sf_types has to be a list!")
         check_type_or_throw_except(
             sf_order, 1, int, "sf_order has to be an int!")
 
-        sf = analyze.calc_structurefactor(
-            analyze.partCfg(), sf_types, sf_order)
+        cdef vector[double] wavevectors
+        cdef vector[double] intensities
+        analyze.calc_structurefactor(
+            analyze.partCfg(), sf_types, sf_order, wavevectors, intensities)
 
-        return np.transpose(analyze.modify_stucturefactor(sf_order, sf.data()))
+        return np.vstack([wavevectors, intensities])
 
     #
     # distribution

testsuite/python/CMakeLists.txt

@@ -156,6 +156,7 @@ python_test(FILE hat.py MAX_NUM_PROC 4)
 python_test(FILE analyze_energy.py MAX_NUM_PROC 2)
 python_test(FILE analyze_mass_related.py MAX_NUM_PROC 4)
 python_test(FILE rdf.py MAX_NUM_PROC 1)
+python_test(FILE sf_simple_lattice.py MAX_NUM_PROC 1)
 python_test(FILE coulomb_mixed_periodicity.py MAX_NUM_PROC 4)
 python_test(FILE coulomb_cloud_wall_duplicated.py MAX_NUM_PROC 4 LABELS gpu)
 python_test(FILE collision_detection.py MAX_NUM_PROC 4)

testsuite/python/sf_simple_lattice.py

@@ -0,0 +1,181 @@
+#
+# Copyright (C) 2017-2021 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published byss
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+import unittest as ut
+import espressomd
+import numpy as np
+import itertools
+
+
+class StructureFactorTest(ut.TestCase):
+    '''
+    Test structure factor analysis against rectangular lattices.
+    We do not check the wavevectors directly, but rather the
+    corresponding SF order, which is more readable (integer value).
+    '''
+
+    box_l = 16
+    part_ty = 0
+    sf_order = 16
+    system = espressomd.System(box_l=[box_l, box_l, box_l])
+
+    def tearDown(self):
+        self.system.part.clear()
+
+    def order(self, wavevectors, a):
+        """
+        Square and rescale wavevectors to recover the corresponding
+        SF order, which is an integer between 0 and ``self.sf_order``.
+        """
+        return (wavevectors * a / (2 * np.pi))**2
+
+    def generate_peaks(self, a, b, c, conditions):
+        '''
+        Generate the main diffraction peaks for crystal structures.
+
+        Parameters
+        ----------
+        a: :obj:`float`
+            Length of the unit cell on the x-axis.
+        b: :obj:`float`
+            Length of the unit cell on the y-axis.
+        c: :obj:`float`
+            Length of the unit cell on the z-axis.
+        conditions: :obj:`function`
+            Reflection conditions for the crystal lattice.
+        '''
+        hkl_ranges = [
+            range(0, self.sf_order + 1),
+            range(-self.sf_order, self.sf_order + 1),
+            range(-self.sf_order, self.sf_order + 1),
+        ]
+        reflections = [np.linalg.norm([h / a, k / b, l / c]) * 2 * np.pi
+                       for (h, k, l) in itertools.product(*hkl_ranges)
+                       if conditions(h, k, l) and (h + k + l != 0) and
+                       (h**2 + k**2 + l**2) <= self.sf_order**2]
+        return np.unique(reflections)
+
+    def test_tetragonal(self):
+        """Check tetragonal lattice."""
+        a = 2
+        b = 4
+        c = 8
+        xen = range(0, self.box_l, a)
+        yen = range(0, self.box_l, b)
+        zen = range(0, self.box_l, c)
+        for i, j, k in itertools.product(xen, yen, zen):
+            self.system.part.add(type=self.part_ty, pos=(i, j, k))
+        wavevectors, intensities = self.system.analysis.structure_factor(
+            sf_types=[self.part_ty], sf_order=self.sf_order)
+        intensities = np.around(intensities, 8)
+        # no reflection conditions on (h,k,l)
+        peaks_ref = self.generate_peaks(a, b, c, lambda h, k, l: True)
+        peaks = wavevectors[np.nonzero(intensities)]
+        np.testing.assert_array_almost_equal(
+            self.order(peaks, a), self.order(peaks_ref[:len(peaks)], a))
+
+    def test_sc(self):
+        """Check simple cubic lattice."""
+        l0 = 4
+        xen = range(0, self.box_l, l0)
+        for i, j, k in itertools.product(xen, repeat=3):
+            self.system.part.add(type=self.part_ty, pos=(i, j, k))
+        wavevectors, intensities = self.system.analysis.structure_factor(
+            sf_types=[self.part_ty], sf_order=self.sf_order)
+        intensities = np.around(intensities, 8)
+        np.testing.assert_array_equal(
+            intensities[np.nonzero(intensities)], len(self.system.part))
+        # no reflection conditions on (h,k,l)
+        peaks = wavevectors[np.nonzero(intensities)]
+        peaks_ref = self.generate_peaks(l0, l0, l0, lambda h, k, l: True)
+        np.testing.assert_array_almost_equal(
+            self.order(peaks, l0), self.order(peaks_ref[:len(peaks)], l0))
+
+    def test_bcc(self):
+        """Check body-centered cubic lattice."""
+        l0 = 4
+        m = l0 / 2
+        xen = range(0, self.box_l, l0)
+        for i, j, k in itertools.product(xen, repeat=3):
+            self.system.part.add(type=self.part_ty, pos=(i, j, k))
+            self.system.part.add(type=self.part_ty, pos=(i + m, j + m, k + m))
+        wavevectors, intensities = self.system.analysis.structure_factor(
+            sf_types=[self.part_ty], sf_order=self.sf_order)
+        intensities = np.around(intensities, 8)
+        np.testing.assert_array_equal(
+            intensities[np.nonzero(intensities)], len(self.system.part))
+        # reflection conditions
+        # (h+k+l) even => F = 2f, otherwise F = 0
+        peaks_ref = self.generate_peaks(
+            l0, l0, l0, lambda h, k, l: (h + k + l) % 2 == 0)
+        peaks = wavevectors[np.nonzero(intensities)]
+        np.testing.assert_array_almost_equal(
+            self.order(peaks, l0), self.order(peaks_ref[:len(peaks)], l0))
+
+    def test_fcc(self):
+        """Check face-centered cubic lattice."""
+        l0 = 4
+        m = l0 / 2
+        xen = range(0, self.box_l, l0)
+        for i, j, k in itertools.product(xen, repeat=3):
+            self.system.part.add(type=self.part_ty, pos=(i, j, k))
+            self.system.part.add(type=self.part_ty, pos=(i + m, j + m, k))
+            self.system.part.add(type=self.part_ty, pos=(i + m, j, k + m))
+            self.system.part.add(type=self.part_ty, pos=(i, j + m, k + m))
+        wavevectors, intensities = self.system.analysis.structure_factor(
+            sf_types=[self.part_ty], sf_order=self.sf_order)
+        intensities = np.around(intensities, 8)
+        np.testing.assert_array_equal(
+            intensities[np.nonzero(intensities)], len(self.system.part))
+        # reflection conditions
+        # (h,k,l) all even or odd => F = 4f, otherwise F = 0
+        peaks_ref = self.generate_peaks(
+            l0, l0, l0, lambda h, k, l:
+            h % 2 == 0 and k % 2 == 0 and l % 2 == 0 or
+            h % 2 == 1 and k % 2 == 1 and l % 2 == 1)
+        peaks = wavevectors[np.nonzero(intensities)]
+        np.testing.assert_array_almost_equal(
+            self.order(peaks, l0), self.order(peaks_ref[:len(peaks)], l0))
+
+    def test_cco(self):
+        """Check c-centered orthorhombic lattice."""
+        l0 = 4
+        m = l0 / 2
+        xen = range(0, self.box_l, l0)
+        for i, j, k in itertools.product(xen, repeat=3):
+            self.system.part.add(type=self.part_ty, pos=(i, j, k))
+            self.system.part.add(type=self.part_ty, pos=(i + m, j + m, k))
+        wavevectors, intensities = self.system.analysis.structure_factor(
+            sf_types=[self.part_ty], sf_order=self.sf_order)
+        intensities = np.around(intensities, 8)
+        # reflection conditions
+        # (h+k) even => F = 2f, otherwise F = 0
+        peaks_ref = self.generate_peaks(
+            l0, l0, l0, lambda h, k, l: (h + k) % 2 == 0)
+        peaks = wavevectors[np.nonzero(intensities)]
+        np.testing.assert_array_almost_equal(
+            self.order(peaks, l0), self.order(peaks_ref[:len(peaks)], l0))
+
+    def test_exceptions(self):
+        with self.assertRaisesRegex(ValueError, 'order has to be a strictly positive number'):
+            self.system.analysis.structure_factor(sf_types=[0], sf_order=0)
+
+
+if __name__ == "__main__":
+    ut.main()