36#include "communication.hpp"
38#include "system/System.hpp"
42#include <boost/mpi/collectives/all_reduce.hpp>
43#include <boost/mpi/collectives/reduce.hpp>
44#include <boost/mpi/operations.hpp>
59void DipolarLayerCorrection::check_gap(
Particle const &p)
const {
61 auto const z = p.
pos()[2];
64 <<
"by " << z - ((z < 0.) ? 0. :
dlc.box_h);
71 auto const local_particles = system.
cell_structure->local_particles();
72 auto const mu_max_local = std::accumulate(
73 local_particles.begin(), local_particles.end(), 0.,
74 [](
double mu,
Particle const &p) { return std::max(mu, p.dipm()); });
76 return boost::mpi::all_reduce(
comm_cart, mu_max_local,
77 boost::mpi::maximum<double>());
82 auto const cc2 = c * c;
83 auto const x3 = x * x * x;
84 auto const a = g * g * g / x + 1.5 * cc2 + 1.5 * g / x3 + 0.75 / (x3 * x);
89 auto const x2 = x * x;
90 auto const a = g * g / x + 2.0 * g / x2 + 2.0 / (x2 * x);
98 for (
auto const &p : particles) {
104 return boost::mpi::all_reduce(
comm_cart, local_dip, std::plus<>());
112 std::vector<Utils::Vector3d> &fs,
113 std::vector<Utils::Vector3d> &ts,
116 auto const facux = 2. * std::numbers::pi * box_geo.
length_inv()[0];
117 auto const facuy = 2. * std::numbers::pi * box_geo.
length_inv()[1];
119 auto const n_local_particles = particles.
size();
120 std::vector<double> ReSjp(n_local_particles);
121 std::vector<double> ReSjm(n_local_particles);
122 std::vector<double> ImSjp(n_local_particles);
123 std::vector<double> ImSjm(n_local_particles);
124 std::vector<double> ReGrad_Mup(n_local_particles);
125 std::vector<double> ImGrad_Mup(n_local_particles);
126 std::vector<double> ReGrad_Mum(n_local_particles);
127 std::vector<double> ImGrad_Mum(n_local_particles);
129 for (
int ix = -kcut; ix <= +kcut; ix++) {
130 for (
int iy = -kcut; iy <= +kcut; iy++) {
131 if (ix == 0 and iy == 0) {
134 auto const gx =
static_cast<double>(ix) * facux;
135 auto const gy =
static_cast<double>(iy) * facuy;
136 auto const gr = sqrt(gx * gx + gy * gy);
139 auto const fa1 = 1. / (gr * (exp(gr * box_geo.
length()[2]) - 1.));
144 double S[4] = {0., 0., 0., 0.};
146 for (
auto const &p : particles) {
147 if (p.
dipm() != 0.) {
148 auto const &pos = p.
pos();
151 auto const a = gx * dip[0] + gy * dip[1];
152 auto const b = gr * dip[2];
153 auto const er = gx * pos[0] + gy * pos[1];
154 auto const ez = gr * pos[2];
155 auto const c = cos(er);
156 auto const d = sin(er);
157 auto const f = exp(ez);
159 ReSjp[ip] = (b * c - a * d) * f;
160 ImSjp[ip] = (c * a + b * d) * f;
161 ReSjm[ip] = (-b * c - a * d) / f;
162 ImSjm[ip] = (c * a - b * d) / f;
163 ReGrad_Mup[ip] = c * f;
164 ReGrad_Mum[ip] = c / f;
165 ImGrad_Mup[ip] = d * f;
166 ImGrad_Mum[ip] = d / f;
176 MPI_Allreduce(MPI_IN_PLACE, S, 4, MPI_DOUBLE, MPI_SUM,
comm_cart);
181 for (
auto &p : particles) {
182 if (p.
dipm() != 0.) {
185 auto const s1 = -(-ReSjp[ip] * S[3] + ImSjp[ip] * S[2]);
186 auto const s2 = +(ReSjm[ip] * S[1] - ImSjm[ip] * S[0]);
187 auto const s3 = -(-ReSjm[ip] * S[1] + ImSjm[ip] * S[0]);
188 auto const s4 = +(ReSjp[ip] * S[3] - ImSjp[ip] * S[2]);
190 auto const s1z = +(ReSjp[ip] * S[2] + ImSjp[ip] * S[3]);
191 auto const s2z = -(ReSjm[ip] * S[0] + ImSjm[ip] * S[1]);
192 auto const s3z = -(ReSjm[ip] * S[0] + ImSjm[ip] * S[1]);
193 auto const s4z = +(ReSjp[ip] * S[2] + ImSjp[ip] * S[3]);
195 auto const ss = s1 + s2 + s3 + s4;
196 fs[ip][0] += fa1 * gx * ss;
197 fs[ip][1] += fa1 * gy * ss;
198 fs[ip][2] += fa1 * gr * (s1z + s2z + s3z + s4z);
202 auto const s1 = -(-ReGrad_Mup[ip] * S[3] + ImGrad_Mup[ip] * S[2]);
203 auto const s2 = +(ReGrad_Mum[ip] * S[1] - ImGrad_Mum[ip] * S[0]);
204 auto const s3 = -(-ReGrad_Mum[ip] * S[1] + ImGrad_Mum[ip] * S[0]);
205 auto const s4 = +(ReGrad_Mup[ip] * S[3] - ImGrad_Mup[ip] * S[2]);
207 auto const s1z = +(ReGrad_Mup[ip] * S[2] + ImGrad_Mup[ip] * S[3]);
208 auto const s2z = -(ReGrad_Mum[ip] * S[0] + ImGrad_Mum[ip] * S[1]);
209 auto const s3z = -(ReGrad_Mum[ip] * S[0] + ImGrad_Mum[ip] * S[1]);
210 auto const s4z = +(ReGrad_Mup[ip] * S[2] + ImGrad_Mup[ip] * S[3]);
212 auto const ss = s1 + s2 + s3 + s4;
213 ts[ip][0] += fa1 * gx * ss;
214 ts[ip][1] += fa1 * gy * ss;
215 ts[ip][2] += fa1 * gr * (s1z + s2z + s3z + s4z);
226 for (
auto const &p : particles) {
227 if (p.
dipm() != 0.) {
237 for (std::size_t j = 0; j < n_local_particles; ++j) {
250 auto const facux = 2. * std::numbers::pi * box_geo.
length_inv()[0];
251 auto const facuy = 2. * std::numbers::pi * box_geo.
length_inv()[1];
254 double sum_S[4] = {0., 0., 0., 0.};
255 for (
int ix = -kcut; ix <= +kcut; ix++) {
256 for (
int iy = -kcut; iy <= +kcut; iy++) {
257 if ((ix == 0 && iy == 0)) {
260 auto const gx =
static_cast<double>(ix) * facux;
261 auto const gy =
static_cast<double>(iy) * facuy;
262 auto const gr = sqrt(gx * gx + gy * gy);
265 auto const fa1 = 1. / (gr * (exp(gr * box_geo.
length()[2]) - 1.));
269 double S[4] = {0., 0., 0., 0.};
271 for (
auto const &p : particles) {
272 if (p.
dipm() != 0.) {
273 auto const &pos = p.
pos();
276 auto const a = gx * dip[0] + gy * dip[1];
277 auto const b = gr * dip[2];
278 auto const er = gx * pos[0] + gy * pos[1];
279 auto const ez = gr * pos[2];
280 auto const c = cos(er);
281 auto const d = sin(er);
282 auto const f = exp(ez);
284 S[0] += (b * c - a * d) * f;
285 S[1] += (c * a + b * d) * f;
286 S[2] += (-b * c - a * d) / f;
287 S[3] += (c * a - b * d) / f;
290 boost::mpi::reduce(
comm_cart, S, 4, sum_S, std::plus<>(), 0);
294 energy += fa1 * 2. * (sum_S[0] * sum_S[2] + sum_S[1] * sum_S[3]);
308 auto const volume = box_geo.volume();
309 auto const correc = 4. * std::numbers::pi / volume;
313 std::vector<Utils::Vector3d> dip_DLC_f(particles.
size());
314 std::vector<Utils::Vector3d> dip_DLC_t(particles.
size());
319 auto const kcut =
static_cast<int>(std::round(
dlc.
far_cut));
335 for (
auto &p : particles) {
338 if (p.
dipm() != 0.) {
360 auto const volume = box_geo.volume();
361 auto const pref =
prefactor * 2. * std::numbers::pi / volume;
366 for (
auto const &p : particles) {
373 auto const k_cut =
static_cast<int>(std::round(
dlc.
far_cut));
374 auto dip_DLC_energy =
387 dip_DLC_energy += pref *
Utils::sqr(box_dip[2]);
393 return dip_DLC_energy;
402 if (p.
dipm() != 0.) {
407 return boost::mpi::all_reduce(
comm_cart, local_n, std::plus<>());
416double DipolarLayerCorrection::tune_far_cut()
const {
421 auto const &box_geo = *system.box_geo;
422 auto const lx = box_geo.length()[0];
423 auto const ly = box_geo.length()[1];
424 auto const lz = box_geo.length()[2];
426 if (std::abs(lx - ly) > 0.001) {
427 throw std::runtime_error(
"DLC tuning: box size in x direction is "
428 "different from y direction. The tuning "
429 "formula requires both to be equal.");
432 auto constexpr limitkc = 200;
433 auto const piarea = std::numbers::pi / (lx * ly);
437 for (
int kc = 1; kc < limitkc; kc++) {
438 auto const gc = kc * 2. * std::numbers::pi / lx;
439 auto const fa0 =
sqrt(9. * exp(+2. * gc * h) *
g1_DLC_dip(gc, lz - h) +
440 9. * exp(-2. * gc * h) *
g1_DLC_dip(gc, lz + h) +
442 auto const fa1 =
sqrt(0.125 * piarea) * fa0;
444 auto const de = nmp * mu_max_sq / (4. * (exp(gc * lz) - 1.)) * (fa1 + fa2);
446 far_cut =
static_cast<double>(kc);
451 throw std::runtime_error(
"DLC tuning failed: maxPWerror too small");
460 void operator()(std::shared_ptr<DipolarDirectSum>
const &solver) {
468 m_actor->
epsilon = solver->dp3m_params.epsilon;
483void DipolarLayerCorrection::sanity_checks_node_grid()
const {
485 if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) {
486 throw std::runtime_error(
"DLC: requires periodicity (True, True, True)");
491 : maxPWerror{maxPWerror}, gap_size{gap_size}, box_h{-1.}, far_cut{far_cut},
492 far_calculated{far_cut == -1.} {
494 throw std::domain_error(
"Parameter 'far_cut' must be > 0");
497 throw std::domain_error(
"Parameter 'maxPWerror' must be > 0");
500 throw std::domain_error(
"Parameter 'gap_size' must be > 0");
505 auto const lz =
get_system().box_geo->length()[2];
507 if (new_box_h < 0.) {
508 throw std::runtime_error(
"DLC gap size (" + std::to_string(
dlc.
gap_size) +
509 ") larger than box length in z-direction (" +
510 std::to_string(lz) +
")");
517 : dlc{parameters}, base_solver{solver} {
Utils::Vector3d const & length() const
Box length.
Utils::Vector3d const & length_inv() const
Inverse box length.
base_type::size_type size() const
std::shared_ptr< CellStructure > cell_structure
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
This file contains the defaults for ESPResSo.
__device__ void vector_product(float const *a, float const *b, float *out)
static void dipolar_force_corrections(int kcut, std::vector< Utils::Vector3d > &fs, std::vector< Utils::Vector3d > &ts, ParticleRange const &particles, BoxGeometry const &box_geo)
Compute the dipolar force and torque corrections.
static double g1_DLC_dip(double g, double x)
static double dipolar_energy_correction(int kcut, ParticleRange const &particles, BoxGeometry const &box_geo)
Compute the dipolar DLC energy correction.
static double g2_DLC_dip(double g, double x)
static Utils::Vector3d calc_slab_dipole(ParticleRange const &particles)
Compute the box magnetic dipole.
static double calc_mu_max(System::System const &system)
Calculate the maximal dipole moment in the system.
static int count_magnetic_particles(System::System const &system)
P3M algorithm for long-range magnetic dipole-dipole interaction.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Vector< T, N > sqrt(Vector< T, N > const &a)
Common functions for dipolar and charge P3M.
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
Lock an actor and modify its parameters.
AdaptSolver(DipolarLayerCorrection *this_ptr)
void operator()(std::shared_ptr< DipolarDirectSum > const &solver)
void operator()(std::shared_ptr< DipolarP3M > const &solver)
Adapt a magnetostatics solver to remove contributions from the z-direction.
double epsilon_correction
double energy_correction(ParticleRange const &particles) const
Calculate the dipolar energy correction.
BaseSolver base_solver
Magnetostatics solver that is adapted.
void add_force_corrections(ParticleRange const &particles) const
Add the dipolar force and torque corrections.
DipolarLayerCorrection(dlc_data &¶meters, BaseSolver &&solver)
std::variant< std::shared_ptr< DipolarP3M >, std::shared_ptr< DipolarDirectSum > > BaseSolver
Struct holding all information for one particle.
auto const & torque() const
auto const & dipm() const
auto const & force() const
Parameters for the DLC method.
bool far_calculated
Flag whether far_cut was set by the user, or calculated by ESPResSo.
dlc_data(double maxPWerror, double gap_size, double far_cut)
double maxPWerror
maximal pairwise error of the potential and force
double far_cut
Cutoff of the exponential sum.
double gap_size
Size of the empty gap.
double box_h
Up to where particles can be found.