ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
coulomb.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#include "config/config.hpp"
21
23
24#ifdef ELECTROSTATICS
25
27
28#include "ParticleRange.hpp"
30#include "actor/visitors.hpp"
32#include "communication.hpp"
34#include "errorhandling.hpp"
35#include "system/System.hpp"
36
37#include <utils/Vector.hpp>
38#include <utils/demangle.hpp>
39
40#include <boost/accumulators/accumulators.hpp>
41#include <boost/accumulators/statistics/sum_kahan.hpp>
42#include <boost/mpi/collectives/broadcast.hpp>
43#include <boost/mpi/collectives/gather.hpp>
44
45#include <algorithm>
46#include <cassert>
47#include <cmath>
48#include <cstdio>
49#include <iomanip>
50#include <limits>
51#include <memory>
52#include <optional>
53#include <sstream>
54#include <stdexcept>
55#include <type_traits>
56#include <variant>
57#include <vector>
58
59namespace Coulomb {
60
62 impl = std::make_unique<Implementation>();
64}
65
67
69 if (impl->solver) {
70 std::visit([](auto const &ptr) { ptr->sanity_checks(); }, *impl->solver);
71 }
72}
73
76 if (impl->solver) {
77 visit_try_catch([](auto &ptr) { ptr->init(); }, *impl->solver);
78 }
79}
80
82 if (impl->solver) {
83 visit_try_catch([](auto &ptr) { ptr->on_boxl_change(); }, *impl->solver);
84 }
85}
86
88 if (impl->solver) {
89 std::visit([](auto &ptr) { ptr->on_node_grid_change(); }, *impl->solver);
90 }
91}
92
94 if (impl->solver) {
95 visit_try_catch([](auto &ptr) { ptr->on_periodicity_change(); },
96 *impl->solver);
97 }
98}
99
101 if (impl->solver) {
102 visit_try_catch([](auto &ptr) { ptr->on_cell_structure_change(); },
103 *impl->solver);
104 }
105}
106
108 explicit LongRangePressure(ParticleRange const &particles)
109 : m_particles{particles} {}
110
111#ifdef P3M
112 auto operator()(std::shared_ptr<CoulombP3M> const &actor) const {
113 return actor->long_range_pressure(m_particles);
114 }
115#endif // P3M
116
117 auto operator()(std::shared_ptr<DebyeHueckel> const &) const {
118 return Utils::Vector9d{};
119 }
120
121 auto operator()(std::shared_ptr<ReactionField> const &) const {
122 return Utils::Vector9d{};
123 }
124
125 template <typename T,
126 std::enable_if_t<!traits::has_pressure<T>::value> * = nullptr>
127 auto operator()(std::shared_ptr<T> const &) const {
128 runtimeWarningMsg() << "Pressure calculation not implemented by "
129 << "electrostatics method " << Utils::demangle<T>();
130 return Utils::Vector9d{};
131 }
132
133private:
134 ParticleRange const &m_particles;
135};
136
139 if (impl->solver) {
140 return std::visit(LongRangePressure(particles), *impl->solver);
141 }
142 return {};
143}
144
146#ifdef P3M
147 auto operator()(std::shared_ptr<CoulombP3M> const &actor) const {
148 return actor->p3m_params.r_cut;
149 }
150 auto
151 operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &actor) const {
152 return std::max(actor->elc.space_layer,
153 std::visit(*this, actor->base_solver));
154 }
155#endif // P3M
156 auto operator()(std::shared_ptr<CoulombMMM1D> const &) const {
157 return std::numeric_limits<double>::infinity();
158 }
159#ifdef SCAFACOS
160 auto operator()(std::shared_ptr<CoulombScafacos> const &actor) const {
161 return actor->get_r_cut();
162 }
163#endif // SCAFACOS
164 auto operator()(std::shared_ptr<ReactionField> const &actor) const {
165 return actor->r_cut;
166 }
167 auto operator()(std::shared_ptr<DebyeHueckel> const &actor) const {
168 return actor->r_cut;
169 }
170};
171
172double Solver::cutoff() const {
173 if (impl->solver) {
174 return std::visit(ShortRangeCutoff(), *impl->solver);
175 }
176 return -1.0;
177}
178
180 template <typename T> void operator()(std::shared_ptr<T> const &) const {}
181
182#ifdef P3M
183 void operator()(std::shared_ptr<CoulombP3M> const &actor) const {
184 actor->count_charged_particles();
185 }
186 void
187 operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &actor) const {
188 std::visit(*this, actor->base_solver);
189 }
190#endif // P3M
191};
192
195 if (impl->solver) {
196 std::visit(EventOnObservableCalc(), *impl->solver);
197 }
199 }
200}
201
203 explicit LongRangeForce(ParticleRange const &particles)
204 : m_particles(particles) {}
205
206#ifdef P3M
207 void operator()(std::shared_ptr<CoulombP3M> const &actor) const {
208 actor->add_long_range_forces(m_particles);
209 }
210 void
211 operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &actor) const {
212 actor->add_long_range_forces(m_particles);
213 }
214#endif // P3M
215#ifdef SCAFACOS
216 void operator()(std::shared_ptr<CoulombScafacos> const &actor) const {
217 actor->add_long_range_forces();
218 }
219#endif
220 /* Several algorithms only provide near-field kernels */
221 void operator()(std::shared_ptr<CoulombMMM1D> const &) const {}
222 void operator()(std::shared_ptr<DebyeHueckel> const &) const {}
223 void operator()(std::shared_ptr<ReactionField> const &) const {}
224
225private:
226 ParticleRange const &m_particles;
227};
228
230 explicit LongRangeEnergy(ParticleRange const &particles)
231 : m_particles(particles) {}
232
233#ifdef P3M
234 auto operator()(std::shared_ptr<CoulombP3M> const &actor) const {
235 return actor->long_range_energy(m_particles);
236 }
237 auto
238 operator()(std::shared_ptr<ElectrostaticLayerCorrection> const &actor) const {
239 return actor->long_range_energy(m_particles);
240 }
241#endif // P3M
242#ifdef SCAFACOS
243 auto operator()(std::shared_ptr<CoulombScafacos> const &actor) const {
244 return actor->long_range_energy();
245 }
246#endif
247 /* Several algorithms only provide near-field kernels */
248 auto operator()(std::shared_ptr<CoulombMMM1D> const &) const { return 0.; }
249 auto operator()(std::shared_ptr<DebyeHueckel> const &) const { return 0.; }
250 auto operator()(std::shared_ptr<ReactionField> const &) const { return 0.; }
251
252private:
253 ParticleRange const &m_particles;
254};
255
256void Solver::calc_long_range_force(ParticleRange const &particles) const {
257 if (impl->solver) {
258 std::visit(LongRangeForce(particles), *impl->solver);
259 }
260}
261
262double Solver::calc_energy_long_range(ParticleRange const &particles) const {
263 if (impl->solver) {
264 return std::visit(LongRangeEnergy(particles), *impl->solver);
265 }
266 return 0.;
267}
268
269/** @brief Compute the net charge rescaled by the smallest non-zero charge. */
270static auto calc_charge_excess_ratio(std::vector<double> const &charges) {
271 using namespace boost::accumulators;
272 using KahanSum = accumulator_set<double, features<tag::sum_kahan>>;
273
274 KahanSum q_sum;
275 auto q_min = std::numeric_limits<double>::infinity();
276
277 for (auto const q : charges) {
278 if (q != 0.) {
279 q_sum(q);
280 q_min = std::min(q_min, std::abs(q));
281 }
282 }
283
284 return std::abs(sum_kahan(q_sum)) / q_min;
285}
286
288 double relative_tolerance) {
289 // collect non-zero particle charges from all nodes
290 std::vector<double> local_charges;
291 for (auto const &p : system.cell_structure->local_particles()) {
292 local_charges.push_back(p.q());
293 }
294 std::vector<std::vector<double>> node_charges;
295 boost::mpi::gather(comm_cart, local_charges, node_charges, 0);
296
297 // run Kahan sum on charges
298 auto excess_ratio = 0.;
299 if (this_node == 0) {
300 auto charges = std::move(local_charges);
301 for (auto it = ++node_charges.begin(); it != node_charges.end(); ++it) {
302 charges.insert(charges.end(), it->begin(), it->end());
303 }
304 excess_ratio = calc_charge_excess_ratio(charges);
305 }
306 boost::mpi::broadcast(comm_cart, excess_ratio, 0);
307
308 if (excess_ratio >= relative_tolerance) {
309 std::ostringstream serializer;
310 serializer << std::scientific << std::setprecision(4);
311 serializer << excess_ratio;
312 throw std::runtime_error(
313 "The system is not charge neutral. Please neutralize the system "
314 "before adding a new actor by adding the corresponding counterions "
315 "to the system. Alternatively you can turn off the electroneutrality "
316 "check by supplying check_neutrality=False when creating the actor. "
317 "In this case you may be simulating a non-neutral system which will "
318 "affect physical observables like e.g. the pressure, the chemical "
319 "potentials of charged species or potential energies of the system. "
320 "Since simulations of non charge neutral systems are special please "
321 "make sure you know what you are doing. The relative charge excess "
322 "was " +
323 serializer.str());
324 }
325}
326
327} // namespace Coulomb
328#endif // ELECTROSTATICS
Vector implementation and trait types for boost qvm interoperability.
A range of particles.
Main system class.
std::shared_ptr< CellStructure > cell_structure
Coulomb::Solver coulomb
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeWarningMsg()
ICC is a method that allows to take into account the influence of arbitrarily shaped dielectric inter...
void check_charge_neutrality(System::System const &system, double relative_tolerance)
Check if the system is charge-neutral.
Definition coulomb.cpp:287
static auto calc_charge_excess_ratio(std::vector< double > const &charges)
Compute the net charge rescaled by the smallest non-zero charge.
Definition coulomb.cpp:270
Solver const & get_coulomb()
Definition coulomb.cpp:66
System & get_system()
void operator()(std::shared_ptr< T > const &) const
Definition coulomb.cpp:180
void operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
Definition coulomb.cpp:187
void operator()(std::shared_ptr< CoulombP3M > const &actor) const
Definition coulomb.cpp:183
auto operator()(std::shared_ptr< ReactionField > const &) const
Definition coulomb.cpp:250
LongRangeEnergy(ParticleRange const &particles)
Definition coulomb.cpp:230
auto operator()(std::shared_ptr< CoulombP3M > const &actor) const
Definition coulomb.cpp:234
auto operator()(std::shared_ptr< CoulombMMM1D > const &) const
Definition coulomb.cpp:248
auto operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
Definition coulomb.cpp:238
auto operator()(std::shared_ptr< CoulombScafacos > const &actor) const
Definition coulomb.cpp:243
auto operator()(std::shared_ptr< DebyeHueckel > const &) const
Definition coulomb.cpp:249
void operator()(std::shared_ptr< ReactionField > const &) const
Definition coulomb.cpp:223
LongRangeForce(ParticleRange const &particles)
Definition coulomb.cpp:203
void operator()(std::shared_ptr< CoulombScafacos > const &actor) const
Definition coulomb.cpp:216
void operator()(std::shared_ptr< CoulombMMM1D > const &) const
Definition coulomb.cpp:221
void operator()(std::shared_ptr< CoulombP3M > const &actor) const
Definition coulomb.cpp:207
void operator()(std::shared_ptr< DebyeHueckel > const &) const
Definition coulomb.cpp:222
void operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
Definition coulomb.cpp:211
LongRangePressure(ParticleRange const &particles)
Definition coulomb.cpp:108
auto operator()(std::shared_ptr< ReactionField > const &) const
Definition coulomb.cpp:121
auto operator()(std::shared_ptr< DebyeHueckel > const &) const
Definition coulomb.cpp:117
auto operator()(std::shared_ptr< T > const &) const
Definition coulomb.cpp:127
auto operator()(std::shared_ptr< CoulombP3M > const &actor) const
Definition coulomb.cpp:112
auto operator()(std::shared_ptr< CoulombP3M > const &actor) const
Definition coulomb.cpp:147
auto operator()(std::shared_ptr< CoulombMMM1D > const &) const
Definition coulomb.cpp:156
auto operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
Definition coulomb.cpp:151
auto operator()(std::shared_ptr< ReactionField > const &actor) const
Definition coulomb.cpp:164
auto operator()(std::shared_ptr< DebyeHueckel > const &actor) const
Definition coulomb.cpp:167
auto operator()(std::shared_ptr< CoulombScafacos > const &actor) const
Definition coulomb.cpp:160
void on_observable_calc()
Definition coulomb.cpp:193
double cutoff() const
Definition coulomb.cpp:172
void sanity_checks() const
Definition coulomb.cpp:68
void on_cell_structure_change()
Definition coulomb.cpp:100
void on_periodicity_change()
Definition coulomb.cpp:93
void calc_long_range_force(ParticleRange const &particles) const
Definition coulomb.cpp:256
void on_node_grid_change()
Definition coulomb.cpp:87
Utils::Vector9d calc_pressure_long_range(ParticleRange const &particles) const
Definition coulomb.cpp:138
void on_coulomb_change()
Definition coulomb.cpp:74
std::unique_ptr< Implementation > impl
Pointer-to-implementation.
void on_boxl_change()
Definition coulomb.cpp:81
bool reinit_on_observable_calc
Whether to reinitialize the solver on observable calculation.
double calc_energy_long_range(ParticleRange const &particles) const
Definition coulomb.cpp:262
void visit_try_catch(Visitor &&visitor, Variant &actor)
Run a kernel on a variant and queue errors.