ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
elc.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 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#include "config/config.hpp"
23
24#ifdef P3M
25
27
30
31#include "BoxGeometry.hpp"
32#include "Particle.hpp"
34#include "ParticleRange.hpp"
36#include "communication.hpp"
37#include "errorhandling.hpp"
38#include "system/System.hpp"
39
40#include <utils/Vector.hpp>
41#include <utils/math/sqr.hpp>
42
43#include <boost/mpi/collectives/all_reduce.hpp>
44#include <boost/range/combine.hpp>
45
46#include <algorithm>
47#include <cassert>
48#include <cmath>
49#include <cstddef>
50#include <functional>
51#include <numbers>
52#include <variant>
53#include <vector>
54
55/** \name Product decomposition data organization
56 * For the cell blocks it is assumed that the lower blocks part is in the
57 * lower half. This has to have positive sign, so that has to be first.
58 */
59/**@{*/
60#define POQESP 0
61#define POQECP 1
62#define POQESM 2
63#define POQECM 3
64
65#define PQESSP 0
66#define PQESCP 1
67#define PQECSP 2
68#define PQECCP 3
69#define PQESSM 4
70#define PQESCM 5
71#define PQECSM 6
72#define PQECCM 7
73/**@}*/
74
75/** ELC axes (x and y directions)*/
76enum class PoQ : int { P, Q };
77/** ELC charge sum/assign protocol: real charges, image charges, or both. */
78enum class ChargeProtocol : int { REAL, IMAGE, BOTH };
79
80/** temporary buffers for product decomposition */
81static std::vector<double> partblk;
82/** collected data from the other cells */
83static double gblcblk[8];
84
85/** structure for caching sin and cos values */
86struct SCCache {
87 double s, c;
88};
89
90/** Cached sin/cos values along the x-axis and y-axis */
91/**@{*/
92static std::vector<SCCache> scxcache;
93static std::vector<SCCache> scycache;
94/**@}*/
95
96/**
97 * @brief Calculate cached sin/cos values for one direction.
98 *
99 * @tparam dir Index of the dimension to consider (e.g. 0 for x ...).
100 *
101 * @param particles Particle to calculate values for
102 * @param n_freq Number of frequencies to calculate per particle
103 * @param u Inverse box length
104 * @return Calculated values.
105 */
106template <std::size_t dir>
107static std::vector<SCCache> calc_sc_cache(ParticleRange const &particles,
108 std::size_t n_freq, double u) {
109 auto constexpr c_2pi = 2. * std::numbers::pi;
110 auto const n_part = particles.size();
111 std::vector<SCCache> ret(n_freq * n_part);
112
113 for (std::size_t freq = 1; freq <= n_freq; freq++) {
114 auto const pref = c_2pi * u * static_cast<double>(freq);
115
116 std::size_t o = (freq - 1) * n_part;
117 for (auto const &p : particles) {
118 auto const arg = pref * p.pos()[dir];
119 ret[o++] = {sin(arg), cos(arg)};
120 }
121 }
122
123 return ret;
124}
125
126static std::pair<std::size_t, std::size_t>
127prepare_sc_cache(ParticleRange const &particles, BoxGeometry const &box_geo,
128 double far_cut) {
129 assert(far_cut >= 0.);
130 auto const n_freq_x =
131 static_cast<std::size_t>(std::ceil(far_cut * box_geo.length()[0]) + 1.);
132 auto const n_freq_y =
133 static_cast<std::size_t>(std::ceil(far_cut * box_geo.length()[1]) + 1.);
134 auto const u_x = box_geo.length_inv()[0];
135 auto const u_y = box_geo.length_inv()[1];
136 scxcache = calc_sc_cache<0>(particles, n_freq_x, u_x);
137 scycache = calc_sc_cache<1>(particles, n_freq_y, u_y);
138 return {n_freq_x, n_freq_y};
139}
140
141/*****************************************************************/
142/* data distribution */
143/*****************************************************************/
144
145static void clear_vec(double *pdc, std::size_t size) {
146 for (std::size_t i = 0; i < size; i++)
147 pdc[i] = 0.;
148}
149
150static void copy_vec(double *pdc_d, double const *pdc_s, std::size_t size) {
151 for (std::size_t i = 0; i < size; i++)
152 pdc_d[i] = pdc_s[i];
153}
154
155static void add_vec(double *pdc_d, double const *pdc_s1, double const *pdc_s2,
156 std::size_t size) {
157 for (std::size_t i = 0; i < size; i++)
158 pdc_d[i] = pdc_s1[i] + pdc_s2[i];
159}
160
161static void addscale_vec(double *pdc_d, double scale, double const *pdc_s1,
162 double const *pdc_s2, std::size_t size) {
163 for (std::size_t i = 0; i < size; i++)
164 pdc_d[i] = scale * pdc_s1[i] + pdc_s2[i];
165}
166
167static void scale_vec(double scale, double *pdc, std::size_t size) {
168 for (std::size_t i = 0; i < size; i++)
169 pdc[i] *= scale;
170}
171
172static double *block(double *p, std::size_t index, std::size_t size) {
173 return &p[index * size];
174}
175
176static void distribute(std::size_t size) {
177 assert(size <= 8);
178 double send_buf[8];
179 copy_vec(send_buf, gblcblk, size);
180 boost::mpi::all_reduce(comm_cart, send_buf, static_cast<int>(size), gblcblk,
181 std::plus<>());
182}
183
184void ElectrostaticLayerCorrection::check_gap(Particle const &p) const {
185 if (p.q() != 0.) {
186 auto const z = p.pos()[2];
187 if (z < 0. or z > elc.box_h) {
188 runtimeErrorMsg() << "Particle " << p.id() << " entered ELC gap "
189 << "region by " << ((z < 0.) ? z : z - elc.box_h);
190 }
191 }
192}
193
194/*****************************************************************/
195/* dipole terms */
196/*****************************************************************/
197
198/** Calculate the dipole force.
199 * See @cite yeh99a.
200 */
201void ElectrostaticLayerCorrection::add_dipole_force(
202 ParticleRange const &particles) const {
203 constexpr std::size_t size = 3;
204 auto const &box_geo = *get_system().box_geo;
205 auto const pref = prefactor * 4. * std::numbers::pi / box_geo.volume();
206
207 /* for non-neutral systems, this shift gives the background contribution
208 * (rsp. for this shift, the DM of the background is zero) */
209 auto const shift = box_geo.length_half()[2];
210
211 // collect moments
212
213 gblcblk[0] = 0.; // sum q_i (z_i - L/2)
214 gblcblk[1] = 0.; // sum q_i z_i
215 gblcblk[2] = 0.; // sum q_i
216
217 for (auto const &p : particles) {
218 check_gap(p);
219 auto const q = p.q();
220 auto const z = p.pos()[2];
221
222 gblcblk[0] += q * (z - shift);
223 gblcblk[1] += q * z;
224 gblcblk[2] += q;
225
227 if (z < elc.space_layer) {
228 gblcblk[0] += elc.delta_mid_bot * q * (-z - shift);
229 gblcblk[2] += elc.delta_mid_bot * q;
230 }
231 if (z > (elc.box_h - elc.space_layer)) {
232 gblcblk[0] += elc.delta_mid_top * q * (2. * elc.box_h - z - shift);
233 gblcblk[2] += elc.delta_mid_top * q;
234 }
235 }
236 }
237
238 gblcblk[0] *= pref;
239 gblcblk[1] *= pref / elc.box_h * box_geo.length()[2];
240 gblcblk[2] *= pref;
241
242 distribute(size);
243
244 // Yeh + Berkowitz dipole term @cite yeh99a
245 auto field_tot = gblcblk[0];
246
247 // Constand potential contribution
248 if (elc.const_pot) {
249 auto const field_induced = gblcblk[1];
250 auto const field_applied = elc.pot_diff / elc.box_h;
251 field_tot -= field_applied + field_induced;
252 }
253
254 for (auto &p : particles) {
255 p.force()[2] -= field_tot * p.q();
256
257 if (!elc.neutralize) {
258 // SUBTRACT the forces of the P3M homogeneous neutralizing background
259 p.force()[2] += gblcblk[2] * p.q() * (p.pos()[2] - shift);
260 }
261 }
262}
263
264/** Calculate the dipole energy.
265 * See @cite yeh99a.
266 */
267double ElectrostaticLayerCorrection::dipole_energy(
268 ParticleRange const &particles) const {
269 constexpr std::size_t size = 7;
270 auto const &box_geo = *get_system().box_geo;
271 auto const pref = prefactor * 2. * std::numbers::pi / box_geo.volume();
272 auto const lz = box_geo.length()[2];
273 /* for nonneutral systems, this shift gives the background contribution
274 (rsp. for this shift, the DM of the background is zero) */
275 auto const shift = box_geo.length_half()[2];
276
277 // collect moments
278
279 gblcblk[0] = 0.; // sum q_i primary box
280 gblcblk[1] = 0.; // sum q_i boundary layers
281 gblcblk[2] = 0.; // sum q_i (z_i - L/2) primary box
282 gblcblk[3] = 0.; // sum q_i (z_i - L/2) boundary layers
283 gblcblk[4] = 0.; // sum q_i (z_i - L/2)^2 primary box
284 gblcblk[5] = 0.; // sum q_i (z_i - L/2)^2 boundary layers
285 gblcblk[6] = 0.; // sum q_i z_i primary box
286
287 for (auto const &p : particles) {
288 check_gap(p);
289 auto const q = p.q();
290 auto const z = p.pos()[2];
291
292 gblcblk[0] += q;
293 gblcblk[2] += q * (z - shift);
294 gblcblk[4] += q * (Utils::sqr(z - shift));
295 gblcblk[6] += q * z;
296
298 if (z < elc.space_layer) {
299 gblcblk[1] += elc.delta_mid_bot * q;
300 gblcblk[3] += elc.delta_mid_bot * q * (-z - shift);
301 gblcblk[5] += elc.delta_mid_bot * q * (Utils::sqr(-z - shift));
302 }
303 if (z > (elc.box_h - elc.space_layer)) {
304 gblcblk[1] += elc.delta_mid_top * q;
305 gblcblk[3] += elc.delta_mid_top * q * (2. * elc.box_h - z - shift);
306 gblcblk[5] +=
307 elc.delta_mid_top * q * (Utils::sqr(2. * elc.box_h - z - shift));
308 }
309 }
310 }
311
312 distribute(size);
313
314 // Yeh + Berkowitz term @cite yeh99a
315 auto energy = 2. * pref * (Utils::sqr(gblcblk[2]) + gblcblk[2] * gblcblk[3]);
316
317 if (!elc.neutralize) {
318 // SUBTRACT the energy of the P3M homogeneous neutralizing background
319 energy += 2. * pref *
320 (-gblcblk[0] * gblcblk[4] -
321 (.25 - .5 / 3.) * Utils::sqr(gblcblk[0] * lz));
322 }
323
325 if (elc.const_pot) {
326 // zero potential difference contribution
327 energy += pref / elc.box_h * lz * Utils::sqr(gblcblk[6]);
328 // external potential shift contribution
329 energy -= 2. * elc.pot_diff / elc.box_h * gblcblk[6];
330 }
331
332 /* counter the P3M homogeneous background contribution to the
333 boundaries. We never need that, since a homogeneous background
334 spanning the artificial boundary layers is aphysical. */
335 energy +=
336 pref * (-(gblcblk[1] * gblcblk[4] + gblcblk[0] * gblcblk[5]) -
337 (1. - 2. / 3.) * gblcblk[0] * gblcblk[1] * Utils::sqr(lz));
338 }
339
340 return this_node == 0 ? energy : 0.;
341}
342
343/*****************************************************************/
344
345struct ImageSum {
346 double delta;
347 double shift;
348 double lz;
349 double dci; // delta complement inverse
350
351 ImageSum(double delta, double shift, double lz)
352 : delta{delta}, shift{shift}, lz{lz}, dci{1. / (1. - delta)} {}
353
354 /** @brief Image sum from the bottom layer. */
355 double b(double q, double z) const {
356 return q * dci * (z - 2. * delta * lz * dci) - q * dci * shift;
357 }
358
359 /** @brief Image sum from the top layer. */
360 double t(double q, double z) const {
361 return q * dci * (z + 2. * delta * lz * dci) - q * dci * shift;
362 }
363};
364
365double
366ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const {
367 constexpr std::size_t size = 4;
368 auto const &box_geo = *get_system().box_geo;
369 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
370 auto const pref = prefactor * 2. * std::numbers::pi * xy_area_inv;
371 auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
372 auto const fac_delta_mid_bot = elc.delta_mid_bot / (1. - delta);
373 auto const fac_delta_mid_top = elc.delta_mid_top / (1. - delta);
374 auto const fac_delta = delta / (1. - delta);
375
376 /* for non-neutral systems, this shift gives the background contribution
377 * (rsp. for this shift, the DM of the background is zero) */
378 auto const shift = box_geo.length_half()[2];
379 auto const lz = box_geo.length()[2];
380
382 if (elc.const_pot) {
383 clear_vec(gblcblk, size);
384 for (auto const &p : particles) {
385 auto const z = p.pos()[2];
386 auto const q = p.q();
387 gblcblk[0] += q;
388 gblcblk[1] += q * (z - shift);
389 if (z < elc.space_layer) {
390 gblcblk[2] -= elc.delta_mid_bot * q;
391 gblcblk[3] -= elc.delta_mid_bot * q * (-z - shift);
392 }
393 if (z > (elc.box_h - elc.space_layer)) {
394 gblcblk[2] += elc.delta_mid_top * q;
395 gblcblk[3] += elc.delta_mid_top * q * (2. * elc.box_h - z - shift);
396 }
397 }
398 } else {
399 // metallic boundaries
400 clear_vec(gblcblk, size);
401 auto const h = elc.box_h;
402 ImageSum const image_sum{delta, shift, lz};
403 for (auto const &p : particles) {
404 auto const z = p.pos()[2];
405 auto const q = p.q();
406 gblcblk[0] += q;
407 gblcblk[1] += q * (z - shift);
409 if (z < elc.space_layer) {
410 gblcblk[2] += fac_delta * (elc.delta_mid_bot + 1.) * q;
411 gblcblk[3] +=
412 q * (image_sum.b(elc.delta_mid_bot * delta, -(2. * h + z)) +
413 image_sum.b(delta, -(2. * h - z)));
414 } else {
415 gblcblk[2] += fac_delta_mid_bot * (1. + elc.delta_mid_top) * q;
416 gblcblk[3] += q * (image_sum.b(elc.delta_mid_bot, -z) +
417 image_sum.b(delta, -(2. * h - z)));
418 }
419 if (z > (h - elc.space_layer)) {
420 // note the minus sign here which is required due to |z_i-z_j|
421 gblcblk[2] -= fac_delta * (elc.delta_mid_top + 1.) * q;
422 gblcblk[3] -=
423 q * (image_sum.t(elc.delta_mid_top * delta, 4. * h - z) +
424 image_sum.t(delta, 2. * h + z));
425 } else {
426 // note the minus sign here which is required due to |z_i-z_j|
427 gblcblk[2] -= fac_delta_mid_top * (1. + elc.delta_mid_bot) * q;
428 gblcblk[3] -= q * (image_sum.t(elc.delta_mid_top, 2. * h - z) +
429 image_sum.t(delta, 2. * h + z));
430 }
431 }
432 }
433 }
434 }
435 distribute(size);
436
437 auto const energy = gblcblk[1] * gblcblk[2] - gblcblk[0] * gblcblk[3];
438 return (this_node == 0) ? -pref * energy : 0.;
439}
440
441void ElectrostaticLayerCorrection::add_z_force(
442 ParticleRange const &particles) const {
443 constexpr std::size_t size = 1;
444 auto const &box_geo = *get_system().box_geo;
445 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
446 auto const pref = prefactor * 2. * std::numbers::pi * xy_area_inv;
447 auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
448 auto const fac_delta_mid_bot = elc.delta_mid_bot / (1. - delta);
449 auto const fac_delta_mid_top = elc.delta_mid_top / (1. - delta);
450 auto const fac_delta = delta / (1. - delta);
451
453 if (elc.const_pot) {
454 clear_vec(gblcblk, size);
455 /* just counter the 2 pi |z| contribution stemming from P3M */
456 for (auto const &p : particles) {
457 auto const z = p.pos()[2];
458 auto const q = p.q();
459 if (z < elc.space_layer)
460 gblcblk[0] -= elc.delta_mid_bot * q;
461 if (z > (elc.box_h - elc.space_layer))
462 gblcblk[0] += elc.delta_mid_top * q;
463 }
464 } else {
465 clear_vec(gblcblk, size);
466 for (auto const &p : particles) {
467 auto const z = p.pos()[2];
468 auto const q = p.q();
469 if (z < elc.space_layer) {
470 gblcblk[0] += fac_delta * (elc.delta_mid_bot + 1.) * q;
471 } else {
472 gblcblk[0] += fac_delta_mid_bot * (elc.delta_mid_top + 1.) * q;
473 }
474 if (z > (elc.box_h - elc.space_layer)) {
475 // note the minus sign here which is required due to |z_i-z_j|
476 gblcblk[0] -= fac_delta * (elc.delta_mid_top + 1.) * q;
477 } else {
478 // note the minus sign here which is required due to |z_i-z_j|
479 gblcblk[0] -= fac_delta_mid_top * (elc.delta_mid_bot + 1.) * q;
480 }
481 }
482 }
483
484 gblcblk[0] *= pref;
485
486 distribute(size);
487
488 for (auto &p : particles) {
489 p.force()[2] += gblcblk[0] * p.q();
490 }
491 }
492}
493
494/*****************************************************************/
495/* PoQ exp sum */
496/*****************************************************************/
497
498/** \name q=0 or p=0 per frequency code */
499/**@{*/
500template <PoQ axis>
501void setup_PoQ(elc_data const &elc, double prefactor, std::size_t index,
502 double omega, ParticleRange const &particles,
503 BoxGeometry const &box_geo) {
504 assert(index >= 1);
505 constexpr std::size_t size = 4;
506 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
507 auto const pref_di = prefactor * 4. * std::numbers::pi * xy_area_inv;
508 auto const pref = -pref_di / expm1(omega * box_geo.length()[2]);
509 double lclimgebot[4], lclimgetop[4], lclimge[4];
510 double fac_delta_mid_bot = 1., fac_delta_mid_top = 1., fac_delta = 1.;
511
512 if (elc.dielectric_contrast_on) {
513 auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
514 auto const fac_elc = 1. / (1. - delta * exp(-omega * 2. * elc.box_h));
515 fac_delta_mid_bot = elc.delta_mid_bot * fac_elc;
516 fac_delta_mid_top = elc.delta_mid_top * fac_elc;
517 fac_delta = fac_delta_mid_bot * elc.delta_mid_top;
518 }
519
520 clear_vec(lclimge, size);
521 clear_vec(gblcblk, size);
522 auto const &sc_cache = (axis == PoQ::P) ? scxcache : scycache;
523
524 std::size_t ic = 0;
525 auto const o = (index - 1) * particles.size();
526 for (auto const &p : particles) {
527 auto const z = p.pos()[2];
528 auto const q = p.q();
529 auto e = exp(omega * z);
530
531 partblk[size * ic + POQESM] = q * sc_cache[o + ic].s / e;
532 partblk[size * ic + POQESP] = q * sc_cache[o + ic].s * e;
533 partblk[size * ic + POQECM] = q * sc_cache[o + ic].c / e;
534 partblk[size * ic + POQECP] = q * sc_cache[o + ic].c * e;
535
536 add_vec(gblcblk, gblcblk, block(partblk.data(), ic, size), size);
537
538 if (elc.dielectric_contrast_on) {
539 if (z < elc.space_layer) { // handle the lower case first
540 // negative sign is okay here as the image is located at -z
541
542 e = exp(-omega * z);
543
544 auto const scale = q * elc.delta_mid_bot;
545
546 lclimgebot[POQESM] = sc_cache[o + ic].s / e;
547 lclimgebot[POQESP] = sc_cache[o + ic].s * e;
548 lclimgebot[POQECM] = sc_cache[o + ic].c / e;
549 lclimgebot[POQECP] = sc_cache[o + ic].c * e;
550
551 addscale_vec(gblcblk, scale, lclimgebot, gblcblk, size);
552
553 e = (exp(omega * (-z - 2. * elc.box_h)) * elc.delta_mid_bot +
554 exp(omega * (+z - 2. * elc.box_h))) *
555 fac_delta;
556 } else {
557 e = (exp(-omega * z) +
558 exp(omega * (z - 2. * elc.box_h)) * elc.delta_mid_top) *
559 fac_delta_mid_bot;
560 }
561
562 lclimge[POQESP] += q * sc_cache[o + ic].s * e;
563 lclimge[POQECP] += q * sc_cache[o + ic].c * e;
564
565 if (z > (elc.box_h - elc.space_layer)) { // handle the upper case now
566 e = exp(omega * (2. * elc.box_h - z));
567
568 auto const scale = q * elc.delta_mid_top;
569
570 lclimgetop[POQESM] = sc_cache[o + ic].s / e;
571 lclimgetop[POQESP] = sc_cache[o + ic].s * e;
572 lclimgetop[POQECM] = sc_cache[o + ic].c / e;
573 lclimgetop[POQECP] = sc_cache[o + ic].c * e;
574
575 addscale_vec(gblcblk, scale, lclimgetop, gblcblk, size);
576
577 e = (exp(omega * (+z - 4. * elc.box_h)) * elc.delta_mid_top +
578 exp(omega * (-z - 2. * elc.box_h))) *
579 fac_delta;
580 } else {
581 e = (exp(omega * (+z - 2. * elc.box_h)) +
582 exp(omega * (-z - 2. * elc.box_h)) * elc.delta_mid_bot) *
583 fac_delta_mid_top;
584 }
585
586 lclimge[POQESM] += q * sc_cache[o + ic].s * e;
587 lclimge[POQECM] += q * sc_cache[o + ic].c * e;
588 }
589
590 ++ic;
591 }
592
593 scale_vec(pref, gblcblk, size);
594
595 if (elc.dielectric_contrast_on) {
596 scale_vec(pref_di, lclimge, size);
597 add_vec(gblcblk, gblcblk, lclimge, size);
598 }
599}
600
601template <PoQ axis> void add_PoQ_force(ParticleRange const &particles) {
602 constexpr auto i = static_cast<int>(axis);
603 constexpr std::size_t size = 4;
604
605 std::size_t ic = 0;
606 for (auto &p : particles) {
607 auto &force = p.force();
608 force[i] += partblk[size * ic + POQESM] * gblcblk[POQECP] -
609 partblk[size * ic + POQECM] * gblcblk[POQESP] +
610 partblk[size * ic + POQESP] * gblcblk[POQECM] -
611 partblk[size * ic + POQECP] * gblcblk[POQESM];
612 force[2] += partblk[size * ic + POQECM] * gblcblk[POQECP] +
613 partblk[size * ic + POQESM] * gblcblk[POQESP] -
614 partblk[size * ic + POQECP] * gblcblk[POQECM] -
615 partblk[size * ic + POQESP] * gblcblk[POQESM];
616 ++ic;
617 }
618}
619
620static double PoQ_energy(double omega, std::size_t n_part) {
621 constexpr std::size_t size = 4;
622
623 auto energy = 0.;
624 for (std::size_t ic = 0; ic < n_part; ic++) {
625 energy += partblk[size * ic + POQECM] * gblcblk[POQECP] +
626 partblk[size * ic + POQESM] * gblcblk[POQESP] +
627 partblk[size * ic + POQECP] * gblcblk[POQECM] +
628 partblk[size * ic + POQESP] * gblcblk[POQESM];
629 }
630
631 return energy / omega;
632}
633/**@}*/
634
635/*****************************************************************/
636/* PQ particle blocks */
637/*****************************************************************/
638
639/** \name p,q <> 0 per frequency code */
640/**@{*/
641static void setup_PQ(elc_data const &elc, double prefactor, std::size_t index_p,
642 std::size_t index_q, double omega,
643 ParticleRange const &particles,
644 BoxGeometry const &box_geo) {
645 assert(index_p >= 1);
646 assert(index_q >= 1);
647 constexpr std::size_t size = 8;
648 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
649 auto const pref_di = prefactor * 8. * std::numbers::pi * xy_area_inv;
650 auto const pref = -pref_di / expm1(omega * box_geo.length()[2]);
651 double lclimgebot[8], lclimgetop[8], lclimge[8];
652 double fac_delta_mid_bot = 1., fac_delta_mid_top = 1., fac_delta = 1.;
653 if (elc.dielectric_contrast_on) {
654 auto const delta = elc.delta_mid_top * elc.delta_mid_bot;
655 auto const fac_elc = 1. / (1. - delta * exp(-omega * 2. * elc.box_h));
656 fac_delta_mid_bot = elc.delta_mid_bot * fac_elc;
657 fac_delta_mid_top = elc.delta_mid_top * fac_elc;
658 fac_delta = fac_delta_mid_bot * elc.delta_mid_top;
659 }
660
661 clear_vec(lclimge, size);
662 clear_vec(gblcblk, size);
663
664 std::size_t ic = 0;
665 auto const ox = (index_p - 1) * particles.size();
666 auto const oy = (index_q - 1) * particles.size();
667 for (auto const &p : particles) {
668 auto const z = p.pos()[2];
669 auto const q = p.q();
670 auto e = exp(omega * z);
671
672 partblk[size * ic + PQESSM] =
673 scxcache[ox + ic].s * scycache[oy + ic].s * q / e;
674 partblk[size * ic + PQESCM] =
675 scxcache[ox + ic].s * scycache[oy + ic].c * q / e;
676 partblk[size * ic + PQECSM] =
677 scxcache[ox + ic].c * scycache[oy + ic].s * q / e;
678 partblk[size * ic + PQECCM] =
679 scxcache[ox + ic].c * scycache[oy + ic].c * q / e;
680
681 partblk[size * ic + PQESSP] =
682 scxcache[ox + ic].s * scycache[oy + ic].s * q * e;
683 partblk[size * ic + PQESCP] =
684 scxcache[ox + ic].s * scycache[oy + ic].c * q * e;
685 partblk[size * ic + PQECSP] =
686 scxcache[ox + ic].c * scycache[oy + ic].s * q * e;
687 partblk[size * ic + PQECCP] =
688 scxcache[ox + ic].c * scycache[oy + ic].c * q * e;
689
690 add_vec(gblcblk, gblcblk, block(partblk.data(), ic, size), size);
691
692 if (elc.dielectric_contrast_on) {
693 if (z < elc.space_layer) { // handle the lower case first
694 // change e to take into account the z position of the images
695
696 e = exp(-omega * z);
697 auto const scale = q * elc.delta_mid_bot;
698
699 lclimgebot[PQESSM] = scxcache[ox + ic].s * scycache[oy + ic].s / e;
700 lclimgebot[PQESCM] = scxcache[ox + ic].s * scycache[oy + ic].c / e;
701 lclimgebot[PQECSM] = scxcache[ox + ic].c * scycache[oy + ic].s / e;
702 lclimgebot[PQECCM] = scxcache[ox + ic].c * scycache[oy + ic].c / e;
703
704 lclimgebot[PQESSP] = scxcache[ox + ic].s * scycache[oy + ic].s * e;
705 lclimgebot[PQESCP] = scxcache[ox + ic].s * scycache[oy + ic].c * e;
706 lclimgebot[PQECSP] = scxcache[ox + ic].c * scycache[oy + ic].s * e;
707 lclimgebot[PQECCP] = scxcache[ox + ic].c * scycache[oy + ic].c * e;
708
709 addscale_vec(gblcblk, scale, lclimgebot, gblcblk, size);
710
711 e = (exp(omega * (-z - 2. * elc.box_h)) * elc.delta_mid_bot +
712 exp(omega * (+z - 2. * elc.box_h))) *
713 fac_delta * q;
714
715 } else {
716
717 e = (exp(-omega * z) +
718 exp(omega * (z - 2. * elc.box_h)) * elc.delta_mid_top) *
719 fac_delta_mid_bot * q;
720 }
721
722 lclimge[PQESSP] += scxcache[ox + ic].s * scycache[oy + ic].s * e;
723 lclimge[PQESCP] += scxcache[ox + ic].s * scycache[oy + ic].c * e;
724 lclimge[PQECSP] += scxcache[ox + ic].c * scycache[oy + ic].s * e;
725 lclimge[PQECCP] += scxcache[ox + ic].c * scycache[oy + ic].c * e;
726
727 if (z > (elc.box_h - elc.space_layer)) { // handle the upper case now
728
729 e = exp(omega * (2. * elc.box_h - z));
730 auto const scale = q * elc.delta_mid_top;
731
732 lclimgetop[PQESSM] = scxcache[ox + ic].s * scycache[oy + ic].s / e;
733 lclimgetop[PQESCM] = scxcache[ox + ic].s * scycache[oy + ic].c / e;
734 lclimgetop[PQECSM] = scxcache[ox + ic].c * scycache[oy + ic].s / e;
735 lclimgetop[PQECCM] = scxcache[ox + ic].c * scycache[oy + ic].c / e;
736
737 lclimgetop[PQESSP] = scxcache[ox + ic].s * scycache[oy + ic].s * e;
738 lclimgetop[PQESCP] = scxcache[ox + ic].s * scycache[oy + ic].c * e;
739 lclimgetop[PQECSP] = scxcache[ox + ic].c * scycache[oy + ic].s * e;
740 lclimgetop[PQECCP] = scxcache[ox + ic].c * scycache[oy + ic].c * e;
741
742 addscale_vec(gblcblk, scale, lclimgetop, gblcblk, size);
743
744 e = (exp(omega * (+z - 4. * elc.box_h)) * elc.delta_mid_top +
745 exp(omega * (-z - 2. * elc.box_h))) *
746 fac_delta * q;
747
748 } else {
749
750 e = (exp(omega * (+z - 2. * elc.box_h)) +
751 exp(omega * (-z - 2. * elc.box_h)) * elc.delta_mid_bot) *
752 fac_delta_mid_top * q;
753 }
754
755 lclimge[PQESSM] += scxcache[ox + ic].s * scycache[oy + ic].s * e;
756 lclimge[PQESCM] += scxcache[ox + ic].s * scycache[oy + ic].c * e;
757 lclimge[PQECSM] += scxcache[ox + ic].c * scycache[oy + ic].s * e;
758 lclimge[PQECCM] += scxcache[ox + ic].c * scycache[oy + ic].c * e;
759 }
760
761 ic++;
762 }
763
764 scale_vec(pref, gblcblk, size);
765 if (elc.dielectric_contrast_on) {
766 scale_vec(pref_di, lclimge, size);
767 add_vec(gblcblk, gblcblk, lclimge, size);
768 }
769}
770
771static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega,
772 ParticleRange const &particles,
773 BoxGeometry const &box_geo) {
774 auto constexpr c_2pi = 2. * std::numbers::pi;
775 auto const pref_x =
776 c_2pi * box_geo.length_inv()[0] * static_cast<double>(index_p) / omega;
777 auto const pref_y =
778 c_2pi * box_geo.length_inv()[1] * static_cast<double>(index_q) / omega;
779 constexpr std::size_t size = 8;
780
781 std::size_t ic = 0;
782 for (auto &p : particles) {
783 auto &force = p.force();
784 force[0] += pref_x * (partblk[size * ic + PQESCM] * gblcblk[PQECCP] +
785 partblk[size * ic + PQESSM] * gblcblk[PQECSP] -
786 partblk[size * ic + PQECCM] * gblcblk[PQESCP] -
787 partblk[size * ic + PQECSM] * gblcblk[PQESSP] +
788 partblk[size * ic + PQESCP] * gblcblk[PQECCM] +
789 partblk[size * ic + PQESSP] * gblcblk[PQECSM] -
790 partblk[size * ic + PQECCP] * gblcblk[PQESCM] -
791 partblk[size * ic + PQECSP] * gblcblk[PQESSM]);
792 force[1] += pref_y * (partblk[size * ic + PQECSM] * gblcblk[PQECCP] +
793 partblk[size * ic + PQESSM] * gblcblk[PQESCP] -
794 partblk[size * ic + PQECCM] * gblcblk[PQECSP] -
795 partblk[size * ic + PQESCM] * gblcblk[PQESSP] +
796 partblk[size * ic + PQECSP] * gblcblk[PQECCM] +
797 partblk[size * ic + PQESSP] * gblcblk[PQESCM] -
798 partblk[size * ic + PQECCP] * gblcblk[PQECSM] -
799 partblk[size * ic + PQESCP] * gblcblk[PQESSM]);
800 force[2] += (partblk[size * ic + PQECCM] * gblcblk[PQECCP] +
801 partblk[size * ic + PQECSM] * gblcblk[PQECSP] +
802 partblk[size * ic + PQESCM] * gblcblk[PQESCP] +
803 partblk[size * ic + PQESSM] * gblcblk[PQESSP] -
804 partblk[size * ic + PQECCP] * gblcblk[PQECCM] -
805 partblk[size * ic + PQECSP] * gblcblk[PQECSM] -
806 partblk[size * ic + PQESCP] * gblcblk[PQESCM] -
807 partblk[size * ic + PQESSP] * gblcblk[PQESSM]);
808 ic++;
809 }
810}
811
812static double PQ_energy(double omega, std::size_t n_part) {
813 constexpr std::size_t size = 8;
814
815 auto energy = 0.;
816 for (std::size_t ic = 0; ic < n_part; ic++) {
817 energy += partblk[size * ic + PQECCM] * gblcblk[PQECCP] +
818 partblk[size * ic + PQECSM] * gblcblk[PQECSP] +
819 partblk[size * ic + PQESCM] * gblcblk[PQESCP] +
820 partblk[size * ic + PQESSM] * gblcblk[PQESSP] +
821 partblk[size * ic + PQECCP] * gblcblk[PQECCM] +
822 partblk[size * ic + PQECSP] * gblcblk[PQECSM] +
823 partblk[size * ic + PQESCP] * gblcblk[PQESCM] +
824 partblk[size * ic + PQESSP] * gblcblk[PQESSM];
825 }
826 return energy / omega;
827}
828/**@}*/
829
830void ElectrostaticLayerCorrection::add_force(
831 ParticleRange const &particles) const {
832 auto constexpr c_2pi = 2. * std::numbers::pi;
833 auto const &box_geo = *get_system().box_geo;
834 auto const n_freqs = prepare_sc_cache(particles, box_geo, elc.far_cut);
835 auto const n_scxcache = std::get<0>(n_freqs);
836 auto const n_scycache = std::get<1>(n_freqs);
837 partblk.resize(particles.size() * 8);
838
839 add_dipole_force(particles);
840 add_z_force(particles);
841
842 /* the second condition is just for the case of numerical accident */
843 for (std::size_t p = 1;
844 box_geo.length_inv()[0] * static_cast<double>(p - 1) < elc.far_cut &&
845 p <= n_scxcache;
846 p++) {
847 auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast<double>(p);
848 setup_PoQ<PoQ::P>(elc, prefactor, p, omega, particles, box_geo);
849 distribute(4);
850 add_PoQ_force<PoQ::P>(particles);
851 }
852
853 for (std::size_t q = 1;
854 box_geo.length_inv()[1] * static_cast<double>(q - 1) < elc.far_cut &&
855 q <= n_scycache;
856 q++) {
857 auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast<double>(q);
858 setup_PoQ<PoQ::Q>(elc, prefactor, q, omega, particles, box_geo);
859 distribute(4);
860 add_PoQ_force<PoQ::Q>(particles);
861 }
862
863 for (std::size_t p = 1;
864 box_geo.length_inv()[0] * static_cast<double>(p - 1) < elc.far_cut &&
865 p <= n_scxcache;
866 p++) {
867 for (std::size_t q = 1;
868 Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p - 1)) +
869 Utils::sqr(box_geo.length_inv()[1] *
870 static_cast<double>(q - 1)) <
871 elc.far_cut2 &&
872 q <= n_scycache;
873 q++) {
874 auto const omega =
875 c_2pi *
876 sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p)) +
877 Utils::sqr(box_geo.length_inv()[1] * static_cast<double>(q)));
878 setup_PQ(elc, prefactor, p, q, omega, particles, box_geo);
879 distribute(8);
880 add_PQ_force(p, q, omega, particles, box_geo);
881 }
882 }
883}
884
885double ElectrostaticLayerCorrection::calc_energy(
886 ParticleRange const &particles) const {
887 auto constexpr c_2pi = 2. * std::numbers::pi;
888 auto const &box_geo = *get_system().box_geo;
889 auto energy = dipole_energy(particles) + z_energy(particles);
890 auto const n_freqs = prepare_sc_cache(particles, box_geo, elc.far_cut);
891 auto const n_scxcache = std::get<0>(n_freqs);
892 auto const n_scycache = std::get<1>(n_freqs);
893
894 auto const n_localpart = particles.size();
895 partblk.resize(n_localpart * 8);
896
897 /* the second condition is just for the case of numerical accident */
898 for (std::size_t p = 1;
899 box_geo.length_inv()[0] * static_cast<double>(p - 1) < elc.far_cut &&
900 p <= n_scxcache;
901 p++) {
902 auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast<double>(p);
903 setup_PoQ<PoQ::P>(elc, prefactor, p, omega, particles, box_geo);
904 distribute(4);
905 energy += PoQ_energy(omega, n_localpart);
906 }
907
908 for (std::size_t q = 1;
909 box_geo.length_inv()[1] * static_cast<double>(q - 1) < elc.far_cut &&
910 q <= n_scycache;
911 q++) {
912 auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast<double>(q);
913 setup_PoQ<PoQ::Q>(elc, prefactor, q, omega, particles, box_geo);
914 distribute(4);
915 energy += PoQ_energy(omega, n_localpart);
916 }
917
918 for (std::size_t p = 1;
919 box_geo.length_inv()[0] * static_cast<double>(p - 1) < elc.far_cut &&
920 p <= n_scxcache;
921 p++) {
922 for (std::size_t q = 1;
923 Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p - 1)) +
924 Utils::sqr(box_geo.length_inv()[1] *
925 static_cast<double>(q - 1)) <
926 elc.far_cut2 &&
927 q <= n_scycache;
928 q++) {
929 auto const omega =
930 c_2pi *
931 sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p)) +
932 Utils::sqr(box_geo.length_inv()[1] * static_cast<double>(q)));
933 setup_PQ(elc, prefactor, p, q, omega, particles, box_geo);
934 distribute(8);
935 energy += PQ_energy(omega, n_localpart);
936 }
937 }
938 /* we count both i<->j and j<->i, so return just half of it */
939 return 0.5 * energy;
940}
941
942double ElectrostaticLayerCorrection::tune_far_cut() const {
943 // Largest reasonable cutoff for far formula
944 auto constexpr maximal_far_cut = 50.;
945 auto const &box_geo = *get_system().box_geo;
946 auto const box_l_x_inv = box_geo.length_inv()[0];
947 auto const box_l_y_inv = box_geo.length_inv()[1];
948 auto const min_inv_boxl = std::min(box_l_x_inv, box_l_y_inv);
949 auto const box_l_z = box_geo.length()[2];
950 // adjust lz according to dielectric layer method
951 auto const lz =
953
954 auto tuned_far_cut = min_inv_boxl;
955 double err;
956 do {
957 auto const pref = 2. * std::numbers::pi * tuned_far_cut;
958 auto const sum = pref + 2. * (box_l_x_inv + box_l_y_inv);
959 auto const den = -expm1(-pref * lz);
960 auto const num1 = exp(pref * (elc.box_h - lz));
961 auto const num2 = exp(-pref * (elc.box_h + lz));
962
963 err = 0.5 / den *
964 (num1 * (sum + 1. / (lz - elc.box_h)) / (lz - elc.box_h) +
965 num2 * (sum + 1. / (lz + elc.box_h)) / (lz + elc.box_h));
966
967 tuned_far_cut += min_inv_boxl;
968 } while (err > elc.maxPWerror and tuned_far_cut < maximal_far_cut);
969 if (tuned_far_cut >= maximal_far_cut) {
970 throw std::runtime_error("ELC tuning failed: maxPWerror too small");
971 }
972 return tuned_far_cut - min_inv_boxl;
973}
974
975static auto calc_total_charge(CellStructure const &cell_structure) {
976 auto local_q = 0.;
977 for (auto const &p : cell_structure.local_particles()) {
978 local_q += p.q();
979 }
980 return boost::mpi::all_reduce(comm_cart, local_q, std::plus<>());
981}
982
983void ElectrostaticLayerCorrection::sanity_checks_periodicity() const {
984 auto const &box_geo = *get_system().box_geo;
985 if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) {
986 throw std::runtime_error("ELC: requires periodicity (True, True, True)");
987 }
988}
989
990void ElectrostaticLayerCorrection::sanity_checks_dielectric_contrasts() const {
992 auto const &cell_structure = *get_system().cell_structure;
993 auto const precision_threshold = std::sqrt(ROUND_ERROR_PREC);
994 auto const total_charge = std::abs(calc_total_charge(cell_structure));
995 if (total_charge >= precision_threshold) {
996 if (elc.const_pot) {
997 // Disable this line to make ELC work again with non-neutral systems
998 // and metallic boundaries
999 throw std::runtime_error("ELC does not currently support non-neutral "
1000 "systems with a dielectric contrast.");
1001 }
1002 // ELC with non-neutral systems and no fully metallic boundaries
1003 // does not work
1004 throw std::runtime_error("ELC does not work for non-neutral systems and "
1005 "non-metallic dielectric contrast.");
1006 }
1007 }
1008}
1009
1010void ElectrostaticLayerCorrection::adapt_solver() {
1011 std::visit(
1012 [this](auto &solver) {
1013 set_prefactor(solver->prefactor);
1014 solver->adapt_epsilon_elc();
1015 assert(solver->p3m_params.epsilon == P3M_EPSILON_METALLIC);
1016 },
1017 base_solver);
1018}
1019
1020void ElectrostaticLayerCorrection::recalc_box_h() {
1021 m_box_geo = get_system().box_geo.get();
1022 auto const box_z = m_box_geo->length()[2];
1023 auto const new_box_h = box_z - elc.gap_size;
1024 if (new_box_h < 0.) {
1025 throw std::runtime_error("ELC gap size (" + std::to_string(elc.gap_size) +
1026 ") larger than box length in z-direction (" +
1027 std::to_string(box_z) + ")");
1028 }
1029 elc.box_h = new_box_h;
1030}
1031
1032void ElectrostaticLayerCorrection::recalc_space_layer() {
1034 auto const p3m_r_cut = std::visit(
1035 [](auto &solver) { return solver->p3m_params.r_cut; }, base_solver);
1036 // recalculate the space layer size:
1037 // 1. set the space_layer to be 1/3 of the gap size, so that box = layer
1038 elc.space_layer = (1. / 3.) * elc.gap_size;
1039 // 2. but make sure we don't overlap with the near-field formula
1040 auto const free_space = elc.gap_size - p3m_r_cut;
1041 // 3. and make sure the space layer is not bigger than half the actual
1042 // simulation box, to avoid overlaps
1043 auto const half_box_h = elc.box_h / 2.;
1044 auto const max_space_layer = std::min(free_space, half_box_h);
1045 if (elc.space_layer > max_space_layer) {
1046 if (max_space_layer <= 0.) {
1047 throw std::runtime_error("P3M real-space cutoff too large for ELC w/ "
1048 "dielectric contrast");
1049 }
1050 elc.space_layer = max_space_layer;
1051 }
1053 }
1054}
1055
1056elc_data::elc_data(double maxPWerror, double gap_size, double far_cut,
1057 bool neutralize, double delta_top, double delta_bot,
1058 bool with_const_pot, double potential_diff)
1059 : maxPWerror{maxPWerror}, gap_size{gap_size}, box_h{-1.}, far_cut{far_cut},
1060 far_cut2{-1.}, far_calculated{far_cut == -1.},
1061 dielectric_contrast_on{delta_top != 0. or delta_bot != 0.},
1062 const_pot{with_const_pot and dielectric_contrast_on},
1063 neutralize{neutralize and !dielectric_contrast_on},
1064 delta_mid_top{std::clamp(delta_top, -1., +1.)},
1065 delta_mid_bot{std::clamp(delta_bot, -1., +1.)},
1066 pot_diff{(with_const_pot) ? potential_diff : 0.},
1067 // initial setup of parameters, may change later when P3M is finally tuned
1068 // set the space_layer to be 1/3 of the gap size, so that box = layer
1069 space_layer{(dielectric_contrast_on) ? gap_size / 3. : 0.},
1070 space_box{gap_size - ((dielectric_contrast_on) ? 2. * space_layer : 0.)} {
1071
1072 auto const delta_range = 1. + std::sqrt(ROUND_ERROR_PREC);
1073 if (far_cut <= 0. and not far_calculated) {
1074 throw std::domain_error("Parameter 'far_cut' must be > 0");
1075 }
1076 if (maxPWerror <= 0.) {
1077 throw std::domain_error("Parameter 'maxPWerror' must be > 0");
1078 }
1079 if (gap_size <= 0.) {
1080 throw std::domain_error("Parameter 'gap_size' must be > 0");
1081 }
1082 if (potential_diff != 0. and not with_const_pot) {
1083 throw std::invalid_argument(
1084 "Parameter 'const_pot' must be True when 'pot_diff' is non-zero");
1085 }
1086 if (delta_top < -delta_range or delta_top > delta_range) {
1087 throw std::domain_error(
1088 "Parameter 'delta_mid_top' must be >= -1 and <= +1");
1089 }
1090 if (delta_bot < -delta_range or delta_bot > delta_range) {
1091 throw std::domain_error(
1092 "Parameter 'delta_mid_bot' must be >= -1 and <= +1");
1093 }
1094 /* Dielectric contrasts: the deltas should be either both -1 or both +1 when
1095 * no constant potential difference is applied. The case of two non-metallic
1096 * parallel boundaries can only be treated with a constant potential. */
1097 if (dielectric_contrast_on and not const_pot and
1098 (std::fabs(1. - delta_mid_top * delta_mid_bot) < ROUND_ERROR_PREC)) {
1099 throw std::domain_error("ELC with two parallel metallic boundaries "
1100 "requires the const_pot option");
1101 }
1102}
1103
1105 elc_data &&parameters, BaseSolver &&solver)
1106 : elc{parameters}, base_solver{solver} {
1107 adapt_solver();
1108}
1109
1110template <ChargeProtocol protocol, typename combined_ranges>
1111void charge_assign(elc_data const &elc, CoulombP3M &solver,
1112 combined_ranges const &p_q_pos_range) {
1113
1114 solver.prepare_fft_mesh(protocol == ChargeProtocol::BOTH or
1115 protocol == ChargeProtocol::IMAGE);
1116
1117 for (auto zipped : p_q_pos_range) {
1118 auto const p_q = boost::get<0>(zipped);
1119 auto const &p_pos = boost::get<1>(zipped);
1120 if (p_q != 0.) {
1121 // assign real charges
1122 if (protocol == ChargeProtocol::BOTH or
1123 protocol == ChargeProtocol::REAL) {
1124 solver.assign_charge(p_q, p_pos, false);
1125 }
1126 // assign image charges
1127 if (protocol == ChargeProtocol::BOTH or
1128 protocol == ChargeProtocol::IMAGE) {
1129 if (p_pos[2] < elc.space_layer) {
1130 auto const q_eff = elc.delta_mid_bot * p_q;
1131 solver.assign_charge(q_eff, {p_pos[0], p_pos[1], -p_pos[2]}, true);
1132 }
1133 if (p_pos[2] > (elc.box_h - elc.space_layer)) {
1134 auto const q_eff = elc.delta_mid_top * p_q;
1135 solver.assign_charge(
1136 q_eff, {p_pos[0], p_pos[1], 2. * elc.box_h - p_pos[2]}, true);
1137 }
1138 }
1139 }
1140 }
1141}
1142
1143template <ChargeProtocol protocol, typename combined_range>
1144void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver,
1145 combined_range const &p_q_pos_range) {
1146
1147 auto local_n = 0;
1148 auto local_q2 = 0.0;
1149 auto local_q = 0.0;
1150 for (auto zipped : p_q_pos_range) {
1151 auto const p_q = boost::get<0>(zipped);
1152 auto const &p_pos = boost::get<1>(zipped);
1153 if (p_q != 0.) {
1154 auto const p_z = p_pos[2];
1155
1156 if (protocol == ChargeProtocol::BOTH or
1157 protocol == ChargeProtocol::REAL) {
1158 local_n++;
1159 local_q2 += Utils::sqr(p_q);
1160 local_q += p_q;
1161 }
1162
1163 if (protocol == ChargeProtocol::BOTH or
1164 protocol == ChargeProtocol::IMAGE) {
1165 if (p_z < elc.space_layer) {
1166 local_n++;
1167 local_q2 += Utils::sqr(elc.delta_mid_bot * p_q);
1168 local_q += elc.delta_mid_bot * p_q;
1169 }
1170
1171 if (p_z > (elc.box_h - elc.space_layer)) {
1172 local_n++;
1173 local_q2 += Utils::sqr(elc.delta_mid_top * p_q);
1174 local_q += elc.delta_mid_top * p_q;
1175 }
1176 }
1177 }
1178 }
1179
1180 auto global_n = 0;
1181 auto global_q2 = 0.;
1182 auto global_q = 0.;
1183 boost::mpi::all_reduce(comm_cart, local_n, global_n, std::plus<>());
1184 boost::mpi::all_reduce(comm_cart, local_q2, global_q2, std::plus<>());
1185 boost::mpi::all_reduce(comm_cart, local_q, global_q, std::plus<>());
1186 solver.count_charged_particles_elc(global_n, global_q2, Utils::sqr(global_q));
1187}
1188
1190 ParticleRange const &particles) const {
1191 auto const energy = std::visit(
1192 [this, &particles](auto const &solver_ptr) {
1193 auto &solver = *solver_ptr;
1194 auto const &box_geo = *get_system().box_geo;
1195
1196 auto p_q_range = ParticlePropertyRange::charge_range(particles);
1197 auto p_pos_range = ParticlePropertyRange::pos_range(particles);
1198 auto p_q_pos_range = boost::combine(p_q_range, p_pos_range);
1199
1200 // assign the original charges (they may not have been assigned yet)
1201 solver.charge_assign(particles);
1202
1204 return solver.long_range_energy(particles);
1205 }
1206
1207 auto energy = 0.;
1208 energy += 0.5 * solver.long_range_energy(particles);
1209 energy +=
1210 0.5 * elc.dielectric_layers_self_energy(solver, box_geo, particles);
1211
1212 // assign both original and image charges
1213 charge_assign<ChargeProtocol::BOTH>(elc, solver, p_q_pos_range);
1214 modify_p3m_sums<ChargeProtocol::BOTH>(elc, solver, p_q_pos_range);
1215 energy += 0.5 * solver.long_range_energy(particles);
1216
1217 // assign only the image charges now
1218 charge_assign<ChargeProtocol::IMAGE>(elc, solver, p_q_pos_range);
1219 modify_p3m_sums<ChargeProtocol::IMAGE>(elc, solver, p_q_pos_range);
1220 energy -= 0.5 * solver.long_range_energy(particles);
1221
1222 // restore modified sums
1223 modify_p3m_sums<ChargeProtocol::REAL>(elc, solver, p_q_pos_range);
1224
1225 return energy;
1226 },
1227 base_solver);
1228 return energy + calc_energy(particles);
1229}
1230
1232 ParticleRange const &particles) const {
1233 std::visit(
1234 [this, &particles](auto const &solver_ptr) {
1235 auto &solver = *solver_ptr;
1236 auto p_q_range = ParticlePropertyRange::charge_range(particles);
1237 auto p_pos_range = ParticlePropertyRange::pos_range(particles);
1238 auto p_q_pos_range = boost::combine(p_q_range, p_pos_range);
1240 auto const &box_geo = *get_system().box_geo;
1241 modify_p3m_sums<ChargeProtocol::BOTH>(elc, solver, p_q_pos_range);
1242 charge_assign<ChargeProtocol::BOTH>(elc, solver, p_q_pos_range);
1243 elc.dielectric_layers_self_forces(solver, box_geo, particles);
1244 } else {
1245 solver.charge_assign(particles);
1246 }
1247 solver.add_long_range_forces(particles);
1249 modify_p3m_sums<ChargeProtocol::REAL>(elc, solver, p_q_pos_range);
1250 }
1251 },
1252 base_solver);
1253 add_force(particles);
1254}
1255
1256#endif
Vector implementation and trait types for boost qvm interoperability.
Utils::Vector3d const & length() const
Box length.
Utils::Vector3d const & length_inv() const
Inverse box length.
A range of particles.
base_type::size_type size() const
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.
Definition config.hpp:66
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.
static void addscale_vec(double *pdc_d, double scale, double const *pdc_s1, double const *pdc_s2, std::size_t size)
Definition elc.cpp:161
#define PQESSM
Definition elc.cpp:69
#define POQESP
Definition elc.cpp:60
static std::pair< std::size_t, std::size_t > prepare_sc_cache(ParticleRange const &particles, BoxGeometry const &box_geo, double far_cut)
Definition elc.cpp:127
#define PQESCP
Definition elc.cpp:66
static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
Definition elc.cpp:771
static auto calc_total_charge(CellStructure const &cell_structure)
Definition elc.cpp:975
static std::vector< double > partblk
temporary buffers for product decomposition
Definition elc.cpp:81
#define PQESCM
Definition elc.cpp:70
static void clear_vec(double *pdc, std::size_t size)
Definition elc.cpp:145
static double * block(double *p, std::size_t index, std::size_t size)
Definition elc.cpp:172
void setup_PoQ(elc_data const &elc, double prefactor, std::size_t index, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
Definition elc.cpp:501
static void distribute(std::size_t size)
Definition elc.cpp:176
void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver, combined_range const &p_q_pos_range)
Definition elc.cpp:1144
static double PoQ_energy(double omega, std::size_t n_part)
Definition elc.cpp:620
static std::vector< SCCache > scxcache
Cached sin/cos values along the x-axis and y-axis.
Definition elc.cpp:92
#define PQECCP
Definition elc.cpp:68
static double PQ_energy(double omega, std::size_t n_part)
Definition elc.cpp:812
void charge_assign(elc_data const &elc, CoulombP3M &solver, combined_ranges const &p_q_pos_range)
Definition elc.cpp:1111
#define POQECP
Definition elc.cpp:61
static std::vector< SCCache > scycache
Definition elc.cpp:93
static void setup_PQ(elc_data const &elc, double prefactor, std::size_t index_p, std::size_t index_q, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
Definition elc.cpp:641
static void copy_vec(double *pdc_d, double const *pdc_s, std::size_t size)
Definition elc.cpp:150
#define PQECCM
Definition elc.cpp:72
#define POQECM
Definition elc.cpp:63
static void add_vec(double *pdc_d, double const *pdc_s1, double const *pdc_s2, std::size_t size)
Definition elc.cpp:155
#define POQESM
Definition elc.cpp:62
static void scale_vec(double scale, double *pdc, std::size_t size)
Definition elc.cpp:167
void add_PoQ_force(ParticleRange const &particles)
Definition elc.cpp:601
ChargeProtocol
ELC charge sum/assign protocol: real charges, image charges, or both.
Definition elc.cpp:78
static std::vector< SCCache > calc_sc_cache(ParticleRange const &particles, std::size_t n_freq, double u)
Calculate cached sin/cos values for one direction.
Definition elc.cpp:107
#define PQESSP
Definition elc.cpp:65
#define PQECSP
Definition elc.cpp:67
static double gblcblk[8]
collected data from the other cells
Definition elc.cpp:83
PoQ
ELC axes (x and y directions)
Definition elc.cpp:76
#define PQECSM
Definition elc.cpp:71
ELC algorithm for long-range Coulomb interactions.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
ParticleRange particles(std::span< Cell *const > cells)
auto charge_range(ParticleRange const &particles)
auto pos_range(ParticleRange const &particles)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Vector< T, N > sqrt(Vector< T, N > const &a)
Definition Vector.hpp:357
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
P3M algorithm for long-range Coulomb interaction.
Describes a cell structure / cell system.
ParticleRange local_particles() const
P3M solver.
Definition p3m.hpp:55
virtual void count_charged_particles_elc(int, double, double)=0
virtual void prepare_fft_mesh(bool reset_weights)=0
virtual void assign_charge(double q, Utils::Vector3d const &real_pos, bool skip_cache)=0
Assign a single charge into the current charge grid.
std::variant< std::shared_ptr< CoulombP3M > > BaseSolver
Definition elc.hpp:169
BaseSolver base_solver
Electrostatics solver that is adapted.
Definition elc.hpp:175
ElectrostaticLayerCorrection(elc_data &&parameters, BaseSolver &&solver)
Definition elc.cpp:1104
BoxGeometry * m_box_geo
Definition elc.hpp:172
void add_long_range_forces(ParticleRange const &particles) const
Definition elc.cpp:1231
double long_range_energy(ParticleRange const &particles) const
Definition elc.cpp:1189
double b(double q, double z) const
Image sum from the bottom layer.
Definition elc.cpp:355
double dci
Definition elc.cpp:349
double shift
Definition elc.cpp:347
ImageSum(double delta, double shift, double lz)
Definition elc.cpp:351
double t(double q, double z) const
Image sum from the top layer.
Definition elc.cpp:360
double lz
Definition elc.cpp:348
double delta
Definition elc.cpp:346
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & q() const
Definition Particle.hpp:508
auto const & pos() const
Definition Particle.hpp:431
auto const & force() const
Definition Particle.hpp:435
auto const & id() const
Definition Particle.hpp:414
structure for caching sin and cos values
Definition elc.cpp:86
double c
Definition elc.cpp:87
double s
Definition elc.cpp:87
Parameters for the ELC method.
Definition elc.hpp:63
double dielectric_layers_self_energy(CoulombP3M const &p3m, BoxGeometry const &box_geo, ParticleRange const &particles) const
self energies of top and bottom layers with their virtual images
Definition elc.hpp:139
double maxPWerror
Maximal allowed pairwise error for the potential and force.
Definition elc.hpp:72
double pot_diff
Constant potential difference.
Definition elc.hpp:111
double box_h
Up to where particles can be found.
Definition elc.hpp:78
bool dielectric_contrast_on
Flag whether there is any dielectric contrast in the system.
Definition elc.hpp:94
elc_data(double maxPWerror, double gap_size, double far_cut, bool neutralize, double delta_top, double delta_bot, bool const_pot, double pot_diff)
Definition elc.cpp:1056
double space_box
The space that is finally left.
Definition elc.hpp:116
bool neutralize
Flag whether the box is neutralized by a homogeneous background.
Definition elc.hpp:104
double far_cut
Cutoff of the exponential sum.
Definition elc.hpp:84
double space_layer
Layer around the dielectric contrast in which we trick around.
Definition elc.hpp:114
bool far_calculated
Flag whether far_cut was set by the user, or calculated by ESPResSo.
Definition elc.hpp:91
double gap_size
Size of the empty gap.
Definition elc.hpp:76
double delta_mid_bot
dielectric contrast in the lower part of the simulation cell.
Definition elc.hpp:109
void dielectric_layers_self_forces(CoulombP3M const &p3m, BoxGeometry const &box_geo, ParticleRange const &particles) const
forces of particles in border layers with themselves
Definition elc.hpp:154
bool const_pot
Flag whether a constant potential difference is applied.
Definition elc.hpp:96
double far_cut2
Squared value of far_cut.
Definition elc.hpp:86
double delta_mid_top
dielectric contrast in the upper part of the simulation cell.
Definition elc.hpp:107