54template <std::
size_t S>
59 auto constexpr exp_limit = 30.;
62 auto const cao =
params.cao;
63 auto const mesh =
params.mesh[0];
71 auto denominator = 0.;
74 m_start, m_stop, indices,
76 auto const norm_sq = nm.norm2();
78 auto const exp_term = f2 * norm_sq;
79 if (exp_term < exp_limit) {
80 auto const f3 = sz * std::exp(-exp_term) / norm_sq;
81 numerator += f3 * Utils::int_pow<S>(d_op * nm);
85 [&](
unsigned dim,
int n) {
86 nm[dim] = shift[dim] + n * mesh;
90 return numerator / (Utils::int_pow<S>(
static_cast<double>(d_op.
norm2())) *
91 Utils::int_pow<2>(denominator));
108template <
typename FloatType, std::
size_t S>
113 auto const size = n_stop - n_start;
116 auto g = std::vector<FloatType>(
Utils::product(size), FloatType(0));
124 auto prefactor = Utils::int_pow<3>(
static_cast<double>(
params.mesh[0])) * 2. *
125 Utils::int_pow<2>(inv_box_l[0]);
127 auto const offset = detail::calc_meshift(
params.mesh,
false)[0];
128 auto const d_op = detail::calc_meshift(
params.mesh,
true)[0];
129 auto const half_mesh =
params.mesh[0] / 2;
133 auto index = std::size_t(0u);
136 n_start, n_stop, indices,
138 if (((indices[0] % half_mesh != 0) or (indices[1] % half_mesh != 0) or
139 (indices[2] % half_mesh != 0))) {
140 g[index] = FloatType(prefactor *
141 G_opt_dipolar<S>(
params, shift_off, d_op_off));
145 [&](
unsigned dim,
int n) {
146 d_op_off[dim] = d_op[n];
147 shift_off[dim] = offset[n];
159 auto const cao =
params.cao;
160 auto const mesh =
params.mesh[0];
168 m_start, m_stop, indices,
170 [&](
unsigned dim,
int n) {
171 fnm[dim] =
math::sinc(offset[dim] + n * mesh);
186template <
typename FloatType>
191 auto const offset = detail::calc_meshift(
params.mesh,
false)[0];
192 auto const d_op = detail::calc_meshift(
params.mesh,
true)[0];
193 auto const half_mesh =
params.mesh[0] / 2;
197 auto index = std::size_t(0u);
201 n_start, n_stop, indices,
203 if (((indices[0] % half_mesh != 0) or (indices[1] % half_mesh != 0) or
204 (indices[2] % half_mesh != 0))) {
206 energy += double(g[index]) * U2 * d_op_off.norm2();
210 [&](
unsigned dim,
int n) {
211 d_op_off[dim] = d_op[n];
212 shift_off[dim] = offset[n];
Vector implementation and trait types for boost qvm interoperability.
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
This file contains the defaults for ESPResSo.
#define P3M_BRILLOUIN
P3M: Number of Brillouin zones taken into account in the calculation of the optimal influence functio...
and std::invocable< Projector, unsigned, int > void for_each_3d(detail::IndexVectorConcept auto &&start, detail::IndexVectorConcept auto &&stop, detail::IndexVectorConcept auto &&counters, Kernel &&kernel, Projector &&projector=detail::noop_projector)
Repeat an operation on every element of a 3D grid.
double G_opt_dipolar(P3MParameters const ¶ms, Utils::Vector3i const &shift, Utils::Vector3i const &d_op)
Calculate the aliasing sums for the optimal influence function.
std::vector< FloatType > grid_influence_function(P3MParameters const ¶ms, Utils::Vector3i const &n_start, Utils::Vector3i const &n_stop, Utils::Vector3d const &inv_box_l)
Map influence function over a grid.
double grid_influence_function_self_energy(P3MParameters const ¶ms, Utils::Vector3i const &n_start, Utils::Vector3i const &n_stop, std::vector< FloatType > const &g)
Calculate self-energy of the influence function.
double G_opt_dipolar_self_energy(P3MParameters const ¶ms, Utils::Vector3i const &shift)
T product(Vector< T, N > const &v)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
Common functions for dipolar and charge P3M.
static SteepestDescentParameters params
Currently active steepest descent instance.
Structure to hold P3M parameters and some dependent variables.