36#include "communication.hpp"
38#include "system/System.hpp"
43#include <boost/mpi/collectives/all_reduce.hpp>
44#include <boost/range/combine.hpp>
106template <std::
size_t dir>
108 std::size_t n_freq,
double u) {
109 auto constexpr c_2pi = 2. * std::numbers::pi;
110 auto const n_part = particles.
size();
111 std::vector<SCCache> ret(n_freq * n_part);
113 for (std::size_t freq = 1; freq <= n_freq; freq++) {
114 auto const pref = c_2pi * u *
static_cast<double>(freq);
116 std::size_t o = (freq - 1) * n_part;
117 for (
auto const &p : particles) {
118 auto const arg = pref * p.pos()[dir];
119 ret[o++] = {sin(arg), cos(arg)};
126static std::pair<std::size_t, std::size_t>
129 assert(far_cut >= 0.);
130 auto const n_freq_x =
131 static_cast<std::size_t
>(std::ceil(far_cut * box_geo.
length()[0]) + 1.);
132 auto const n_freq_y =
133 static_cast<std::size_t
>(std::ceil(far_cut * box_geo.
length()[1]) + 1.);
136 scxcache = calc_sc_cache<0>(particles, n_freq_x, u_x);
137 scycache = calc_sc_cache<1>(particles, n_freq_y, u_y);
138 return {n_freq_x, n_freq_y};
146 for (std::size_t i = 0; i < size; i++)
150static void copy_vec(
double *pdc_d,
double const *pdc_s, std::size_t size) {
151 for (std::size_t i = 0; i < size; i++)
155static void add_vec(
double *pdc_d,
double const *pdc_s1,
double const *pdc_s2,
157 for (std::size_t i = 0; i < size; i++)
158 pdc_d[i] = pdc_s1[i] + pdc_s2[i];
161static void addscale_vec(
double *pdc_d,
double scale,
double const *pdc_s1,
162 double const *pdc_s2, std::size_t size) {
163 for (std::size_t i = 0; i < size; i++)
164 pdc_d[i] = scale * pdc_s1[i] + pdc_s2[i];
167static void scale_vec(
double scale,
double *pdc, std::size_t size) {
168 for (std::size_t i = 0; i < size; i++)
172static double *
block(
double *p, std::size_t index, std::size_t size) {
173 return &p[index * size];
180 boost::mpi::all_reduce(
comm_cart, send_buf,
static_cast<int>(size),
gblcblk,
184void ElectrostaticLayerCorrection::check_gap(
Particle const &p)
const {
186 auto const z = p.
pos()[2];
189 <<
"region by " << ((z < 0.) ? z : z -
elc.box_h);
201void ElectrostaticLayerCorrection::add_dipole_force(
203 constexpr std::size_t size = 3;
205 auto const pref =
prefactor * 4. * std::numbers::pi / box_geo.volume();
209 auto const shift = box_geo.length_half()[2];
219 auto const q = p.
q();
220 auto const z = p.
pos()[2];
249 auto const field_induced =
gblcblk[1];
251 field_tot -= field_applied + field_induced;
255 p.
force()[2] -= field_tot * p.
q();
267double ElectrostaticLayerCorrection::dipole_energy(
269 constexpr std::size_t size = 7;
271 auto const pref =
prefactor * 2. * std::numbers::pi / box_geo.volume();
272 auto const lz = box_geo.length()[2];
275 auto const shift = box_geo.length_half()[2];
289 auto const q = p.
q();
290 auto const z = p.
pos()[2];
319 energy += 2. * pref *
355 double b(
double q,
double z)
const {
360 double t(
double q,
double z)
const {
366ElectrostaticLayerCorrection::z_energy(
ParticleRange const &particles)
const {
367 constexpr std::size_t size = 4;
369 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
370 auto const pref =
prefactor * 2. * std::numbers::pi * xy_area_inv;
374 auto const fac_delta = delta / (1. - delta);
378 auto const shift = box_geo.length_half()[2];
379 auto const lz = box_geo.length()[2];
385 auto const z = p.
pos()[2];
386 auto const q = p.
q();
404 auto const z = p.
pos()[2];
405 auto const q = p.
q();
438 return (
this_node == 0) ? -pref * energy : 0.;
441void ElectrostaticLayerCorrection::add_z_force(
443 constexpr std::size_t size = 1;
445 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
446 auto const pref =
prefactor * 2. * std::numbers::pi * xy_area_inv;
450 auto const fac_delta = delta / (1. - delta);
457 auto const z = p.
pos()[2];
458 auto const q = p.
q();
467 auto const z = p.
pos()[2];
468 auto const q = p.
q();
505 constexpr std::size_t size = 4;
507 auto const pref_di = prefactor * 4. * std::numbers::pi * xy_area_inv;
508 auto const pref = -pref_di / expm1(omega * box_geo.
length()[2]);
509 double lclimgebot[4], lclimgetop[4], lclimge[4];
510 double fac_delta_mid_bot = 1., fac_delta_mid_top = 1., fac_delta = 1.;
514 auto const fac_elc = 1. / (1. - delta * exp(-omega * 2. * elc.
box_h));
525 auto const o = (index - 1) * particles.
size();
526 for (
auto const &p : particles) {
527 auto const z = p.
pos()[2];
528 auto const q = p.
q();
529 auto e = exp(omega * z);
546 lclimgebot[
POQESM] = sc_cache[o + ic].s / e;
547 lclimgebot[
POQESP] = sc_cache[o + ic].s * e;
548 lclimgebot[
POQECM] = sc_cache[o + ic].c / e;
549 lclimgebot[
POQECP] = sc_cache[o + ic].c * e;
554 exp(omega * (+z - 2. * elc.
box_h))) *
557 e = (exp(-omega * z) +
562 lclimge[
POQESP] += q * sc_cache[o + ic].s * e;
563 lclimge[
POQECP] += q * sc_cache[o + ic].c * e;
566 e = exp(omega * (2. * elc.
box_h - z));
570 lclimgetop[
POQESM] = sc_cache[o + ic].s / e;
571 lclimgetop[
POQESP] = sc_cache[o + ic].s * e;
572 lclimgetop[
POQECM] = sc_cache[o + ic].c / e;
573 lclimgetop[
POQECP] = sc_cache[o + ic].c * e;
578 exp(omega * (-z - 2. * elc.
box_h))) *
581 e = (exp(omega * (+z - 2. * elc.
box_h)) +
586 lclimge[
POQESM] += q * sc_cache[o + ic].s * e;
587 lclimge[
POQECM] += q * sc_cache[o + ic].c * e;
602 constexpr auto i =
static_cast<int>(axis);
603 constexpr std::size_t size = 4;
606 for (
auto &p : particles) {
607 auto &force = p.
force();
621 constexpr std::size_t size = 4;
624 for (std::size_t ic = 0; ic < n_part; ic++) {
631 return energy / omega;
642 std::size_t index_q,
double omega,
645 assert(index_p >= 1);
646 assert(index_q >= 1);
647 constexpr std::size_t size = 8;
649 auto const pref_di = prefactor * 8. * std::numbers::pi * xy_area_inv;
650 auto const pref = -pref_di / expm1(omega * box_geo.
length()[2]);
651 double lclimgebot[8], lclimgetop[8], lclimge[8];
652 double fac_delta_mid_bot = 1., fac_delta_mid_top = 1., fac_delta = 1.;
655 auto const fac_elc = 1. / (1. - delta * exp(-omega * 2. * elc.
box_h));
665 auto const ox = (index_p - 1) * particles.
size();
666 auto const oy = (index_q - 1) * particles.
size();
667 for (
auto const &p : particles) {
668 auto const z = p.
pos()[2];
669 auto const q = p.
q();
670 auto e = exp(omega * z);
712 exp(omega * (+z - 2. * elc.
box_h))) *
717 e = (exp(-omega * z) +
719 fac_delta_mid_bot * q;
729 e = exp(omega * (2. * elc.
box_h - z));
745 exp(omega * (-z - 2. * elc.
box_h))) *
750 e = (exp(omega * (+z - 2. * elc.
box_h)) +
752 fac_delta_mid_top * q;
771static void add_PQ_force(std::size_t index_p, std::size_t index_q,
double omega,
774 auto constexpr c_2pi = 2. * std::numbers::pi;
776 c_2pi * box_geo.
length_inv()[0] *
static_cast<double>(index_p) / omega;
778 c_2pi * box_geo.
length_inv()[1] *
static_cast<double>(index_q) / omega;
779 constexpr std::size_t size = 8;
782 for (
auto &p : particles) {
783 auto &force = p.
force();
812static double PQ_energy(
double omega, std::size_t n_part) {
813 constexpr std::size_t size = 8;
816 for (std::size_t ic = 0; ic < n_part; ic++) {
826 return energy / omega;
830void ElectrostaticLayerCorrection::add_force(
832 auto constexpr c_2pi = 2. * std::numbers::pi;
835 auto const n_scxcache = std::get<0>(n_freqs);
836 auto const n_scycache = std::get<1>(n_freqs);
839 add_dipole_force(particles);
840 add_z_force(particles);
843 for (std::size_t p = 1;
844 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
847 auto const omega = c_2pi * box_geo.length_inv()[0] *
static_cast<double>(p);
848 setup_PoQ<PoQ::P>(
elc,
prefactor, p, omega, particles, box_geo);
850 add_PoQ_force<PoQ::P>(particles);
853 for (std::size_t q = 1;
854 box_geo.length_inv()[1] *
static_cast<double>(q - 1) <
elc.
far_cut &&
857 auto const omega = c_2pi * box_geo.length_inv()[1] *
static_cast<double>(q);
858 setup_PoQ<PoQ::Q>(
elc,
prefactor, q, omega, particles, box_geo);
860 add_PoQ_force<PoQ::Q>(particles);
863 for (std::size_t p = 1;
864 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
867 for (std::size_t q = 1;
868 Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p - 1)) +
870 static_cast<double>(q - 1)) <
876 sqrt(
Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p)) +
877 Utils::sqr(box_geo.length_inv()[1] *
static_cast<double>(q)));
885double ElectrostaticLayerCorrection::calc_energy(
887 auto constexpr c_2pi = 2. * std::numbers::pi;
889 auto energy = dipole_energy(particles) + z_energy(particles);
891 auto const n_scxcache = std::get<0>(n_freqs);
892 auto const n_scycache = std::get<1>(n_freqs);
895 partblk.resize(n_localpart * 8);
898 for (std::size_t p = 1;
899 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
902 auto const omega = c_2pi * box_geo.length_inv()[0] *
static_cast<double>(p);
903 setup_PoQ<PoQ::P>(
elc,
prefactor, p, omega, particles, box_geo);
908 for (std::size_t q = 1;
909 box_geo.length_inv()[1] *
static_cast<double>(q - 1) <
elc.
far_cut &&
912 auto const omega = c_2pi * box_geo.length_inv()[1] *
static_cast<double>(q);
913 setup_PoQ<PoQ::Q>(
elc,
prefactor, q, omega, particles, box_geo);
918 for (std::size_t p = 1;
919 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
922 for (std::size_t q = 1;
923 Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p - 1)) +
925 static_cast<double>(q - 1)) <
931 sqrt(
Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p)) +
932 Utils::sqr(box_geo.length_inv()[1] *
static_cast<double>(q)));
942double ElectrostaticLayerCorrection::tune_far_cut()
const {
944 auto constexpr maximal_far_cut = 50.;
946 auto const box_l_x_inv = box_geo.length_inv()[0];
947 auto const box_l_y_inv = box_geo.length_inv()[1];
948 auto const min_inv_boxl = std::min(box_l_x_inv, box_l_y_inv);
949 auto const box_l_z = box_geo.length()[2];
954 auto tuned_far_cut = min_inv_boxl;
957 auto const pref = 2. * std::numbers::pi * tuned_far_cut;
958 auto const sum = pref + 2. * (box_l_x_inv + box_l_y_inv);
959 auto const den = -expm1(-pref * lz);
960 auto const num1 = exp(pref * (
elc.
box_h - lz));
961 auto const num2 = exp(-pref * (
elc.
box_h + lz));
967 tuned_far_cut += min_inv_boxl;
968 }
while (err >
elc.
maxPWerror and tuned_far_cut < maximal_far_cut);
969 if (tuned_far_cut >= maximal_far_cut) {
970 throw std::runtime_error(
"ELC tuning failed: maxPWerror too small");
972 return tuned_far_cut - min_inv_boxl;
980 return boost::mpi::all_reduce(
comm_cart, local_q, std::plus<>());
983void ElectrostaticLayerCorrection::sanity_checks_periodicity()
const {
985 if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) {
986 throw std::runtime_error(
"ELC: requires periodicity (True, True, True)");
990void ElectrostaticLayerCorrection::sanity_checks_dielectric_contrasts()
const {
992 auto const &cell_structure = *
get_system().cell_structure;
995 if (total_charge >= precision_threshold) {
999 throw std::runtime_error(
"ELC does not currently support non-neutral "
1000 "systems with a dielectric contrast.");
1004 throw std::runtime_error(
"ELC does not work for non-neutral systems and "
1005 "non-metallic dielectric contrast.");
1010void ElectrostaticLayerCorrection::adapt_solver() {
1012 [
this](
auto &solver) {
1014 solver->adapt_epsilon_elc();
1020void ElectrostaticLayerCorrection::recalc_box_h() {
1024 if (new_box_h < 0.) {
1025 throw std::runtime_error(
"ELC gap size (" + std::to_string(
elc.
gap_size) +
1026 ") larger than box length in z-direction (" +
1027 std::to_string(box_z) +
")");
1032void ElectrostaticLayerCorrection::recalc_space_layer() {
1034 auto const p3m_r_cut = std::visit(
1035 [](
auto &solver) {
return solver->p3m_params.r_cut; },
base_solver);
1043 auto const half_box_h =
elc.
box_h / 2.;
1044 auto const max_space_layer = std::min(free_space, half_box_h);
1046 if (max_space_layer <= 0.) {
1047 throw std::runtime_error(
"P3M real-space cutoff too large for ELC w/ "
1048 "dielectric contrast");
1057 bool neutralize,
double delta_top,
double delta_bot,
1058 bool with_const_pot,
double potential_diff)
1059 : maxPWerror{maxPWerror}, gap_size{gap_size}, box_h{-1.}, far_cut{far_cut},
1060 far_cut2{-1.}, far_calculated{far_cut == -1.},
1061 dielectric_contrast_on{delta_top != 0. or delta_bot != 0.},
1062 const_pot{with_const_pot and dielectric_contrast_on},
1063 neutralize{neutralize and !dielectric_contrast_on},
1064 delta_mid_top{std::clamp(delta_top, -1., +1.)},
1065 delta_mid_bot{std::clamp(delta_bot, -1., +1.)},
1066 pot_diff{(with_const_pot) ? potential_diff : 0.},
1069 space_layer{(dielectric_contrast_on) ? gap_size / 3. : 0.},
1070 space_box{gap_size - ((dielectric_contrast_on) ? 2. * space_layer : 0.)} {
1074 throw std::domain_error(
"Parameter 'far_cut' must be > 0");
1077 throw std::domain_error(
"Parameter 'maxPWerror' must be > 0");
1080 throw std::domain_error(
"Parameter 'gap_size' must be > 0");
1082 if (potential_diff != 0. and not with_const_pot) {
1083 throw std::invalid_argument(
1084 "Parameter 'const_pot' must be True when 'pot_diff' is non-zero");
1086 if (delta_top < -delta_range or delta_top > delta_range) {
1087 throw std::domain_error(
1088 "Parameter 'delta_mid_top' must be >= -1 and <= +1");
1090 if (delta_bot < -delta_range or delta_bot > delta_range) {
1091 throw std::domain_error(
1092 "Parameter 'delta_mid_bot' must be >= -1 and <= +1");
1099 throw std::domain_error(
"ELC with two parallel metallic boundaries "
1100 "requires the const_pot option");
1106 : elc{parameters}, base_solver{solver} {
1110template <ChargeProtocol protocol,
typename combined_ranges>
1112 combined_ranges
const &p_q_pos_range) {
1117 for (
auto zipped : p_q_pos_range) {
1118 auto const p_q = boost::get<0>(zipped);
1119 auto const &p_pos = boost::get<1>(zipped);
1131 solver.
assign_charge(q_eff, {p_pos[0], p_pos[1], -p_pos[2]},
true);
1136 q_eff, {p_pos[0], p_pos[1], 2. * elc.
box_h - p_pos[2]},
true);
1143template <ChargeProtocol protocol,
typename combined_range>
1145 combined_range
const &p_q_pos_range) {
1148 auto local_q2 = 0.0;
1150 for (
auto zipped : p_q_pos_range) {
1151 auto const p_q = boost::get<0>(zipped);
1152 auto const &p_pos = boost::get<1>(zipped);
1154 auto const p_z = p_pos[2];
1181 auto global_q2 = 0.;
1183 boost::mpi::all_reduce(
comm_cart, local_n, global_n, std::plus<>());
1184 boost::mpi::all_reduce(
comm_cart, local_q2, global_q2, std::plus<>());
1185 boost::mpi::all_reduce(
comm_cart, local_q, global_q, std::plus<>());
1191 auto const energy = std::visit(
1192 [
this, &particles](
auto const &solver_ptr) {
1193 auto &solver = *solver_ptr;
1198 auto p_q_pos_range = boost::combine(p_q_range, p_pos_range);
1201 solver.charge_assign(particles);
1204 return solver.long_range_energy(particles);
1208 energy += 0.5 * solver.long_range_energy(particles);
1213 charge_assign<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1214 modify_p3m_sums<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1215 energy += 0.5 * solver.long_range_energy(particles);
1218 charge_assign<ChargeProtocol::IMAGE>(
elc, solver, p_q_pos_range);
1219 modify_p3m_sums<ChargeProtocol::IMAGE>(
elc, solver, p_q_pos_range);
1220 energy -= 0.5 * solver.long_range_energy(particles);
1223 modify_p3m_sums<ChargeProtocol::REAL>(
elc, solver, p_q_pos_range);
1228 return energy + calc_energy(particles);
1234 [
this, &particles](
auto const &solver_ptr) {
1235 auto &solver = *solver_ptr;
1238 auto p_q_pos_range = boost::combine(p_q_range, p_pos_range);
1241 modify_p3m_sums<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1242 charge_assign<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1245 solver.charge_assign(particles);
1247 solver.add_long_range_forces(particles);
1249 modify_p3m_sums<ChargeProtocol::REAL>(
elc, solver, p_q_pos_range);
1253 add_force(particles);
Vector implementation and trait types for boost qvm interoperability.
Utils::Vector3d const & length() const
Box length.
Utils::Vector3d const & length_inv() const
Inverse box length.
void set_prefactor(double new_prefactor)
double prefactor
Electrostatics prefactor.
base_type::size_type size() const
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
T image_sum(InputIterator begin, InputIterator end, InputIterator it, bool with_replicas, Utils::Vector3i const &ncut, BoxGeometry const &box_geo, T init, F f)
Sum over all pairs with periodic images.
static void addscale_vec(double *pdc_d, double scale, double const *pdc_s1, double const *pdc_s2, std::size_t size)
static std::pair< std::size_t, std::size_t > prepare_sc_cache(ParticleRange const &particles, BoxGeometry const &box_geo, double far_cut)
static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
static auto calc_total_charge(CellStructure const &cell_structure)
static std::vector< double > partblk
temporary buffers for product decomposition
static void clear_vec(double *pdc, std::size_t size)
static double * block(double *p, std::size_t index, std::size_t size)
void setup_PoQ(elc_data const &elc, double prefactor, std::size_t index, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
static void distribute(std::size_t size)
void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver, combined_range const &p_q_pos_range)
static double PoQ_energy(double omega, std::size_t n_part)
static std::vector< SCCache > scxcache
Cached sin/cos values along the x-axis and y-axis.
static double PQ_energy(double omega, std::size_t n_part)
void charge_assign(elc_data const &elc, CoulombP3M &solver, combined_ranges const &p_q_pos_range)
static std::vector< SCCache > scycache
static void setup_PQ(elc_data const &elc, double prefactor, std::size_t index_p, std::size_t index_q, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
static void copy_vec(double *pdc_d, double const *pdc_s, std::size_t size)
static void add_vec(double *pdc_d, double const *pdc_s1, double const *pdc_s2, std::size_t size)
static void scale_vec(double scale, double *pdc, std::size_t size)
void add_PoQ_force(ParticleRange const &particles)
ChargeProtocol
ELC charge sum/assign protocol: real charges, image charges, or both.
static std::vector< SCCache > calc_sc_cache(ParticleRange const &particles, std::size_t n_freq, double u)
Calculate cached sin/cos values for one direction.
static double gblcblk[8]
collected data from the other cells
PoQ
ELC axes (x and y directions)
ELC algorithm for long-range Coulomb interactions.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
ParticleRange particles(std::span< Cell *const > cells)
auto charge_range(ParticleRange const &particles)
auto pos_range(ParticleRange const &particles)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Vector< T, N > sqrt(Vector< T, N > const &a)
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
P3M algorithm for long-range Coulomb interaction.
Describes a cell structure / cell system.
ParticleRange local_particles() const
virtual void count_charged_particles_elc(int, double, double)=0
virtual void prepare_fft_mesh(bool reset_weights)=0
virtual void assign_charge(double q, Utils::Vector3d const &real_pos, bool skip_cache)=0
Assign a single charge into the current charge grid.
std::variant< std::shared_ptr< CoulombP3M > > BaseSolver
BaseSolver base_solver
Electrostatics solver that is adapted.
ElectrostaticLayerCorrection(elc_data &¶meters, BaseSolver &&solver)
void add_long_range_forces(ParticleRange const &particles) const
double long_range_energy(ParticleRange const &particles) const
double b(double q, double z) const
Image sum from the bottom layer.
ImageSum(double delta, double shift, double lz)
double t(double q, double z) const
Image sum from the top layer.
Struct holding all information for one particle.
auto const & force() const
structure for caching sin and cos values
Parameters for the ELC method.
double dielectric_layers_self_energy(CoulombP3M const &p3m, BoxGeometry const &box_geo, ParticleRange const &particles) const
self energies of top and bottom layers with their virtual images
double maxPWerror
Maximal allowed pairwise error for the potential and force.
double pot_diff
Constant potential difference.
double box_h
Up to where particles can be found.
bool dielectric_contrast_on
Flag whether there is any dielectric contrast in the system.
elc_data(double maxPWerror, double gap_size, double far_cut, bool neutralize, double delta_top, double delta_bot, bool const_pot, double pot_diff)
double space_box
The space that is finally left.
bool neutralize
Flag whether the box is neutralized by a homogeneous background.
double far_cut
Cutoff of the exponential sum.
double space_layer
Layer around the dielectric contrast in which we trick around.
bool far_calculated
Flag whether far_cut was set by the user, or calculated by ESPResSo.
double gap_size
Size of the empty gap.
double delta_mid_bot
dielectric contrast in the lower part of the simulation cell.
void dielectric_layers_self_forces(CoulombP3M const &p3m, BoxGeometry const &box_geo, ParticleRange const &particles) const
forces of particles in border layers with themselves
bool const_pot
Flag whether a constant potential difference is applied.
double far_cut2
Squared value of far_cut.
double delta_mid_top
dielectric contrast in the upper part of the simulation cell.