32#include "communication.hpp"
35#include "system/System.hpp"
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>
62 impl = std::make_unique<Implementation>();
70 std::visit([](
auto const &ptr) { ptr->sanity_checks(); }, *
impl->solver);
89 std::visit([](
auto &ptr) { ptr->on_node_grid_change(); }, *
impl->solver);
109 : m_particles{particles} {}
112 auto operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
113 return actor->long_range_pressure(m_particles);
117 auto operator()(std::shared_ptr<DebyeHueckel>
const &)
const {
121 auto operator()(std::shared_ptr<ReactionField>
const &)
const {
125 template <
typename T,
126 std::enable_if_t<!traits::has_pressure<T>::value> * =
nullptr>
129 <<
"electrostatics method " << Utils::demangle<T>();
147 auto operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
148 return actor->p3m_params.r_cut;
151 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
152 return std::max(actor->elc.space_layer,
153 std::visit(*
this, actor->base_solver));
156 auto operator()(std::shared_ptr<CoulombMMM1D>
const &)
const {
157 return std::numeric_limits<double>::infinity();
160 auto operator()(std::shared_ptr<CoulombScafacos>
const &actor)
const {
161 return actor->get_r_cut();
164 auto operator()(std::shared_ptr<ReactionField>
const &actor)
const {
167 auto operator()(std::shared_ptr<DebyeHueckel>
const &actor)
const {
180 template <
typename T>
void operator()(std::shared_ptr<T>
const &)
const {}
183 void operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
184 actor->count_charged_particles();
187 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
188 std::visit(*
this, actor->base_solver);
204 : m_particles(particles) {}
207 void operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
208 actor->add_long_range_forces(m_particles);
211 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
212 actor->add_long_range_forces(m_particles);
216 void operator()(std::shared_ptr<CoulombScafacos>
const &actor)
const {
217 actor->add_long_range_forces();
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 {}
231 : m_particles(particles) {}
234 auto operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
235 return actor->long_range_energy(m_particles);
238 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
239 return actor->long_range_energy(m_particles);
243 auto operator()(std::shared_ptr<CoulombScafacos>
const &actor)
const {
244 return actor->long_range_energy();
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.; }
271 using namespace boost::accumulators;
272 using KahanSum = accumulator_set<double, features<tag::sum_kahan>>;
275 auto q_min = std::numeric_limits<double>::infinity();
277 for (
auto const q : charges) {
280 q_min = std::min(q_min, std::abs(q));
284 return std::abs(sum_kahan(q_sum)) / q_min;
288 double relative_tolerance) {
290 std::vector<double> local_charges;
292 local_charges.push_back(p.q());
294 std::vector<std::vector<double>> node_charges;
295 boost::mpi::gather(
comm_cart, local_charges, node_charges, 0);
298 auto excess_ratio = 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());
306 boost::mpi::broadcast(
comm_cart, excess_ratio, 0);
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 "
Vector implementation and trait types for boost qvm interoperability.
std::shared_ptr< CellStructure > cell_structure
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.
static auto calc_charge_excess_ratio(std::vector< double > const &charges)
Compute the net charge rescaled by the smallest non-zero charge.
Solver const & get_coulomb()
void operator()(std::shared_ptr< T > const &) const
void operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
void operator()(std::shared_ptr< CoulombP3M > const &actor) const
auto operator()(std::shared_ptr< ReactionField > const &) const
LongRangeEnergy(ParticleRange const &particles)
auto operator()(std::shared_ptr< CoulombP3M > const &actor) const
auto operator()(std::shared_ptr< CoulombMMM1D > const &) const
auto operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
auto operator()(std::shared_ptr< CoulombScafacos > const &actor) const
auto operator()(std::shared_ptr< DebyeHueckel > const &) const
void operator()(std::shared_ptr< ReactionField > const &) const
LongRangeForce(ParticleRange const &particles)
void operator()(std::shared_ptr< CoulombScafacos > const &actor) const
void operator()(std::shared_ptr< CoulombMMM1D > const &) const
void operator()(std::shared_ptr< CoulombP3M > const &actor) const
void operator()(std::shared_ptr< DebyeHueckel > const &) const
void operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
LongRangePressure(ParticleRange const &particles)
auto operator()(std::shared_ptr< ReactionField > const &) const
auto operator()(std::shared_ptr< DebyeHueckel > const &) const
auto operator()(std::shared_ptr< T > const &) const
auto operator()(std::shared_ptr< CoulombP3M > const &actor) const
auto operator()(std::shared_ptr< CoulombP3M > const &actor) const
auto operator()(std::shared_ptr< CoulombMMM1D > const &) const
auto operator()(std::shared_ptr< ElectrostaticLayerCorrection > const &actor) const
auto operator()(std::shared_ptr< ReactionField > const &actor) const
auto operator()(std::shared_ptr< DebyeHueckel > const &actor) const
auto operator()(std::shared_ptr< CoulombScafacos > const &actor) const
void on_observable_calc()
void sanity_checks() const
void on_cell_structure_change()
void on_periodicity_change()
void calc_long_range_force(ParticleRange const &particles) const
void on_node_grid_change()
Utils::Vector9d calc_pressure_long_range(ParticleRange const &particles) const
std::unique_ptr< Implementation > impl
Pointer-to-implementation.
bool reinit_on_observable_calc
Whether to reinitialize the solver on observable calculation.
double calc_energy_long_range(ParticleRange const &particles) const
void visit_try_catch(Visitor &&visitor, Variant &actor)
Run a kernel on a variant and queue errors.