ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
dp3m.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2024 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/** @file
23 * P3M algorithm for long-range magnetic dipole-dipole interaction.
24 *
25 * We use here a P3M (Particle-Particle Particle-Mesh) method based
26 * on the dipolar Ewald summation. Details of the used method can be
27 * found in @cite hockney88a and @cite deserno98a @cite deserno98b.
28 *
29 * Further reading: @cite cerda08d
30 */
31
32#pragma once
33
34#include "config/config.hpp"
35
36#ifdef DP3M
37
39
40#include "p3m/common.hpp"
41#include "p3m/data_struct.hpp"
42
43#include "Particle.hpp"
44#include "ParticleRange.hpp"
45
46#include <utils/Vector.hpp>
48#include <utils/math/sqr.hpp>
49
50#include <cmath>
51#include <numbers>
52
53#ifdef NPT
54/** Update the NpT virial */
55void npt_add_virial_magnetic_contribution(double energy);
56#endif
57
58/** @brief Dipolar P3M solver. */
59struct DipolarP3M : public Dipoles::Actor<DipolarP3M> {
61
62public:
64
65 virtual ~DipolarP3M() = default;
66
67 [[nodiscard]] virtual bool is_tuned() const noexcept = 0;
68 [[nodiscard]] virtual bool is_gpu() const noexcept = 0;
69 [[nodiscard]] virtual bool is_double_precision() const noexcept = 0;
70
71 virtual void on_activation() = 0;
72 /** @brief Recalculate all box-length-dependent parameters. */
80 /** @brief Recalculate all derived parameters. */
81 virtual void init() = 0;
82
89
90 /**
91 * @brief Count the number of magnetic particles and calculate
92 * the sum of the squared dipole moments.
93 */
94 virtual void count_magnetic_particles() = 0;
95
96 /** Assign the physical dipoles using the tabulated assignment function. */
97 virtual void dipole_assign(ParticleRange const &particles) = 0;
98
99 /**
100 * @brief Tune dipolar P3M parameters to desired accuracy.
101 *
102 * The parameters
103 * @ref P3MParameters::mesh "mesh",
104 * @ref P3MParameters::cao "cao",
105 * @ref P3MParameters::r_cut_iL "r_cut_iL" and
106 * @ref P3MParameters::alpha_L "alpha_L" are tuned to obtain the target
107 * @ref P3MParameters::accuracy "accuracy" in optimal time.
108 * These parameters are stored in the @ref dp3m_params object.
109 *
110 * The function utilizes the analytic expression of the error estimate
111 * for the dipolar P3M method in the paper of @cite cerda08d in
112 * order to obtain the rms error in the force for a system of N randomly
113 * distributed particles in a cubic box. For the real space error, the
114 * estimate in @cite kolafa92a is used.
115 *
116 * Parameter ranges if not given explicitly in the constructor:
117 * - @p mesh is set up such that the number of mesh points is equal to the
118 * number of magnetic dipolar particles
119 * - @p cao explores all possible values
120 * - @p alpha_L is tuned for each tuple (@p r_cut_iL, @p mesh, @p cao) and
121 * calculated assuming that the error contributions of real and reciprocal
122 * space should be equal
123 *
124 * After checking if the total error lies below the target accuracy,
125 * the time needed for one force calculation is measured. Parameters
126 * that minimize the runtime are kept.
127 *
128 * The function is based on routines of the program HE_Q.cpp for charges
129 * written by M. Deserno.
130 */
131 virtual void tune() = 0;
132
133 /** Compute the k-space part of energies. */
134 virtual double long_range_energy(ParticleRange const &particles) = 0;
135
136 /** Compute the k-space part of forces. */
137 virtual void add_long_range_forces(ParticleRange const &particles) = 0;
138
139 /** Calculate real-space contribution of p3m dipolar pair forces and torques.
140 * If NPT is compiled in, update the NpT virial.
141 */
142 inline ParticleForce pair_force(Particle const &p1, Particle const &p2,
143 Utils::Vector3d const &d, double dist2,
144 double dist) const {
145 if ((p1.dipm() == 0.) || (p2.dipm() == 0.) || dist >= dp3m_params.r_cut ||
146 dist <= 0.)
147 return {};
148
149 auto const dip1 = p1.calc_dip();
150 auto const dip2 = p2.calc_dip();
151 auto const alpsq = dp3m_params.alpha * dp3m_params.alpha;
152 auto const adist = dp3m_params.alpha * dist;
153#if USE_ERFC_APPROXIMATION
154 auto const erfc_part_ri = Utils::AS_erfc_part(adist) / dist;
155#else
156 auto const erfc_part_ri = erfc(adist) / dist;
157#endif
158
159 // Calculate scalar multiplications for vectors mi, mj, rij
160 auto const mimj = dip1 * dip2;
161
162 auto const mir = dip1 * d;
163 auto const mjr = dip2 * d;
164
165 auto const coeff = 2. * dp3m_params.alpha * std::numbers::inv_sqrtpi;
166 auto const dist2i = 1. / dist2;
167 auto const exp_adist2 = exp(-Utils::sqr(adist));
168
169 auto const B_r = (dp3m_params.accuracy > 5e-06)
170 ? (erfc_part_ri + coeff) * exp_adist2 * dist2i
171 : (erfc(adist) / dist + coeff * exp_adist2) * dist2i;
172
173 auto const common_term = alpsq * coeff * exp_adist2;
174 auto const C_r = dist2i * (3. * B_r + 2. * common_term);
175 auto const D_r = dist2i * (5. * C_r + 4. * common_term * alpsq);
176
177 // Calculate real-space forces
178 auto const force = prefactor * ((mimj * d + dip1 * mjr + dip2 * mir) * C_r -
179 mir * mjr * D_r * d);
180
181 // Calculate vector multiplications for vectors mi, mj, rij
182 auto const mixmj = vector_product(dip1, dip2);
183 auto const mixr = vector_product(dip1, d);
184
185 // Calculate real-space torques
186 auto const torque = prefactor * (-mixmj * B_r + mixr * (mjr * C_r));
187#ifdef NPT
188#if USE_ERFC_APPROXIMATION
189 auto const fac = prefactor * p1.dipm() * p2.dipm() * exp_adist2;
190#else
191 auto const fac = prefactor * p1.dipm() * p2.dipm();
192#endif
193 auto const energy = fac * (mimj * B_r - mir * mjr * C_r);
195#endif // NPT
196 return ParticleForce{force, torque};
197 }
198
199 /** Calculate real-space contribution of dipolar pair energy. */
200 inline double pair_energy(Particle const &p1, Particle const &p2,
201 Utils::Vector3d const &d, double dist2,
202 double dist) const {
203 if ((p1.dipm() == 0.) || (p2.dipm() == 0.) || dist >= dp3m_params.r_cut ||
204 dist <= 0.)
205 return {};
206
207 auto const dip1 = p1.calc_dip();
208 auto const dip2 = p2.calc_dip();
209
210 auto const alpsq = dp3m_params.alpha * dp3m_params.alpha;
211 auto const adist = dp3m_params.alpha * dist;
212
213#if USE_ERFC_APPROXIMATION
214 auto const erfc_part_ri = Utils::AS_erfc_part(adist) / dist;
215#else
216 auto const erfc_part_ri = erfc(adist) / dist;
217#endif
218
219 // Calculate scalar multiplications for vectors mi, mj, rij
220 auto const mimj = dip1 * dip2;
221 auto const mir = dip1 * d;
222 auto const mjr = dip2 * d;
223
224 auto const coeff = 2. * dp3m_params.alpha * std::numbers::inv_sqrtpi;
225 auto const dist2i = 1. / dist2;
226 auto const exp_adist2 = exp(-Utils::sqr(adist));
227
228 auto const B_r = (dp3m_params.accuracy > 5e-06)
229 ? dist2i * (erfc_part_ri + coeff) * exp_adist2
230 : dist2i * (erfc(adist) / dist + coeff * exp_adist2);
231 auto const C_r = (3. * B_r + 2. * alpsq * coeff * exp_adist2) * dist2i;
232
233 return prefactor * (mimj * B_r - mir * mjr * C_r);
234 }
235
236protected:
237 /** Calculate self-energy in k-space. */
238 virtual double calc_average_self_energy_k_space() const = 0;
239
240 /** Calculate energy correction that minimizes the error.
241 * This quantity is only calculated once and is cached.
242 */
243 virtual void calc_energy_correction() = 0;
244
245 /** Calculate the influence function for the dipolar forces and torques. */
247
248 /** Calculate the influence function for the dipolar energy. */
250
251 /** Compute the dipolar surface terms */
252 virtual double calc_surface_term(bool force_flag, bool energy_flag,
253 ParticleRange const &particles) = 0;
254
255 /** Checks for correctness of the k-space cutoff. */
256 void sanity_checks_boxl() const;
257 void sanity_checks_node_grid() const;
258 void sanity_checks_periodicity() const;
259 void sanity_checks_cell_structure() const;
260
261 virtual void scaleby_box_l() = 0;
262};
263
264#endif // DP3M
Vector implementation and trait types for boost qvm interoperability.
double prefactor
Magnetostatics prefactor.
A range of particles.
This file contains the defaults for ESPResSo.
__device__ void vector_product(float const *a, float const *b, float *out)
void npt_add_virial_magnetic_contribution(double energy)
Update the NpT virial.
Definition dp3m.cpp:903
constexpr T AS_erfc_part(T d)
Approximate by applying a formula from chapter 7.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Common functions for dipolar and charge P3M.
Dipolar P3M solver.
Definition dp3m.hpp:59
void sanity_checks_boxl() const
Checks for correctness of the k-space cutoff.
Definition dp3m.cpp:822
virtual void init()=0
Recalculate all derived parameters.
virtual bool is_gpu() const noexcept=0
void sanity_checks() const
Definition dp3m.hpp:83
virtual void scaleby_box_l()=0
virtual void tune()=0
Tune dipolar P3M parameters to desired accuracy.
virtual void calc_energy_correction()=0
Calculate energy correction that minimizes the error.
void on_node_grid_change() const
Definition dp3m.hpp:74
virtual bool is_double_precision() const noexcept=0
void on_boxl_change()
Recalculate all box-length-dependent parameters.
Definition dp3m.hpp:73
void sanity_checks_cell_structure() const
Definition dp3m.cpp:856
virtual double long_range_energy(ParticleRange const &particles)=0
Compute the k-space part of energies.
virtual ~DipolarP3M()=default
virtual void dipole_assign(ParticleRange const &particles)=0
Assign the physical dipoles using the tabulated assignment function.
P3MParameters const & dp3m_params
Definition dp3m.hpp:60
virtual void on_activation()=0
virtual void calc_influence_function_force()=0
Calculate the influence function for the dipolar forces and torques.
virtual double calc_average_self_energy_k_space() const =0
Calculate self-energy in k-space.
virtual bool is_tuned() const noexcept=0
void sanity_checks_periodicity() const
Definition dp3m.cpp:848
DipolarP3M(P3MParameters const &dp3m_params)
Definition dp3m.hpp:63
void sanity_checks_node_grid() const
Definition dp3m.cpp:871
double pair_energy(Particle const &p1, Particle const &p2, Utils::Vector3d const &d, double dist2, double dist) const
Calculate real-space contribution of dipolar pair energy.
Definition dp3m.hpp:200
void on_periodicity_change() const
Definition dp3m.hpp:75
virtual void count_magnetic_particles()=0
Count the number of magnetic particles and calculate the sum of the squared dipole moments.
virtual double calc_surface_term(bool force_flag, bool energy_flag, ParticleRange const &particles)=0
Compute the dipolar surface terms.
void on_cell_structure_change()
Definition dp3m.hpp:76
virtual void calc_influence_function_energy()=0
Calculate the influence function for the dipolar energy.
virtual void add_long_range_forces(ParticleRange const &particles)=0
Compute the k-space part of forces.
ParticleForce pair_force(Particle const &p1, Particle const &p2, Utils::Vector3d const &d, double dist2, double dist) const
Calculate real-space contribution of p3m dipolar pair forces and torques.
Definition dp3m.hpp:142
Structure to hold P3M parameters and some dependent variables.
double alpha
unscaled alpha_L for use with fast inline functions only
double accuracy
accuracy of the actual parameter set.
double r_cut
unscaled r_cut_iL for use with fast inline functions only
Force information on a particle.
Definition Particle.hpp:290
Struct holding all information for one particle.
Definition Particle.hpp:395
auto calc_dip() const
Definition Particle.hpp:495
auto const & dipm() const
Definition Particle.hpp:493