ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
dipolar_direct_sum.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#include "config/config.hpp"
23
24#ifdef DIPOLES
25
27
28#include "BoxGeometry.hpp"
29#include "cells.hpp"
30#include "communication.hpp"
31#include "errorhandling.hpp"
32#include "system/System.hpp"
33
35#include <utils/constants.hpp>
36#include <utils/math/sqr.hpp>
38
39#include <boost/mpi/collectives.hpp>
40#include <boost/mpi/communicator.hpp>
41#include <boost/range/counting_range.hpp>
42
43#include <mpi.h>
44
45#include <algorithm>
46#include <cassert>
47#include <cmath>
48#include <iterator>
49#include <stdexcept>
50#include <tuple>
51#include <utility>
52#include <vector>
53
54/**
55 * @brief Pair force of two interacting dipoles.
56 *
57 * @param d Distance vector.
58 * @param m1 Dipole moment of one particle.
59 * @param m2 Dipole moment of the other particle.
60 *
61 * @return Resulting force.
62 */
63static auto pair_force(Utils::Vector3d const &d, Utils::Vector3d const &m1,
64 Utils::Vector3d const &m2) {
65 auto const pe2 = m1 * d;
66 auto const pe3 = m2 * d;
67
68 auto const r2 = d.norm2();
69 auto const r = std::sqrt(r2);
70 auto const r5 = r2 * r2 * r;
71 auto const r7 = r5 * r2;
72
73 auto const a = 3.0 * (m1 * m2) / r5;
74 auto const b = -15.0 * pe2 * pe3 / r7;
75
76 auto const f = (a + b) * d + 3.0 * (pe3 * m1 + pe2 * m2) / r5;
77 auto const r3 = r2 * r;
78 auto const t =
79 -vector_product(m1, m2) / r3 + 3.0 * pe3 * vector_product(m1, d) / r5;
80
81 return ParticleForce{f, t};
82}
83
84/**
85 * @brief Pair potential for two interacting dipoles.
86 *
87 * @param d Distance vector.
88 * @param m1 Dipole moment of one particle.
89 * @param m2 Dipole moment of the other particle.
90 *
91 * @return Interaction energy.
92 */
93static auto pair_potential(Utils::Vector3d const &d, Utils::Vector3d const &m1,
94 Utils::Vector3d const &m2) {
95 auto const r2 = d * d;
96 auto const r = sqrt(r2);
97 auto const r3 = r2 * r;
98 auto const r5 = r3 * r2;
99
100 auto const pe1 = m1 * m2;
101 auto const pe2 = m1 * d;
102 auto const pe3 = m2 * d;
103
104 return pe1 / r3 - 3.0 * pe2 * pe3 / r5;
105}
106
107/**
108 * @brief Dipole field contribution from a particle with dipole moment @c m1
109 * at a distance @c d.
110 *
111 * @param d Distance vector.
112 * @param m1 Dipole moment of one particle.
113 *
114 * @return Utils::Vector3d containing dipole field components.
115 */
116static auto dipole_field(Utils::Vector3d const &d, Utils::Vector3d const &m1) {
117 auto const r2 = d * d;
118 auto const r = sqrt(r2);
119 auto const r3 = r2 * r;
120 auto const r5 = r3 * r2;
121 auto const pe2 = m1 * d;
122
123 return 3.0 * pe2 * d / r5 - m1 / r3;
124}
125
126/**
127 * @brief Call kernel for every 3d index in a sphere around the origin.
128 *
129 * This calls a callable for all index-triples
130 * that are within ball around the origin with
131 * radius |ncut|.
132 *
133 * @tparam F Callable
134 * @param ncut Limits in the three directions,
135 * all non-zero elements have to be
136 * the same number.
137 * @param f will be called for each index triple
138 * within the limits of @p ncut.
139 */
140template <typename F> void for_each_image(Utils::Vector3i const &ncut, F f) {
141 auto const ncut2 = ncut.norm2();
142
143 /* This runs over the index "cube"
144 * [-ncut[0], ncut[0]] x ... x [-ncut[2], ncut[2]]
145 * (inclusive on both sides), and calls f with
146 * all the elements as argument. Counting range
147 * is a range that just enumerates a range.
148 */
150 [&](int nx, int ny, int nz) {
151 if (nx * nx + ny * ny + nz * nz <= ncut2) {
152 f(nx, ny, nz);
153 }
154 },
155 boost::counting_range(-ncut[0], ncut[0] + 1),
156 boost::counting_range(-ncut[1], ncut[1] + 1),
157 boost::counting_range(-ncut[2], ncut[2] + 1));
158}
159
160/**
161 * @brief Position and dipole moment of one particle.
162 */
163struct PosMom {
166
167 template <class Archive> void serialize(Archive &ar, long int) {
168 ar & pos & m;
169 }
170};
171
172/**
173 * @brief Sum over all pairs with periodic images.
174 *
175 * This implements the "primed" pair sum, the sum over all
176 * pairs between one particles and all other particles,
177 * including @p ncut periodic replicas in each direction.
178 * Primed means that in the primary replica the self-interaction
179 * is excluded, but not with the other periodic replicas. E.g.
180 * a particle does not interact with its self, but does with
181 * its periodically shifted versions.
182 *
183 * @param begin Iterator pointing to begin of particle range
184 * @param end Iterator pointing past the end of particle range
185 * @param it Pointer to particle that is considered
186 * @param with_replicas If periodic replicas are to be considered
187 * at all. If false, distances are calculated as Euclidean
188 * distances, and not using minimum image convention.
189 * @param ncut Number of replicas in each direction.
190 * @param box_geo Box geometry.
191 * @param init Initial value of the sum.
192 * @param f Binary operation mapping distance and moment of the
193 * interaction partner to the value to be summed up for this pair.
194 *
195 * @return The total sum.
196 */
197template <class InputIterator, class T, class F>
198T image_sum(InputIterator begin, InputIterator end, InputIterator it,
199 bool with_replicas, Utils::Vector3i const &ncut,
200 BoxGeometry const &box_geo, T init, F f) {
201
202 auto const &box_l = box_geo.length();
203 for (auto jt = begin; jt != end; ++jt) {
204 auto const exclude_primary = (it == jt);
205 auto const primary_distance = (with_replicas)
206 ? (it->pos - jt->pos)
207 : box_geo.get_mi_vector(it->pos, jt->pos);
208
209 for_each_image(ncut, [&](int nx, int ny, int nz) {
210 if (!(exclude_primary && nx == 0 && ny == 0 && nz == 0)) {
211 auto const rn =
212 primary_distance +
213 Utils::Vector3d{nx * box_l[0], ny * box_l[1], nz * box_l[2]};
214 init += f(rn, jt->m);
215 }
216 });
217 }
218
219 return init;
220}
221
222static auto gather_particle_data(BoxGeometry const &box_geo,
223 ParticleRange const &particles) {
224 auto const &comm = ::comm_cart;
225 std::vector<Particle *> local_particles;
226 std::vector<PosMom> local_posmom;
227 std::vector<PosMom> all_posmom;
228 std::vector<boost::mpi::request> reqs;
229
230 local_particles.reserve(particles.size());
231 local_posmom.reserve(particles.size());
232
233 for (auto &p : particles) {
234 if (p.dipm() != 0.0) {
235 local_particles.emplace_back(&p);
236 local_posmom.emplace_back(
237 PosMom{box_geo.folded_position(p.pos()), p.calc_dip()});
238 }
239 }
240
241 auto const local_size = static_cast<int>(local_posmom.size());
242 std::vector<int> all_sizes;
243 boost::mpi::all_gather(comm, local_size, all_sizes);
244
245 auto const offset =
246 std::accumulate(all_sizes.begin(), all_sizes.begin() + comm.rank(), 0);
247 auto const total_size =
248 std::accumulate(all_sizes.begin() + comm.rank(), all_sizes.end(), offset);
249
250 if (comm.size() > 1) {
251 all_posmom.resize(total_size);
252 reqs = Utils::Mpi::iall_gatherv(comm, local_posmom.data(), local_size,
253 all_posmom.data(), all_sizes.data());
254 } else {
255 std::swap(all_posmom, local_posmom);
256 }
257
258 return std::make_tuple(std::move(local_particles), std::move(all_posmom),
259 std::move(reqs), offset);
260}
261
262static auto get_n_cut(BoxGeometry const &box_geo, int n_replicas) {
263 return n_replicas * Utils::Vector3i{static_cast<int>(box_geo.periodic(0)),
264 static_cast<int>(box_geo.periodic(1)),
265 static_cast<int>(box_geo.periodic(2))};
266}
267
268/**
269 * @brief Calculate and add the interaction forces/torques to the particles.
270 *
271 * This employs a parallel N-square loop over all particle pairs.
272 * The computation the partitioned into several steps so that the
273 * communication latency can be hidden behinder some local computation:
274 *
275 * 1. The local particle positions and momenta are packed into
276 * one array.
277 * 2. The asynchronous distribution of the local arrays to all
278 * ranks is started.
279 * 3. The interaction for the local pairs is started, here every
280 * pair is visited only once, and the force is added to both particles.
281 * 4. Wait for the data from the other nodes.
282 * 5. Calculate the interaction with the rest of the particles. Here
283 * every pair is visited twice (not necessarily on the same rank)
284 * so that no reduction of the forces is needed.
285 *
286 * Logically this is equivalent to the potential calculation
287 * in @ref DipolarDirectSum::long_range_energy, which calculates
288 * a naive N-square sum, but has better performance and scaling.
289 */
291 ParticleRange const &particles) const {
292 auto const &box_geo = *get_system().box_geo;
293 auto const &box_l = box_geo.length();
294 auto [local_particles, all_posmom, reqs, offset] =
295 gather_particle_data(box_geo, particles);
296
297 /* Number of image boxes considered */
298 auto const ncut = get_n_cut(box_geo, n_replicas);
299 auto const with_replicas = (ncut.norm2() > 0);
300
301 /* Range of particles we calculate the ia for on this node */
302 auto const local_posmom_begin = all_posmom.begin() + offset;
303 auto const local_posmom_end =
304 local_posmom_begin + static_cast<long>(local_particles.size());
305
306 /* Output iterator for the force */
307 auto p = local_particles.begin();
308
309 /* IA with local particles */
310 for (auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) {
311 /* IA with own images */
312 auto fi = image_sum(
313 it, std::next(it), it, with_replicas, ncut, box_geo, ParticleForce{},
314 [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) {
315 return pair_force(rn, it->m, mj);
316 });
317
318 /* IA with other local particles */
319 auto q = std::next(p);
320 for (auto jt = std::next(it); jt != local_posmom_end; ++jt, ++q) {
321 auto const d = (with_replicas) ? (it->pos - jt->pos)
322 : box_geo.get_mi_vector(it->pos, jt->pos);
323
324 ParticleForce fij{};
325 ParticleForce fji{};
326 for_each_image(ncut, [&](int nx, int ny, int nz) {
327 auto const rn =
328 d + Utils::Vector3d{nx * box_l[0], ny * box_l[1], nz * box_l[2]};
329 auto const pf = pair_force(rn, it->m, jt->m);
330 fij += pf;
331 fji.f -= pf.f;
332 /* Conservation of angular momentum mandates that
333 * 0 = t_i + r_ij x F_ij + t_j */
334 fji.torque += vector_product(pf.f, rn) - pf.torque;
335 });
336
337 fi += fij;
338 (*q)->force() += prefactor * fji.f;
339 (*q)->torque() += prefactor * fji.torque;
340 }
341
342 (*p)->force() += prefactor * fi.f;
343 (*p)->torque() += prefactor * fi.torque;
344 }
345
346 /* Wait for the rest of the data to arrive */
347 boost::mpi::wait_all(reqs.begin(), reqs.end());
348
349 /* Output iterator for the force */
350 p = local_particles.begin();
351
352 /* Interaction with all the other particles */
353 for (auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) {
354 // red particles
355 auto fi =
356 image_sum(all_posmom.begin(), local_posmom_begin, it, with_replicas,
357 ncut, box_geo, ParticleForce{},
358 [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) {
359 return pair_force(rn, it->m, mj);
360 });
361
362 // black particles
363 fi += image_sum(local_posmom_end, all_posmom.end(), it, with_replicas, ncut,
364 box_geo, ParticleForce{},
365 [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) {
366 return pair_force(rn, it->m, mj);
367 });
368
369 (*p)->force() += prefactor * fi.f;
370 (*p)->torque() += prefactor * fi.torque;
371 }
372}
373
374/**
375 * @brief Calculate the interaction potential.
376 *
377 * This employs a parallel N-square loop over all particle pairs.
378 */
379double
381 auto const &box_geo = *get_system().box_geo;
382 auto [local_particles, all_posmom, reqs, offset] =
383 gather_particle_data(box_geo, particles);
384
385 /* Number of image boxes considered */
386 auto const ncut = get_n_cut(box_geo, n_replicas);
387 auto const with_replicas = (ncut.norm2() > 0);
388
389 /* Wait for the rest of the data to arrive */
390 boost::mpi::wait_all(reqs.begin(), reqs.end());
391
392 /* Range of particles we calculate the ia for on this node */
393 auto const local_posmom_begin = all_posmom.begin() + offset;
394 auto const local_posmom_end =
395 local_posmom_begin + static_cast<long>(local_particles.size());
396
397 auto u = 0.;
398 for (auto it = local_posmom_begin; it != local_posmom_end; ++it) {
399 u = image_sum(it, all_posmom.end(), it, with_replicas, ncut, box_geo, u,
400 [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) {
401 return pair_potential(rn, it->m, mj);
402 });
403 }
404
405 return prefactor * u;
406}
407
408/**
409 * @brief Calculate total dipole field of each particle.
410 *
411 * This employs a parallel N-square loop over all particles.
412 * Logically this is equivalent to the potential calculation
413 * in @ref DipolarDirectSum::long_range_energy, which calculates
414 * a naive N-square sum. The difference is summation range,
415 * and the kernel calculates the dipole field rather than the energy.
416 */
417#ifdef DIPOLE_FIELD_TRACKING
419 ParticleRange const &particles) const {
420 auto const &box_geo = *get_system().box_geo;
421 /* collect particle data */
422 auto [local_particles, all_posmom, reqs, offset] =
423 gather_particle_data(box_geo, particles);
424
425 auto const ncut = get_n_cut(box_geo, n_replicas);
426 auto const with_replicas = (ncut.norm2() > 0);
427
428 boost::mpi::wait_all(reqs.begin(), reqs.end());
429
430 auto const local_posmom_begin = all_posmom.begin() + offset;
431 auto const local_posmom_end =
432 local_posmom_begin + static_cast<long>(local_particles.size());
433
434 Utils::Vector3d u_init = {0., 0., 0.};
435 auto p = local_particles.begin();
436 for (auto pi = local_posmom_begin; pi != local_posmom_end; ++pi, ++p) {
437 auto const u = image_sum(
438 all_posmom.begin(), all_posmom.end(), pi, with_replicas, ncut, box_geo,
439 u_init, [](Utils::Vector3d const &rn, Utils::Vector3d const &mj) {
440 return dipole_field(rn, mj);
441 });
442 (*p)->dip_fld() = prefactor * u;
443 }
444}
445#endif
446
447DipolarDirectSum::DipolarDirectSum(double prefactor, int n_replicas) {
449 this->n_replicas = n_replicas;
450 if (n_replicas < 0) {
451 throw std::domain_error("Parameter 'n_replicas' must be >= 0");
452 }
453}
454
455#endif // DIPOLES
float f[3]
float u[3]
This file contains everything related to the global cell structure / cell system.
Utils::Vector3d const & length() const
Box length.
constexpr bool periodic(unsigned coord) const
Check periodicity in direction.
auto folded_position(Utils::Vector3d const &p) const
Calculate coordinates folded to primary simulation box.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
void set_prefactor(double new_prefactor)
double prefactor
Magnetostatics prefactor.
A range of particles.
base_type::size_type size() const
T norm2() const
Definition Vector.hpp:130
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
T image_sum(InputIterator begin, InputIterator end, InputIterator it, bool with_replicas, Utils::Vector3i const &ncut, BoxGeometry const &box_geo, T init, F f)
Sum over all pairs with periodic images.
void for_each_image(Utils::Vector3i const &ncut, F f)
Call kernel for every 3d index in a sphere around the origin.
static auto dipole_field(Utils::Vector3d const &d, Utils::Vector3d const &m1)
Dipole field contribution from a particle with dipole moment m1 at a distance d.
static auto pair_potential(Utils::Vector3d const &d, Utils::Vector3d const &m1, Utils::Vector3d const &m2)
Pair potential for two interacting dipoles.
static auto gather_particle_data(BoxGeometry const &box_geo, ParticleRange const &particles)
static auto get_n_cut(BoxGeometry const &box_geo, int n_replicas)
static auto pair_force(Utils::Vector3d const &d, Utils::Vector3d const &m1, Utils::Vector3d const &m2)
Pair force of two interacting dipoles.
__device__ void vector_product(float const *a, float const *b, float *out)
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
auto iall_gatherv(boost::mpi::communicator const &comm, T const *in_values, int in_size, T *out_values, int const *sizes)
void cartesian_product(const Body &op, const ForwardRange &...rng)
Call op with each element of the cartesian product set of rng.
void dipole_field_at_part(ParticleRange const &particles) const
Calculate total dipole field of each particle.
DipolarDirectSum(double prefactor, int n_replicas)
void add_long_range_forces(ParticleRange const &particles) const
Calculate and add the interaction forces/torques to the particles.
double long_range_energy(ParticleRange const &particles) const
Calculate the interaction potential.
Force information on a particle.
Definition Particle.hpp:290
Utils::Vector3d f
force.
Definition Particle.hpp:314
Position and dipole moment of one particle.
void serialize(Archive &ar, long int)
Utils::Vector3d m
Utils::Vector3d pos
DEVICE_QUALIFIER constexpr iterator begin() noexcept
Definition Array.hpp:128