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