30#include "communication.hpp"
32#include "system/System.hpp"
38#include <boost/mpi/collectives.hpp>
39#include <boost/mpi/communicator.hpp>
64 auto const pe2 = m1 * d;
65 auto const pe3 = m2 * d;
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;
72 auto const a = 3.0 * (m1 * m2) / r5;
73 auto const b = -15.0 * pe2 * pe3 / r7;
75 auto const f = (a + b) * d + 3.0 * (pe3 * m1 + pe2 * m2) / r5;
76 auto const r3 = r2 * r;
94 auto const r2 = d * d;
95 auto const r = sqrt(r2);
96 auto const r3 = r2 * r;
97 auto const r5 = r3 * r2;
99 auto const pe1 = m1 * m2;
100 auto const pe2 = m1 * d;
101 auto const pe3 = m2 * d;
103 return pe1 / r3 - 3.0 * pe2 * pe3 / r5;
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;
122 return 3.0 * pe2 * d / r5 - m1 / r3;
140 auto const ncut2 = ncut.
norm2();
149 [&](
int nx,
int ny,
int nz) {
150 if (nx * nx + ny * ny + nz * nz <= ncut2) {
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));
166 template <
class Archive>
void serialize(Archive &ar,
long int) {
196template <
class InputIterator,
class T,
class F>
197T
image_sum(InputIterator begin, InputIterator end, InputIterator it,
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)
209 if (!(exclude_primary && nx == 0 && ny == 0 && nz == 0)) {
213 init += f(rn, jt->m);
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;
229 local_particles.reserve(particles.
size());
230 local_posmom.reserve(particles.
size());
232 for (
auto &p : particles) {
233 if (p.dipm() != 0.0) {
234 local_particles.emplace_back(&p);
235 local_posmom.emplace_back(
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);
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);
249 if (comm.size() > 1) {
250 all_posmom.resize(total_size);
252 all_posmom.data(), all_sizes.data());
254 std::swap(all_posmom, local_posmom);
257 return std::make_tuple(std::move(local_particles), std::move(all_posmom),
258 std::move(reqs), offset);
263 static_cast<int>(box_geo.
periodic(1)),
264 static_cast<int>(box_geo.
periodic(2))};
292 auto const &box_l = box_geo.length();
293 auto [local_particles, all_posmom, reqs, offset] =
298 auto const with_replicas = (ncut.norm2() > 0);
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());
306 auto p = local_particles.begin();
309 for (
auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) {
312 it, std::next(it), it, with_replicas, ncut, box_geo,
ParticleForce{},
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);
338 (*q)->torque() +=
prefactor * fji.torque;
346 boost::mpi::wait_all(reqs.begin(), reqs.end());
349 p = local_particles.begin();
352 for (
auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) {
355 image_sum(all_posmom.begin(), local_posmom_begin, it, with_replicas,
358 return pair_force(rn, it->m, mj);
362 fi +=
image_sum(local_posmom_end, all_posmom.end(), it, with_replicas, ncut,
365 return pair_force(rn, it->m, mj);
381 auto [local_particles, all_posmom, reqs, offset] =
386 auto const with_replicas = (ncut.norm2() > 0);
389 boost::mpi::wait_all(reqs.begin(), reqs.end());
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());
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,
400 return pair_potential(rn, it->m, mj);
416#ifdef DIPOLE_FIELD_TRACKING
421 auto [local_particles, all_posmom, reqs, offset] =
425 auto const with_replicas = (ncut.norm2() > 0);
427 boost::mpi::wait_all(reqs.begin(), reqs.end());
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());
434 auto p = local_particles.
begin();
435 for (
auto pi = local_posmom_begin; pi != local_posmom_end; ++pi, ++p) {
437 all_posmom.begin(), all_posmom.end(), pi, with_replicas, ncut, box_geo,
439 return dipole_field(rn, mj);
450 throw std::domain_error(
"Parameter 'n_replicas' must be >= 0");
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.
base_type::size_type size() const
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.
Position and dipole moment of one particle.
void serialize(Archive &ar, long int)
DEVICE_QUALIFIER constexpr iterator begin() noexcept