ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
p3m/common.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2024 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22/** \file
23 * Common functions for dipolar and charge P3M.
24 *
25 * We use here a P3M (Particle-Particle Particle-Mesh) method based
26 * on the Ewald summation. Details of the used method can be found in
27 * @cite hockney88a and @cite deserno98a @cite deserno98b. The file p3m
28 * contains only the Particle-Mesh part.
29 *
30 * Further reading: @cite ewald21a, @cite hockney88a, @cite deserno98a,
31 * @cite deserno98b, @cite deserno00e, @cite deserno00b, @cite cerda08d
32 *
33 */
34
35#pragma once
36
37#include "config/config.hpp"
38
39#include <utils/Vector.hpp>
40
41#include <algorithm>
42#include <array>
43#include <vector>
44
45/** This value indicates metallic boundary conditions. */
46auto constexpr P3M_EPSILON_METALLIC = 0.0;
47
48#if defined(P3M) or defined(DP3M)
49
50#include "LocalBox.hpp"
51
52#include <cstddef>
53#include <span>
54#include <stdexcept>
55
56enum class Arch { CPU, GPU };
57
58/** @brief Structure to hold P3M parameters and some dependent variables. */
60 /** tuning or production? */
61 bool tuning;
62 /** Ewald splitting parameter (0<alpha<1), rescaled to
63 * @p alpha_L = @p alpha * @p box_l. */
64 double alpha_L;
65 /** cutoff radius for real space electrostatics (>0), rescaled to
66 * @p r_cut_iL = @p r_cut * @p box_l_i. */
67 double r_cut_iL;
68 /** number of mesh points per coordinate direction (>0). */
70 /** offset of the first mesh point (lower left corner) from the
71 * coordinate origin ([0,1[). */
73 /** charge assignment order ([0,7]). */
74 int cao;
75 /** accuracy of the actual parameter set. */
76 double accuracy;
77
78 /** epsilon of the "surrounding dielectric". */
79 double epsilon;
80 /** cutoff for charge assignment. */
82 /** mesh constant. */
84 /** inverse mesh constant. */
86 /** unscaled @ref P3MParameters::alpha_L "alpha_L" for use with fast
87 * inline functions only */
88 double alpha;
89 /** unscaled @ref P3MParameters::r_cut_iL "r_cut_iL" for use with fast
90 * inline functions only */
91 double r_cut;
92 /** number of points unto which a single charge is interpolated, i.e.
93 * @ref P3MParameters::cao "cao" cubed */
94 int cao3;
95
96 P3MParameters(bool tuning, double epsilon, double r_cut,
98 int cao, double alpha, double accuracy)
99 : tuning{tuning}, alpha_L{0.}, r_cut_iL{0.}, mesh{mesh},
101 cao_cut{}, a{}, ai{}, alpha{alpha}, r_cut{r_cut}, cao3{-1} {
102
103 auto constexpr value_to_tune = -1.;
104
105 if (epsilon < 0.) {
106 throw std::domain_error("Parameter 'epsilon' must be >= 0");
107 }
108
109 if (accuracy <= 0.) {
110 throw std::domain_error("Parameter 'accuracy' must be > 0");
111 }
112
113 if (r_cut <= 0.) {
114 if (tuning and r_cut == value_to_tune) {
115 this->r_cut = 0.;
116 } else {
117 throw std::domain_error("Parameter 'r_cut' must be > 0");
118 }
119 }
120
121 if (alpha <= 0.) {
122 if (tuning and alpha == value_to_tune) {
123 this->alpha = 0.;
124 } else {
125 throw std::domain_error("Parameter 'alpha' must be > 0");
126 }
127 }
128
129 if (not(mesh >= Utils::Vector3i::broadcast(1) or
130 ((mesh[0] >= 1) and (mesh == Utils::Vector3i{{mesh[0], -1, -1}})) or
131 (tuning and mesh == Utils::Vector3i::broadcast(-1)))) {
132 throw std::domain_error("Parameter 'mesh' must be > 0");
133 }
134
135 if (not(mesh_off >= Utils::Vector3d::broadcast(0.) and
138 this->mesh_off = Utils::Vector3d::broadcast(P3M_MESHOFF);
139 } else {
140 throw std::domain_error("Parameter 'mesh_off' must be >= 0 and <= 1");
141 }
142 }
143
144 if ((cao < 1 or cao > 7) and (not tuning or cao != -1)) {
145 throw std::domain_error("Parameter 'cao' must be >= 1 and <= 7");
146 }
147
148 if (not tuning and (Utils::Vector3i::broadcast(cao) > mesh)) {
149 throw std::domain_error("Parameter 'cao' cannot be larger than 'mesh'");
150 }
151 }
152
153 /**
154 * @brief Recalculate quantities derived from the mesh and box length:
155 * @ref P3MParameters::a "a",
156 * @ref P3MParameters::ai "ai" and
157 * @ref P3MParameters::cao_cut "cao_cut".
158 */
162 cao_cut = (static_cast<double>(cao) / 2.) * a;
163 }
164
165 /**
166 * @brief Convert spatial position to grid position.
167 * To get the grid index, round the result to the nearest integer.
168 */
169 auto calc_grid_pos(Utils::Vector3d const &pos) const {
170 return Utils::hadamard_product(pos, ai) - mesh_off;
171 }
172};
173
174/** @brief Properties of the local mesh. */
176 /** dimension (size) of local mesh. */
178 /** number of local mesh points. */
179 int size;
180 /** index of lower left corner of the
181 local mesh in the global mesh. */
182 int ld_ind[3];
183 /** position of the first local mesh point. */
184 double ld_pos[3];
185 /** dimension of mesh inside node domain. */
186 int inner[3];
187 /** inner left down grid point */
188 int in_ld[3];
189 /** inner up right grid point + (1,1,1) */
190 int in_ur[3];
191 /** number of margin mesh points. */
192 int margin[6];
193 /** number of margin mesh points from neighbour nodes */
194 int r_margin[6];
195 /** offset between mesh lines of the last dimension */
197 /** offset between mesh lines of the two last dimensions */
199
200 /**
201 * @brief Recalculate quantities derived from the mesh and box length:
202 * @ref P3MLocalMesh::ld_pos "ld_pos" (position of the left down mesh).
203 */
205 // spatial position of left down mesh point
206 for (auto i = 0u; i < 3u; i++) {
207 ld_pos[i] = (ld_ind[i] + params.mesh_off[i]) * params.a[i];
208 }
209 }
210
211 /**
212 * @brief Calculate properties of the local FFT mesh
213 * for the charge assignment process.
214 */
216 LocalBox const &local_geo, double skin,
217 double space_layer);
218};
219
220/** @brief Local mesh FFT buffers. */
221template <typename FloatType> struct P3MFFTMesh {
222 /** @brief real-space scalar mesh for charge assignment and FFT. */
223 std::span<FloatType> rs_scalar;
224 /** @brief real-space vector meshes for the electric or dipolar field. */
225 std::array<std::span<FloatType>, 3> rs_fields;
226
227 /** @brief Indices of the lower left corner of the local mesh grid. */
229 /** @brief Indices of the upper right corner of the local mesh grid. */
231 /** @brief Extents of the local mesh grid. */
233
234 /** @brief number of permutations in k_space */
235 int ks_pnum = 0;
236};
237
238#endif // defined(P3M) or defined(DP3M)
239
240namespace detail {
241/** Calculate indices that shift @ref P3MParameters::mesh "mesh" by `mesh/2`.
242 * For each mesh size @f$ n @f$ in @c mesh_size, create a sequence of integer
243 * values @f$ \left( 0, \ldots, \lfloor n/2 \rfloor, -\lfloor n/2 \rfloor,
244 * \ldots, -1\right) @f$ if @c zero_out_midpoint is false, otherwise
245 * @f$ \left( 0, \ldots, \lfloor n/2 - 1 \rfloor, 0, -\lfloor n/2 \rfloor,
246 * \ldots, -1\right) @f$.
247 */
248std::array<std::vector<int>, 3> inline calc_meshift(
249 Utils::Vector3i const &mesh_size, bool zero_out_midpoint = false) {
250 std::array<std::vector<int>, 3> ret{};
251
252 for (auto i = 0u; i < 3u; ++i) {
253 ret[i] = std::vector<int>(static_cast<std::size_t>(mesh_size[i]));
254
255 for (int j = 1; j <= mesh_size[i] / 2; j++) {
256 ret[i][j] = j;
257 ret[i][mesh_size[i] - j] = -j;
258 }
259 if (zero_out_midpoint)
260 ret[i][mesh_size[i] / 2] = 0;
261 }
262
263 return ret;
264}
265} // namespace detail
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.
Definition Vector.hpp:110
This file contains the defaults for ESPResSo.
#define P3M_MESHOFF
P3M: Default for offset of first mesh point from the origin (left down corner of the simulation box).
Definition config.hpp:46
auto hadamard_division(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:422
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:379
Arch
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
static SteepestDescentParameters params
Currently active steepest descent instance.
Local mesh FFT buffers.
std::array< std::span< FloatType >, 3 > rs_fields
real-space vector meshes for the electric or dipolar field.
int ks_pnum
number of permutations in k_space
Utils::Vector3i start
Indices of the lower left corner of the local mesh grid.
std::span< FloatType > rs_scalar
real-space scalar mesh for charge assignment and FFT.
Utils::Vector3i stop
Indices of the upper right corner of the local mesh grid.
Utils::Vector3i size
Extents of the local mesh grid.
Properties of the local mesh.
Utils::Vector3i dim
dimension (size) of local mesh.
int in_ur[3]
inner up right grid point + (1,1,1)
int size
number of local mesh points.
void recalc_ld_pos(P3MParameters const &params)
Recalculate quantities derived from the mesh and box length: ld_pos (position of the left down mesh).
void calc_local_ca_mesh(P3MParameters const &params, LocalBox const &local_geo, double skin, double space_layer)
Calculate properties of the local FFT mesh for the charge assignment process.
Definition common.cpp:34
int in_ld[3]
inner left down grid point
int r_margin[6]
number of margin mesh points from neighbour nodes
int margin[6]
number of margin mesh points.
int q_2_off
offset between mesh lines of the last dimension
double ld_pos[3]
position of the first local mesh point.
int inner[3]
dimension of mesh inside node domain.
int ld_ind[3]
index of lower left corner of the local mesh in the global mesh.
int q_21_off
offset between mesh lines of the two last dimensions
Structure to hold P3M parameters and some dependent variables.
auto calc_grid_pos(Utils::Vector3d const &pos) const
Convert spatial position to grid position.
Utils::Vector3d cao_cut
cutoff for charge assignment.
double alpha
unscaled alpha_L for use with fast inline functions only
P3MParameters(bool tuning, double epsilon, double r_cut, Utils::Vector3i const &mesh, Utils::Vector3d const &mesh_off, int cao, double alpha, double accuracy)
double r_cut_iL
cutoff radius for real space electrostatics (>0), rescaled to r_cut_iL = r_cut * box_l_i.
int cao
charge assignment order ([0,7]).
double accuracy
accuracy of the actual parameter set.
double alpha_L
Ewald splitting parameter (0.
int cao3
number of points unto which a single charge is interpolated, i.e.
Utils::Vector3d mesh_off
offset of the first mesh point (lower left corner) from the coordinate origin ([0,...
Utils::Vector3d ai
inverse mesh constant.
double r_cut
unscaled r_cut_iL for use with fast inline functions only
void recalc_a_ai_cao_cut(Utils::Vector3d const &box_l)
Recalculate quantities derived from the mesh and box length: a, ai and cao_cut.
bool tuning
tuning or production?
Utils::Vector3i mesh
number of mesh points per coordinate direction (>0).
double epsilon
epsilon of the "surrounding dielectric".
Utils::Vector3d a
mesh constant.