32#include "communication.hpp"
35#include "system/System.hpp"
42#include <boost/accumulators/accumulators.hpp>
43#include <boost/accumulators/statistics/sum_kahan.hpp>
44#include <boost/mpi/collectives/broadcast.hpp>
45#include <boost/mpi/collectives/gather.hpp>
64 impl = std::make_unique<Implementation>();
72 std::visit([](
auto const &ptr) { ptr->sanity_checks(); }, *
impl->solver);
91 std::visit([](
auto &ptr) { ptr->on_node_grid_change(); }, *
impl->solver);
111 : m_particles{particles} {}
114 auto operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
115 return actor->long_range_pressure(m_particles);
119 auto operator()(std::shared_ptr<DebyeHueckel>
const &)
const {
123 auto operator()(std::shared_ptr<ReactionField>
const &)
const {
127 template <
typename T,
128 std::enable_if_t<!traits::has_pressure<T>::value> * =
nullptr>
131 <<
"electrostatics method " << Utils::demangle<T>();
149 auto operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
150 return actor->p3m.params.r_cut;
153 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
154 return std::max(actor->elc.space_layer,
155 std::visit(*
this, actor->base_solver));
159 auto operator()(std::shared_ptr<CoulombMMM1DGpu>
const &)
const {
160 return std::numeric_limits<double>::infinity();
163 auto operator()(std::shared_ptr<CoulombMMM1D>
const &)
const {
164 return std::numeric_limits<double>::infinity();
167 auto operator()(std::shared_ptr<CoulombScafacos>
const &actor)
const {
168 return actor->get_r_cut();
171 auto operator()(std::shared_ptr<ReactionField>
const &actor)
const {
174 auto operator()(std::shared_ptr<DebyeHueckel>
const &actor)
const {
187 template <
typename T>
void operator()(std::shared_ptr<T>
const &)
const {}
190 void operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
191 actor->count_charged_particles();
194 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
195 std::visit(*
this, actor->base_solver);
211 : m_particles(particles) {}
214 void operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
215 actor->add_long_range_forces(m_particles);
218 void operator()(std::shared_ptr<CoulombP3MGPU>
const &actor)
const {
219 actor->add_long_range_forces(m_particles);
223 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
224 actor->add_long_range_forces(m_particles);
228 void operator()(std::shared_ptr<CoulombMMM1DGpu>
const &actor)
const {
229 actor->add_long_range_forces();
233 void operator()(std::shared_ptr<CoulombScafacos>
const &actor)
const {
234 actor->add_long_range_forces();
238 void operator()(std::shared_ptr<CoulombMMM1D>
const &)
const {}
239 void operator()(std::shared_ptr<DebyeHueckel>
const &)
const {}
240 void operator()(std::shared_ptr<ReactionField>
const &)
const {}
248 : m_particles(particles) {}
251 auto operator()(std::shared_ptr<CoulombP3M>
const &actor)
const {
252 return actor->long_range_energy(m_particles);
255 operator()(std::shared_ptr<ElectrostaticLayerCorrection>
const &actor)
const {
256 return actor->long_range_energy(m_particles);
260 auto operator()(std::shared_ptr<CoulombMMM1DGpu>
const &actor)
const {
261 actor->add_long_range_energy();
266 auto operator()(std::shared_ptr<CoulombScafacos>
const &actor)
const {
267 return actor->long_range_energy();
271 auto operator()(std::shared_ptr<CoulombMMM1D>
const &)
const {
return 0.; }
272 auto operator()(std::shared_ptr<DebyeHueckel>
const &)
const {
return 0.; }
273 auto operator()(std::shared_ptr<ReactionField>
const &)
const {
return 0.; }
294 using namespace boost::accumulators;
295 using KahanSum = accumulator_set<double, features<tag::sum_kahan>>;
298 auto q_min = std::numeric_limits<double>::infinity();
300 for (
auto const q : charges) {
303 q_min = std::min(q_min, std::abs(q));
307 return std::abs(sum_kahan(q_sum)) / q_min;
311 double relative_tolerance) {
313 std::vector<double> local_charges;
315 local_charges.push_back(p.q());
317 std::vector<std::vector<double>> node_charges;
318 boost::mpi::gather(
comm_cart, local_charges, node_charges, 0);
321 auto excess_ratio = 0.;
323 auto charges = std::move(local_charges);
324 for (
auto it = ++node_charges.begin(); it != node_charges.end(); ++it) {
325 charges.insert(charges.end(), it->begin(), it->end());
329 boost::mpi::broadcast(
comm_cart, excess_ratio, 0);
331 if (excess_ratio >= relative_tolerance) {
332 std::ostringstream serializer;
333 serializer << std::scientific << std::setprecision(4);
334 serializer << excess_ratio;
335 throw std::runtime_error(
336 "The system is not charge neutral. Please neutralize the system "
337 "before adding a new actor by adding the corresponding counterions "
338 "to the system. Alternatively you can turn off the electroneutrality "
339 "check by supplying check_neutrality=False when creating the actor. "
340 "In this case you may be simulating a non-neutral system which will "
341 "affect physical observables like e.g. the pressure, the chemical "
342 "potentials of charged species or potential energies of the system. "
343 "Since simulations of non charge neutral systems are special please "
344 "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< CoulombMMM1DGpu > const &actor) 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< CoulombP3MGPU > const &actor) 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< CoulombMMM1DGpu > const &actor) 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< CoulombMMM1DGpu > const &) 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.