ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
icc.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/** \file
23 * Functions to compute the electric field acting on the induced charges,
24 * excluding forces other than the electrostatic ones. Detailed information
25 * about the ICC* method is included in the corresponding header file
26 * \ref icc.hpp.
27 */
28
29#include "config/config.hpp"
30
31#ifdef ELECTROSTATICS
32
33#include "icc.hpp"
34
35#include "Particle.hpp"
36#include "ParticleRange.hpp"
37#include "PropagationMode.hpp"
38#include "actor/visitors.hpp"
40#include "communication.hpp"
44#include "errorhandling.hpp"
46#include "system/System.hpp"
47
48#include <boost/mpi/collectives/all_reduce.hpp>
49#include <boost/mpi/operations.hpp>
50
51#include <algorithm>
52#include <cmath>
53#include <cstddef>
54#include <limits>
55#include <numbers>
56#include <stdexcept>
57#include <variant>
58#include <vector>
59
60/** Calculate the electrostatic forces between source charges (= real charges)
61 * and wall charges. For each electrostatic method, the proper functions
62 * for short- and long-range parts are called. Long-range parts are calculated
63 * directly, short-range parts need helper functions according to the particle
64 * data organisation. This is a modified version of
65 * @ref System::System::calculate_forces.
66 */
67static void force_calc_icc(
68 CellStructure &cell_structure, ParticleRange const &particles,
69 ParticleRange const &ghost_particles,
72 // reset forces
73 for (auto &p : particles) {
74 p.force() = {};
75 }
76 for (auto &p : ghost_particles) {
77 p.force() = {};
78 }
79
80 // calc ICC forces
81 cell_structure.non_bonded_loop(
82 [coulomb_kernel_ptr = get_ptr(coulomb_kernel),
83 elc_kernel_ptr = get_ptr(elc_kernel)](Particle &p1, Particle &p2,
84 Distance const &d) {
85 auto const q1q2 = p1.q() * p2.q();
86 if (q1q2 != 0.) {
87 auto force = (*coulomb_kernel_ptr)(q1q2, d.vec21, std::sqrt(d.dist2));
88 p1.force() += force;
89 p2.force() -= force;
90#ifdef P3M
91 if (elc_kernel_ptr) {
92 (*elc_kernel_ptr)(p1, p2, q1q2);
93 }
94#endif // P3M
95 }
96 });
97}
98
99void ICCStar::iteration(CellStructure &cell_structure,
100 ParticleRange const &particles,
101 ParticleRange const &ghost_particles) {
102
103 try {
104 sanity_check();
105 } catch (std::runtime_error const &err) {
106 runtimeErrorMsg() << err.what();
107 return;
108 }
109
110 auto &system = get_system();
111 auto const &coulomb = system.coulomb;
112 auto const prefactor = std::visit(
113 [](auto const &ptr) { return ptr->prefactor; }, *coulomb.impl->solver);
114 auto const pref = 1. / (prefactor * 2. * std::numbers::pi);
115 auto const kernel = coulomb.pair_force_kernel();
116 auto const elc_kernel = coulomb.pair_force_elc_kernel();
118
119 auto global_max_rel_diff = 0.;
120
121 for (int j = 0; j < icc_cfg.max_iterations; j++) {
122 auto charge_density_max = 0.;
123
124 // calculate electrostatic forces (SR+LR) excluding self-interactions
125 force_calc_icc(cell_structure, particles, ghost_particles, kernel,
126 elc_kernel);
127 system.coulomb.calc_long_range_force(particles);
128 cell_structure.ghosts_reduce_forces();
129
130 auto max_rel_diff = 0.;
131
132 for (auto &p : particles) {
133 auto const pid = p.id();
134 if (pid >= icc_cfg.first_id and pid < icc_cfg.n_icc + icc_cfg.first_id) {
135 if (p.q() == 0.) {
137 << "ICC found zero electric charge on a particle. This must "
138 "never happen";
139 break;
140 }
141 auto const id = p.id() - icc_cfg.first_id;
142 /* the dielectric-related prefactor: */
143 auto const eps_in = icc_cfg.epsilons[id];
144 auto const eps_out = icc_cfg.eps_out;
145 auto const del_eps = (eps_in - eps_out) / (eps_in + eps_out);
146 /* calculate the electric field at the certain position */
147 auto const local_e_field = p.force() / p.q() + icc_cfg.ext_field;
148
149 if (local_e_field.norm2() == 0.) {
151 << "ICC found zero electric field on a charge. This must "
152 "never happen";
153 }
154
155 auto const charge_density_old = p.q() / icc_cfg.areas[id];
156
157 charge_density_max =
158 std::max(charge_density_max, std::abs(charge_density_old));
159
160 auto const charge_density_update =
161 del_eps * pref * (local_e_field * icc_cfg.normals[id]) +
163 icc_cfg.sigmas[id];
164 /* relative variation: never use an estimator which can be negative
165 * here */
166 auto const charge_density_new =
167 (1. - icc_cfg.relaxation) * charge_density_old +
168 (icc_cfg.relaxation) * charge_density_update;
169
170 /* Take the largest error to check for convergence */
171 auto const relative_difference =
172 std::abs((charge_density_new - charge_density_old) /
173 (charge_density_max +
174 std::abs(charge_density_new + charge_density_old)));
175
176 max_rel_diff = std::max(max_rel_diff, relative_difference);
177
178 p.q() = charge_density_new * icc_cfg.areas[id];
179
180 /* check if the charge now is more than 1e6, to determine if ICC still
181 * leads to reasonable results. This is kind of an arbitrary measure
182 * but does a good job of spotting divergence! */
183 if (std::abs(p.q()) > 1e6) {
185 << "Particle with id " << p.id() << " has a charge (q=" << p.q()
186 << ") that is too large for the ICC algorithm";
187
188 max_rel_diff = std::numeric_limits<double>::max();
189 break;
190 }
191 }
192 }
193
194 /* Update charges on ghosts. */
196
198
199 boost::mpi::all_reduce(comm_cart, max_rel_diff, global_max_rel_diff,
200 boost::mpi::maximum<double>());
201
202 if (global_max_rel_diff < icc_cfg.convergence)
203 break;
204 }
205
206 if (global_max_rel_diff > icc_cfg.convergence) {
208 << "ICC failed to converge in the given number of maximal steps.";
209 }
210
211 system.on_particle_charge_change();
212}
213
215 if (n_icc <= 0)
216 throw std::domain_error("Parameter 'n_icc' must be >= 1");
217 if (convergence <= 0.)
218 throw std::domain_error("Parameter 'convergence' must be > 0");
219 if (relaxation < 0. or relaxation > 2.)
220 throw std::domain_error("Parameter 'relaxation' must be >= 0 and <= 2");
221 if (max_iterations <= 0)
222 throw std::domain_error("Parameter 'max_iterations' must be > 0");
223 if (first_id < 0)
224 throw std::domain_error("Parameter 'first_id' must be >= 0");
225 if (eps_out <= 0.)
226 throw std::domain_error("Parameter 'eps_out' must be > 0");
227 if (areas.size() != static_cast<std::size_t>(n_icc))
228 throw std::invalid_argument("Parameter 'areas' has incorrect shape");
229 if (epsilons.size() != static_cast<std::size_t>(n_icc))
230 throw std::invalid_argument("Parameter 'epsilons' has incorrect shape");
231 if (sigmas.size() != static_cast<std::size_t>(n_icc))
232 throw std::invalid_argument("Parameter 'sigmas' has incorrect shape");
233 if (normals.size() != static_cast<std::size_t>(n_icc))
234 throw std::invalid_argument("Parameter 'normals' has incorrect shape");
235}
236
238 data.sanity_checks();
239 icc_cfg = std::move(data);
240}
241
243 sanity_check();
244 auto &system = get_system();
245 system.on_particle_charge_change();
246}
247
249 template <typename T> void operator()(std::shared_ptr<T> const &) const {}
250#ifdef P3M
251#ifdef CUDA
252 void operator()(std::shared_ptr<CoulombP3M> const &p) const {
253 if (p->is_gpu()) {
254 throw std::runtime_error("ICC does not work with P3MGPU");
255 }
256 }
257#endif // CUDA
258 void
259 operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &actor) const {
260 if (actor->elc.dielectric_contrast_on) {
261 throw std::runtime_error("ICC conflicts with ELC dielectric contrast");
262 }
263 std::visit(*this, actor->base_solver);
264 }
265#endif // P3M
266 [[noreturn]] void operator()(std::shared_ptr<DebyeHueckel> const &) const {
267 throw std::runtime_error("ICC does not work with DebyeHueckel.");
268 }
269 [[noreturn]] void operator()(std::shared_ptr<ReactionField> const &) const {
270 throw std::runtime_error("ICC does not work with ReactionField.");
271 }
272};
273
276#ifdef NPT
277 if (get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO) {
278 throw std::runtime_error("ICC does not work in the NPT ensemble");
279 }
280#endif
281}
282
284 auto &system = get_system();
285 if (system.coulomb.impl->solver) {
286 std::visit(SanityChecksICC(), *system.coulomb.impl->solver);
287 } else {
288 throw std::runtime_error("An electrostatics solver is needed by ICC");
289 }
290}
291
293 if (coulomb.impl->extension) {
294 if (auto icc = std::get_if<std::shared_ptr<ICCStar>>(
295 get_ptr(coulomb.impl->extension))) {
296 (**icc).iteration(*cell_structure, cell_structure->local_particles(),
297 cell_structure->ghost_particles());
298 }
299 }
300}
301
302#endif // ELECTROSTATICS
@ INTEG_METHOD_NPT_ISO
A range of particles.
void update_icc_particles()
Definition icc.cpp:292
boost::mpi::communicator comm_cart
The communicator.
This file contains the defaults for ESPResSo.
const T * get_ptr(std::optional< T > const &opt)
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
static void force_calc_icc(CellStructure &cell_structure, ParticleRange const &particles, ParticleRange const &ghost_particles, Coulomb::ShortRangeForceKernel::result_type const &coulomb_kernel, Coulomb::ShortRangeForceCorrectionsKernel::result_type const &elc_kernel)
Calculate the electrostatic forces between source charges (= real charges) and wall charges.
Definition icc.cpp:67
ICC is a method that allows to take into account the influence of arbitrarily shaped dielectric inter...
@ DATA_PART_PROPERTIES
Particle::p.
P3M algorithm for long-range Coulomb interaction.
Describes a cell structure / cell system.
void ghosts_update(unsigned data_parts)
Update ghost particles.
void ghosts_reduce_forces()
Add forces from ghost particles to real particles.
void non_bonded_loop(PairKernel pair_kernel)
Non-bonded pair loop.
std::optional< kernel_type > result_type
std::optional< kernel_type > result_type
Distance vector and length handed to pair kernels.
void sanity_checks_active_solver() const
Definition icc.cpp:283
void iteration(CellStructure &cell_structure, ParticleRange const &particles, ParticleRange const &ghost_particles)
The main iterative scheme, where the surface element charges are calculated self-consistently.
Definition icc.cpp:99
icc_data icc_cfg
ICC parameters.
Definition icc.hpp:92
void on_activation() const
Definition icc.cpp:242
void sanity_check() const
Definition icc.cpp:274
ICCStar(icc_data data)
Definition icc.cpp:237
Struct holding all information for one particle.
Definition Particle.hpp:395
void operator()(std::shared_ptr< T > const &) const
Definition icc.cpp:249
void operator()(std::shared_ptr< ReactionField > const &) const
Definition icc.cpp:269
void operator()(std::shared_ptr< CoulombP3M > const &p) const
Definition icc.cpp:252
void operator()(std::shared_ptr< DebyeHueckel > const &) const
Definition icc.cpp:266
void operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
Definition icc.cpp:259
ICC data structure.
Definition icc.hpp:61
double convergence
convergence criteria
Definition icc.hpp:75
int first_id
first ICC particle id
Definition icc.hpp:85
double relaxation
relaxation parameter
Definition icc.hpp:81
std::vector< Utils::Vector3d > normals
surface normal vectors
Definition icc.hpp:77
int max_iterations
maximum number of iterations
Definition icc.hpp:65
int n_icc
First id of ICC particle.
Definition icc.hpp:63
double eps_out
bulk dielectric constant
Definition icc.hpp:67
void sanity_checks() const
Definition icc.cpp:214
std::vector< double > sigmas
surface charge density of the particles
Definition icc.hpp:73
std::vector< double > epsilons
dielectric constants of the particles
Definition icc.hpp:71
Utils::Vector3d ext_field
external electric field
Definition icc.hpp:79
int citeration
last number of iterations
Definition icc.hpp:83
std::vector< double > areas
areas of the particles
Definition icc.hpp:69