ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
elc.hpp
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/** @file
23 * @brief ELC algorithm for long-range Coulomb interactions.
24 *
25 * Implementation of the ELC method for the calculation of the electrostatic
26 * interaction in two dimensional periodic systems. For details on the method
27 * see MMM in general. The ELC method works together with any three-dimensional
28 * method, for example @ref p3m.hpp "P3M", with metallic boundary conditions.
29 */
30
31#pragma once
32
33#include "config/config.hpp"
34
35#ifdef P3M
36
37#include "actor/traits.hpp"
38
40
41#include "BoxGeometry.hpp"
42#include "Particle.hpp"
43#include "ParticleRange.hpp"
44
45#include <utils/Vector.hpp>
46#include <utils/math/sqr.hpp>
47
48#include <memory>
49#include <optional>
50#include <stdexcept>
51#include <string>
52#include <utility>
53#include <variant>
54
56
57namespace traits {
58template <>
60} // namespace traits
61
62/** @brief Parameters for the ELC method */
63struct elc_data {
64 elc_data(double maxPWerror, double gap_size, double far_cut, bool neutralize,
65 double delta_top, double delta_bot, bool const_pot, double pot_diff);
66
67 /**
68 * @brief Maximal allowed pairwise error for the potential and force.
69 * Note that this counts for the plain 1/r contribution
70 * alone, without the prefactor and the charge prefactor.
71 */
72 double maxPWerror;
73 /** Size of the empty gap. Note that ELC relies on the user to make sure
74 * that this condition is fulfilled.
75 */
76 double gap_size;
77 /** Up to where particles can be found. */
78 double box_h;
79 /**
80 * @brief Cutoff of the exponential sum.
81 * Since in all other MMM methods this is the far formula,
82 * it is given the same name here.
83 */
84 double far_cut;
85 /** Squared value of #far_cut. */
86 double far_cut2;
87 /** Flag whether #far_cut was set by the user, or calculated by ESPResSo.
88 * In the latter case, the cutoff will be adapted if important parameters,
89 * such as the box dimensions, change.
90 */
92
93 /// @brief Flag whether there is any dielectric contrast in the system.
95 /// @brief Flag whether a constant potential difference is applied.
97 /**
98 * @brief Flag whether the box is neutralized by a homogeneous background.
99 * If true, use a homogeneous neutralizing background for non-neutral
100 * systems. Unlike the 3D case, this background adds an additional
101 * force pointing towards the system center, and the gap size
102 * enters into the value of the forces, so be careful with this.
103 */
105
106 /// dielectric contrast in the upper part of the simulation cell.
108 /// dielectric contrast in the lower part of the simulation cell.
110 /// @brief Constant potential difference.
111 double pot_diff;
112
113 /** Layer around the dielectric contrast in which we trick around. */
115 /** The space that is finally left. */
116 double space_box;
117
118 /// pairwise contributions from the lowest and top layers
119 template <typename Kernel>
121 BoxGeometry const &box_geo,
122 Utils::Vector3d const &pos1,
123 Utils::Vector3d const &pos2, double q1q2,
124 Kernel &&kernel) const {
125 if (pos1[2] < space_layer) {
126 auto const q_eff = delta_mid_bot * q1q2;
127 auto const d = box_geo.get_mi_vector(pos2, {pos1[0], pos1[1], -pos1[2]});
128 kernel(q_eff, d);
129 }
130 if (pos1[2] > (box_h - space_layer)) {
131 auto const q_eff = delta_mid_top * q1q2;
132 auto const z = 2. * box_h - pos1[2];
133 auto const d = box_geo.get_mi_vector(pos2, {pos1[0], pos1[1], z});
134 kernel(q_eff, d);
135 }
136 }
137
138 /// self energies of top and bottom layers with their virtual images
140 BoxGeometry const &box_geo,
141 ParticleRange const &particles) const {
142 auto energy = 0.;
143 for (auto const &p : particles) {
145 p3m, box_geo, p.pos(), p.pos(), Utils::sqr(p.q()),
146 [&](double q1q2, Utils::Vector3d const &d) {
147 energy += p3m.pair_energy(q1q2, d.norm());
148 });
149 }
150 return energy;
151 }
152
153 /// forces of particles in border layers with themselves
155 BoxGeometry const &box_geo,
156 ParticleRange const &particles) const {
157 for (auto &p : particles) {
159 p3m, box_geo, p.pos(), p.pos(), Utils::sqr(p.q()),
160 [&](double q1q2, Utils::Vector3d const &d) {
161 p.force() += p3m.pair_force(q1q2, d, d.norm());
162 });
163 }
164 }
165};
166
168 : public Coulomb::Actor<ElectrostaticLayerCorrection> {
169 using BaseSolver = std::variant<std::shared_ptr<CoulombP3M>>;
170
173
174 /** Electrostatics solver that is adapted. */
176
177 ElectrostaticLayerCorrection(elc_data &&parameters, BaseSolver &&solver);
178
180 visit_base_solver(
181 [this](auto &solver) { solver->bind_system(m_system.lock()); });
182 sanity_checks_periodicity();
183 sanity_checks_cell_structure();
185 sanity_checks_dielectric_contrasts();
186 /* Most ELC parameters do not depend on the P3M parameters,
187 * but the P3M parameters depend on the ELC parameters during tuning,
188 * therefore ELC needs to be tuned before P3M. */
189 recalc_box_h();
190 recalc_far_cut();
191 visit_base_solver([](auto &solver) { solver->on_activation(); });
192 /* With dielectric contrasts, the ELC space layer depends
193 * on the P3M real-space cutoff, and the P3M FFT parameters
194 * depend on the ELC space layer */
196 recalc_space_layer();
197 visit_base_solver([](auto &actor) { actor->init(); });
198 }
199 }
200 /** @brief Recalculate all box-length-dependent parameters. */
202 visit_base_solver([](auto &actor) { actor->on_boxl_change(); });
203 recalc_box_h();
204 recalc_far_cut();
205 recalc_space_layer();
206 }
207 void on_node_grid_change() const {
208 visit_base_solver([](auto &solver) { solver->on_node_grid_change(); });
209 }
211 sanity_checks_periodicity();
212 visit_base_solver([](auto &solver) { solver->on_periodicity_change(); });
213 }
215 sanity_checks_cell_structure();
216 visit_base_solver([](auto &solver) { solver->on_cell_structure_change(); });
217 recalc_box_h();
218 recalc_far_cut();
220 recalc_space_layer();
221 visit_base_solver([](auto &actor) { actor->init(); });
222 }
223 }
224 /** @brief Recalculate all derived parameters. */
225 void init() {
226 visit_base_solver([](auto &actor) { actor->init(); });
227 recalc_box_h();
228 recalc_far_cut();
229 recalc_space_layer();
230 }
231
232 void sanity_checks() const {
233 sanity_checks_periodicity();
234 sanity_checks_cell_structure();
236 sanity_checks_dielectric_contrasts();
237 visit_base_solver([](auto &actor) { actor->sanity_checks(); });
238 }
239
240 /**
241 * @brief Veto real-space cutoff values that are incompatible with ELC.
242 * When ELC is used with dielectric contrasts, the short-range cutoff needs
243 * to be smaller than the gap size to allow placement of the image charges.
244 */
245 std::optional<std::string> veto_r_cut(double r_cut) const {
246 if (elc.dielectric_contrast_on and r_cut >= elc.gap_size) {
247 return {std::string("conflict with ELC w/ dielectric contrasts")};
248 }
249 return {};
250 }
251
252 /** @brief Calculate short-range pair energy correction. */
253 double pair_energy_correction(Particle const &p1, Particle const &p2,
254 double q1q2) const {
255 double energy = 0.;
257 energy = std::visit(
258 [this, &p1, &p2, q1q2](auto &p3m_ptr) {
259 auto const &pos1 = p1.pos();
260 auto const &pos2 = p2.pos();
261 auto const &p3m = *p3m_ptr;
262 auto energy = 0.;
264 p3m, *m_box_geo, pos1, pos2, q1q2,
265 [&](double q_eff, Utils::Vector3d const &d) {
266 energy += p3m.pair_energy(q_eff, d.norm());
267 });
269 p3m, *m_box_geo, pos2, pos1, q1q2,
270 [&](double q_eff, Utils::Vector3d const &d) {
271 energy += p3m.pair_energy(q_eff, d.norm());
272 });
273 return energy / 2.;
274 },
276 }
277 return energy;
278 }
279
280 /** @brief Add short-range pair force corrections. */
282 double q1q2) const {
284 std::visit(
285 [this, &p1, &p2, q1q2](auto &p3m_ptr) {
286 auto const &pos1 = p1.pos();
287 auto const &pos2 = p2.pos();
288 auto const &p3m = *p3m_ptr;
290 p3m, *m_box_geo, pos1, pos2, q1q2,
291 [&](double q_eff, Utils::Vector3d const &d) {
292 p1.force() += p3m.pair_force(q_eff, d, d.norm());
293 });
295 p3m, *m_box_geo, pos2, pos1, q1q2,
296 [&](double q_eff, Utils::Vector3d const &d) {
297 p2.force() += p3m.pair_force(q_eff, d, d.norm());
298 });
299 },
301 }
302 }
303
304 double long_range_energy(ParticleRange const &particles) const;
305 void add_long_range_forces(ParticleRange const &particles) const;
306
307private:
308 /** Check if a charged particle is in the gap region. */
309 void check_gap(Particle const &p) const;
310 double tune_far_cut() const;
311 void adapt_solver();
312 /** pairwise contributions from the lowest and top layers to the energy */
313 double dipole_energy(ParticleRange const &particles) const;
314 void add_dipole_force(ParticleRange const &particles) const;
315 double z_energy(ParticleRange const &particles) const;
316 void add_z_force(ParticleRange const &particles) const;
317
318 void recalc_box_h();
319 void recalc_far_cut() {
320 if (elc.far_calculated) {
321 elc.far_cut = tune_far_cut();
322 }
324 }
325 void recalc_space_layer();
326
327 void sanity_checks_cell_structure() const {}
328 void sanity_checks_periodicity() const;
329 void sanity_checks_dielectric_contrasts() const;
330
331 /// the force calculation
332 void add_force(ParticleRange const &particles) const;
333 /// the energy calculation
334 double calc_energy(ParticleRange const &particles) const;
335
336 template <class Visitor> void visit_base_solver(Visitor &&visitor) const {
337 std::visit(visitor, base_solver);
338 }
339};
340
341#endif // P3M
Vector implementation and trait types for boost qvm interoperability.
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
A range of particles.
std::weak_ptr< System > m_system
T norm() const
Definition Vector.hpp:138
This file contains the defaults for ESPResSo.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
P3M algorithm for long-range Coulomb interaction.
P3M solver.
Definition p3m.hpp:55
double pair_energy_correction(Particle const &p1, Particle const &p2, double q1q2) const
Calculate short-range pair energy correction.
Definition elc.hpp:253
void on_boxl_change()
Recalculate all box-length-dependent parameters.
Definition elc.hpp:201
void on_periodicity_change() const
Definition elc.hpp:210
void sanity_checks() const
Definition elc.hpp:232
std::optional< std::string > veto_r_cut(double r_cut) const
Veto real-space cutoff values that are incompatible with ELC.
Definition elc.hpp:245
std::variant< std::shared_ptr< CoulombP3M > > BaseSolver
Definition elc.hpp:169
BaseSolver base_solver
Electrostatics solver that is adapted.
Definition elc.hpp:175
void init()
Recalculate all derived parameters.
Definition elc.hpp:225
void add_pair_force_corrections(Particle &p1, Particle &p2, double q1q2) const
Add short-range pair force corrections.
Definition elc.hpp:281
void on_node_grid_change() const
Definition elc.hpp:207
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
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & pos() const
Definition Particle.hpp:431
auto const & force() const
Definition Particle.hpp:435
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
void dielectric_layers_contribution(CoulombP3M const &, BoxGeometry const &box_geo, Utils::Vector3d const &pos1, Utils::Vector3d const &pos2, double q1q2, Kernel &&kernel) const
pairwise contributions from the lowest and top layers
Definition elc.hpp:120
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
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
Whether an actor is a layer correction method.
Definition traits.hpp:26