30#include "communication.hpp"
32#include "system/System.hpp"
39#include <boost/mpi/collectives.hpp>
40#include <boost/mpi/communicator.hpp>
41#include <boost/range/counting_range.hpp>
65 auto const pe2 = m1 * d;
66 auto const pe3 = m2 * d;
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;
73 auto const a = 3.0 * (m1 * m2) / r5;
74 auto const b = -15.0 * pe2 * pe3 / r7;
76 auto const f = (a + b) * d + 3.0 * (pe3 * m1 + pe2 * m2) / r5;
77 auto const r3 = r2 * r;
95 auto const r2 = d * d;
96 auto const r = sqrt(r2);
97 auto const r3 = r2 * r;
98 auto const r5 = r3 * r2;
100 auto const pe1 = m1 * m2;
101 auto const pe2 = m1 * d;
102 auto const pe3 = m2 * d;
104 return pe1 / r3 - 3.0 * pe2 * pe3 / r5;
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;
123 return 3.0 * pe2 * d / r5 - m1 / r3;
141 auto const ncut2 = ncut.
norm2();
150 [&](
int nx,
int ny,
int nz) {
151 if (nx * nx + ny * ny + nz * nz <= ncut2) {
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));
167 template <
class Archive>
void serialize(Archive &ar,
long int) {
197template <
class InputIterator,
class T,
class F>
198T
image_sum(InputIterator begin, InputIterator end, InputIterator it,
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)
210 if (!(exclude_primary && nx == 0 && ny == 0 && nz == 0)) {
214 init +=
f(rn, jt->m);
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;
230 local_particles.reserve(particles.
size());
231 local_posmom.reserve(particles.
size());
233 for (
auto &p : particles) {
234 if (p.dipm() != 0.0) {
235 local_particles.emplace_back(&p);
236 local_posmom.emplace_back(
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);
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);
250 if (comm.size() > 1) {
251 all_posmom.resize(total_size);
253 all_posmom.data(), all_sizes.data());
255 std::swap(all_posmom, local_posmom);
258 return std::make_tuple(std::move(local_particles), std::move(all_posmom),
259 std::move(reqs), offset);
264 static_cast<int>(box_geo.
periodic(1)),
265 static_cast<int>(box_geo.
periodic(2))};
293 auto const &box_l = box_geo.length();
294 auto [local_particles, all_posmom, reqs, offset] =
299 auto const with_replicas = (ncut.norm2() > 0);
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());
307 auto p = local_particles.begin();
310 for (
auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) {
313 it, std::next(it), it, with_replicas, ncut, box_geo,
ParticleForce{},
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);
339 (*q)->torque() +=
prefactor * fji.torque;
347 boost::mpi::wait_all(reqs.begin(), reqs.end());
350 p = local_particles.begin();
353 for (
auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) {
356 image_sum(all_posmom.begin(), local_posmom_begin, it, with_replicas,
359 return pair_force(rn, it->m, mj);
363 fi +=
image_sum(local_posmom_end, all_posmom.end(), it, with_replicas, ncut,
366 return pair_force(rn, it->m, mj);
382 auto [local_particles, all_posmom, reqs, offset] =
387 auto const with_replicas = (ncut.norm2() > 0);
390 boost::mpi::wait_all(reqs.begin(), reqs.end());
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());
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,
401 return pair_potential(rn, it->m, mj);
417#ifdef DIPOLE_FIELD_TRACKING
422 auto [local_particles, all_posmom, reqs, offset] =
426 auto const with_replicas = (ncut.norm2() > 0);
428 boost::mpi::wait_all(reqs.begin(), reqs.end());
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());
435 auto p = local_particles.
begin();
436 for (
auto pi = local_posmom_begin; pi != local_posmom_end; ++pi, ++p) {
438 all_posmom.begin(), all_posmom.end(), pi, with_replicas, ncut, box_geo,
440 return dipole_field(rn, mj);
451 throw std::domain_error(
"Parameter 'n_replicas' must be >= 0");
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.
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