50#include "communication.hpp"
54#include "system/System.hpp"
62#include <boost/mpi/collectives/all_reduce.hpp>
63#include <boost/mpi/collectives/reduce.hpp>
80#error "The FFTW3 library shouldn't be visible in this translation unit"
83template <
typename FloatType, Arch Architecture>
86 double local_mu2 = 0.;
88 for (
auto const &p : get_system().
cell_structure->local_particles()) {
90 local_mu2 += p.calc_dip().norm2();
95 boost::mpi::all_reduce(
comm_cart, local_mu2, dp3m.sum_mu2, std::plus<>());
96 boost::mpi::all_reduce(
comm_cart, local_n, dp3m.sum_dip_part, std::plus<>());
100 int n_c_part,
double sum_q2,
double alpha_L);
103 int n_c_part,
double sum_q2,
110 double sum_q2,
double x1,
double x2,
double xacc,
111 double tuned_accuracy);
113template <
typename FloatType, Arch Architecture>
117 auto const &box_geo = *get_system().
box_geo;
119 dp3m.params, dp3m.mesh.start, dp3m.mesh.stop, dp3m.g_energy);
122 boost::mpi::reduce(
comm_cart, node_phi, phi, std::plus<>(), 0);
123 phi /= 3. * box_geo.length()[0] * Utils::int_pow<3>(dp3m.params.mesh[0]);
124 return phi * std::numbers::pi;
127template <
typename FloatType, Arch Architecture>
130 assert(dp3m.params.cao >= 1 and dp3m.params.cao <= 7);
131 assert(dp3m.params.alpha > 0.);
133 auto const &system = get_system();
134 auto const &box_geo = *system.box_geo;
135 auto const &local_geo = *system.local_geo;
136 auto const verlet_skin = system.cell_structure->get_verlet_skin();
138 dp3m.params.cao3 = Utils::int_pow<3>(dp3m.params.cao);
139 dp3m.params.recalc_a_ai_cao_cut(box_geo.length());
142 dp3m.local_mesh.calc_local_ca_mesh(dp3m.params, local_geo, verlet_skin, 0.);
143 dp3m.fft_buffers->init_halo();
144 dp3m.fft->init(dp3m.params);
145 dp3m.mesh.ks_pnum = dp3m.fft->get_ks_pnum();
146 dp3m.fft_buffers->init_meshes(dp3m.fft->get_ca_mesh_size());
147 dp3m.update_mesh_views();
148 dp3m.calc_differential_operator();
161 typename std::remove_reference_t<
decltype(dp3m)>::value_type;
162 auto const weights = p3m_calculate_interpolation_weights<cao>(
163 real_pos, dp3m.params.ai, dp3m.local_mesh);
164 p3m_interpolate<cao>(
165 dp3m.local_mesh, weights, [&dip, &dp3m](
int ind,
double w) {
166 dp3m.mesh.rs_fields[0u][ind] += value_type(w * dip[0u]);
167 dp3m.mesh.rs_fields[1u][ind] += value_type(w * dip[1u]);
168 dp3m.mesh.rs_fields[2u][ind] += value_type(w * dip[2u]);
171 dp3m.inter_weights.template store<cao>(weights);
176template <
typename FloatType, Arch Architecture>
179 dp3m.inter_weights.reset(dp3m.params.cao);
182 for (
auto &rs_mesh_field : dp3m.mesh.rs_fields)
183 for (
int j = 0; j < dp3m.local_mesh.size; j++)
184 rs_mesh_field[j] = 0.;
186 for (
auto const &p : particles) {
187 if (p.dipm() != 0.) {
188 Utils::integral_parameter<int, AssignDipole, 1, 7>(dp3m.params.cao, dp3m,
189 p.pos(), p.calc_dip());
200 auto p_index = std::size_t{0ul};
202 for (
auto &p : particles) {
203 if (p.dipm() != 0.) {
204 auto const weights = dp3m.inter_weights.template load<cao>(p_index);
208 [&E, &dp3m, d_rs](
int ind,
double w) {
209 E[d_rs] += w * double(dp3m.mesh.rs_scalar[ind]);
224 auto p_index = std::size_t{0ul};
226 for (
auto &p : particles) {
227 if (p.dipm() != 0.) {
228 auto const weights = dp3m.inter_weights.template load<cao>(p_index);
232 [&E, &dp3m](
int ind,
double w) {
233 E[0u] += w * double(dp3m.mesh.rs_fields[0u][ind]);
234 E[1u] += w * double(dp3m.mesh.rs_fields[1u][ind]);
235 E[2u] += w * double(dp3m.mesh.rs_fields[2u][ind]);
238 p.force()[d_rs] += p.calc_dip() * prefac * E;
246template <
typename FloatType, Arch Architecture>
248 bool force_flag,
bool energy_flag,
ParticleRange const &particles) {
251 auto const &system = get_system();
252 auto const &box_geo = *system.box_geo;
253 auto const dipole_prefac = prefactor / Utils::int_pow<3>(dp3m.params.mesh[0]);
255 auto const npt_flag =
258 auto constexpr npt_flag =
false;
261 if (dp3m.sum_mu2 > 0.) {
262 dipole_assign(particles);
263 dp3m.fft_buffers->perform_vector_halo_gather();
264 for (
auto &rs_mesh : dp3m.fft_buffers->get_vector_mesh()) {
265 dp3m.fft->forward_fft(rs_mesh);
267 dp3m.update_mesh_views();
271 if (energy_flag or npt_flag) {
275 if (dp3m.sum_mu2 > 0.) {
280 auto const &offset = dp3m.mesh.start;
281 auto const &d_op = dp3m.d_op[0u];
282 auto const &mesh_dip = dp3m.mesh.rs_fields;
283 auto const [KX, KY, KZ] = dp3m.fft->get_permutations();
285 auto index = std::size_t(0u);
286 auto it_energy = dp3m.g_energy.begin();
287 auto node_energy = 0.;
288 for_each_3d(mesh_start, dp3m.mesh.size, indices, [&]() {
289 auto const shift = indices + offset;
291 auto const re = mesh_dip[0u][index] * FloatType(d_op[shift[KX]]) +
292 mesh_dip[1u][index] * FloatType(d_op[shift[KY]]) +
293 mesh_dip[2u][index] * FloatType(d_op[shift[KZ]]);
296 auto const im = mesh_dip[0u][index] * FloatType(d_op[shift[KX]]) +
297 mesh_dip[1u][index] * FloatType(d_op[shift[KY]]) +
298 mesh_dip[2u][index] * FloatType(d_op[shift[KZ]]);
300 node_energy += *it_energy * (Utils::sqr(re) + Utils::sqr(im));
301 std::advance(it_energy, 1);
304 node_energy *= dipole_prefac * std::numbers::pi * box_geo.length_inv()[0];
305 boost::mpi::reduce(
comm_cart, node_energy, energy, std::plus<>(), 0);
307 if (dp3m.energy_correction == 0.)
308 calc_energy_correction();
312 energy -= prefactor * dp3m.sum_mu2 * std::numbers::inv_sqrtpi *
313 (2. / 3.) * Utils::int_pow<3>(dp3m.params.alpha);
316 energy += prefactor * dp3m.energy_correction / box_geo.volume();
326 if (dp3m.sum_mu2 > 0.) {
327 auto const wavenumber = 2. * std::numbers::pi * box_geo.length_inv()[0u];
329 auto const &mesh_stop = dp3m.mesh.size;
330 auto const offset = dp3m.mesh.start;
331 auto const &d_op = dp3m.d_op[0u];
332 auto const [KX, KY, KZ] = dp3m.fft->get_permutations();
333 auto &mesh_dip = dp3m.mesh.rs_fields;
335 auto index = std::size_t(0u);
336 dp3m.ks_scalar.resize(dp3m.mesh.rs_scalar.size());
339 auto it_energy = dp3m.g_energy.begin();
341 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
342 auto const shift = indices + offset;
344 auto const re = mesh_dip[0u][index] * FloatType(d_op[shift[KX]]) +
345 mesh_dip[1u][index] * FloatType(d_op[shift[KY]]) +
346 mesh_dip[2u][index] * FloatType(d_op[shift[KZ]]);
347 dp3m.ks_scalar[index] = *it_energy * re;
350 auto const im = mesh_dip[0u][index] * FloatType(d_op[shift[KX]]) +
351 mesh_dip[1u][index] * FloatType(d_op[shift[KY]]) +
352 mesh_dip[2u][index] * FloatType(d_op[shift[KZ]]);
353 dp3m.ks_scalar[index] = *it_energy * im;
355 std::advance(it_energy, 1);
359 for (
int d = 0; d < 3; d++) {
361 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
362 auto const d_op_val = FloatType(d_op[indices[d] + offset[d]]);
363 dp3m.mesh.rs_scalar[index] = d_op_val * dp3m.ks_scalar[index];
365 dp3m.mesh.rs_scalar[index] = d_op_val * dp3m.ks_scalar[index];
368 dp3m.fft->backward_fft(dp3m.fft_buffers->get_scalar_mesh());
369 dp3m.fft_buffers->perform_scalar_halo_spread();
371 auto const d_rs = (d + dp3m.mesh.ks_pnum) % 3;
372 Utils::integral_parameter<int, AssignTorques, 1, 7>(
373 dp3m.params.cao, dp3m, dipole_prefac * wavenumber, d_rs, particles);
384 auto it_force = dp3m.g_force.begin();
386 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
387 auto const shift = indices + offset;
389 auto const re = mesh_dip[0u][index] * FloatType(d_op[shift[KX]]) +
390 mesh_dip[1u][index] * FloatType(d_op[shift[KY]]) +
391 mesh_dip[2u][index] * FloatType(d_op[shift[KZ]]);
394 auto const im = mesh_dip[0u][index] * FloatType(d_op[shift[KX]]) +
395 mesh_dip[1u][index] * FloatType(d_op[shift[KY]]) +
396 mesh_dip[2u][index] * FloatType(d_op[shift[KZ]]);
398 dp3m.ks_scalar[index - 2] = *it_force * im;
399 dp3m.ks_scalar[index - 1] = *it_force * (-re);
400 std::advance(it_force, 1);
404 for (
int d = 0; d < 3; d++) {
406 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
407 auto const d_op_val = FloatType(d_op[indices[d] + offset[d]]);
408 auto const shift = indices + offset;
409 auto const f1 = d_op_val * dp3m.ks_scalar[index];
410 mesh_dip[0u][index] = FloatType(d_op[shift[KX]]) * f1;
411 mesh_dip[1u][index] = FloatType(d_op[shift[KY]]) * f1;
412 mesh_dip[2u][index] = FloatType(d_op[shift[KZ]]) * f1;
414 auto const f2 = d_op_val * dp3m.ks_scalar[index];
415 mesh_dip[0u][index] = FloatType(d_op[shift[KX]]) * f2;
416 mesh_dip[1u][index] = FloatType(d_op[shift[KY]]) * f2;
417 mesh_dip[2u][index] = FloatType(d_op[shift[KZ]]) * f2;
420 for (
auto &rs_mesh : dp3m.fft_buffers->get_vector_mesh()) {
421 dp3m.fft->backward_fft(rs_mesh);
423 dp3m.fft_buffers->perform_vector_halo_spread();
425 auto const d_rs = (d + dp3m.mesh.ks_pnum) % 3;
426 Utils::integral_parameter<int, AssignForces, 1, 7>(
427 dp3m.params.cao, dp3m, dipole_prefac *
Utils::sqr(wavenumber), d_rs,
434 auto const surface_term =
435 calc_surface_term(force_flag, energy_flag or npt_flag, particles);
437 energy += surface_term;
442 fprintf(stderr,
"dipolar_P3M at this moment is added to p_vir[0]\n");
444 if (not energy_flag) {
451template <
typename FloatType, Arch Architecture>
453 bool force_flag,
bool energy_flag,
ParticleRange const &particles) {
454 auto const &box_geo = *get_system().
box_geo;
455 auto const pref = prefactor * 4. * std::numbers::pi / box_geo.volume() /
456 (2. * dp3m.params.epsilon + 1.);
457 auto const n_local_part = particles.
size();
461 std::vector<double> mx(n_local_part);
462 std::vector<double> my(n_local_part);
463 std::vector<double> mz(n_local_part);
466 for (
auto const &p : particles) {
467 auto const dip = p.calc_dip();
476 for (std::size_t i = 0u; i < n_local_part; i++) {
477 local_dip[0u] += mx[i];
478 local_dip[1u] += my[i];
479 local_dip[2u] += mz[i];
482 boost::mpi::all_reduce(
comm_cart, local_dip, std::plus<>());
487 for (std::size_t i = 0u; i < n_local_part; i++) {
488 sum_e += mx[i] * box_dip[0] + my[i] * box_dip[1] + mz[i] * box_dip[2];
491 0.5 * pref * boost::mpi::all_reduce(
comm_cart, sum_e, std::plus<>());
496 std::vector<double> sumix(n_local_part);
497 std::vector<double> sumiy(n_local_part);
498 std::vector<double> sumiz(n_local_part);
500 for (std::size_t i = 0u; i < n_local_part; i++) {
501 sumix[i] = my[i] * box_dip[2u] - mz[i] * box_dip[1u];
502 sumiy[i] = mz[i] * box_dip[0u] - mx[i] * box_dip[2u];
503 sumiz[i] = mx[i] * box_dip[1u] - my[i] * box_dip[0u];
507 for (
auto &p : particles) {
508 auto &torque = p.torque();
509 torque[0u] -= pref * sumix[ip];
510 torque[1u] -= pref * sumiy[ip];
511 torque[2u] -= pref * sumiz[ip];
519template <
typename FloatType, Arch Architecture>
521 dp3m.g_force = grid_influence_function<FloatType, 3>(
522 dp3m.params, dp3m.mesh.start, dp3m.mesh.stop,
523 get_system().box_geo->length_inv());
526template <
typename FloatType, Arch Architecture>
528 dp3m.g_energy = grid_influence_function<FloatType, 2>(
529 dp3m.params, dp3m.mesh.start, dp3m.mesh.stop,
530 get_system().box_geo->length_inv());
533template <
typename FloatType, Arch Architecture>
536 int m_mesh_max = -1, m_mesh_min = -1;
540 double prefactor,
int timings)
547 std::optional<std::string>
554 m_logger = std::make_unique<TuningLogger>(
562 std::tuple<double, double, double, double>
564 double r_cut_iL)
const override {
566 double alpha_L, rs_err, ks_err;
578 0.0001 * box_geo.length()[0], 5. * box_geo.length()[0], 0.0001,
601 auto const expo = std::log(std::cbrt(dp3m.
sum_dip_part)) / std::log(2.);
603 m_mesh_min =
static_cast<int>(std::round(std::pow(2., std::floor(expo))));
607 m_mesh_min = m_mesh_max = dp3m.
params.
mesh[0];
615 for (
auto tmp_mesh = m_mesh_min; tmp_mesh <= m_mesh_max; tmp_mesh += 2) {
620 auto const trial_time =
621 get_m_time(trial_params.mesh, trial_params.cao, trial_params.r_cut_iL,
622 trial_params.alpha_L, trial_params.accuracy);
632 if (trial_time < time_best) {
635 tuned_params = trial_params;
636 time_best = tuned_params.time = trial_time;
647template <
typename FloatType, Arch Architecture>
649 auto &system = get_system();
650 auto const &box_geo = *system.box_geo;
657 if (not is_tuned()) {
660 throw std::runtime_error(
661 "DipolarP3M: no dipolar particles in the system");
665 system, dp3m, prefactor, tune_timings);
674 system.on_dipoles_change();
685 double mesh_i,
int cao,
double alpha_L_i) {
689 auto const factor1 =
Utils::sqr(std::numbers::pi * alpha_L_i);
697 mesh_start, mesh_stop, indices,
699 auto const norm_sq = nm.norm2();
700 auto const ex = std::exp(-factor1 * norm_sq);
703 alias2 += U2 * ex * std::pow(shift * nm, 3) / norm_sq;
705 [&](
unsigned dim,
int n) {
706 nm[dim] = shift[dim] + n * mesh;
710 return std::make_pair(alias1, alias2);
715 int n_c_part,
double sum_q2,
double alpha_L) {
718 auto const mesh_i = 1. /
static_cast<double>(mesh);
719 auto const alpha_L_i = 1. / alpha_L;
721 auto const mesh_start = -mesh_stop;
727 mesh_start, mesh_stop, indices,
729 if ((indices[0] != 0) or (indices[1] != 0) or (indices[2] != 0)) {
730 auto const n2 = indices.norm2();
732 auto const [alias1, alias2] =
736 Utils::int_pow<3>(
static_cast<double>(n2));
743 [&values, &mesh_i, cotangent_sum](
unsigned dim,
int n) {
744 values[dim] = cotangent_sum(n, mesh_i);
747 return 8. *
Utils::sqr(std::numbers::pi) / 3. * sum_q2 *
748 sqrt(he_q / n_c_part) / Utils::int_pow<4>(box_size);
758 int n_c_part,
double sum_q2,
760 double d_error_f, d_cc, d_dc, d_con;
762 auto const d_rcut = r_cut_iL * box_size;
768 auto const d_c = sum_q2 * exp(-d_a2 * d_rcut2);
772 d_dc = 8. * Utils::int_pow<3>(d_a2) * Utils::int_pow<3>(d_rcut2) +
773 20. *
Utils::sqr(d_a2) * d_rcut4 + 30. * d_a2 * d_rcut2 + 15.;
775 d_con = 1. / sqrt(Utils::int_pow<3>(box_size) *
Utils::sqr(d_a2) * d_rcut *
776 Utils::sqr(d_rcut4) *
static_cast<double>(n_c_part));
778 d_error_f = d_c * d_con *
780 (2. / 15.) *
Utils::sqr(d_dc) - (13. / 15.) * d_cc * d_dc);
790 double sum_q2,
double x1,
double x2,
double xacc,
791 double tuned_accuracy) {
792 constexpr int JJ_RTBIS_MAX = 40;
794 auto const constant = tuned_accuracy / std::numbers::sqrt2;
802 if (f1 * f2 >= 0.0) {
803 throw std::runtime_error(
804 "Root must be bracketed for bisection in dp3m_rtbisection");
808 double rtb = f1 < 0.0 ? (dx = x2 - x1, x1) : (dx = x1 - x2, x2);
809 for (
int j = 1; j <= JJ_RTBIS_MAX; j++) {
810 auto const xmid = rtb + (dx *= 0.5);
816 if (fabs(dx) < xacc || fmid == 0.0)
819 throw std::runtime_error(
"Too many bisections in dp3m_rtbisection");
824 auto const &box_geo = *system.box_geo;
825 auto const &local_geo = *system.local_geo;
826 for (
auto i = 0u; i < 3u; i++) {
829 std::stringstream msg;
831 <<
" is larger than half of box dimension " << box_geo.length()[i];
832 throw std::runtime_error(msg.str());
835 std::stringstream msg;
837 <<
" is larger than local box dimension " << local_geo.length()[i];
838 throw std::runtime_error(msg.str());
842 if ((box_geo.length()[0] != box_geo.length()[1]) or
843 (box_geo.length()[1] != box_geo.length()[2])) {
844 throw std::runtime_error(
"DipolarP3M: requires a cubic box");
850 if (!box_geo.periodic(0) or !box_geo.periodic(1) or !box_geo.periodic(2)) {
851 throw std::runtime_error(
852 "DipolarP3M: requires periodicity (True, True, True)");
857 auto const &local_geo = *
get_system().local_geo;
860 throw std::runtime_error(
861 "DipolarP3M: requires the regular or hybrid decomposition cell system");
865 throw std::runtime_error(
866 "DipolarP3M: does not work with the hybrid decomposition cell system, "
867 "if using more than one MPI node");
873 if (node_grid[0] < node_grid[1] or node_grid[1] < node_grid[2]) {
874 throw std::runtime_error(
875 "DipolarP3M: node grid must be sorted, largest first");
879template <
typename FloatType, Arch Architecture>
881 auto const &box_geo = *get_system().
box_geo;
886 sanity_checks_boxl();
887 calc_influence_function_force();
888 calc_influence_function_energy();
892template <
typename FloatType, Arch Architecture>
894 auto const &box_geo = *get_system().
box_geo;
895 auto const Ukp3m = calc_average_self_energy_k_space() * box_geo.volume();
896 auto const Ewald_volume = Utils::int_pow<3>(dp3m.
params.
alpha_L);
897 auto const Eself = -2. * Ewald_volume * std::numbers::inv_sqrtpi / 3.;
899 -dp3m.
sum_mu2 * (Ukp3m + Eself + 2. * std::numbers::pi / 3.);
@ HYBRID
Hybrid decomposition.
@ REGULAR
Regular decomposition.
Vector implementation and trait types for boost qvm interoperability.
void on_solver_change() const override
DipolarTuningAlgorithm(System::System &system, decltype(dp3m) &input_dp3m, double prefactor, int timings)
void setup_logger(bool verbose) override
TuningAlgorithm::Parameters get_time() override
P3MParameters & get_params() override
std::tuple< double, double, double, double > calculate_accuracy(Utils::Vector3i const &mesh, int cao, double r_cut_iL) const override
std::optional< std::string > layer_correction_veto_r_cut(double) const override
void determine_mesh_limits() override
base_type::size_type size() const
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
Tuning algorithm for P3M.
double get_m_time(Utils::Vector3i const &mesh, int &tuned_cao, double &tuned_r_cut_iL, double &tuned_alpha_L, double &tuned_accuracy)
Get the optimal alpha and the corresponding computation time for a fixed mesh.
static auto constexpr time_sentinel
Value for invalid time measurements.
static auto constexpr max_n_consecutive_trials
Maximal number of consecutive trials that don't improve runtime.
System::System & m_system
void determine_cao_limits(int initial_cao)
Determine a sensible range for the charge assignment order.
void determine_r_cut_limits()
Determine a sensible range for the real-space cutoff.
std::unique_ptr< TuningLogger > m_logger
static auto constexpr time_granularity
Granularity of the time measurement (milliseconds).
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
Communicator communicator
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
#define P3M_BRILLOUIN
P3M: Number of Brillouin zones taken into account in the calculation of the optimal influence functio...
__device__ void vector_product(float const *a, float const *b, float *out)
static int count_magnetic_particles(System::System const &system)
static auto dp3m_tune_aliasing_sums(Utils::Vector3i const &shift, int mesh, double mesh_i, int cao, double alpha_L_i)
Tuning dipolar-P3M.
static double dp3m_k_space_error(double box_size, int mesh, int cao, int n_c_part, double sum_q2, double alpha_L)
Calculate the k-space error of dipolar-P3M.
void npt_add_virial_magnetic_contribution(double energy)
Update the NpT virial.
double dp3m_rtbisection(double box_size, double r_cut_iL, int n_c_part, double sum_q2, double x1, double x2, double xacc, double tuned_accuracy)
Compute the value of alpha through a bisection method.
static double dp3m_real_space_error(double box_size, double r_cut_iL, int n_c_part, double sum_q2, double alpha_L)
Calculate the value of the errors for the REAL part of the force in terms of the splitting parameter ...
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 ...
Routines, row decomposition, data structures and communication for the 3D-FFT.
and std::invocable< Projector, unsigned, int > void for_each_3d(detail::IndexVectorConcept auto &&start, detail::IndexVectorConcept auto &&stop, detail::IndexVectorConcept auto &&counters, Kernel &&kernel, Projector &&projector=detail::noop_projector)
Repeat an operation on every element of a 3D grid.
double grid_influence_function_self_energy(P3MParameters const ¶ms, Utils::Vector3i const &n_start, Utils::Vector3i const &n_stop, std::vector< FloatType > const &g)
Calculate self-energy of the influence function.
void p3m_interpolate(P3MLocalMesh const &local_mesh, InterpolationWeights< cao > const &weights, Kernel kernel)
P3M grid interpolation.
T product(Vector< T, N > const &v)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
auto get_analytic_cotangent_sum_kernel(int cao)
void npt_add_virial_contribution(double energy)
Exports for the NpT code.
Common functions for dipolar and charge P3M.
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
Utils::Vector3i node_grid
void calc_energy_correction() override
double long_range_kernel(bool force_flag, bool energy_flag, ParticleRange const &particles)
Compute the k-space part of forces and energies.
double calc_surface_term(bool force_flag, bool energy_flag, ParticleRange const &particles) override
void dipole_assign(ParticleRange const &particles) override
void calc_influence_function_energy() override
void count_magnetic_particles() override
void calc_influence_function_force() override
void scaleby_box_l() override
double calc_average_self_energy_k_space() const override
void sanity_checks_boxl() const
Checks for correctness of the k-space cutoff.
void sanity_checks_cell_structure() const
P3MParameters const & dp3m_params
void sanity_checks_periodicity() const
void sanity_checks_node_grid() const
void recalc_ld_pos(P3MParameters const ¶ms)
Recalculate quantities derived from the mesh and box length: ld_pos (position of the left down mesh).
Structure to hold P3M parameters and some dependent variables.
Utils::Vector3d cao_cut
cutoff for charge assignment.
double alpha
unscaled alpha_L for use with fast inline functions only
double r_cut_iL
cutoff radius for real space electrostatics (>0), rescaled to r_cut_iL = r_cut * box_l_i.
double accuracy
accuracy of the actual parameter set.
double alpha_L
Ewald splitting parameter (0.
double r_cut
unscaled r_cut_iL for use with fast inline functions only
void recalc_a_ai_cao_cut(Utils::Vector3d const &box_l)
Recalculate quantities derived from the mesh and box length: a, ai and cao_cut.
bool tuning
tuning or production?
Utils::Vector3i mesh
number of mesh points per coordinate direction (>0).
void operator()(auto &dp3m, Utils::Vector3d const &real_pos, Utils::Vector3d const &dip) const
void operator()(auto &dp3m, double prefac, int d_rs, ParticleRange const &particles) const
void operator()(auto &dp3m, double prefac, int d_rs, ParticleRange const &particles) const
double energy_correction
cached k-space self-energy correction
int sum_dip_part
number of dipolar particles (only on head node).
double sum_mu2
Sum of square of magnetic dipoles (only on head node).
P3MLocalMesh local_mesh
Local mesh properties.
P3MParameters params
P3M base parameters.