From 37351cb8acaf3353e08409343e9b6049c744d7a9 Mon Sep 17 00:00:00 2001 From: stekajack Date: Fri, 26 Mar 2021 15:11:07 +0100 Subject: [PATCH 1/6] tests: Add testcase for structure factor --- testsuite/python/CMakeLists.txt | 1 + testsuite/python/sf_simple_lattice.py | 63 +++++++++++++++++++++++++++ 2 files changed, 64 insertions(+) create mode 100644 testsuite/python/sf_simple_lattice.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index a94ba79bc5d..59288b00a36 100644 --- a/testsuite/python/CMakeLists.txt +++ b/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) diff --git a/testsuite/python/sf_simple_lattice.py b/testsuite/python/sf_simple_lattice.py new file mode 100644 index 00000000000..f094a666a5a --- /dev/null +++ b/testsuite/python/sf_simple_lattice.py @@ -0,0 +1,63 @@ +# +# Copyright (C) 2017-2019 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 . +# + +import unittest as ut +import espressomd +import numpy as np +from itertools import product + + +class sfTest(ut.TestCase): + + box_l = 20 + part_ty = 0 + sf_order = 20 + s = espressomd.System(box_l=[box_l, box_l, box_l]) + + def make_lattice(self): + xen, yen, zen = [ + x for x in range(0, self.box_l, 1)], [ + x for x in range(0, self.box_l, 2)], [ + x for x in range(0, self.box_l, 4)] + for i, j, k in product(xen, yen, zen): + self.s.part.add(type=self.part_ty, pos=(i, j, k)) + + def peaks(self): + a, b, c = [0, 2 * np.pi], [0, np.pi], [0, np.pi / 2] + self.diags = [np.sqrt(a**2 + b**2 + c**2) + for (a, b, c) in product(a, b, c)] + self.diags.remove(0) + + def scatter_lattice(self): + sf_observable_x, sf_observable_y = self.s.analysis.structure_factor( + sf_types=[self.part_ty], sf_order=self.sf_order) + self.sf_data = [(x, y) + for (x, y) in zip(sf_observable_x, sf_observable_y)] + self.sf_data.sort(key=lambda tup: tup[1]) + + def test(self): + self.make_lattice() + self.peaks() + self.scatter_lattice() + self.assertTrue([np.any(np.isclose(element, self.sf_data, rtol=1e-02)) + for element in self.diags]) + + +if __name__ == "__main__": + ut.main() From 1c2ddf716efb9bd63252a8c859dde42ffe4d4944 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 26 Mar 2021 18:25:06 +0100 Subject: [PATCH 2/6] core: Fix undefined behavior --- src/core/statistics.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index 318ebff65b8..5b3b9c9570a 100644 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -267,9 +267,7 @@ std::vector calc_structurefactor(PartCfg &partCfg, std::vector const &p_types, int order) { auto const order2 = order * order; - std::vector ff; - ff.resize(2 * order2); - ff[2 * order2] = 0; + std::vector ff(2 * order2 + 1); auto const twoPI_L = 2 * Utils::pi() / box_geo.length()[0]; if (order < 1) { From 310c90faa62afc2d7ff740def0877e66f861c8db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 6 Apr 2021 17:55:52 +0100 Subject: [PATCH 3/6] tests: Test multiple tetragonal lattices Rewrite test to check multiple variants of the tetragonal lattice. Analytical solutions are derived from reflection conditions. --- testsuite/python/sf_simple_lattice.py | 170 +++++++++++++++++++++----- 1 file changed, 137 insertions(+), 33 deletions(-) diff --git a/testsuite/python/sf_simple_lattice.py b/testsuite/python/sf_simple_lattice.py index f094a666a5a..afc67e0d6b8 100644 --- a/testsuite/python/sf_simple_lattice.py +++ b/testsuite/python/sf_simple_lattice.py @@ -1,5 +1,5 @@ # -# Copyright (C) 2017-2019 The ESPResSo project +# Copyright (C) 2017-2021 The ESPResSo project # # This file is part of ESPResSo. # @@ -20,43 +20,147 @@ import unittest as ut import espressomd import numpy as np -from itertools import product +import itertools -class sfTest(ut.TestCase): +class StructureFactorTest(ut.TestCase): + ''' + Test structure factor analysis against rectangular lattices. + We do not check the wavevectors directly, but rather the + quantity (wavevectors * 2 / pi)^2, which is more readable. + ''' - box_l = 20 + box_l = 16 part_ty = 0 - sf_order = 20 - s = espressomd.System(box_l=[box_l, box_l, box_l]) - - def make_lattice(self): - xen, yen, zen = [ - x for x in range(0, self.box_l, 1)], [ - x for x in range(0, self.box_l, 2)], [ - x for x in range(0, self.box_l, 4)] - for i, j, k in product(xen, yen, zen): - self.s.part.add(type=self.part_ty, pos=(i, j, k)) - - def peaks(self): - a, b, c = [0, 2 * np.pi], [0, np.pi], [0, np.pi / 2] - self.diags = [np.sqrt(a**2 + b**2 + c**2) - for (a, b, c) in product(a, b, c)] - self.diags.remove(0) - - def scatter_lattice(self): - sf_observable_x, sf_observable_y = self.s.analysis.structure_factor( + sf_order = 16 + system = espressomd.System(box_l=[box_l, box_l, box_l]) + + def tearDown(self): + self.system.part.clear() + + def generate_peaks(self, a, b, c, condition, cutoff=box_l): + ''' + 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. + condition: :obj:`function` + Reflection conditions for the crystal lattice. + cutoff: :obj:`float` (optional) + Cutoff value for the generation of reflections. + ''' + reflections = [np.linalg.norm([h * a, k * b, l * c]) + for (h, k, l) in itertools.product(np.arange(0, 10), repeat=3) + if condition(h, k, l)] + return (np.array([x for x in sorted(set(reflections)) + if 0. < x < 1.001 * cutoff]) / 4) ** 2 + + def test_tetragonal(self): + """Check tetragonal lattice.""" + xen = range(0, self.box_l, 2) + yen = range(0, self.box_l, 4) + zen = range(0, self.box_l, 8) + 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_int = np.around(intensities).astype(int) + # check all reflections + # no condition on (h,k,l) + peaks_ref = self.generate_peaks(2, 4, 8, lambda h, k, l: True) + peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 + np.testing.assert_array_almost_equal(peaks, peaks_ref) + + def test_sc(self): + """Check simple cubic lattice.""" + xen = range(0, self.box_l, 4) + 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_int = np.around(intensities).astype(int) + np.testing.assert_array_almost_equal(intensities, intensities_int) + intensities_ref = np.zeros(intensities.shape) + intensities_ref[np.nonzero(intensities_int)] = len(self.system.part) + np.testing.assert_array_equal(intensities_int, intensities_ref) + # check all reflections + # no condition on (h,k,l) + peaks_ref = self.generate_peaks(4, 4, 4, lambda h, k, l: True) + peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 + np.testing.assert_array_almost_equal(peaks, peaks_ref) + + def test_bcc(self): + """Check body-centered cubic lattice.""" + xen = range(0, self.box_l, 4) + 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 + 2, j + 2, k + 2)) + wavevectors, intensities = self.system.analysis.structure_factor( + sf_types=[self.part_ty], sf_order=self.sf_order) + intensities_int = np.around(intensities).astype(int) + np.testing.assert_array_almost_equal(intensities, intensities_int) + intensities_ref = np.zeros(intensities.shape) + intensities_ref[np.nonzero(intensities_int)] = len(self.system.part) + np.testing.assert_array_equal(intensities_int, intensities_ref) + # check all reflections + # (h+k+l) even => F = 2f, otherwise F = 0 + peaks_ref = self.generate_peaks( + 4, 4, 4, lambda h, k, l: (h + k + l) % 2 == 0) + peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 + np.testing.assert_array_almost_equal(peaks, peaks_ref) + + def test_fcc(self): + """Check face-centered cubic lattice.""" + xen = range(0, self.box_l, 4) + 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 + 2, j + 2, k)) + self.system.part.add(type=self.part_ty, pos=(i + 2, j, k + 2)) + self.system.part.add(type=self.part_ty, pos=(i, j + 2, k + 2)) + wavevectors, intensities = self.system.analysis.structure_factor( + sf_types=[self.part_ty], sf_order=self.sf_order) + intensities_int = np.around(intensities).astype(int) + np.testing.assert_array_almost_equal(intensities, intensities_int) + intensities_ref = np.zeros(intensities.shape) + intensities_ref[np.nonzero(intensities_int)] = len(self.system.part) + np.testing.assert_array_equal(intensities_int, intensities_ref) + # check all reflections + # (h,k,l) all even or odd => F = 4f, otherwise F = 0 + peaks_ref = self.generate_peaks( + 4, 4, 4, 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_int)] * 2 / np.pi)**2 + np.testing.assert_array_almost_equal(peaks, peaks_ref) + + def test_cco(self): + """Check c-centered orthorhombic lattice.""" + xen = range(0, self.box_l, 4) + 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 + 2, j + 2, k)) + wavevectors, intensities = self.system.analysis.structure_factor( sf_types=[self.part_ty], sf_order=self.sf_order) - self.sf_data = [(x, y) - for (x, y) in zip(sf_observable_x, sf_observable_y)] - self.sf_data.sort(key=lambda tup: tup[1]) - - def test(self): - self.make_lattice() - self.peaks() - self.scatter_lattice() - self.assertTrue([np.any(np.isclose(element, self.sf_data, rtol=1e-02)) - for element in self.diags]) + intensities_int = np.around(intensities).astype(int) + # check all reflections + peaks_ref = self.generate_peaks(4, 4, 4, lambda h, k, l: True) + peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 + np.testing.assert_array_almost_equal(peaks, peaks_ref) + # check main reflections + # (h+k) even => F = 2f, otherwise F = 0 + main_peaks_ref = 4 * self.generate_peaks( + 4, 4, 4, lambda h, k, l: (h + k) % 2 == 0, self.box_l / 2) + peaks = (wavevectors * 2 / np.pi)**2 + for main_peak in main_peaks_ref: + idx = (np.abs(peaks - main_peak)).argmin() + self.assertAlmostEqual(peaks[idx], main_peak) + self.assertAlmostEqual(intensities[idx], len(self.system.part)) if __name__ == "__main__": From 0eb09baad64438e29bde9069dc8d9c264ced7488 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 6 Apr 2021 18:44:43 +0200 Subject: [PATCH 4/6] core: Refactor structure factor code --- src/core/statistics.cpp | 90 +++++++++++++-------------- src/core/statistics.hpp | 4 +- src/python/espressomd/analyze.pxd | 4 +- src/python/espressomd/analyze.pyx | 4 +- testsuite/python/sf_simple_lattice.py | 4 ++ 5 files changed, 55 insertions(+), 51 deletions(-) diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index 5b3b9c9570a..b5f848cb413 100644 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -39,8 +39,10 @@ #include #include +#include #include #include +#include /**************************************************************************************** * basic observables calculation @@ -266,73 +268,71 @@ void calc_part_distribution(PartCfg &partCfg, std::vector const &p1_types, std::vector calc_structurefactor(PartCfg &partCfg, std::vector const &p_types, int order) { - auto const order2 = order * order; - std::vector ff(2 * order2 + 1); + if (order < 1) + throw std::domain_error("order has to be a strictly positive number"); + + auto const order_sq = order * order; + std::vector 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 qi = 0; qi < 2 * order_sq; 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 <= 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++; - } + } + int n = 0; + for (auto const &p : partCfg) { + for (int t : p_types) { + if (p.p.type == t) + n++; } - for (int qi = 0; qi < order2; qi++) - if (ff[2 * qi + 1] != 0) - ff[2 * qi] /= n * ff[2 * qi + 1]; } + for (int qi = 0; qi < order_sq; qi++) + if (ff[2 * qi + 1] != 0) + ff[2 * qi] /= n * ff[2 * qi + 1]; + return ff; } -std::vector> modify_stucturefactor(int order, - double const *sf) { - int length = 0; +std::vector> +modify_stucturefactor(int order, std::vector const &sf) { - for (int i = 0; i < order * order; i++) { + assert(order > 0); + + auto const order_sq = order * order; + + int length = 0; + for (int i = 0; i < order_sq; i++) { if (sf[2 * i + 1] > 0) { length++; } } auto const qfak = 2.0 * Utils::pi() / box_geo.length()[0]; - std::vector intern; - intern.assign(2, 0.0); - std::vector> structure_factor; - structure_factor.assign(length, intern); + std::vector> structure_factor(length, + std::vector(2)); int cnt = 0; - for (int i = 0; i < order * order; i++) { + for (int i = 0; i < order_sq; i++) { if (sf[2 * i + 1] > 0) { structure_factor[cnt][0] = qfak * sqrt(i + 1); structure_factor[cnt][1] = sf[2 * i]; diff --git a/src/core/statistics.hpp b/src/core/statistics.hpp index d24c9a50f5a..7ce5fd24477 100644 --- a/src/core/statistics.hpp +++ b/src/core/statistics.hpp @@ -106,8 +106,8 @@ std::vector calc_structurefactor(PartCfg &partCfg, std::vector const &p_types, int order); -std::vector> modify_stucturefactor(int order, - double const *sf); +std::vector> +modify_stucturefactor(int order, std::vector const &sf); /** Calculate the center of mass of a special type of the current configuration. * \param part_type type of the particle diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index 80966c4c022..42c1fb8c3a6 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -61,8 +61,8 @@ 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 vector[double] calc_structurefactor(PartCfg & , const vector[int] & p_types, int order) except + + cdef vector[vector[double]] modify_stucturefactor(int order, const vector[double] & sf) 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) diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 8dac5d94b6e..7f07c2f2b45 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -629,7 +629,7 @@ 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!") @@ -637,7 +637,7 @@ class Analysis: sf = analyze.calc_structurefactor( analyze.partCfg(), sf_types, sf_order) - return np.transpose(analyze.modify_stucturefactor(sf_order, sf.data())) + return np.transpose(analyze.modify_stucturefactor(sf_order, sf)) # # distribution diff --git a/testsuite/python/sf_simple_lattice.py b/testsuite/python/sf_simple_lattice.py index afc67e0d6b8..cd2574f399c 100644 --- a/testsuite/python/sf_simple_lattice.py +++ b/testsuite/python/sf_simple_lattice.py @@ -162,6 +162,10 @@ def test_cco(self): self.assertAlmostEqual(peaks[idx], main_peak) self.assertAlmostEqual(intensities[idx], len(self.system.part)) + 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() From a5d4223516f904327f92c7c46d236dd23a4183d7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 6 Apr 2021 19:16:02 +0200 Subject: [PATCH 5/6] core: Combine the two Structure Factor functions --- src/core/statistics.cpp | 46 ++++++++++--------------------- src/core/statistics.hpp | 17 ++++++------ src/python/espressomd/analyze.pxd | 3 +- src/python/espressomd/analyze.pyx | 8 ++++-- 4 files changed, 29 insertions(+), 45 deletions(-) diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp index b5f848cb413..54cc6d0c9e8 100644 --- a/src/core/statistics.cpp +++ b/src/core/statistics.cpp @@ -265,9 +265,10 @@ void calc_part_distribution(PartCfg &partCfg, std::vector const &p1_types, dist[i] /= (double)cnt; } -std::vector calc_structurefactor(PartCfg &partCfg, - std::vector const &p_types, - int order) { +void calc_structurefactor(PartCfg &partCfg, std::vector const &p_types, + int order, std::vector &wavevectors, + std::vector &intensities) { + if (order < 1) throw std::domain_error("order has to be a strictly positive number"); @@ -275,9 +276,6 @@ std::vector calc_structurefactor(PartCfg &partCfg, std::vector ff(2 * order_sq + 1); auto const twoPI_L = 2 * Utils::pi() / box_geo.length()[0]; - for (int qi = 0; qi < 2 * order_sq; 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++) { @@ -299,46 +297,32 @@ std::vector calc_structurefactor(PartCfg &partCfg, } } } - int n = 0; + + long n_particles = 0l; for (auto const &p : partCfg) { for (int t : p_types) { if (p.p.type == t) - n++; + n_particles++; } } - for (int qi = 0; qi < order_sq; qi++) - if (ff[2 * qi + 1] != 0) - ff[2 * qi] /= n * ff[2 * qi + 1]; - - return ff; -} - -std::vector> -modify_stucturefactor(int order, std::vector const &sf) { - - assert(order > 0); - - auto const order_sq = order * order; int length = 0; - for (int i = 0; i < order_sq; 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(n_particles) * ff[2 * qi + 1]; length++; } } - auto const qfak = 2.0 * Utils::pi() / box_geo.length()[0]; - std::vector> structure_factor(length, - std::vector(2)); + wavevectors.resize(length); + intensities.resize(length); int cnt = 0; for (int i = 0; i < order_sq; i++) { - if (sf[2 * i + 1] > 0) { - structure_factor[cnt][0] = qfak * sqrt(i + 1); - structure_factor[cnt][1] = sf[2 * i]; + if (ff[2 * i + 1] != 0) { + wavevectors[cnt] = twoPI_L * sqrt(i + 1); + intensities[cnt] = ff[2 * i]; cnt++; } } - - return structure_factor; } diff --git a/src/core/statistics.hpp b/src/core/statistics.hpp index 7ce5fd24477..5b7f95e2ffa 100644 --- a/src/core/statistics.hpp +++ b/src/core/statistics.hpp @@ -98,16 +98,15 @@ void calc_part_distribution(PartCfg &partCfg, std::vector 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 calc_structurefactor(PartCfg &partCfg, - std::vector const &p_types, - int order); - -std::vector> -modify_stucturefactor(int order, std::vector const &sf); +void calc_structurefactor(PartCfg &partCfg, std::vector const &p_types, + int order, std::vector &wavevectors, + std::vector &intensities); /** Calculate the center of mass of a special type of the current configuration. * \param part_type type of the particle diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index 42c1fb8c3a6..e1dc4317e03 100644 --- a/src/python/espressomd/analyze.pxd +++ b/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) except + - cdef vector[vector[double]] modify_stucturefactor(int order, const vector[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) diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 7f07c2f2b45..28b235575c3 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -634,10 +634,12 @@ class Analysis: 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)) + return np.vstack([wavevectors, intensities]) # # distribution From 1c3be47d15d81c98683b04487fd3ee519a9011b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 7 Apr 2021 18:30:46 +0200 Subject: [PATCH 6/6] tests: Make SF test independent of parameters Replace magic values by unit cell parameters. Fix a bug in the reflection peaks generator. The test is now independent of the box length, sf_order and unit cell parameters. --- testsuite/python/sf_simple_lattice.py | 148 ++++++++++++++------------ 1 file changed, 79 insertions(+), 69 deletions(-) diff --git a/testsuite/python/sf_simple_lattice.py b/testsuite/python/sf_simple_lattice.py index cd2574f399c..d1f010a8794 100644 --- a/testsuite/python/sf_simple_lattice.py +++ b/testsuite/python/sf_simple_lattice.py @@ -27,7 +27,7 @@ class StructureFactorTest(ut.TestCase): ''' Test structure factor analysis against rectangular lattices. We do not check the wavevectors directly, but rather the - quantity (wavevectors * 2 / pi)^2, which is more readable. + corresponding SF order, which is more readable (integer value). ''' box_l = 16 @@ -38,7 +38,14 @@ class StructureFactorTest(ut.TestCase): def tearDown(self): self.system.part.clear() - def generate_peaks(self, a, b, c, condition, cutoff=box_l): + 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. @@ -50,117 +57,120 @@ def generate_peaks(self, a, b, c, condition, cutoff=box_l): Length of the unit cell on the y-axis. c: :obj:`float` Length of the unit cell on the z-axis. - condition: :obj:`function` + conditions: :obj:`function` Reflection conditions for the crystal lattice. - cutoff: :obj:`float` (optional) - Cutoff value for the generation of reflections. ''' - reflections = [np.linalg.norm([h * a, k * b, l * c]) - for (h, k, l) in itertools.product(np.arange(0, 10), repeat=3) - if condition(h, k, l)] - return (np.array([x for x in sorted(set(reflections)) - if 0. < x < 1.001 * cutoff]) / 4) ** 2 + 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.""" - xen = range(0, self.box_l, 2) - yen = range(0, self.box_l, 4) - zen = range(0, self.box_l, 8) + 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_int = np.around(intensities).astype(int) - # check all reflections - # no condition on (h,k,l) - peaks_ref = self.generate_peaks(2, 4, 8, lambda h, k, l: True) - peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 - np.testing.assert_array_almost_equal(peaks, peaks_ref) + 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.""" - xen = range(0, self.box_l, 4) + 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_int = np.around(intensities).astype(int) - np.testing.assert_array_almost_equal(intensities, intensities_int) - intensities_ref = np.zeros(intensities.shape) - intensities_ref[np.nonzero(intensities_int)] = len(self.system.part) - np.testing.assert_array_equal(intensities_int, intensities_ref) - # check all reflections - # no condition on (h,k,l) - peaks_ref = self.generate_peaks(4, 4, 4, lambda h, k, l: True) - peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 - np.testing.assert_array_almost_equal(peaks, peaks_ref) + 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.""" - xen = range(0, self.box_l, 4) + 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 + 2, j + 2, k + 2)) + 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_int = np.around(intensities).astype(int) - np.testing.assert_array_almost_equal(intensities, intensities_int) - intensities_ref = np.zeros(intensities.shape) - intensities_ref[np.nonzero(intensities_int)] = len(self.system.part) - np.testing.assert_array_equal(intensities_int, intensities_ref) - # check all reflections + 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( - 4, 4, 4, lambda h, k, l: (h + k + l) % 2 == 0) - peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 - np.testing.assert_array_almost_equal(peaks, peaks_ref) + 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.""" - xen = range(0, self.box_l, 4) + 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 + 2, j + 2, k)) - self.system.part.add(type=self.part_ty, pos=(i + 2, j, k + 2)) - self.system.part.add(type=self.part_ty, pos=(i, j + 2, k + 2)) + 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_int = np.around(intensities).astype(int) - np.testing.assert_array_almost_equal(intensities, intensities_int) - intensities_ref = np.zeros(intensities.shape) - intensities_ref[np.nonzero(intensities_int)] = len(self.system.part) - np.testing.assert_array_equal(intensities_int, intensities_ref) - # check all reflections + 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( - 4, 4, 4, lambda h, k, l: + 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_int)] * 2 / np.pi)**2 - np.testing.assert_array_almost_equal(peaks, peaks_ref) + 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.""" - xen = range(0, self.box_l, 4) + 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 + 2, j + 2, 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_int = np.around(intensities).astype(int) - # check all reflections - peaks_ref = self.generate_peaks(4, 4, 4, lambda h, k, l: True) - peaks = (wavevectors[np.nonzero(intensities_int)] * 2 / np.pi)**2 - np.testing.assert_array_almost_equal(peaks, peaks_ref) - # check main reflections + intensities = np.around(intensities, 8) + # reflection conditions # (h+k) even => F = 2f, otherwise F = 0 - main_peaks_ref = 4 * self.generate_peaks( - 4, 4, 4, lambda h, k, l: (h + k) % 2 == 0, self.box_l / 2) - peaks = (wavevectors * 2 / np.pi)**2 - for main_peak in main_peaks_ref: - idx = (np.abs(peaks - main_peak)).argmin() - self.assertAlmostEqual(peaks[idx], main_peak) - self.assertAlmostEqual(intensities[idx], len(self.system.part)) + 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'):