From 6328f651549efd72e6eabfe6c65e44a587562059 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 8 Dec 2020 19:11:54 +0100 Subject: [PATCH 01/14] tests: Fix broken ScaFaCoS test system.actors.clear() doesn't clear the DLC actor on the first pass. In that situation, if the DLC check fails, it exits early before the calls to `del system.actors[0]` and triggers the tearDown() method, which cannot remove DLC. The system is then in an undetermined state for the ScaFaCoS check, which throws an exception. --- testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py index 4d160e2b09f..56ed5071b03 100644 --- a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py +++ b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py @@ -47,6 +47,9 @@ class Dipolar_p3m_mdlc_p2nfft(ut.TestCase): def tearDown(self): self.system.part.clear() self.system.actors.clear() + if self.system.actors: + # workaround for DLC actor bug + self.system.actors.clear() def test_mdlc(self): s = self.system @@ -94,9 +97,6 @@ def test_mdlc(self): self.assertLessEqual(abs(err_t), tol_t, "Torque difference too large") self.assertLessEqual(abs(err_f), tol_f, "Force difference too large") - del s.actors[0] - del s.actors[0] - def test_p3m(self): s = self.system rho = 0.09 From 5ddc3f0178878566c8d59d1034f888dc269de618 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 9 Dec 2020 23:12:10 +0100 Subject: [PATCH 02/14] tests: Cleanup and remove code duplication --- .../dipolar_mdlc_p3m_scafacos_p2nfft.py | 58 ++++-------- testsuite/python/scafacos_dipoles_1d_2d.py | 92 +++++++++---------- 2 files changed, 62 insertions(+), 88 deletions(-) diff --git a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py index 56ed5071b03..2783350e1db 100644 --- a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py +++ b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py @@ -1,5 +1,4 @@ # Copyright (C) 2010-2019 The ESPResSo project - # # This file is part of ESPResSo. # @@ -42,7 +41,6 @@ class Dipolar_p3m_mdlc_p2nfft(ut.TestCase): system.time_step = 0.01 system.cell_system.skin = .4 system.periodicity = [1, 1, 1] - system.thermostat.turn_off() def tearDown(self): self.system.part.clear() @@ -51,6 +49,9 @@ def tearDown(self): # workaround for DLC actor bug self.system.actors.clear() + def vector_error(self, a, b): + return np.sum(np.linalg.norm(a - b, axis=1)) / np.sqrt(a.shape[0]) + def test_mdlc(self): s = self.system rho = 0.3 @@ -60,34 +61,25 @@ def test_mdlc(self): n_particle = 100 particle_radius = 0.5 - box_l = pow(((4 * n_particle * np.pi) / (3 * rho)), - 1.0 / 3.0) * particle_radius + box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius s.box_l = 3 * [box_l] - f = open(abspath("data/mdlc_reference_data_energy.dat")) - ref_E = float(f.readline()) - f.close() + ref_E_path = abspath("data/mdlc_reference_data_energy.dat") + ref_E = float(np.genfromtxt(ref_E_path)) # Particles data = np.genfromtxt( abspath("data/mdlc_reference_data_forces_torques.dat")) - for p in data[:, :]: - s.part.add(id=int(p[0]), pos=p[1:4], dip=p[4:7]) + s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) s.part[:].rotation = (1, 1, 1) p3m = magnetostatics.DipolarP3M(prefactor=1, mesh=32, accuracy=1E-4) dlc = magnetostatic_extensions.DLC(maxPWerror=1E-5, gap_size=2.) s.actors.add(p3m) s.actors.add(dlc) - s.thermostat.turn_off() s.integrator.run(0) - err_f = np.sum(np.linalg.norm( - s.part[:].f - data[:, 7:10], axis=1)) / np.sqrt(data.shape[0]) - err_t = np.sum(np.linalg.norm( - s.part[:].torque_lab - data[:, 10:13], axis=1)) / np.sqrt(data.shape[0]) + err_f = self.vector_error(s.part[:].f, data[:, 7:10]) + err_t = self.vector_error(s.part[:].torque_lab, data[:, 10:13]) err_e = s.analysis.energy()["dipolar"] - ref_E - print("Energy difference", err_e) - print("Force difference", err_f) - print("Torque difference", err_t) tol_f = 2E-3 tol_t = 2E-3 @@ -106,14 +98,12 @@ def test_p3m(self): n_particle = 1000 particle_radius = 1 - box_l = pow(((4 * n_particle * np.pi) / (3 * rho)), - 1.0 / 3.0) * particle_radius + box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius s.box_l = 3 * [box_l] # Particles data = np.genfromtxt(abspath("data/p3m_magnetostatics_system.data")) - for p in data[:, :]: - s.part.add(id=int(p[0]), pos=p[1:4], dip=p[4:7]) + s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) s.part[:].rotation = (1, 1, 1) p3m = magnetostatics.DipolarP3M( @@ -122,15 +112,10 @@ def test_p3m(self): s.integrator.run(0) expected = np.genfromtxt( abspath("data/p3m_magnetostatics_expected.data"))[:, 1:] - err_f = np.sum(np.linalg.norm( - s.part[:].f - expected[:, 0:3], axis=1)) / np.sqrt(data.shape[0]) - err_t = np.sum(np.linalg.norm( - s.part[:].torque_lab - expected[:, 3:6], axis=1)) / np.sqrt(data.shape[0]) + err_f = self.vector_error(s.part[:].f, expected[:, 0:3]) + err_t = self.vector_error(s.part[:].torque_lab, expected[:, 3:6]) ref_E = 5.570 err_e = s.analysis.energy()["dipolar"] - ref_E - print("Energy difference", err_e) - print("Force difference", err_f) - print("Torque difference", err_t) tol_f = 2E-3 tol_t = 2E-3 @@ -150,15 +135,13 @@ def test_scafacos_dipoles(self): n_particle = 1000 particle_radius = 1 - box_l = pow(((4 * n_particle * np.pi) / (3 * rho)), - 1.0 / 3.0) * particle_radius + box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius s.box_l = 3 * [box_l] # Particles data = np.genfromtxt(abspath("data/p3m_magnetostatics_system.data")) - for p in data[:, :]: - s.part.add(id=int(p[0]), pos=p[1:4], - dip=p[4:7], rotation=(1, 1, 1)) + s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) + s.part[:].rotation = (1, 1, 1) scafacos = magnetostatics.Scafacos( prefactor=1, @@ -177,15 +160,10 @@ def test_scafacos_dipoles(self): s.integrator.run(0) expected = np.genfromtxt( abspath("data/p3m_magnetostatics_expected.data"))[:, 1:] - err_f = np.sum(np.linalg.norm( - s.part[:].f - expected[:, 0:3], axis=1)) / np.sqrt(data.shape[0]) - err_t = np.sum(np.linalg.norm( - s.part[:].torque_lab - expected[:, 3:6], axis=1)) / np.sqrt(data.shape[0]) + err_f = self.vector_error(s.part[:].f, expected[:, 0:3]) + err_t = self.vector_error(s.part[:].torque_lab, expected[:, 3:6]) ref_E = 5.570 err_e = s.analysis.energy()["dipolar"] - ref_E - print("Energy difference", err_e) - print("Force difference", err_f) - print("Torque difference", err_t) tol_f = 2E-3 tol_t = 2E-3 diff --git a/testsuite/python/scafacos_dipoles_1d_2d.py b/testsuite/python/scafacos_dipoles_1d_2d.py index 4442236e607..fb9c2e9f5c5 100644 --- a/testsuite/python/scafacos_dipoles_1d_2d.py +++ b/testsuite/python/scafacos_dipoles_1d_2d.py @@ -31,7 +31,21 @@ @utx.skipIfMissingFeatures(["SCAFACOS_DIPOLES"]) class Scafacos1d2d(ut.TestCase): + system = espressomd.System(box_l=[1.0, 1.0, 1.0]) + system.time_step = 0.01 + system.cell_system.skin = 0.5 + system.periodicity = [1, 1, 1] + + def tearDown(self): + self.system.part.clear() + self.system.actors.clear() + self.system.periodicity = [1, 1, 1] + + def vector_error(self, a, b): + return np.sum(np.linalg.norm(a - b, axis=1)) / np.sqrt(a.shape[0]) + def test_scafacos(self): + s = self.system rho = 0.3 # This is only for box size calculation. The actual particle number is @@ -40,18 +54,10 @@ def test_scafacos(self): particle_radius = 0.5 - ################################################# - - box_l = pow(((4 * n_particle * np.pi) / (3 * rho)), - 1.0 / 3.0) * particle_radius - skin = 0.5 - - s = espressomd.System(box_l=[1.0, 1.0, 1.0]) - # give Espresso some parameters - s.time_step = 0.01 - s.cell_system.skin = skin + box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius s.box_l = 3 * [box_l] - for dim in 2, 1: + + for dim in (2, 1): print("Dimension", dim) # Read reference data @@ -62,15 +68,14 @@ def test_scafacos(self): s.periodicity = [1, 0, 0] file_prefix = "data/scafacos_dipoles_1d" - with open(abspath(file_prefix + "_reference_data_energy.dat")) as f: - ref_E = float(f.readline()) + ref_E_path = abspath(file_prefix + "_reference_data_energy.dat") + ref_E = float(np.genfromtxt(ref_E_path)) # Particles data = np.genfromtxt(abspath( file_prefix + "_reference_data_forces_torques.dat")) - for p in data[:, :]: - s.part.add( - id=int(p[0]), pos=p[1:4], dip=p[4:7], rotation=(1, 1, 1)) + s.part.add(pos=data[:, 1:4], dip=data[:, 4:7]) + s.part[:].rotation = (1, 1, 1) if dim == 2: scafacos = magnetostatics.Scafacos( @@ -92,42 +97,33 @@ def test_scafacos(self): s.box_l = np.array((1, 1, 1.3)) * box_l else: - if dim == 1: - # 1d periodic in x - scafacos = magnetostatics.Scafacos( - prefactor=1, - method_name="p2nfft", - method_params={ - "p2nfft_verbose_tuning": 1, - "pnfft_N": "32,128,128", - "pnfft_direct": 0, - "p2nfft_r_cut": 2.855, - "p2nfft_alpha": "1.5", - "p2nfft_intpol_order": "-1", - "p2nfft_reg_kernel_name": "ewald", - "p2nfft_p": "16", - "p2nfft_ignore_tolerance": "1", - "pnfft_window_name": "bspline", - "pnfft_m": "8", - "pnfft_diff_ik": "1", - "p2nfft_epsB": "0.125"}) - s.box_l = np.array((1, 1, 1)) * box_l - s.actors.add(scafacos) - else: - raise Exception("This shouldn't happen.") - s.thermostat.turn_off() + # 1d periodic in x + scafacos = magnetostatics.Scafacos( + prefactor=1, + method_name="p2nfft", + method_params={ + "p2nfft_verbose_tuning": 1, + "pnfft_N": "32,128,128", + "pnfft_direct": 0, + "p2nfft_r_cut": 2.855, + "p2nfft_alpha": "1.5", + "p2nfft_intpol_order": "-1", + "p2nfft_reg_kernel_name": "ewald", + "p2nfft_p": "16", + "p2nfft_ignore_tolerance": "1", + "pnfft_window_name": "bspline", + "pnfft_m": "8", + "pnfft_diff_ik": "1", + "p2nfft_epsB": "0.125"}) + s.box_l = np.array((1, 1, 1)) * box_l + s.actors.add(scafacos) s.integrator.run(0) # Calculate errors - err_f = np.sum(np.linalg.norm( - s.part[:].f - data[:, 7:10], axis=1)) / np.sqrt(data.shape[0]) - err_t = np.sum(np.linalg.norm( - s.part[:].torque_lab - data[:, 10:13], axis=1)) / np.sqrt(data.shape[0]) + err_f = self.vector_error(s.part[:].f, data[:, 7:10]) + err_t = self.vector_error(s.part[:].torque_lab, data[:, 10:13]) err_e = s.analysis.energy()["dipolar"] - ref_E - print("Energy difference", err_e) - print("Force difference", err_f) - print("Torque difference", err_t) tol_f = 2E-3 tol_t = 2E-3 @@ -141,7 +137,7 @@ def test_scafacos(self): abs(err_f), tol_f, "Force difference too large") s.part.clear() - del s.actors[0] + s.actors.clear() if __name__ == "__main__": From 67a6fcfd73ecd4c0ca7fa804cd2cc141b3f7d386 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 8 Dec 2020 19:45:44 +0100 Subject: [PATCH 03/14] core: Simplify ScaFaCoS interface --- .../coulomb_inline.hpp | 2 +- .../scafacos.cpp | 60 ++++++++----------- .../scafacos.hpp | 5 +- src/scafacos/Scafacos.cpp | 13 ++-- 4 files changed, 33 insertions(+), 47 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp index a9b58edad4b..8c629727e21 100644 --- a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp @@ -58,7 +58,7 @@ inline Utils::Vector3d central_force(double const q1q2, break; #ifdef SCAFACOS case COULOMB_SCAFACOS: - Scafacos::add_pair_force(q1q2, d.data(), dist, f.data()); + Scafacos::add_pair_force(q1q2, d, dist, f); break; #endif default: diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index 2c8fc72da1d..bc6bc188111 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -44,6 +44,7 @@ #include #include +#include #include #include @@ -51,7 +52,6 @@ #include #include #include -#include #include #include #include @@ -63,14 +63,18 @@ namespace Scafacos { -/** Get available scafacos methods */ +/** Get available ScaFaCoS methods */ std::list available_methods() { return Scafacos::available_methods(); } -/** Encapsulation for the particle data needed by scafacos */ +/** Encapsulation for the particle data needed by ScaFaCoS */ struct ScafacosData { - int update_particle_data(); + /** @brief Collect particle data in continuous arrays as required + * by ScaFaCoS. + */ + void update_particle_data(); + /** @brief Write forces back to particles. */ void update_particle_forces() const; /** Inputs */ std::vector charges, positions, dipoles; @@ -81,8 +85,7 @@ struct ScafacosData { static Scafacos *scafacos = nullptr; static ScafacosData particles; -/** \brief Collect particle data in continuous arrays as required by fcs */ -int ScafacosData::update_particle_data() { +void ScafacosData::update_particle_data() { positions.clear(); if (!dipolar()) { @@ -109,16 +112,13 @@ int ScafacosData::update_particle_data() { #endif } } - - return static_cast(positions.size() / 3); } -/** \brief Write forces back to particles */ void ScafacosData::update_particle_forces() const { - int it = 0; if (positions.empty()) return; + int it = 0; for (auto &p : cell_structure.local_particles()) { if (!dipolar()) { p.f.f += coulomb.prefactor * p.p.q * @@ -167,17 +167,15 @@ void ScafacosData::update_particle_forces() const { } } -void add_pair_force(double q1q2, const double *d, double dist, double *force) { +void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, + Utils::Vector3d &force) { if (dist > get_r_cut()) return; assert(scafacos); - const double field = scafacos->pair_force(dist); - const double fak = q1q2 * field / dist; - - for (int i = 0; i < 3; i++) { - force[i] -= fak * d[i]; - } + auto const field = scafacos->pair_force(dist); + auto const fak = q1q2 * field / dist; + force -= fak * d; } double pair_energy(double q1q2, double dist) { @@ -200,8 +198,7 @@ void add_long_range_force() { #endif } } else - throw std::runtime_error( - "Scafacos internal error. Instance pointer is not valid."); + throw std::runtime_error("Scafacos not initialized"); particles.update_particle_forces(); } @@ -213,17 +210,13 @@ double long_range_energy() { scafacos->run(particles.charges, particles.positions, particles.fields, particles.potentials); return 0.5 * coulomb.prefactor * - std::inner_product(particles.charges.begin(), - particles.charges.end(), - particles.potentials.begin(), 0.0); + boost::inner_product(particles.charges, particles.potentials, 0.0); } #ifdef SCAFACOS_DIPOLES scafacos->run_dipolar(particles.dipoles, particles.positions, particles.fields, particles.potentials); return -0.5 * dipole.prefactor * - std::inner_product(particles.dipoles.begin(), - particles.dipoles.end(), - particles.potentials.begin(), 0.0); + boost::inner_product(particles.dipoles, particles.potentials, 0.0); #endif } @@ -241,8 +234,7 @@ REGISTER_CALLBACK(set_r_cut_and_tune_local) /** Determine runtime for a specific cutoff */ double time_r_cut(double r_cut) { /* Set cutoff to time */ - mpi_call(set_r_cut_and_tune_local, r_cut); - set_r_cut_and_tune_local(r_cut); + mpi_call_all(set_r_cut_and_tune_local, r_cut); return time_force_calc(10); } @@ -301,10 +293,8 @@ void tune() { static void set_params_safe(const std::string &method, const std::string ¶ms, bool dipolar_ia, int n_part) { - if (scafacos) { - delete scafacos; - scafacos = nullptr; - } + delete scafacos; + scafacos = nullptr; scafacos = new Scafacos(method, comm_cart, params); @@ -376,10 +366,8 @@ void set_dipolar(bool d) { void free_handle() { if (this_node == 0) mpi_call(free_handle); - if (scafacos) { - delete scafacos; - scafacos = nullptr; - } + delete scafacos; + scafacos = nullptr; } REGISTER_CALLBACK(free_handle) @@ -387,7 +375,7 @@ REGISTER_CALLBACK(free_handle) void update_system_params() { // If scafacos is not active, do nothing if (!scafacos) { - throw std::runtime_error("Scafacos object not there"); + throw std::runtime_error("Scafacos not initialized"); } int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0, diff --git a/src/core/electrostatics_magnetostatics/scafacos.hpp b/src/core/electrostatics_magnetostatics/scafacos.hpp index a1272eeb64e..6a30fcbe61e 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.hpp +++ b/src/core/electrostatics_magnetostatics/scafacos.hpp @@ -28,13 +28,16 @@ #include "config.hpp" +#include + #include #include namespace Scafacos { #if defined(SCAFACOS) /** Near-field pair force */ -void add_pair_force(double q1q2, const double *d, double dist, double *force); +void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, + Utils::Vector3d &force); /** Near-field pair energy */ double pair_energy(double q1q2, double dist); /** Long range part */ diff --git a/src/scafacos/Scafacos.cpp b/src/scafacos/Scafacos.cpp index a9b2711c8d5..2bd1da3fafb 100644 --- a/src/scafacos/Scafacos.cpp +++ b/src/scafacos/Scafacos.cpp @@ -162,15 +162,10 @@ void Scafacos::set_common_parameters(const double *box_l, boxa[0] = box_l[0]; boxb[1] = box_l[1]; boxc[2] = box_l[2]; - // Does scafacos calculate the near field part - // For charges, if the method supports it, Es calculates near field - int sr = 0; - if (!dipolar() && has_near) { - sr = 0; - } else { - // Scafacos does near field calc - sr = 1; - } + // Does scafacos calculate the near field part? + // For charges, ESPResSo calculates near field if the method supports it. + // For dipoles, ScaFaCoS calculates near field. + int const sr = dipolar() || !has_near; handle_error(fcs_set_common(handle, sr, boxa, boxb, boxc, off, periodicity, total_particles)); #ifdef FCS_ENABLE_DIPOLES From b3641c4d1558ab12e1ec2e31251ab12a51805ddf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 8 Dec 2020 22:13:23 +0100 Subject: [PATCH 04/14] core: Split ScaFaCoS Coulomb/Dipoles --- .../scafacos.cpp | 249 +++++++++++------- src/python/espressomd/scafacos.pxd | 2 +- 2 files changed, 155 insertions(+), 96 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index bc6bc188111..5daba687789 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -70,102 +70,144 @@ std::list available_methods() { /** Encapsulation for the particle data needed by ScaFaCoS */ struct ScafacosData { + virtual ~ScafacosData() {} /** @brief Collect particle data in continuous arrays as required * by ScaFaCoS. */ - void update_particle_data(); + virtual void update_particle_data() = 0; /** @brief Write forces back to particles. */ - void update_particle_forces() const; - /** Inputs */ - std::vector charges, positions, dipoles; + virtual void update_particle_forces() const = 0; + virtual double long_range_energy() = 0; + virtual void add_long_range_force() = 0; + virtual void tune() = 0; + +protected: /** Outputs */ std::vector fields, potentials; }; -static Scafacos *scafacos = nullptr; -static ScafacosData particles; +struct ScafacosDataCoulomb : ScafacosData { + ScafacosDataCoulomb() {} + void update_particle_data() override; + void update_particle_forces() const override; + double long_range_energy() override; + void add_long_range_force() override; + void tune() override; -void ScafacosData::update_particle_data() { - positions.clear(); +private: + /** Inputs */ + std::vector positions, charges; +}; - if (!dipolar()) { - charges.clear(); - } else { #ifdef SCAFACOS_DIPOLES - dipoles.clear(); -#endif +struct ScafacosDataDipoles : ScafacosData { + ScafacosDataDipoles() {} + void update_particle_data() override; + void update_particle_forces() const override; + double long_range_energy() override; + void add_long_range_force() override; + void tune() override { + runtimeErrorMsg() << "Tuning unavailable for ScaFaCoS dipoles."; } +private: + /** Inputs */ + std::vector positions, dipoles; +}; +#endif + +static Scafacos *scafacos = nullptr; +static ScafacosData *scafacos_data = nullptr; + +void ScafacosDataCoulomb::tune() { scafacos->tune(charges, positions); } + +void ScafacosDataCoulomb::update_particle_data() { + positions.clear(); + charges.clear(); + for (auto const &p : cell_structure.local_particles()) { auto const pos = folded_position(p.r.p, box_geo); positions.push_back(pos[0]); positions.push_back(pos[1]); positions.push_back(pos[2]); - if (!dipolar()) { - charges.push_back(p.p.q); - } else { + charges.push_back(p.p.q); + } +} + #ifdef SCAFACOS_DIPOLES - auto const dip = p.calc_dip(); - dipoles.push_back(dip[0]); - dipoles.push_back(dip[1]); - dipoles.push_back(dip[2]); -#endif - } +void ScafacosDataDipoles::update_particle_data() { + positions.clear(); + dipoles.clear(); + + for (auto const &p : cell_structure.local_particles()) { + auto const pos = folded_position(p.r.p, box_geo); + positions.push_back(pos[0]); + positions.push_back(pos[1]); + positions.push_back(pos[2]); + auto const dip = p.calc_dip(); + dipoles.push_back(dip[0]); + dipoles.push_back(dip[1]); + dipoles.push_back(dip[2]); } } +#endif -void ScafacosData::update_particle_forces() const { +void ScafacosDataCoulomb::update_particle_forces() const { if (positions.empty()) return; int it = 0; for (auto &p : cell_structure.local_particles()) { - if (!dipolar()) { - p.f.f += coulomb.prefactor * p.p.q * - Utils::Vector3d(Utils::Span(&(fields[it]), 3)); - it += 3; - } else { -#ifdef SCAFACOS_DIPOLES - // Indices - // 3 "potential" values per particles (see below) - int const it_p = 3 * it; - // 6 "field" values per particles (see below) - int const it_f = 6 * it; - - // The scafacos term "potential" here in fact refers to the magnetic - // field - // So, the torques are given by m \times B - auto const dip = p.calc_dip(); - auto const t = vector_product( - dip, - Utils::Vector3d(Utils::Span(&(potentials[it_p]), 3))); - // The force is given by G m, where G is a matrix - // which comes from the "fields" output of scafacos like this - // 0 1 2 - // 1 3 4 - // 2 4 5 - // where the numbers refer to indices in the "field" output from scafacos - auto const G = Utils::Vector{ - {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]}, - {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]}, - {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}}; - auto const f = G * dip; - - // Add to particles - p.f.f += dipole.prefactor * f; - p.f.torque += dipole.prefactor * t; - it++; -#endif - } + p.f.f += coulomb.prefactor * p.p.q * + Utils::Vector3d(Utils::Span(&(fields[it]), 3)); + it += 3; } /* Check that the particle number did not change */ - if (!dipolar()) { - assert(it == fields.size()); - } else { - assert(it == positions.size() / 3); + assert(it == fields.size()); +} + +#ifdef SCAFACOS_DIPOLES +void ScafacosDataDipoles::update_particle_forces() const { + if (positions.empty()) + return; + + int it = 0; + for (auto &p : cell_structure.local_particles()) { + // Indices + // 3 "potential" values per particles (see below) + int const it_p = 3 * it; + // 6 "field" values per particles (see below) + int const it_f = 6 * it; + + // The scafacos term "potential" here in fact refers to the magnetic + // field. So, the torques are given by m \times B + auto const dip = p.calc_dip(); + auto const t = vector_product( + dip, + Utils::Vector3d(Utils::Span(&(potentials[it_p]), 3))); + // The force is given by G m, where G is a matrix + // which comes from the "fields" output of scafacos like this + // 0 1 2 + // 1 3 4 + // 2 4 5 + // where the numbers refer to indices in the "field" output from scafacos + auto const G = Utils::Vector{ + {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]}, + {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]}, + {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}}; + auto const f = G * dip; + + // Add to particles + p.f.f += dipole.prefactor * f; + p.f.torque += dipole.prefactor * t; + it++; } + + /* Check that the particle number did not change */ + assert(it == positions.size() / 3); } +#endif void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, Utils::Vector3d &force) { @@ -173,6 +215,7 @@ void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, return; assert(scafacos); + assert(scafacos_data); auto const field = scafacos->pair_force(dist); auto const fak = q1q2 * field / dist; force -= fak * d; @@ -184,49 +227,55 @@ double pair_energy(double q1q2, double dist) { return 0.; } -void add_long_range_force() { - particles.update_particle_data(); +void ScafacosDataCoulomb::add_long_range_force() { + update_particle_data(); + scafacos->run(charges, positions, fields, potentials); + update_particle_forces(); +} - if (scafacos) { - if (!dipolar()) { - scafacos->run(particles.charges, particles.positions, particles.fields, - particles.potentials); - } else { #ifdef SCAFACOS_DIPOLES - scafacos->run_dipolar(particles.dipoles, particles.positions, - particles.fields, particles.potentials); +void ScafacosDataDipoles::add_long_range_force() { + update_particle_data(); + scafacos->run_dipolar(dipoles, positions, fields, potentials); + update_particle_forces(); +} #endif - } - } else + +void add_long_range_force() { + if (!scafacos) { throw std::runtime_error("Scafacos not initialized"); + } + scafacos_data->add_long_range_force(); +} - particles.update_particle_forces(); +double ScafacosDataCoulomb::long_range_energy() { + update_particle_data(); + scafacos->run(charges, positions, fields, potentials); + return 0.5 * coulomb.prefactor * + boost::inner_product(charges, potentials, 0.0); } -double long_range_energy() { - if (scafacos) { - particles.update_particle_data(); - if (!dipolar()) { - scafacos->run(particles.charges, particles.positions, particles.fields, - particles.potentials); - return 0.5 * coulomb.prefactor * - boost::inner_product(particles.charges, particles.potentials, 0.0); - } #ifdef SCAFACOS_DIPOLES - scafacos->run_dipolar(particles.dipoles, particles.positions, - particles.fields, particles.potentials); - return -0.5 * dipole.prefactor * - boost::inner_product(particles.dipoles, particles.potentials, 0.0); +double ScafacosDataDipoles::long_range_energy() { + update_particle_data(); + scafacos->run_dipolar(dipoles, positions, fields, potentials); + return -0.5 * dipole.prefactor * + boost::inner_product(dipoles, potentials, 0.0); +} #endif - } +double long_range_energy() { + if (scafacos) { + return scafacos_data->long_range_energy(); + } return 0.0; } + void set_r_cut_and_tune_local(double r_cut) { - particles.update_particle_data(); + scafacos_data->update_particle_data(); scafacos->set_r_cut(r_cut); - scafacos->tune(particles.charges, particles.positions); + scafacos_data->tune(); } REGISTER_CALLBACK(set_r_cut_and_tune_local) @@ -272,7 +321,7 @@ void tune_r_cut() { } void tune() { - particles.update_particle_data(); + scafacos_data->update_particle_data(); /* Check whether we have to do a bisection for the short range cutoff */ /* Check if there is a user supplied cutoff */ @@ -286,7 +335,7 @@ void tune() { } } else { // ESPResSo is not affected by a short range cutoff. Tune in parallel - scafacos->tune(particles.charges, particles.positions); + scafacos_data->tune(); } } @@ -294,7 +343,9 @@ static void set_params_safe(const std::string &method, const std::string ¶ms, bool dipolar_ia, int n_part) { delete scafacos; + delete scafacos_data; scafacos = nullptr; + scafacos_data = nullptr; scafacos = new Scafacos(method, comm_cart, params); @@ -302,14 +353,16 @@ static void set_params_safe(const std::string &method, box_geo.periodic(2) != 0}; scafacos->set_dipolar(dipolar_ia); -#ifdef DIPOLES +#ifdef SCAFACOS_DIPOLES if (dipolar_ia) { dipole.method = DIPOLAR_SCAFACOS; + scafacos_data = new ScafacosDataDipoles(); } #endif #ifdef ELECTROSTATICS if (!dipolar_ia) { coulomb.method = COULOMB_SCAFACOS; + scafacos_data = new ScafacosDataCoulomb(); } #endif scafacos->set_common_parameters(box_geo.length().data(), per, n_part); @@ -348,6 +401,10 @@ double get_r_cut() { void set_parameters(const std::string &method, const std::string ¶ms, bool dipolar_ia) { mpi_call_all(set_params_safe, method, params, dipolar_ia, get_n_part()); + if (!scafacos) + throw std::runtime_error("Scafacos not initialized"); + if (!scafacos_data) + throw std::runtime_error("Scafacos data object not initialized"); } bool dipolar() { @@ -367,7 +424,9 @@ void free_handle() { if (this_node == 0) mpi_call(free_handle); delete scafacos; + delete scafacos_data; scafacos = nullptr; + scafacos_data = nullptr; } REGISTER_CALLBACK(free_handle) diff --git a/src/python/espressomd/scafacos.pxd b/src/python/espressomd/scafacos.pxd index d8c4d8465bc..666255e6c3a 100644 --- a/src/python/espressomd/scafacos.pxd +++ b/src/python/espressomd/scafacos.pxd @@ -24,7 +24,7 @@ from libcpp cimport bool from libcpp.list cimport list IF SCAFACOS: cdef extern from "electrostatics_magnetostatics/scafacos.hpp" namespace "Scafacos": - void set_parameters(string & method_name, string & params, bool dipolar) + void set_parameters(string & method_name, string & params, bool dipolar) except+ string get_method_and_parameters() list[string] available_methods_core "Scafacos::available_methods" () void free_handle() From 6cf160a8ef2a47ab3c5bf655b144bb81fdbd6bd6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 9 Dec 2020 16:21:05 +0100 Subject: [PATCH 05/14] core: Couple ScaFaCoS interface to its context --- .../scafacos.cpp | 122 +++++++----------- .../scafacos.hpp | 3 - 2 files changed, 48 insertions(+), 77 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index 5daba687789..783335aae22 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -69,8 +69,9 @@ std::list available_methods() { } /** Encapsulation for the particle data needed by ScaFaCoS */ -struct ScafacosData { - virtual ~ScafacosData() {} +struct ScafacosContext : Scafacos { + using Scafacos::Scafacos; + virtual ~ScafacosContext() {} /** @brief Collect particle data in continuous arrays as required * by ScaFaCoS. */ @@ -86,13 +87,22 @@ struct ScafacosData { std::vector fields, potentials; }; -struct ScafacosDataCoulomb : ScafacosData { - ScafacosDataCoulomb() {} +struct ScafacosContextCoulomb : ScafacosContext { + using ScafacosContext::ScafacosContext; void update_particle_data() override; void update_particle_forces() const override; - double long_range_energy() override; - void add_long_range_force() override; - void tune() override; + double long_range_energy() override { + update_particle_data(); + run(charges, positions, fields, potentials); + return 0.5 * coulomb.prefactor * + boost::inner_product(charges, potentials, 0.0); + } + void add_long_range_force() override { + update_particle_data(); + run(charges, positions, fields, potentials); + update_particle_forces(); + } + void tune() override { Scafacos::tune(charges, positions); } private: /** Inputs */ @@ -100,12 +110,21 @@ struct ScafacosDataCoulomb : ScafacosData { }; #ifdef SCAFACOS_DIPOLES -struct ScafacosDataDipoles : ScafacosData { - ScafacosDataDipoles() {} +struct ScafacosContextDipoles : ScafacosContext { + using ScafacosContext::ScafacosContext; void update_particle_data() override; void update_particle_forces() const override; - double long_range_energy() override; - void add_long_range_force() override; + double long_range_energy() { + update_particle_data(); + run_dipolar(dipoles, positions, fields, potentials); + return -0.5 * dipole.prefactor * + boost::inner_product(dipoles, potentials, 0.0); + } + void add_long_range_force() override { + update_particle_data(); + run_dipolar(dipoles, positions, fields, potentials); + update_particle_forces(); + } void tune() override { runtimeErrorMsg() << "Tuning unavailable for ScaFaCoS dipoles."; } @@ -116,12 +135,9 @@ struct ScafacosDataDipoles : ScafacosData { }; #endif -static Scafacos *scafacos = nullptr; -static ScafacosData *scafacos_data = nullptr; +static ScafacosContext *scafacos = nullptr; -void ScafacosDataCoulomb::tune() { scafacos->tune(charges, positions); } - -void ScafacosDataCoulomb::update_particle_data() { +void ScafacosContextCoulomb::update_particle_data() { positions.clear(); charges.clear(); @@ -135,7 +151,7 @@ void ScafacosDataCoulomb::update_particle_data() { } #ifdef SCAFACOS_DIPOLES -void ScafacosDataDipoles::update_particle_data() { +void ScafacosContextDipoles::update_particle_data() { positions.clear(); dipoles.clear(); @@ -152,7 +168,7 @@ void ScafacosDataDipoles::update_particle_data() { } #endif -void ScafacosDataCoulomb::update_particle_forces() const { +void ScafacosContextCoulomb::update_particle_forces() const { if (positions.empty()) return; @@ -168,7 +184,7 @@ void ScafacosDataCoulomb::update_particle_forces() const { } #ifdef SCAFACOS_DIPOLES -void ScafacosDataDipoles::update_particle_forces() const { +void ScafacosContextDipoles::update_particle_forces() const { if (positions.empty()) return; @@ -215,7 +231,6 @@ void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, return; assert(scafacos); - assert(scafacos_data); auto const field = scafacos->pair_force(dist); auto const fak = q1q2 * field / dist; force -= fak * d; @@ -227,55 +242,25 @@ double pair_energy(double q1q2, double dist) { return 0.; } -void ScafacosDataCoulomb::add_long_range_force() { - update_particle_data(); - scafacos->run(charges, positions, fields, potentials); - update_particle_forces(); -} - -#ifdef SCAFACOS_DIPOLES -void ScafacosDataDipoles::add_long_range_force() { - update_particle_data(); - scafacos->run_dipolar(dipoles, positions, fields, potentials); - update_particle_forces(); -} -#endif - void add_long_range_force() { if (!scafacos) { throw std::runtime_error("Scafacos not initialized"); } - scafacos_data->add_long_range_force(); + scafacos->add_long_range_force(); } -double ScafacosDataCoulomb::long_range_energy() { - update_particle_data(); - scafacos->run(charges, positions, fields, potentials); - return 0.5 * coulomb.prefactor * - boost::inner_product(charges, potentials, 0.0); -} - -#ifdef SCAFACOS_DIPOLES -double ScafacosDataDipoles::long_range_energy() { - update_particle_data(); - scafacos->run_dipolar(dipoles, positions, fields, potentials); - return -0.5 * dipole.prefactor * - boost::inner_product(dipoles, potentials, 0.0); -} -#endif - double long_range_energy() { if (scafacos) { - return scafacos_data->long_range_energy(); + return scafacos->long_range_energy(); } return 0.0; } void set_r_cut_and_tune_local(double r_cut) { - scafacos_data->update_particle_data(); + scafacos->update_particle_data(); scafacos->set_r_cut(r_cut); - scafacos_data->tune(); + scafacos->tune(); } REGISTER_CALLBACK(set_r_cut_and_tune_local) @@ -321,7 +306,7 @@ void tune_r_cut() { } void tune() { - scafacos_data->update_particle_data(); + scafacos->update_particle_data(); /* Check whether we have to do a bisection for the short range cutoff */ /* Check if there is a user supplied cutoff */ @@ -335,7 +320,7 @@ void tune() { } } else { // ESPResSo is not affected by a short range cutoff. Tune in parallel - scafacos_data->tune(); + scafacos->tune(); } } @@ -343,28 +328,28 @@ static void set_params_safe(const std::string &method, const std::string ¶ms, bool dipolar_ia, int n_part) { delete scafacos; - delete scafacos_data; scafacos = nullptr; - scafacos_data = nullptr; - - scafacos = new Scafacos(method, comm_cart, params); int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0, box_geo.periodic(2) != 0}; - scafacos->set_dipolar(dipolar_ia); #ifdef SCAFACOS_DIPOLES if (dipolar_ia) { dipole.method = DIPOLAR_SCAFACOS; - scafacos_data = new ScafacosDataDipoles(); + scafacos = new ScafacosContextDipoles(method, comm_cart, params); } #endif #ifdef ELECTROSTATICS if (!dipolar_ia) { coulomb.method = COULOMB_SCAFACOS; - scafacos_data = new ScafacosDataCoulomb(); + scafacos = new ScafacosContextCoulomb(method, comm_cart, params); } #endif + if (!scafacos) { + return; + } + + scafacos->set_dipolar(dipolar_ia); scafacos->set_common_parameters(box_geo.length().data(), per, n_part); on_coulomb_change(); @@ -403,8 +388,6 @@ void set_parameters(const std::string &method, const std::string ¶ms, mpi_call_all(set_params_safe, method, params, dipolar_ia, get_n_part()); if (!scafacos) throw std::runtime_error("Scafacos not initialized"); - if (!scafacos_data) - throw std::runtime_error("Scafacos data object not initialized"); } bool dipolar() { @@ -413,20 +396,11 @@ bool dipolar() { throw std::runtime_error("Scafacos not initialized"); } -void set_dipolar(bool d) { - if (scafacos) { - scafacos->set_dipolar(d); - } else - runtimeErrorMsg() << "Scafacos not initialized."; -} - void free_handle() { if (this_node == 0) mpi_call(free_handle); delete scafacos; - delete scafacos_data; scafacos = nullptr; - scafacos_data = nullptr; } REGISTER_CALLBACK(free_handle) diff --git a/src/core/electrostatics_magnetostatics/scafacos.hpp b/src/core/electrostatics_magnetostatics/scafacos.hpp index 6a30fcbe61e..a5cd0590c25 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.hpp +++ b/src/core/electrostatics_magnetostatics/scafacos.hpp @@ -54,9 +54,6 @@ double get_r_cut(); /** Is scafacos used for dipolar interactions */ bool dipolar(); -/** Choose whether scafacos is used for dipolar interactions */ -void set_dipolar(bool d); - /** Reinit scafacos number of particles, box shape and periodicity */ void update_system_params(); From ea35259f567a380dff14f97adee7d099bda8e9af Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 9 Dec 2020 20:39:14 +0100 Subject: [PATCH 06/14] core: Completely separate ScaFaCoS Coulomb/Dipoles There are now two static globals: one for ScaFaCoS Coulomb and one for ScaFaCoS dipoles. We can now in principle run a magnetostatics method and an electrostatics method simultaneously. --- .../CMakeLists.txt | 1 + .../ScafacosContext.cpp | 167 ++++++++ .../ScafacosContext.hpp | 118 ++++++ .../electrostatics_magnetostatics/coulomb.cpp | 10 +- .../coulomb_inline.hpp | 4 +- .../electrostatics_magnetostatics/dipole.cpp | 9 +- .../scafacos.cpp | 376 ++++++------------ .../scafacos.hpp | 28 +- src/core/event.cpp | 4 +- src/python/espressomd/electrostatics.pyx | 2 +- src/python/espressomd/magnetostatics.pyx | 2 +- src/python/espressomd/scafacos.pxd | 4 +- src/python/espressomd/scafacos.pyx | 2 +- 13 files changed, 452 insertions(+), 275 deletions(-) create mode 100644 src/core/electrostatics_magnetostatics/ScafacosContext.cpp create mode 100644 src/core/electrostatics_magnetostatics/ScafacosContext.hpp diff --git a/src/core/electrostatics_magnetostatics/CMakeLists.txt b/src/core/electrostatics_magnetostatics/CMakeLists.txt index 498dcf79e30..fcd3037dfbf 100644 --- a/src/core/electrostatics_magnetostatics/CMakeLists.txt +++ b/src/core/electrostatics_magnetostatics/CMakeLists.txt @@ -14,6 +14,7 @@ target_sources( ${CMAKE_CURRENT_SOURCE_DIR}/p3m-dipolar.cpp ${CMAKE_CURRENT_SOURCE_DIR}/p3m_gpu.cpp ${CMAKE_CURRENT_SOURCE_DIR}/scafacos.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/ScafacosContext.cpp ${CMAKE_CURRENT_SOURCE_DIR}/fft.cpp ${CMAKE_CURRENT_SOURCE_DIR}/coulomb.cpp ${CMAKE_CURRENT_SOURCE_DIR}/dipole.cpp diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.cpp b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp new file mode 100644 index 00000000000..6c9fd15f850 --- /dev/null +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp @@ -0,0 +1,167 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 by + * 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 . + */ +/** @file + * Provide a C-like interface for ScaFaCoS. + */ + +#include "config.hpp" + +#if defined(SCAFACOS) + +#include "electrostatics_magnetostatics/ScafacosContext.hpp" + +#include "cells.hpp" +#include "communication.hpp" +#include "electrostatics_magnetostatics/coulomb.hpp" +#include "electrostatics_magnetostatics/dipole.hpp" +#include "grid.hpp" +#include "particle_data.hpp" + +#include + +#include +#include + +#include + +namespace Scafacos { + +std::string ScafacosContext::get_method_and_parameters() { + std::string representation = get_method() + " " + get_parameters(); + std::replace(representation.begin(), representation.end(), ',', ' '); + return representation; +} + +void ScafacosContext::update_system_params() { + int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0, + box_geo.periodic(2) != 0}; + + auto const n_part = boost::mpi::all_reduce( + comm_cart, cell_structure.local_particles().size(), std::plus<>()); + set_common_parameters(box_geo.length().data(), per, n_part); +} + +void ScafacosContextCoulomb::update_particle_data() { + positions.clear(); + charges.clear(); + + for (auto const &p : cell_structure.local_particles()) { + auto const pos = folded_position(p.r.p, box_geo); + positions.push_back(pos[0]); + positions.push_back(pos[1]); + positions.push_back(pos[2]); + charges.push_back(p.p.q); + } +} + +#ifdef SCAFACOS_DIPOLES +void ScafacosContextDipoles::update_particle_data() { + positions.clear(); + dipoles.clear(); + + for (auto const &p : cell_structure.local_particles()) { + auto const pos = folded_position(p.r.p, box_geo); + positions.push_back(pos[0]); + positions.push_back(pos[1]); + positions.push_back(pos[2]); + auto const dip = p.calc_dip(); + dipoles.push_back(dip[0]); + dipoles.push_back(dip[1]); + dipoles.push_back(dip[2]); + } +} +#endif + +void ScafacosContextCoulomb::update_particle_forces() const { + if (positions.empty()) + return; + + int it = 0; + for (auto &p : cell_structure.local_particles()) { + p.f.f += coulomb.prefactor * p.p.q * + Utils::Vector3d(Utils::Span(&(fields[it]), 3)); + it += 3; + } + + /* Check that the particle number did not change */ + assert(it == fields.size()); +} + +#ifdef SCAFACOS_DIPOLES +void ScafacosContextDipoles::update_particle_forces() const { + if (positions.empty()) + return; + + int it = 0; + for (auto &p : cell_structure.local_particles()) { + // Indices + // 3 "potential" values per particles (see below) + int const it_p = 3 * it; + // 6 "field" values per particles (see below) + int const it_f = 6 * it; + + // The scafacos term "potential" here in fact refers to the magnetic + // field. So, the torques are given by m \times B + auto const dip = p.calc_dip(); + auto const t = vector_product( + dip, + Utils::Vector3d(Utils::Span(&(potentials[it_p]), 3))); + // The force is given by G m, where G is a matrix + // which comes from the "fields" output of scafacos like this + // 0 1 2 + // 1 3 4 + // 2 4 5 + // where the numbers refer to indices in the "field" output from scafacos + auto const G = Utils::Vector{ + {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]}, + {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]}, + {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}}; + auto const f = G * dip; + + // Add to particles + p.f.f += dipole.prefactor * f; + p.f.torque += dipole.prefactor * t; + it++; + } + + /* Check that the particle number did not change */ + assert(it == positions.size() / 3); +} +#endif + +double ScafacosContextCoulomb::long_range_energy() { + update_particle_data(); + run(charges, positions, fields, potentials); + return 0.5 * coulomb.prefactor * + boost::inner_product(charges, potentials, 0.0); +} + +#ifdef SCAFACOS_DIPOLES +double ScafacosContextDipoles::long_range_energy() { + update_particle_data(); + run_dipolar(dipoles, positions, fields, potentials); + return -0.5 * dipole.prefactor * + boost::inner_product(dipoles, potentials, 0.0); +} +#endif + +} // namespace Scafacos +#endif /* SCAFACOS */ diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp new file mode 100644 index 00000000000..a27a6e1310a --- /dev/null +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp @@ -0,0 +1,118 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 by + * 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 . + */ +#ifndef SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP +#define SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP +/** @file + * Provide a C-like interface for ScaFaCoS. + */ + +#include "config.hpp" + +#if defined(SCAFACOS) + +#include "Scafacos.hpp" + +#include + +#include +#include + +#if defined(SCAFACOS_DIPOLES) && !defined(FCS_ENABLE_DIPOLES) +#error \ + "SCAFACOS_DIPOLES requires dipoles support in scafacos library (FCS_ENABLE_DIPOLES)." +#endif + +namespace Scafacos { + +/** Encapsulation for the particle data needed by ScaFaCoS */ +struct ScafacosContext : Scafacos { + using Scafacos::Scafacos; + virtual ~ScafacosContext() {} + /** @brief Collect particle data in continuous arrays as required + * by ScaFaCoS. + */ + virtual void update_particle_data() = 0; + /** @brief Write forces back to particles. */ + virtual void update_particle_forces() const = 0; + virtual double long_range_energy() = 0; + virtual void add_long_range_force() = 0; + /** @brief Reinitialize number of particles, box shape and periodicity. */ + void update_system_params(); + void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, + Utils::Vector3d &force) { + if (dist > r_cut()) + return; + + auto const field = Scafacos::pair_force(dist); + auto const fak = q1q2 * field / dist; + force -= fak * d; + } + double pair_energy(double q1q2, double dist) { + if (dist > r_cut()) + return 0.; + + return q1q2 * Scafacos::pair_energy(dist); + } + std::string get_method_and_parameters(); + +protected: + /** Outputs */ + std::vector fields, potentials; +}; + +struct ScafacosContextCoulomb : ScafacosContext { + using ScafacosContext::ScafacosContext; + void update_particle_data() override; + void update_particle_forces() const override; + double long_range_energy() override; + void add_long_range_force() override { + update_particle_data(); + run(charges, positions, fields, potentials); + update_particle_forces(); + } + void tune() { Scafacos::tune(charges, positions); } + +private: + /** Inputs */ + std::vector positions, charges; +}; + +#ifdef SCAFACOS_DIPOLES +struct ScafacosContextDipoles : ScafacosContext { + using ScafacosContext::ScafacosContext; + void update_particle_data() override; + void update_particle_forces() const override; + double long_range_energy() override; + void add_long_range_force() override { + update_particle_data(); + run_dipolar(dipoles, positions, fields, potentials); + update_particle_forces(); + } + +private: + /** Inputs */ + std::vector positions, dipoles; +}; +#endif + +} // namespace Scafacos +#endif // SCAFACOS +#endif // SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index af3a2720a37..10ddd0f6831 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -122,7 +122,7 @@ double cutoff(const Utils::Vector3d &box_l) { return rf_params.r_cut; #ifdef SCAFACOS case COULOMB_SCAFACOS: - return Scafacos::get_r_cut(); + return Scafacos::Coulomb::get_r_cut(); #endif default: return -1.0; @@ -217,7 +217,7 @@ void on_boxl_change() { break; #ifdef SCAFACOS case COULOMB_SCAFACOS: - Scafacos::update_system_params(); + Scafacos::Coulomb::update_system_params(); break; #endif default: @@ -291,8 +291,7 @@ void calc_long_range_force(const ParticleRange &particles) { #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: - assert(!Scafacos::dipolar()); - Scafacos::add_long_range_force(); + Scafacos::Coulomb::add_long_range_force(); break; #endif default: @@ -350,8 +349,7 @@ double calc_energy_long_range(const ParticleRange &particles) { #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: - assert(!Scafacos::dipolar()); - energy += Scafacos::long_range_energy(); + energy += Scafacos::Coulomb::long_range_energy(); break; #endif default: diff --git a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp index 8c629727e21..76308293eab 100644 --- a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp @@ -58,7 +58,7 @@ inline Utils::Vector3d central_force(double const q1q2, break; #ifdef SCAFACOS case COULOMB_SCAFACOS: - Scafacos::add_pair_force(q1q2, d, dist, f); + Scafacos::Coulomb::add_pair_force(q1q2, d, dist, f); break; #endif default: @@ -157,7 +157,7 @@ inline double pair_energy(Particle const &p1, Particle const &p2, #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: - return Scafacos::pair_energy(q1q2, dist); + return Scafacos::Coulomb::pair_energy(q1q2, dist); #endif case COULOMB_DH: return dh_coulomb_pair_energy(q1q2, dist); diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index f188320be3a..3d3f495122d 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -140,7 +140,7 @@ void on_boxl_change() { #endif #ifdef SCAFACOS case DIPOLAR_SCAFACOS: - Scafacos::update_system_params(); + Scafacos::Dipoles::update_system_params(); break; #endif default: @@ -202,8 +202,7 @@ void calc_long_range_force(const ParticleRange &particles) { #endif #ifdef SCAFACOS_DIPOLES case DIPOLAR_SCAFACOS: - assert(Scafacos::dipolar()); - Scafacos::add_long_range_force(); + Scafacos::Dipoles::add_long_range_force(); #endif case DIPOLAR_NONE: break; @@ -249,8 +248,8 @@ double calc_energy_long_range(const ParticleRange &particles) { #endif #ifdef SCAFACOS_DIPOLES case DIPOLAR_SCAFACOS: - assert(Scafacos::dipolar()); - energy = Scafacos::long_range_energy(); + energy = Scafacos::Dipoles::long_range_energy(); + break; #endif case DIPOLAR_NONE: break; diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index 783335aae22..f08828da053 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -31,6 +31,7 @@ #include "Scafacos.hpp" #include "cells.hpp" #include "communication.hpp" +#include "electrostatics_magnetostatics/ScafacosContext.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/dipole.hpp" #include "errorhandling.hpp" @@ -44,7 +45,6 @@ #include #include -#include #include #include @@ -56,11 +56,6 @@ #include #include -#if defined(SCAFACOS_DIPOLES) && !defined(FCS_ENABLE_DIPOLES) -#error \ - "SCAFACOS_DIPOLES requires dipoles support in scafacos library (FCS_ENABLE_DIPOLES)." -#endif - namespace Scafacos { /** Get available ScaFaCoS methods */ @@ -68,197 +63,135 @@ std::list available_methods() { return Scafacos::available_methods(); } -/** Encapsulation for the particle data needed by ScaFaCoS */ -struct ScafacosContext : Scafacos { - using Scafacos::Scafacos; - virtual ~ScafacosContext() {} - /** @brief Collect particle data in continuous arrays as required - * by ScaFaCoS. - */ - virtual void update_particle_data() = 0; - /** @brief Write forces back to particles. */ - virtual void update_particle_forces() const = 0; - virtual double long_range_energy() = 0; - virtual void add_long_range_force() = 0; - virtual void tune() = 0; - -protected: - /** Outputs */ - std::vector fields, potentials; -}; - -struct ScafacosContextCoulomb : ScafacosContext { - using ScafacosContext::ScafacosContext; - void update_particle_data() override; - void update_particle_forces() const override; - double long_range_energy() override { - update_particle_data(); - run(charges, positions, fields, potentials); - return 0.5 * coulomb.prefactor * - boost::inner_product(charges, potentials, 0.0); - } - void add_long_range_force() override { - update_particle_data(); - run(charges, positions, fields, potentials); - update_particle_forces(); - } - void tune() override { Scafacos::tune(charges, positions); } - -private: - /** Inputs */ - std::vector positions, charges; -}; - -#ifdef SCAFACOS_DIPOLES -struct ScafacosContextDipoles : ScafacosContext { - using ScafacosContext::ScafacosContext; - void update_particle_data() override; - void update_particle_forces() const override; - double long_range_energy() { - update_particle_data(); - run_dipolar(dipoles, positions, fields, potentials); - return -0.5 * dipole.prefactor * - boost::inner_product(dipoles, potentials, 0.0); - } - void add_long_range_force() override { - update_particle_data(); - run_dipolar(dipoles, positions, fields, potentials); - update_particle_forces(); - } - void tune() override { - runtimeErrorMsg() << "Tuning unavailable for ScaFaCoS dipoles."; - } - -private: - /** Inputs */ - std::vector positions, dipoles; -}; -#endif +namespace Dipoles { -static ScafacosContext *scafacos = nullptr; +static ScafacosContextDipoles *scafacos = nullptr; -void ScafacosContextCoulomb::update_particle_data() { - positions.clear(); - charges.clear(); +void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, + Utils::Vector3d &force) { + if (!scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + return scafacos->add_pair_force(q1q2, d, dist, force); +} - for (auto const &p : cell_structure.local_particles()) { - auto const pos = folded_position(p.r.p, box_geo); - positions.push_back(pos[0]); - positions.push_back(pos[1]); - positions.push_back(pos[2]); - charges.push_back(p.p.q); - } +double pair_energy(double q1q2, double dist) { + if (!scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + return scafacos->pair_energy(q1q2, dist); } -#ifdef SCAFACOS_DIPOLES -void ScafacosContextDipoles::update_particle_data() { - positions.clear(); - dipoles.clear(); - - for (auto const &p : cell_structure.local_particles()) { - auto const pos = folded_position(p.r.p, box_geo); - positions.push_back(pos[0]); - positions.push_back(pos[1]); - positions.push_back(pos[2]); - auto const dip = p.calc_dip(); - dipoles.push_back(dip[0]); - dipoles.push_back(dip[1]); - dipoles.push_back(dip[2]); - } +void add_long_range_force() { + if (!scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + scafacos->add_long_range_force(); } -#endif -void ScafacosContextCoulomb::update_particle_forces() const { - if (positions.empty()) - return; +double long_range_energy() { + if (!scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + return scafacos->long_range_energy(); +} - int it = 0; - for (auto &p : cell_structure.local_particles()) { - p.f.f += coulomb.prefactor * p.p.q * - Utils::Vector3d(Utils::Span(&(fields[it]), 3)); - it += 3; - } +double get_r_cut() { + if (!scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + return scafacos->r_cut(); +} - /* Check that the particle number did not change */ - assert(it == fields.size()); +void update_system_params() { + if (!scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + scafacos->update_system_params(); } -#ifdef SCAFACOS_DIPOLES -void ScafacosContextDipoles::update_particle_forces() const { - if (positions.empty()) - return; +static void set_parameters_worker(const std::string &method, + const std::string ¶ms) { + delete scafacos; + scafacos = nullptr; - int it = 0; - for (auto &p : cell_structure.local_particles()) { - // Indices - // 3 "potential" values per particles (see below) - int const it_p = 3 * it; - // 6 "field" values per particles (see below) - int const it_f = 6 * it; - - // The scafacos term "potential" here in fact refers to the magnetic - // field. So, the torques are given by m \times B - auto const dip = p.calc_dip(); - auto const t = vector_product( - dip, - Utils::Vector3d(Utils::Span(&(potentials[it_p]), 3))); - // The force is given by G m, where G is a matrix - // which comes from the "fields" output of scafacos like this - // 0 1 2 - // 1 3 4 - // 2 4 5 - // where the numbers refer to indices in the "field" output from scafacos - auto const G = Utils::Vector{ - {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]}, - {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]}, - {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}}; - auto const f = G * dip; - - // Add to particles - p.f.f += dipole.prefactor * f; - p.f.torque += dipole.prefactor * t; - it++; + scafacos = new ScafacosContextDipoles(method, comm_cart, params); + if (!scafacos) { + runtimeErrorMsg() << "Scafacos Dipoles failed to initialize"; + return; } - /* Check that the particle number did not change */ - assert(it == positions.size() / 3); + scafacos->set_dipolar(true); + scafacos->update_system_params(); + + dipole.method = DIPOLAR_SCAFACOS; + on_coulomb_change(); } -#endif + +REGISTER_CALLBACK(set_parameters_worker) + +} // namespace Dipoles + +namespace Coulomb { + +static ScafacosContextCoulomb *scafacos = nullptr; void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, Utils::Vector3d &force) { - if (dist > get_r_cut()) - return; - - assert(scafacos); - auto const field = scafacos->pair_force(dist); - auto const fak = q1q2 * field / dist; - force -= fak * d; + if (!scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + return scafacos->add_pair_force(q1q2, d, dist, force); } double pair_energy(double q1q2, double dist) { - if (dist <= get_r_cut()) - return q1q2 * scafacos->pair_energy(dist); - return 0.; + if (!scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + return scafacos->pair_energy(q1q2, dist); } void add_long_range_force() { - if (!scafacos) { - throw std::runtime_error("Scafacos not initialized"); - } + if (!scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); scafacos->add_long_range_force(); } double long_range_energy() { - if (scafacos) { - return scafacos->long_range_energy(); + if (!scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + return scafacos->long_range_energy(); +} + +double get_r_cut() { + if (!scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + return scafacos->r_cut(); +} + +void update_system_params() { + if (!scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + scafacos->update_system_params(); +} + +void tune(); + +static void set_parameters_worker(const std::string &method, + const std::string ¶ms) { + delete scafacos; + scafacos = nullptr; + + scafacos = new ScafacosContextCoulomb(method, comm_cart, params); + if (!scafacos) { + runtimeErrorMsg() << "Scafacos Coulomb failed to initialize"; + return; } - return 0.0; + + scafacos->set_dipolar(false); + scafacos->update_system_params(); + + coulomb.method = COULOMB_SCAFACOS; + on_coulomb_change(); + tune(); } +REGISTER_CALLBACK(set_parameters_worker) + void set_r_cut_and_tune_local(double r_cut) { scafacos->update_particle_data(); - scafacos->set_r_cut(r_cut); scafacos->tune(); } @@ -323,100 +256,45 @@ void tune() { scafacos->tune(); } } +} // namespace Coulomb -static void set_params_safe(const std::string &method, - const std::string ¶ms, bool dipolar_ia, - int n_part) { - delete scafacos; - scafacos = nullptr; - - int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0, - box_geo.periodic(2) != 0}; - -#ifdef SCAFACOS_DIPOLES - if (dipolar_ia) { - dipole.method = DIPOLAR_SCAFACOS; - scafacos = new ScafacosContextDipoles(method, comm_cart, params); - } -#endif -#ifdef ELECTROSTATICS - if (!dipolar_ia) { - coulomb.method = COULOMB_SCAFACOS; - scafacos = new ScafacosContextCoulomb(method, comm_cart, params); - } -#endif - if (!scafacos) { - return; - } - - scafacos->set_dipolar(dipolar_ia); - scafacos->set_common_parameters(box_geo.length().data(), per, n_part); - - on_coulomb_change(); - - if (!dipolar_ia) { - tune(); - } -} - -REGISTER_CALLBACK(set_params_safe) - -/** Bend result from scafacos back to original format */ -std::string get_method_and_parameters() { - if (!scafacos) { - return std::string(); +void free_handle(bool dipolar) { + if (this_node == 0) + mpi_call(free_handle, dipolar); + if (dipolar) { + delete Dipoles::scafacos; + Dipoles::scafacos = nullptr; + } else { + delete Coulomb::scafacos; + Coulomb::scafacos = nullptr; } - - std::string p = scafacos->get_method() + " " + scafacos->get_parameters(); - - std::replace(p.begin(), p.end(), ',', ' '); - - return p; } -double get_r_cut() { - if (scafacos) { - if (!scafacos->has_near) - return 0; - return scafacos->r_cut(); - } - return 0.0; -} +REGISTER_CALLBACK(free_handle) void set_parameters(const std::string &method, const std::string ¶ms, - bool dipolar_ia) { - mpi_call_all(set_params_safe, method, params, dipolar_ia, get_n_part()); - if (!scafacos) - throw std::runtime_error("Scafacos not initialized"); -} - -bool dipolar() { - if (scafacos) - return scafacos->dipolar(); - throw std::runtime_error("Scafacos not initialized"); -} - -void free_handle() { - if (this_node == 0) - mpi_call(free_handle); - delete scafacos; - scafacos = nullptr; + bool dipolar) { + if (dipolar) { + mpi_call_all(Dipoles::set_parameters_worker, method, params); + if (!Dipoles::scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + } else { + mpi_call_all(Coulomb::set_parameters_worker, method, params); + if (!Coulomb::scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + } } -REGISTER_CALLBACK(free_handle) - -void update_system_params() { - // If scafacos is not active, do nothing - if (!scafacos) { - throw std::runtime_error("Scafacos not initialized"); +std::string get_method_and_parameters(bool dipolar) { + if (dipolar) { + if (!Dipoles::scafacos) + throw std::runtime_error("Scafacos Dipoles not initialized"); + return Dipoles::scafacos->get_method_and_parameters(); + } else { + if (!Coulomb::scafacos) + throw std::runtime_error("Scafacos Coulomb not initialized"); + return Coulomb::scafacos->get_method_and_parameters(); } - - int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0, - box_geo.periodic(2) != 0}; - - auto const n_part = boost::mpi::all_reduce( - comm_cart, cell_structure.local_particles().size(), std::plus<>()); - scafacos->set_common_parameters(box_geo.length().data(), per, n_part); } } // namespace Scafacos diff --git a/src/core/electrostatics_magnetostatics/scafacos.hpp b/src/core/electrostatics_magnetostatics/scafacos.hpp index a5cd0590c25..8db8e456b9f 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.hpp +++ b/src/core/electrostatics_magnetostatics/scafacos.hpp @@ -35,6 +35,8 @@ namespace Scafacos { #if defined(SCAFACOS) + +namespace Dipoles { /** Near-field pair force */ void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, Utils::Vector3d &force); @@ -47,20 +49,34 @@ double long_range_energy(); /** Get parameters */ std::string get_method_and_parameters(); /** Set parameters */ -void set_parameters(const std::string &method, const std::string ¶ms, - bool dipolar); +void set_parameters(const std::string &method, const std::string ¶ms); double get_r_cut(); -/** Is scafacos used for dipolar interactions */ -bool dipolar(); +void update_system_params(); +} // namespace Dipoles + +namespace Coulomb { +/** Near-field pair force */ +void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, + Utils::Vector3d &force); +/** Near-field pair energy */ +double pair_energy(double q1q2, double dist); +/** Long range part */ +void add_long_range_force(); +/** Calculate long range energy contribution */ +double long_range_energy(); +double get_r_cut(); -/** Reinit scafacos number of particles, box shape and periodicity */ void update_system_params(); +} // namespace Coulomb #endif /* SCAFACOS */ std::list available_methods(); -void free_handle(); +std::string get_method_and_parameters(bool dipolar); +void set_parameters(const std::string &method, const std::string ¶ms, + bool dipolar); +void free_handle(bool dipolar); } // namespace Scafacos #endif diff --git a/src/core/event.cpp b/src/core/event.cpp index 9050f774e9c..ebe0b29e9e8 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -280,12 +280,12 @@ void on_parameter_change(int field) { #ifdef SCAFACOS #ifdef ELECTROSTATICS if (coulomb.method == COULOMB_SCAFACOS) { - Scafacos::update_system_params(); + Scafacos::Coulomb::update_system_params(); } #endif #ifdef DIPOLES if (dipole.method == DIPOLAR_SCAFACOS) { - Scafacos::update_system_params(); + Scafacos::Dipoles::update_system_params(); } #endif #endif diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx index 86f157782a7..7f64ebad395 100644 --- a/src/python/espressomd/electrostatics.pyx +++ b/src/python/espressomd/electrostatics.pyx @@ -691,5 +691,5 @@ IF ELECTROSTATICS: def _deactivate_method(self): super()._deactivate_method() - scafacos.free_handle() + scafacos.free_handle(self.dipolar) mpi_bcast_coulomb_params() diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx index 85d2d33c4f7..7e3f29f138f 100644 --- a/src/python/espressomd/magnetostatics.pyx +++ b/src/python/espressomd/magnetostatics.pyx @@ -345,7 +345,7 @@ IF DIPOLES == 1: def _deactivate_method(self): dipole.method = DIPOLAR_NONE - scafacos.free_handle() + scafacos.free_handle(self.dipolar) mpi_bcast_coulomb_params() def default_params(self): diff --git a/src/python/espressomd/scafacos.pxd b/src/python/espressomd/scafacos.pxd index 666255e6c3a..a14e1bc4656 100644 --- a/src/python/espressomd/scafacos.pxd +++ b/src/python/espressomd/scafacos.pxd @@ -25,6 +25,6 @@ from libcpp.list cimport list IF SCAFACOS: cdef extern from "electrostatics_magnetostatics/scafacos.hpp" namespace "Scafacos": void set_parameters(string & method_name, string & params, bool dipolar) except+ - string get_method_and_parameters() + string get_method_and_parameters(bool dipolar) except+ + void free_handle(bool dipolar) list[string] available_methods_core "Scafacos::available_methods" () - void free_handle() diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx index da6133a81e0..99105219625 100644 --- a/src/python/espressomd/scafacos.pyx +++ b/src/python/espressomd/scafacos.pyx @@ -60,7 +60,7 @@ IF SCAFACOS == 1: # Parameters are returned as strings # First word is method name, rest are key value pairs # which we convert to a dict - p = to_str(get_method_and_parameters()).split(" ") + p = to_str(get_method_and_parameters(self.dipolar)).split(" ") res = {} res["method_name"] = p[0] From 3656ec3e9101f6ad7ff886c28de0d92b170888c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 9 Dec 2020 20:41:49 +0100 Subject: [PATCH 07/14] core: Pull out ScaFaCoS tuning --- .../ScafacosContext.cpp | 80 +++++++++++++++++ .../ScafacosContext.hpp | 3 +- .../scafacos.cpp | 85 +------------------ .../scafacos.hpp | 2 + src/python/espressomd/scafacos.pxd | 2 +- 5 files changed, 87 insertions(+), 85 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.cpp b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp index 6c9fd15f850..3e01c1c5792 100644 --- a/src/core/electrostatics_magnetostatics/ScafacosContext.cpp +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp @@ -32,15 +32,27 @@ #include "communication.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/dipole.hpp" +#include "electrostatics_magnetostatics/scafacos.hpp" #include "grid.hpp" +#include "integrate.hpp" #include "particle_data.hpp" +#include "tuning.hpp" #include #include +#include #include +#include #include +#include +#include +#include +#include +#include +#include +#include namespace Scafacos { @@ -163,5 +175,73 @@ double ScafacosContextDipoles::long_range_energy() { } #endif +static void set_r_cut_and_tune_local(double r_cut) { + Coulomb::set_r_cut_and_tune_local(r_cut); +} + +REGISTER_CALLBACK(set_r_cut_and_tune_local) + +/** Determine runtime for a specific cutoff */ +static double time_r_cut(double r_cut) { + /* Set cutoff to time */ + mpi_call_all(set_r_cut_and_tune_local, r_cut); + + return time_force_calc(10); +} + +/** Determine the optimal cutoff by bisection */ +static void tune_r_cut() { + const double tune_limit = 1e-3; + + auto const min_box_l = *boost::min_element(box_geo.length()); + auto const min_local_box_l = *boost::min_element(local_geo.length()); + + /* scafacos p3m and Ewald do not accept r_cut 0 for no good reason */ + double r_min = 1.0; + double r_max = std::min(min_local_box_l, min_box_l / 2.0) - skin; + double t_min = 0; + double t_max = std::numeric_limits::max(); + + /* Run bisection */ + while (std::fabs(r_min - r_max) > tune_limit) { + const double dr = 0.5 * (r_max - r_min); + const double t_mid = time_r_cut(r_min + dr); + t_min = time_r_cut(r_min); + t_max = time_r_cut(r_max); + + if (t_min <= 0.0 or t_max <= 0.0) { + break; + } + + if (t_mid > t_min) { + r_max = r_min += dr; + } else { + r_min += dr; + } + } +} + +void ScafacosContextCoulomb::tune() { + update_particle_data(); + + // Check whether we have to do a bisection for the short range cutoff + // Check if there is a user-supplied cutoff + if (has_near && r_cut() <= 0.0) { + // run tuning on the head node + if (this_node == 0) { + tune_r_cut(); + } + } else { + // ESPResSo is not affected by a short range cutoff -> tune in parallel + Scafacos::tune(charges, positions); + } +} + +void ScafacosContextCoulomb::set_r_cut_and_tune_local(double r_cut) { + update_particle_data(); + set_r_cut(r_cut); + Scafacos::tune(charges, positions); +} + } // namespace Scafacos #endif /* SCAFACOS */ diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp index a27a6e1310a..598e40787c5 100644 --- a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp @@ -88,7 +88,8 @@ struct ScafacosContextCoulomb : ScafacosContext { run(charges, positions, fields, potentials); update_particle_forces(); } - void tune() { Scafacos::tune(charges, positions); } + void tune(); + void set_r_cut_and_tune_local(double r_cut); private: /** Inputs */ diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index f08828da053..3d2579b1783 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -29,7 +29,6 @@ #include "electrostatics_magnetostatics/scafacos.hpp" #include "Scafacos.hpp" -#include "cells.hpp" #include "communication.hpp" #include "electrostatics_magnetostatics/ScafacosContext.hpp" #include "electrostatics_magnetostatics/coulomb.hpp" @@ -37,24 +36,12 @@ #include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" -#include "integrate.hpp" -#include "particle_data.hpp" -#include "tuning.hpp" #include -#include -#include - -#include -#include -#include -#include -#include #include #include #include -#include namespace Scafacos { @@ -167,8 +154,6 @@ void update_system_params() { scafacos->update_system_params(); } -void tune(); - static void set_parameters_worker(const std::string &method, const std::string ¶ms) { delete scafacos; @@ -185,77 +170,15 @@ static void set_parameters_worker(const std::string &method, coulomb.method = COULOMB_SCAFACOS; on_coulomb_change(); - tune(); + scafacos->tune(); } REGISTER_CALLBACK(set_parameters_worker) void set_r_cut_and_tune_local(double r_cut) { - scafacos->update_particle_data(); - scafacos->set_r_cut(r_cut); - scafacos->tune(); -} - -REGISTER_CALLBACK(set_r_cut_and_tune_local) - -/** Determine runtime for a specific cutoff */ -double time_r_cut(double r_cut) { - /* Set cutoff to time */ - mpi_call_all(set_r_cut_and_tune_local, r_cut); - - return time_force_calc(10); -} - -/** Determine the optimal cutoff by bisection */ -void tune_r_cut() { - const double tune_limit = 1e-3; - - auto const min_box_l = *boost::min_element(box_geo.length()); - auto const min_local_box_l = *boost::min_element(local_geo.length()); - - /* scafacos p3m and Ewald do not accept r_cut 0 for no good reason */ - double r_min = 1.0; - double r_max = std::min(min_local_box_l, min_box_l / 2.0) - skin; - double t_min = 0; - double t_max = std::numeric_limits::max(); - - /* Run bisection */ - while (std::fabs(r_min - r_max) > tune_limit) { - const double dr = 0.5 * (r_max - r_min); - const double t_mid = time_r_cut(r_min + dr); - t_min = time_r_cut(r_min); - t_max = time_r_cut(r_max); - - if (t_min <= 0.0 or t_max <= 0.0) { - break; - } - - if (t_mid > t_min) { - r_max = r_min += dr; - } else { - r_min += dr; - } - } + scafacos->set_r_cut_and_tune_local(r_cut); } -void tune() { - scafacos->update_particle_data(); - - /* Check whether we have to do a bisection for the short range cutoff */ - /* Check if there is a user supplied cutoff */ - if ((scafacos->has_near) && (scafacos->r_cut() <= 0.0)) { - // Tuning of r_cut needs to run on the master node because it relies on - // master-slave mode communication - if (this_node == 0) { - tune_r_cut(); - } else { - return; // Tune on the master node will issue mpi calls - } - } else { - // ESPResSo is not affected by a short range cutoff. Tune in parallel - scafacos->tune(); - } -} } // namespace Coulomb void free_handle(bool dipolar) { @@ -276,12 +199,8 @@ void set_parameters(const std::string &method, const std::string ¶ms, bool dipolar) { if (dipolar) { mpi_call_all(Dipoles::set_parameters_worker, method, params); - if (!Dipoles::scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); } else { mpi_call_all(Coulomb::set_parameters_worker, method, params); - if (!Coulomb::scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); } } diff --git a/src/core/electrostatics_magnetostatics/scafacos.hpp b/src/core/electrostatics_magnetostatics/scafacos.hpp index 8db8e456b9f..c77d335040b 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.hpp +++ b/src/core/electrostatics_magnetostatics/scafacos.hpp @@ -68,6 +68,8 @@ double long_range_energy(); double get_r_cut(); void update_system_params(); + +void set_r_cut_and_tune_local(double r_cut); } // namespace Coulomb #endif /* SCAFACOS */ diff --git a/src/python/espressomd/scafacos.pxd b/src/python/espressomd/scafacos.pxd index a14e1bc4656..d215ddea7fa 100644 --- a/src/python/espressomd/scafacos.pxd +++ b/src/python/espressomd/scafacos.pxd @@ -24,7 +24,7 @@ from libcpp cimport bool from libcpp.list cimport list IF SCAFACOS: cdef extern from "electrostatics_magnetostatics/scafacos.hpp" namespace "Scafacos": - void set_parameters(string & method_name, string & params, bool dipolar) except+ + void set_parameters(string & method_name, string & params, bool dipolar) string get_method_and_parameters(bool dipolar) except+ void free_handle(bool dipolar) list[string] available_methods_core "Scafacos::available_methods" () From 66daee4023b7afb0a39e91e156012dbc52fe644b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 9 Dec 2020 22:46:57 +0100 Subject: [PATCH 08/14] core: Use type erasure Simplify the ScaFaCoS interface by hiding implementation details behind a pointer to a class that doesn't depend on Scafacos.hpp. --- .../ScafacosContext.cpp | 8 +- .../ScafacosContext.hpp | 31 ++- .../ScafacosContextBase.hpp | 70 +++++++ .../electrostatics_magnetostatics/coulomb.cpp | 8 +- .../coulomb_inline.hpp | 4 +- .../electrostatics_magnetostatics/dipole.cpp | 8 +- .../scafacos.cpp | 177 ++++++------------ .../scafacos.hpp | 57 ++---- src/core/event.cpp | 6 +- 9 files changed, 178 insertions(+), 191 deletions(-) create mode 100644 src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.cpp b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp index 3e01c1c5792..eca126e5d94 100644 --- a/src/core/electrostatics_magnetostatics/ScafacosContext.cpp +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp @@ -175,16 +175,16 @@ double ScafacosContextDipoles::long_range_energy() { } #endif -static void set_r_cut_and_tune_local(double r_cut) { - Coulomb::set_r_cut_and_tune_local(r_cut); +static void set_r_cut_and_tune_local_worker(double r_cut) { + set_r_cut_and_tune_local(r_cut); } -REGISTER_CALLBACK(set_r_cut_and_tune_local) +REGISTER_CALLBACK(set_r_cut_and_tune_local_worker) /** Determine runtime for a specific cutoff */ static double time_r_cut(double r_cut) { /* Set cutoff to time */ - mpi_call_all(set_r_cut_and_tune_local, r_cut); + mpi_call_all(set_r_cut_and_tune_local_worker, r_cut); return time_force_calc(10); } diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp index 598e40787c5..8d28f12681f 100644 --- a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp @@ -20,8 +20,11 @@ */ #ifndef SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP #define SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP -/** @file - * Provide a C-like interface for ScaFaCoS. +/** + * @file + * @ref Scafacos::ScafacosContextBase implements the interface of the + * ScaFaCoS bridge. It is further derived for the coulombic and dipolar + * versions of ScaFaCoS. */ #include "config.hpp" @@ -29,6 +32,7 @@ #if defined(SCAFACOS) #include "Scafacos.hpp" +#include "electrostatics_magnetostatics/ScafacosContextBase.hpp" #include @@ -43,21 +47,13 @@ namespace Scafacos { /** Encapsulation for the particle data needed by ScaFaCoS */ -struct ScafacosContext : Scafacos { +struct ScafacosContext : ScafacosContextBase, Scafacos { using Scafacos::Scafacos; - virtual ~ScafacosContext() {} - /** @brief Collect particle data in continuous arrays as required - * by ScaFaCoS. - */ - virtual void update_particle_data() = 0; - /** @brief Write forces back to particles. */ - virtual void update_particle_forces() const = 0; - virtual double long_range_energy() = 0; - virtual void add_long_range_force() = 0; - /** @brief Reinitialize number of particles, box shape and periodicity. */ - void update_system_params(); + using ScafacosContextBase::ScafacosContextBase; + ~ScafacosContext() override = default; + void update_system_params() override; void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, - Utils::Vector3d &force) { + Utils::Vector3d &force) override { if (dist > r_cut()) return; @@ -65,13 +61,14 @@ struct ScafacosContext : Scafacos { auto const fak = q1q2 * field / dist; force -= fak * d; } - double pair_energy(double q1q2, double dist) { + double pair_energy(double q1q2, double dist) override { if (dist > r_cut()) return 0.; return q1q2 * Scafacos::pair_energy(dist); } - std::string get_method_and_parameters(); + std::string get_method_and_parameters() override; + double get_r_cut() const override { return Scafacos::r_cut(); } protected: /** Outputs */ diff --git a/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp b/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp new file mode 100644 index 00000000000..7f8e2635cb6 --- /dev/null +++ b/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp @@ -0,0 +1,70 @@ +/* + * Copyright (C) 2010-2020 The ESPResSo project + * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 + * Max-Planck-Institute for Polymer Research, Theory Group + * + * 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 by + * 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 . + */ +#ifndef SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXTBASE_HPP +#define SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXTBASE_HPP +/** + * @file + * @ref Scafacos::ScafacosContextBase provides the public interface + * of the ScaFaCoS bridge. It relies on type erasure to hide the + * ScaFaCoS implementation details from the ESPResSo core. It is + * implemented by @ref Scafacos::ScafacosContext. + */ + +#include "config.hpp" + +#if defined(SCAFACOS) + +#include + +#include + +namespace Scafacos { + +/** + * @brief Public interface to the ScaFaCoS context. + * Implementations of this class will store particle positions and + * charges/dipoles in flat arrays, as required by ScaFaCoS, as well + * as the output arrays. + */ +struct ScafacosContextBase { + virtual ~ScafacosContextBase() = default; + /** @brief Collect particle data in continuous arrays. */ + virtual void update_particle_data() = 0; + /** @brief Write forces back to particles. */ + virtual void update_particle_forces() const = 0; + /** @brief Calculate long-range part of the energy. */ + virtual double long_range_energy() = 0; + /** @brief Add long-range part of the forces to particles. */ + virtual void add_long_range_force() = 0; + /** @brief Add near-field pair force. */ + virtual void add_pair_force(double q1q2, Utils::Vector3d const &d, + double dist, Utils::Vector3d &force) = 0; + /** @brief Calculate near-field pair energy. */ + virtual double pair_energy(double q1q2, double dist) = 0; + /** @brief Reinitialize number of particles, box shape and periodicity. */ + virtual void update_system_params() = 0; + virtual double get_r_cut() const = 0; + virtual std::string get_method_and_parameters() = 0; +}; + +} // namespace Scafacos +#endif // SCAFACOS +#endif // SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXTBASE_HPP diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 10ddd0f6831..1fdf502f41d 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -122,7 +122,7 @@ double cutoff(const Utils::Vector3d &box_l) { return rf_params.r_cut; #ifdef SCAFACOS case COULOMB_SCAFACOS: - return Scafacos::Coulomb::get_r_cut(); + return Scafacos::fcs_coulomb()->get_r_cut(); #endif default: return -1.0; @@ -217,7 +217,7 @@ void on_boxl_change() { break; #ifdef SCAFACOS case COULOMB_SCAFACOS: - Scafacos::Coulomb::update_system_params(); + Scafacos::fcs_coulomb()->update_system_params(); break; #endif default: @@ -291,7 +291,7 @@ void calc_long_range_force(const ParticleRange &particles) { #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: - Scafacos::Coulomb::add_long_range_force(); + Scafacos::fcs_coulomb()->add_long_range_force(); break; #endif default: @@ -349,7 +349,7 @@ double calc_energy_long_range(const ParticleRange &particles) { #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: - energy += Scafacos::Coulomb::long_range_energy(); + energy += Scafacos::fcs_coulomb()->long_range_energy(); break; #endif default: diff --git a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp index 76308293eab..d291eb7f095 100644 --- a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp @@ -58,7 +58,7 @@ inline Utils::Vector3d central_force(double const q1q2, break; #ifdef SCAFACOS case COULOMB_SCAFACOS: - Scafacos::Coulomb::add_pair_force(q1q2, d, dist, f); + Scafacos::fcs_coulomb()->add_pair_force(q1q2, d, dist, f); break; #endif default: @@ -157,7 +157,7 @@ inline double pair_energy(Particle const &p1, Particle const &p2, #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: - return Scafacos::Coulomb::pair_energy(q1q2, dist); + return Scafacos::fcs_coulomb()->pair_energy(q1q2, dist); #endif case COULOMB_DH: return dh_coulomb_pair_energy(q1q2, dist); diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index 3d3f495122d..15716439d6a 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -138,9 +138,9 @@ void on_boxl_change() { dp3m_scaleby_box_l(); break; #endif -#ifdef SCAFACOS +#ifdef SCAFACOS_DIPOLES case DIPOLAR_SCAFACOS: - Scafacos::Dipoles::update_system_params(); + Scafacos::fcs_dipoles()->update_system_params(); break; #endif default: @@ -202,7 +202,7 @@ void calc_long_range_force(const ParticleRange &particles) { #endif #ifdef SCAFACOS_DIPOLES case DIPOLAR_SCAFACOS: - Scafacos::Dipoles::add_long_range_force(); + Scafacos::fcs_dipoles()->add_long_range_force(); #endif case DIPOLAR_NONE: break; @@ -248,7 +248,7 @@ double calc_energy_long_range(const ParticleRange &particles) { #endif #ifdef SCAFACOS_DIPOLES case DIPOLAR_SCAFACOS: - energy = Scafacos::Dipoles::long_range_energy(); + energy = Scafacos::fcs_dipoles()->long_range_energy(); break; #endif case DIPOLAR_NONE: diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp index 3d2579b1783..81c09e77823 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.cpp +++ b/src/core/electrostatics_magnetostatics/scafacos.cpp @@ -18,9 +18,6 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -/** @file - * Provide a C-like interface for ScaFaCoS. - */ #include "config.hpp" @@ -50,146 +47,91 @@ std::list available_methods() { return Scafacos::available_methods(); } -namespace Dipoles { - -static ScafacosContextDipoles *scafacos = nullptr; - -void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, - Utils::Vector3d &force) { - if (!scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - return scafacos->add_pair_force(q1q2, d, dist, force); -} - -double pair_energy(double q1q2, double dist) { - if (!scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - return scafacos->pair_energy(q1q2, dist); -} - -void add_long_range_force() { - if (!scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - scafacos->add_long_range_force(); -} +namespace { +#ifdef SCAFACOS_DIPOLES +ScafacosContextDipoles *dipoles_instance = nullptr; +#endif +ScafacosContextCoulomb *coulomb_instance = nullptr; +} // namespace -double long_range_energy() { - if (!scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - return scafacos->long_range_energy(); -} - -double get_r_cut() { - if (!scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - return scafacos->r_cut(); +#ifdef SCAFACOS_DIPOLES +ScafacosContextBase *fcs_dipoles() { + if (!dipoles_instance) { + throw std::runtime_error( + "Attempted access to uninitialized Scafacos Dipoles instance."); + } + return dipoles_instance; } +#endif -void update_system_params() { - if (!scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - scafacos->update_system_params(); +ScafacosContextBase *fcs_coulomb() { + if (!coulomb_instance) { + throw std::runtime_error( + "Attempted access to uninitialized Scafacos Coulomb instance."); + } + return coulomb_instance; } -static void set_parameters_worker(const std::string &method, - const std::string ¶ms) { - delete scafacos; - scafacos = nullptr; +#ifdef SCAFACOS_DIPOLES +static void set_parameters_dipoles_worker(const std::string &method, + const std::string ¶ms) { + delete dipoles_instance; + dipoles_instance = nullptr; - scafacos = new ScafacosContextDipoles(method, comm_cart, params); - if (!scafacos) { + auto *instance = new ScafacosContextDipoles(method, comm_cart, params); + if (!instance) { runtimeErrorMsg() << "Scafacos Dipoles failed to initialize"; return; } + dipoles_instance = instance; - scafacos->set_dipolar(true); - scafacos->update_system_params(); + instance->set_dipolar(true); + instance->update_system_params(); dipole.method = DIPOLAR_SCAFACOS; on_coulomb_change(); } -REGISTER_CALLBACK(set_parameters_worker) +REGISTER_CALLBACK(set_parameters_dipoles_worker) +#endif -} // namespace Dipoles +static void set_parameters_coulomb_worker(const std::string &method, + const std::string ¶ms) { + delete coulomb_instance; + coulomb_instance = nullptr; -namespace Coulomb { - -static ScafacosContextCoulomb *scafacos = nullptr; - -void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, - Utils::Vector3d &force) { - if (!scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - return scafacos->add_pair_force(q1q2, d, dist, force); -} - -double pair_energy(double q1q2, double dist) { - if (!scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - return scafacos->pair_energy(q1q2, dist); -} - -void add_long_range_force() { - if (!scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - scafacos->add_long_range_force(); -} - -double long_range_energy() { - if (!scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - return scafacos->long_range_energy(); -} - -double get_r_cut() { - if (!scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - return scafacos->r_cut(); -} - -void update_system_params() { - if (!scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - scafacos->update_system_params(); -} - -static void set_parameters_worker(const std::string &method, - const std::string ¶ms) { - delete scafacos; - scafacos = nullptr; - - scafacos = new ScafacosContextCoulomb(method, comm_cart, params); - if (!scafacos) { + auto *instance = new ScafacosContextCoulomb(method, comm_cart, params); + if (!instance) { runtimeErrorMsg() << "Scafacos Coulomb failed to initialize"; return; } + coulomb_instance = instance; - scafacos->set_dipolar(false); - scafacos->update_system_params(); + instance->set_dipolar(false); + instance->update_system_params(); coulomb.method = COULOMB_SCAFACOS; on_coulomb_change(); - scafacos->tune(); + instance->tune(); } -REGISTER_CALLBACK(set_parameters_worker) +REGISTER_CALLBACK(set_parameters_coulomb_worker) void set_r_cut_and_tune_local(double r_cut) { - scafacos->set_r_cut_and_tune_local(r_cut); + coulomb_instance->set_r_cut_and_tune_local(r_cut); } -} // namespace Coulomb - void free_handle(bool dipolar) { if (this_node == 0) mpi_call(free_handle, dipolar); if (dipolar) { - delete Dipoles::scafacos; - Dipoles::scafacos = nullptr; +#ifdef SCAFACOS_DIPOLES + delete dipoles_instance; + dipoles_instance = nullptr; +#endif } else { - delete Coulomb::scafacos; - Coulomb::scafacos = nullptr; + delete coulomb_instance; + coulomb_instance = nullptr; } } @@ -198,22 +140,23 @@ REGISTER_CALLBACK(free_handle) void set_parameters(const std::string &method, const std::string ¶ms, bool dipolar) { if (dipolar) { - mpi_call_all(Dipoles::set_parameters_worker, method, params); +#ifdef SCAFACOS_DIPOLES + mpi_call_all(set_parameters_dipoles_worker, method, params); +#endif } else { - mpi_call_all(Coulomb::set_parameters_worker, method, params); + mpi_call_all(set_parameters_coulomb_worker, method, params); } } std::string get_method_and_parameters(bool dipolar) { if (dipolar) { - if (!Dipoles::scafacos) - throw std::runtime_error("Scafacos Dipoles not initialized"); - return Dipoles::scafacos->get_method_and_parameters(); - } else { - if (!Coulomb::scafacos) - throw std::runtime_error("Scafacos Coulomb not initialized"); - return Coulomb::scafacos->get_method_and_parameters(); +#ifdef SCAFACOS_DIPOLES + return fcs_dipoles()->get_method_and_parameters(); +#else + return std::string(); +#endif } + return fcs_coulomb()->get_method_and_parameters(); } } // namespace Scafacos diff --git a/src/core/electrostatics_magnetostatics/scafacos.hpp b/src/core/electrostatics_magnetostatics/scafacos.hpp index c77d335040b..4587623e96c 100644 --- a/src/core/electrostatics_magnetostatics/scafacos.hpp +++ b/src/core/electrostatics_magnetostatics/scafacos.hpp @@ -21,64 +21,41 @@ /** \file * This file contains the c-type wrapper interface to the (oop-) scafacos - * interface. */ + * interface. + */ #ifndef ES_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOS_HPP #define ES_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOS_HPP #include "config.hpp" +#if defined(SCAFACOS) + +#include "electrostatics_magnetostatics/ScafacosContextBase.hpp" + #include #include #include namespace Scafacos { -#if defined(SCAFACOS) -namespace Dipoles { -/** Near-field pair force */ -void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, - Utils::Vector3d &force); -/** Near-field pair energy */ -double pair_energy(double q1q2, double dist); -/** Long range part */ -void add_long_range_force(); -/** Calculate long range energy contribution */ -double long_range_energy(); -/** Get parameters */ -std::string get_method_and_parameters(); -/** Set parameters */ -void set_parameters(const std::string &method, const std::string ¶ms); -double get_r_cut(); - -void update_system_params(); -} // namespace Dipoles - -namespace Coulomb { -/** Near-field pair force */ -void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, - Utils::Vector3d &force); -/** Near-field pair energy */ -double pair_energy(double q1q2, double dist); -/** Long range part */ -void add_long_range_force(); -/** Calculate long range energy contribution */ -double long_range_energy(); -double get_r_cut(); +/** @brief Access the per-MPI-node ScaFaCoS Coulomb instance */ +ScafacosContextBase *fcs_coulomb(); +#ifdef SCAFACOS_DIPOLES +/** @brief Access the per-MPI-node ScaFaCoS dipoles instance */ +ScafacosContextBase *fcs_dipoles(); +#endif -void update_system_params(); +std::string get_method_and_parameters(bool dipolar); +void set_parameters(const std::string &method, const std::string ¶ms, + bool dipolar); +void free_handle(bool dipolar); void set_r_cut_and_tune_local(double r_cut); -} // namespace Coulomb - -#endif /* SCAFACOS */ std::list available_methods(); -std::string get_method_and_parameters(bool dipolar); -void set_parameters(const std::string &method, const std::string ¶ms, - bool dipolar); -void free_handle(bool dipolar); } // namespace Scafacos +#endif // SCAFACOS #endif diff --git a/src/core/event.cpp b/src/core/event.cpp index ebe0b29e9e8..08f5af1372d 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -280,12 +280,12 @@ void on_parameter_change(int field) { #ifdef SCAFACOS #ifdef ELECTROSTATICS if (coulomb.method == COULOMB_SCAFACOS) { - Scafacos::Coulomb::update_system_params(); + Scafacos::fcs_coulomb()->update_system_params(); } #endif -#ifdef DIPOLES +#ifdef SCAFACOS_DIPOLES if (dipole.method == DIPOLAR_SCAFACOS) { - Scafacos::Dipoles::update_system_params(); + Scafacos::fcs_dipoles()->update_system_params(); } #endif #endif From c7e7df72ca7ebbda7f60425fa346cb85f97c925c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 10 Dec 2020 16:00:49 +0100 Subject: [PATCH 09/14] python: Allow two instances of ScaFaCoS --- doc/sphinx/electrostatics.rst | 4 +--- doc/sphinx/installation.rst | 2 +- doc/sphinx/magnetostatics.rst | 4 +--- src/python/espressomd/scafacos.pyx | 12 ------------ 4 files changed, 3 insertions(+), 19 deletions(-) diff --git a/doc/sphinx/electrostatics.rst b/doc/sphinx/electrostatics.rst index e31caff4bd7..df447c6d5d8 100644 --- a/doc/sphinx/electrostatics.rst +++ b/doc/sphinx/electrostatics.rst @@ -386,6 +386,4 @@ for an accuracy of :math:`10^{-3}`:: For details of the various methods and their parameters please refer to the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a -shared library. ScaFaCoS can be used only once, either for Coulomb or for -dipolar interactions. - +shared library. diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst index de728c46179..e21ebaf6ad6 100644 --- a/doc/sphinx/installation.rst +++ b/doc/sphinx/installation.rst @@ -303,7 +303,7 @@ General features :ref:`Magnetostatics / Dipolar interactions` :ref:`Electrostatics` -- ``SCAFACOS_DIPOLES`` +- ``SCAFACOS_DIPOLES`` This activates magnetostatics methods of ScaFaCoS. - ``ROTATION`` Switch on rotational degrees of freedom for the particles, as well as the corresponding quaternion integrator. diff --git a/doc/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst index 52fc41df91d..21fb78f8194 100644 --- a/doc/sphinx/magnetostatics.rst +++ b/doc/sphinx/magnetostatics.rst @@ -211,6 +211,4 @@ root mean square accuracy for the magnetic field. For details of the various methods and their parameters please refer to the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a -shared library. ScaFaCoS can be used only once, either for Coulomb or for -dipolar interactions. - +shared library. diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx index 99105219625..fa15288e371 100644 --- a/src/python/espressomd/scafacos.pyx +++ b/src/python/espressomd/scafacos.pyx @@ -108,18 +108,6 @@ IF SCAFACOS == 1: return res def _set_params_in_es_core(self): - # Verify that scafacos is not used for electrostatics and dipoles - # at the same time - IF ELECTROSTATICS == 1: - if self.dipolar and < int > electrostatics.coulomb.method == electrostatics.COULOMB_SCAFACOS: - raise Exception( - "Scafacos cannot be used for dipoles and charges at the same time") - - IF DIPOLES == 1: - if not self.dipolar and < int > magnetostatics.dipole.method == magnetostatics.DIPOLAR_SCAFACOS: - raise Exception( - "Scafacos cannot be used for dipoles and charges at the same time") - # Convert the parameter dictionary to a list of strings method_params = self._params["method_params"] param_string = "" From 38c73b5af389f1704ae5e2790a8ed0a7d391445b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 10 Dec 2020 16:48:46 +0100 Subject: [PATCH 10/14] python: Fix bug in getter --- src/python/espressomd/scafacos.pyx | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx index fa15288e371..d258c267e81 100644 --- a/src/python/espressomd/scafacos.pyx +++ b/src/python/espressomd/scafacos.pyx @@ -62,25 +62,22 @@ IF SCAFACOS == 1: # which we convert to a dict p = to_str(get_method_and_parameters(self.dipolar)).split(" ") res = {} - res["method_name"] = p[0] + res["method_name"] = p.pop(0) method_params = {} - i = 1 - while i < len(p): - pname = p[i] - i += 1 + while p: + pname = p.pop(0) # The first item after the parameter name is always a value # But in the case of array-like properties, there can also # follow several values. Therefore, we treat the next # words as part of the value, if they begin with a digit - pvalues = [p[i]] - i += 1 - if i >= len(p): - break - while p[i][:1] in "-0123456789": - pvalues.append(p[i]) - i += 1 + pvalues = [p.pop(0)] + while p: + if p[0][0] in "-0123456789": + pvalues.append(p.pop(0)) + else: + break # If there is one value, cast away the list if len(pvalues) == 1: @@ -89,8 +86,6 @@ IF SCAFACOS == 1: # Cast array elements to strings and join them by commas # to achieve consistency with setting array-likes # such as "pnfft_n":"128,128,128" - for j in range(len(pvalues)): - pvalues[j] = str(pvalues[j]) pvalues = ",".join(pvalues) method_params[pname] = pvalues From 9b65a905559f29c24bfca3768b60213f1c2b1b41 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 10 Dec 2020 16:56:47 +0100 Subject: [PATCH 11/14] tests: Add testcase for the ScaFaCoS interface --- testsuite/python/CMakeLists.txt | 1 + testsuite/python/scafacos_interface.py | 200 +++++++++++++++++++++++++ 2 files changed, 201 insertions(+) create mode 100644 testsuite/python/scafacos_interface.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index a40807f3964..42e659a6e68 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -109,6 +109,7 @@ python_test(FILE p3m_gpu.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE particle.py MAX_NUM_PROC 4) python_test(FILE pressure.py MAX_NUM_PROC 4) python_test(FILE scafacos_dipoles_1d_2d.py MAX_NUM_PROC 4) +python_test(FILE scafacos_interface.py MAX_NUM_PROC 4) python_test(FILE tabulated.py MAX_NUM_PROC 2) python_test(FILE particle_slice.py MAX_NUM_PROC 4) python_test(FILE rigid_bond.py MAX_NUM_PROC 4) diff --git a/testsuite/python/scafacos_interface.py b/testsuite/python/scafacos_interface.py new file mode 100644 index 00000000000..6b6c8a5b0c0 --- /dev/null +++ b/testsuite/python/scafacos_interface.py @@ -0,0 +1,200 @@ +# +# Copyright (C) 2020 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 by +# 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 numpy as np +import unittest as ut +import unittest_decorators as utx +import espressomd +import espressomd.electrostatics +import espressomd.magnetostatics +import espressomd.scafacos + + +@utx.skipIfMissingFeatures(["SCAFACOS"]) +class ScafacosInterface(ut.TestCase): + + system = espressomd.System(box_l=3 * [5]) + system.time_step = 0.01 + system.cell_system.skin = 0.5 + system.periodicity = [1, 1, 1] + + def tearDown(self): + self.system.part.clear() + self.system.actors.clear() + if self.system.actors: + # workaround for ScaFaCoS actor bug + self.system.actors.clear() + + def test_available_methods(self): + # all methods that are accessible from the ScaFaCoS bridge in ESPResSo + scafacos_methods = { + "direct", "ewald", "fmm", "memd", "mmm1d", "mmm2d", + "p2nfft", "p3m", "pepc", "pp3mg", "vmg", "wolf"} + + # all methods that were compiled in when building ScaFaCoS + available_methods = espressomd.scafacos.available_methods() + self.assertGreaterEqual(len(available_methods), 1) + for method in available_methods: + self.assertIn(method, scafacos_methods) + + def test_actor_coulomb(self): + system = self.system + + system.actors.add(espressomd.electrostatics.Scafacos( + prefactor=0.5, + method_name="p3m", + method_params={ + "p3m_r_cut": 1.0, + "p3m_grid": 32, + "p3m_cao": 7})) + actor = system.actors[0] + params = actor.get_params() + self.assertEqual(params["prefactor"], 0.5) + self.assertEqual(params["method_name"], "p3m") + self.assertEqual(params["method_params"], + {'p3m_cao': '7', 'p3m_r_cut': '1.0', 'p3m_grid': '32'}) + + @utx.skipIfMissingFeatures(["SCAFACOS_DIPOLES"]) + def test_actor_dipoles(self): + system = self.system + + method_params = { + "p2nfft_verbose_tuning": "0", + "pnfft_N": "32,32,32", + "pnfft_n": "32,32,32", + "pnfft_window_name": "bspline", + "pnfft_m": "4", + "p2nfft_ignore_tolerance": "1", + "pnfft_diff_ik": "0", + "p2nfft_r_cut": "11", + "p2nfft_alpha": "0.37"} + system.actors.add(espressomd.magnetostatics.Scafacos( + prefactor=1.2, + method_name="p2nfft", + method_params=method_params)) + actor = system.actors[0] + params = actor.get_params() + self.assertEqual(params["prefactor"], 1.2) + self.assertEqual(params["method_name"], "p2nfft") + self.assertEqual(params["method_params"], method_params) + + def p3m_data(self): + system = self.system + + p3m = espressomd.electrostatics.P3M( + prefactor=0.5, + accuracy=5e-4, + mesh=32, + cao=7, + r_cut=1.0) + system.actors.add(p3m) + + dp3m = espressomd.magnetostatics.DipolarP3M( + prefactor=1.0, + accuracy=1e-5, + cao=7, + mesh=48, + epsilon="metallic") + system.actors.add(dp3m) + + system.integrator.run(0, recalc_forces=True) + ref_E_coulomb = system.analysis.energy()["coulomb"] + ref_E_dipoles = system.analysis.energy()["dipolar"] + ref_forces = np.copy(system.part[:].f) + ref_torques = np.copy(system.part[:].torque_lab) + + system.actors.clear() + system.actors.clear() + + return (ref_E_coulomb, ref_E_dipoles, ref_forces, ref_torques) + + def fcs_data(self): + system = self.system + + scafacos_coulomb = espressomd.electrostatics.Scafacos( + prefactor=0.5, + method_name="p3m", + method_params={ + "p3m_r_cut": 1.0, + "p3m_grid": 32, + "p3m_cao": 7}) + system.actors.add(scafacos_coulomb) + + scafacos_dipoles = espressomd.magnetostatics.Scafacos( + prefactor=1.0, + method_name="p2nfft", + method_params={ + "p2nfft_verbose_tuning": 0, + "pnfft_N": "32,32,32", + "pnfft_n": "32,32,32", + "pnfft_window_name": "bspline", + "pnfft_m": "4", + "p2nfft_ignore_tolerance": "1", + "pnfft_diff_ik": "0", + "p2nfft_r_cut": "11", + "p2nfft_alpha": "0.37"}) + system.actors.add(scafacos_dipoles) + + system.integrator.run(0, recalc_forces=True) + ref_E_coulomb = system.analysis.energy()["coulomb"] + ref_E_dipoles = system.analysis.energy()["dipolar"] + ref_forces = np.copy(system.part[:].f) + ref_torques = np.copy(system.part[:].torque_lab) + + system.actors.clear() + if system.actors: + # workaround for ScaFaCoS actor bug + system.actors.clear() + + return (ref_E_coulomb, ref_E_dipoles, ref_forces, ref_torques) + + @utx.skipIfMissingFeatures(["SCAFACOS_DIPOLES", "LENNARD_JONES"]) + def test_electrostatics_plus_magnetostatics(self): + # check that two instances of ScaFaCoS can be used + system = self.system + + # add particles + N = 100 + np.random.seed(42) + system.part.add(pos=np.random.uniform(0, system.box_l[0], (N, 3)), + dip=np.random.uniform(0, 1, (N, 3)), + q=np.sign((np.arange(N) % 2) * 2 - 1), + rotation=N * [(1, 1, 1)]) + + # minimize system + system.non_bonded_inter[0, 0].lennard_jones.set_params( + epsilon=1.0, sigma=1.0, cutoff=2**(1.0 / 6.0), shift="auto") + system.integrator.set_steepest_descent( + f_max=10, + gamma=0.001, + max_displacement=0.01) + system.integrator.run(100) + + # compute forces and energies + p3m_E_coulomb, p3m_E_dipoles, p3m_forces, p3m_torques = self.p3m_data() + fcs_E_coulomb, fcs_E_dipoles, fcs_forces, fcs_torques = self.fcs_data() + + self.assertAlmostEqual(fcs_E_coulomb, p3m_E_coulomb, delta=1e-4) + self.assertAlmostEqual(fcs_E_dipoles, p3m_E_dipoles, delta=1e-4) + np.testing.assert_allclose(fcs_forces, p3m_forces, rtol=1e-3) + np.testing.assert_allclose(fcs_torques, p3m_torques, rtol=1e-3) + + +if __name__ == "__main__": + ut.main() From c6a15ae2954ad017502a54735465c369d6e7f9c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 10 Dec 2020 20:19:58 +0100 Subject: [PATCH 12/14] core: Inline the short-range energy/force kernels --- .../ScafacosContext.hpp | 17 ++++------------ .../ScafacosContextBase.hpp | 20 ++++++++++++++++--- 2 files changed, 21 insertions(+), 16 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp index 8d28f12681f..e22d1039b4a 100644 --- a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp +++ b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp @@ -52,20 +52,11 @@ struct ScafacosContext : ScafacosContextBase, Scafacos { using ScafacosContextBase::ScafacosContextBase; ~ScafacosContext() override = default; void update_system_params() override; - void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, - Utils::Vector3d &force) override { - if (dist > r_cut()) - return; - - auto const field = Scafacos::pair_force(dist); - auto const fak = q1q2 * field / dist; - force -= fak * d; + double get_pair_force(double dist) const override { + return Scafacos::pair_force(dist); } - double pair_energy(double q1q2, double dist) override { - if (dist > r_cut()) - return 0.; - - return q1q2 * Scafacos::pair_energy(dist); + double get_pair_energy(double dist) const override { + return Scafacos::pair_energy(dist); } std::string get_method_and_parameters() override; double get_r_cut() const override { return Scafacos::r_cut(); } diff --git a/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp b/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp index 7f8e2635cb6..b3f6b32a219 100644 --- a/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp +++ b/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp @@ -55,13 +55,27 @@ struct ScafacosContextBase { /** @brief Add long-range part of the forces to particles. */ virtual void add_long_range_force() = 0; /** @brief Add near-field pair force. */ - virtual void add_pair_force(double q1q2, Utils::Vector3d const &d, - double dist, Utils::Vector3d &force) = 0; + inline void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist, + Utils::Vector3d &force) { + if (dist > get_r_cut()) + return; + + auto const field = get_pair_force(dist); + auto const fak = q1q2 * field / dist; + force -= fak * d; + } /** @brief Calculate near-field pair energy. */ - virtual double pair_energy(double q1q2, double dist) = 0; + inline double pair_energy(double q1q2, double dist) { + if (dist > get_r_cut()) + return 0.; + + return q1q2 * get_pair_energy(dist); + } /** @brief Reinitialize number of particles, box shape and periodicity. */ virtual void update_system_params() = 0; virtual double get_r_cut() const = 0; + virtual double get_pair_force(double dist) const = 0; + virtual double get_pair_energy(double dist) const = 0; virtual std::string get_method_and_parameters() = 0; }; From 1f9bc11d0231ed10a250f2a97e5bd3f2d7ad7135 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 10 Dec 2020 21:49:01 +0100 Subject: [PATCH 13/14] tests: Avoid DP3M bisection error on 3 MPI ranks --- testsuite/python/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 42e659a6e68..2f9aba81870 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -109,7 +109,7 @@ python_test(FILE p3m_gpu.py MAX_NUM_PROC 4 LABELS gpu) python_test(FILE particle.py MAX_NUM_PROC 4) python_test(FILE pressure.py MAX_NUM_PROC 4) python_test(FILE scafacos_dipoles_1d_2d.py MAX_NUM_PROC 4) -python_test(FILE scafacos_interface.py MAX_NUM_PROC 4) +python_test(FILE scafacos_interface.py MAX_NUM_PROC 2) python_test(FILE tabulated.py MAX_NUM_PROC 2) python_test(FILE particle_slice.py MAX_NUM_PROC 4) python_test(FILE rigid_bond.py MAX_NUM_PROC 4) From 691847ae7c5b774a982ee03d133b3bf7a61b5732 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 10 Dec 2020 21:50:20 +0100 Subject: [PATCH 14/14] scafacos: Fix undefined behavior UBSAN report: "reference binding to null pointer of type 'double'". When the system has no particles, the ScaFaCoS input vectors are empty and cannot be dereferenced. --- src/scafacos/Scafacos.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/scafacos/Scafacos.cpp b/src/scafacos/Scafacos.cpp index 2bd1da3fafb..881a7518be4 100644 --- a/src/scafacos/Scafacos.cpp +++ b/src/scafacos/Scafacos.cpp @@ -124,10 +124,11 @@ void Scafacos::run(std::vector &charges, std::vector &positions, fields.resize(3 * local_n_part); potentials.resize(local_n_part); - handle_error(fcs_tune(handle, local_n_part, &(positions[0]), &(charges[0]))); + handle_error( + fcs_tune(handle, local_n_part, positions.data(), charges.data())); - handle_error(fcs_run(handle, local_n_part, &(positions[0]), &(charges[0]), - &(fields[0]), &(potentials[0]))); + handle_error(fcs_run(handle, local_n_part, positions.data(), charges.data(), + fields.data(), potentials.data())); } #ifdef FCS_ENABLE_DIPOLES @@ -142,8 +143,8 @@ void Scafacos::run_dipolar(std::vector &dipoles, potentials.resize(3 * local_n_part); handle_error(fcs_set_dipole_particles(handle, static_cast(local_n_part), - &(positions[0]), &(dipoles[0]), - &(fields[0]), &(potentials[0]))); + positions.data(), dipoles.data(), + fields.data(), potentials.data())); handle_error(fcs_run(handle, 0, nullptr, nullptr, nullptr, nullptr)); } #endif @@ -151,7 +152,7 @@ void Scafacos::run_dipolar(std::vector &dipoles, void Scafacos::tune(std::vector &charges, std::vector &positions) { handle_error( - fcs_tune(handle, charges.size(), &(positions[0]), &(charges[0]))); + fcs_tune(handle, charges.size(), positions.data(), charges.data())); } void Scafacos::set_common_parameters(const double *box_l,