ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
influence_function_dipolar.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002-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#pragma once
23
24#include "config/config.hpp"
25
26#if defined(DP3M)
27
28#include "p3m/common.hpp"
29#include "p3m/for_each_3d.hpp"
30#include "p3m/math.hpp"
31
32#include <utils/Vector.hpp>
34#include <utils/math/sqr.hpp>
35
36#include <cmath>
37#include <cstddef>
38#include <functional>
39#include <numbers>
40#include <vector>
41
42/** Calculate the aliasing sums for the optimal influence function.
43 *
44 * Calculates the aliasing sums in the numerator and denominator of
45 * the expression for the optimal influence function (see
46 * @cite hockney88a : 8-22, p. 275).
47 *
48 * \tparam S order (2 for energy, 3 for forces)
49 * \param params DP3M parameters
50 * \param shift shift for a given n-vector
51 * \param d_op differential operator for a given n-vector
52 * \return The result of the fraction.
53 */
54template <std::size_t S>
56 Utils::Vector3i const &d_op) {
57
58 auto constexpr limit = P3M_BRILLOUIN;
59 auto constexpr exp_limit = 30.;
60 auto constexpr m_start = Utils::Vector3i::broadcast(-limit);
61 auto constexpr m_stop = Utils::Vector3i::broadcast(limit + 1);
62 auto const cao = params.cao;
63 auto const mesh = params.mesh[0];
64 auto const offset =
65 static_cast<Utils::Vector3d>(shift) / static_cast<double>(mesh);
66 auto const f2 = Utils::sqr(std::numbers::pi / params.alpha_L);
67 auto indices = Utils::Vector3i{};
68 auto nm = Utils::Vector3i{};
69 auto fnm = Utils::Vector3d{};
70 auto numerator = 0.;
71 auto denominator = 0.;
72
74 m_start, m_stop, indices,
75 [&]() {
76 auto const norm_sq = nm.norm2();
77 auto const sz = std::pow(Utils::product(fnm), 2 * cao);
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);
82 }
83 denominator += sz;
84 },
85 [&](unsigned dim, int n) {
86 nm[dim] = shift[dim] + n * mesh;
87 fnm[dim] = math::sinc(offset[dim] + n * mesh);
88 });
89
90 return numerator / (Utils::int_pow<S>(static_cast<double>(d_op.norm2())) *
91 Utils::int_pow<2>(denominator));
92}
93
94/**
95 * @brief Map influence function over a grid.
96 *
97 * This evaluates the optimal influence function @ref G_opt_dipolar
98 * over a regular grid of k vectors, and returns the values as a vector.
99 *
100 * @tparam S Order of the differential operator, e.g. 2 for energy, 3 for force
101 *
102 * @param params DP3M parameters
103 * @param n_start Lower left corner of the grid
104 * @param n_stop Upper right corner of the grid.
105 * @param inv_box_l Inverse box length
106 * @return Values of the influence function at regular grid points.
107 */
108template <typename FloatType, std::size_t S>
109std::vector<FloatType> grid_influence_function(
110 P3MParameters const &params, Utils::Vector3i const &n_start,
111 Utils::Vector3i const &n_stop, Utils::Vector3d const &inv_box_l) {
112
113 auto const size = n_stop - n_start;
114
115 /* The influence function grid */
116 auto g = std::vector<FloatType>(Utils::product(size), FloatType(0));
117
118 /* Skip influence function calculation in tuning mode,
119 the results need not be correct for timing. */
120 if (params.tuning) {
121 return g;
122 }
123
124 auto prefactor = Utils::int_pow<3>(static_cast<double>(params.mesh[0])) * 2. *
125 Utils::int_pow<2>(inv_box_l[0]);
126
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;
130 auto indices = Utils::Vector3i{};
131 auto shift_off = Utils::Vector3i{};
132 auto d_op_off = Utils::Vector3i{};
133 auto index = std::size_t(0u);
134
136 n_start, n_stop, indices,
137 [&]() {
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));
142 }
143 ++index;
144 },
145 [&](unsigned dim, int n) {
146 d_op_off[dim] = d_op[n];
147 shift_off[dim] = offset[n];
148 });
149
150 return g;
151}
152
154 Utils::Vector3i const &shift) {
155
156 auto constexpr limit = P3M_BRILLOUIN + 1;
157 auto constexpr m_start = Utils::Vector3i::broadcast(-limit);
158 auto constexpr m_stop = Utils::Vector3i::broadcast(limit + 1);
159 auto const cao = params.cao;
160 auto const mesh = params.mesh[0];
161 auto const offset =
162 static_cast<Utils::Vector3d>(shift) / static_cast<double>(mesh);
163 auto indices = Utils::Vector3i{};
164 auto fnm = Utils::Vector3d{};
165 auto energy = 0.;
166
168 m_start, m_stop, indices,
169 [&]() { energy += std::pow(Utils::product(fnm), 2 * cao); },
170 [&](unsigned dim, int n) {
171 fnm[dim] = math::sinc(offset[dim] + n * mesh);
172 });
173
174 return energy;
175}
176
177/**
178 * @brief Calculate self-energy of the influence function.
179 *
180 * @param params DP3M parameters
181 * @param n_start Lower left corner of the grid
182 * @param n_stop Upper right corner of the grid.
183 * @param g Energies on the grid.
184 * @return Total self-energy.
185 */
186template <typename FloatType>
188 P3MParameters const &params, Utils::Vector3i const &n_start,
189 Utils::Vector3i const &n_stop, std::vector<FloatType> const &g) {
190
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;
194 auto indices = Utils::Vector3i{};
195 auto shift_off = Utils::Vector3i{};
196 auto d_op_off = Utils::Vector3i{};
197 auto index = std::size_t(0u);
198 auto energy = 0.;
199
201 n_start, n_stop, indices,
202 [&]() {
203 if (((indices[0] % half_mesh != 0) or (indices[1] % half_mesh != 0) or
204 (indices[2] % half_mesh != 0))) {
205 auto const U2 = G_opt_dipolar_self_energy(params, shift_off);
206 energy += double(g[index]) * U2 * d_op_off.norm2();
207 }
208 ++index;
209 },
210 [&](unsigned dim, int n) {
211 d_op_off[dim] = d_op[n];
212 shift_off[dim] = offset[n];
213 });
214
215 return energy;
216}
217
218#endif // defined(DP3M)
Vector implementation and trait types for boost qvm interoperability.
T norm2() const
Definition Vector.hpp:137
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_BRILLOUIN
P3M: Number of Brillouin zones taken into account in the calculation of the optimal influence functio...
Definition config.hpp:53
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 &params, 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 &params, 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 &params, 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 &params, Utils::Vector3i const &shift)
T product(Vector< T, N > const &v)
Definition Vector.hpp:374
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
Definition math.hpp:53
Common functions for dipolar and charge P3M.
static SteepestDescentParameters params
Currently active steepest descent instance.
Structure to hold P3M parameters and some dependent variables.