40#include "communication.hpp"
46#include "system/System.hpp"
48#include <boost/mpi/collectives/all_reduce.hpp>
49#include <boost/mpi/operations.hpp>
73 for (
auto &p : particles) {
76 for (
auto &p : ghost_particles) {
82 [coulomb_kernel_ptr =
get_ptr(coulomb_kernel),
85 auto const q1q2 = p1.q() * p2.q();
87 auto force = (*coulomb_kernel_ptr)(q1q2, d.vec21, std::sqrt(d.dist2));
92 (*elc_kernel_ptr)(p1, p2, q1q2);
105 }
catch (std::runtime_error
const &err) {
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();
119 auto global_max_rel_diff = 0.;
122 auto charge_density_max = 0.;
125 force_calc_icc(cell_structure, particles, ghost_particles, kernel,
127 system.coulomb.calc_long_range_force(particles);
130 auto max_rel_diff = 0.;
132 for (
auto &p : particles) {
133 auto const pid = p.id();
137 <<
"ICC found zero electric charge on a particle. This must "
145 auto const del_eps = (eps_in - eps_out) / (eps_in + eps_out);
149 if (local_e_field.norm2() == 0.) {
151 <<
"ICC found zero electric field on a charge. This must "
155 auto const charge_density_old = p.q() /
icc_cfg.
areas[id];
158 std::max(charge_density_max, std::abs(charge_density_old));
160 auto const charge_density_update =
166 auto const charge_density_new =
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)));
176 max_rel_diff = std::max(max_rel_diff, relative_difference);
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";
188 max_rel_diff = std::numeric_limits<double>::max();
199 boost::mpi::all_reduce(
comm_cart, max_rel_diff, global_max_rel_diff,
200 boost::mpi::maximum<double>());
208 <<
"ICC failed to converge in the given number of maximal steps.";
211 system.on_particle_charge_change();
216 throw std::domain_error(
"Parameter 'n_icc' must be >= 1");
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");
222 throw std::domain_error(
"Parameter 'max_iterations' must be > 0");
224 throw std::domain_error(
"Parameter 'first_id' must be >= 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");
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");
234 throw std::invalid_argument(
"Parameter 'normals' has incorrect shape");
245 system.on_particle_charge_change();
249 template <
typename T>
void operator()(std::shared_ptr<T>
const &)
const {}
252 void operator()(std::shared_ptr<CoulombP3M>
const &p)
const {
254 throw std::runtime_error(
"ICC does not work with P3MGPU");
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");
263 std::visit(*
this, actor->base_solver);
266 [[noreturn]]
void operator()(std::shared_ptr<DebyeHueckel>
const &)
const {
267 throw std::runtime_error(
"ICC does not work with DebyeHueckel.");
269 [[noreturn]]
void operator()(std::shared_ptr<ReactionField>
const &)
const {
270 throw std::runtime_error(
"ICC does not work with ReactionField.");
278 throw std::runtime_error(
"ICC does not work in the NPT ensemble");
285 if (system.coulomb.impl->solver) {
288 throw std::runtime_error(
"An electrostatics solver is needed by ICC");
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());
void update_icc_particles()
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.
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
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.
icc_data icc_cfg
ICC parameters.
void on_activation() const
void sanity_check() const
Struct holding all information for one particle.
void operator()(std::shared_ptr< T > const &) const
void operator()(std::shared_ptr< ReactionField > const &) const
void operator()(std::shared_ptr< CoulombP3M > const &p) const
void operator()(std::shared_ptr< DebyeHueckel > const &) const
void operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
double convergence
convergence criteria
int first_id
first ICC particle id
double relaxation
relaxation parameter
std::vector< Utils::Vector3d > normals
surface normal vectors
int max_iterations
maximum number of iterations
int n_icc
First id of ICC particle.
double eps_out
bulk dielectric constant
void sanity_checks() const
std::vector< double > sigmas
surface charge density of the particles
std::vector< double > epsilons
dielectric constants of the particles
Utils::Vector3d ext_field
external electric field
int citeration
last number of iterations
std::vector< double > areas
areas of the particles