ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
p3m.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 * P3M algorithm for long-range Coulomb interaction.
24 *
25 * We use a P3M (Particle-Particle Particle-Mesh) method based on the
26 * Ewald summation. Details of the used method can be found in
27 * @cite hockney88a and @cite deserno98a @cite deserno98b.
28 *
29 * Further reading: @cite ewald21a, @cite hockney88a, @cite deserno98a,
30 * @cite deserno98b, @cite deserno00e, @cite deserno00b, @cite cerda08d.
31 *
32 * Implementation in p3m.cpp.
33 */
34
35#pragma once
36
37#include "config/config.hpp"
38
39#ifdef P3M
40
42
43#include "p3m/common.hpp"
44#include "p3m/data_struct.hpp"
45
46#include "ParticleRange.hpp"
47
48#include <utils/Vector.hpp>
50
51#include <cmath>
52#include <numbers>
53
54/** @brief P3M solver. */
55struct CoulombP3M : public Coulomb::Actor<CoulombP3M> {
57
58public:
60
61 virtual ~CoulombP3M() = default;
62
63 [[nodiscard]] virtual bool is_tuned() const noexcept = 0;
64 [[nodiscard]] virtual bool is_gpu() const noexcept = 0;
65 [[nodiscard]] virtual bool is_double_precision() const noexcept = 0;
66
67 /** @brief Recalculate all derived parameters. */
68 virtual void init() = 0;
69 virtual void on_activation() = 0;
70
71 /** @brief Recalculate all box-length-dependent parameters. */
86
87 /**
88 * Count the number of charged particles and calculate
89 * the sum of the squared charges.
90 */
91 virtual void count_charged_particles() = 0;
92 virtual void count_charged_particles_elc(int, double, double) = 0;
93 virtual void adapt_epsilon_elc() = 0;
94
95 /**
96 * @brief Tune P3M parameters to desired accuracy.
97 *
98 * The parameters
99 * @ref P3MParameters::mesh "mesh",
100 * @ref P3MParameters::cao "cao",
101 * @ref P3MParameters::r_cut_iL "r_cut_iL" and
102 * @ref P3MParameters::alpha_L "alpha_L" are tuned to obtain the target
103 * @ref P3MParameters::accuracy "accuracy" in optimal time.
104 * These parameters are stored in the @ref p3m_params object.
105 *
106 * The function utilizes the analytic expression of the error estimate
107 * for the P3M method in @cite hockney88a (eq. (8.23)) in
108 * order to obtain the rms error in the force for a system of N randomly
109 * distributed particles in a cubic box.
110 * For the real space error the estimate of Kolafa/Perram is used.
111 *
112 * Parameter ranges if not given explicitly in the constructor:
113 * - @p mesh explores the range 0-512 (biased toward values less than 128)
114 * - @p cao explores all possible values
115 * - @p alpha_L is tuned for each tuple (@p r_cut_iL, @p mesh, @p cao) and
116 * calculated assuming that the error contributions of real and reciprocal
117 * space should be equal
118 *
119 * After checking if the total error lies below the target accuracy,
120 * the time needed for one force calculation is measured. Parameters
121 * that minimize the runtime are kept.
122 *
123 * The function is based on routines of the program HE_Q.cpp written by M.
124 * Deserno.
125 */
126 virtual void tune() = 0;
127
128 /** Assign the physical charges using the tabulated charge assignment
129 * function.
130 */
131 virtual void charge_assign(ParticleRange const &particles) = 0;
132
133 /**
134 * @brief Assign a single charge into the current charge grid.
135 *
136 * @param[in] q Particle charge
137 * @param[in] real_pos Particle position in real space
138 * @param[in] skip_cache Skip interpolation weights caching.
139 */
140 virtual void assign_charge(double q, Utils::Vector3d const &real_pos,
141 bool skip_cache) = 0;
142
143 virtual void prepare_fft_mesh(bool reset_weights) = 0;
144
145 /** Calculate real-space contribution of p3m Coulomb pair forces. */
147 double dist) const {
148 if ((q1q2 == 0.) || dist >= p3m_params.r_cut || dist <= 0.) {
149 return {};
150 }
151 auto const alpha = p3m_params.alpha;
152 auto const adist = alpha * dist;
153 auto const exp_adist_sq = exp(-adist * adist);
154 auto const dist_sq = dist * dist;
155 auto const two_a_sqrt_pi_i = 2. * alpha * std::numbers::inv_sqrtpi;
156#if USE_ERFC_APPROXIMATION
157 auto const erfc_part_ri = Utils::AS_erfc_part(adist) / dist;
158 auto const fac = exp_adist_sq * (erfc_part_ri + two_a_sqrt_pi_i) / dist_sq;
159#else
160 auto const erfc_part_ri = erfc(adist) / dist;
161 auto const fac = (erfc_part_ri + two_a_sqrt_pi_i * exp_adist_sq) / dist_sq;
162#endif
163 return (fac * prefactor * q1q2) * d;
164 }
165
166 /** Calculate real-space contribution of Coulomb pair energy. */
167 // Eq. (3.6) @cite deserno00b
168 double pair_energy(double q1q2, double dist) const {
169 if ((q1q2 == 0.) || dist >= p3m_params.r_cut || dist <= 0.) {
170 return {};
171 }
172 auto const adist = p3m_params.alpha * dist;
173#if USE_ERFC_APPROXIMATION
174 auto const erfc_part_ri = Utils::AS_erfc_part(adist) / dist;
175 return prefactor * q1q2 * erfc_part_ri * exp(-adist * adist);
176#else
177 auto const erfc_part_ri = erfc(adist) / dist;
178 return prefactor * q1q2 * erfc_part_ri;
179#endif
180 }
181
182 /** Compute the k-space part of the pressure tensor */
184
185 /** Compute the k-space part of energies. */
186 virtual double long_range_energy(ParticleRange const &) = 0;
187
188 /** Compute the k-space part of forces. */
189 virtual void add_long_range_forces(ParticleRange const &) = 0;
190
191protected:
194
195 /** Checks for correctness of the k-space cutoff. */
196 void sanity_checks_boxl() const;
197 void sanity_checks_node_grid() const;
198 void sanity_checks_periodicity() const;
199 void sanity_checks_cell_structure() const;
200
201 virtual void scaleby_box_l() = 0;
202};
203
204#endif // P3M
Vector implementation and trait types for boost qvm interoperability.
double prefactor
Electrostatics prefactor.
A range of particles.
This file contains the defaults for ESPResSo.
constexpr T AS_erfc_part(T d)
Approximate by applying a formula from chapter 7.
Common functions for dipolar and charge P3M.
P3M solver.
Definition p3m.hpp:55
virtual ~CoulombP3M()=default
virtual bool is_gpu() const noexcept=0
void sanity_checks_periodicity() const
Definition p3m.cpp:824
void on_cell_structure_change()
Definition p3m.hpp:75
virtual void tune()=0
Tune P3M parameters to desired accuracy.
virtual void add_long_range_forces(ParticleRange const &)=0
Compute the k-space part of forces.
virtual void count_charged_particles_elc(int, double, double)=0
CoulombP3M(P3MParameters const &p3m_params)
Definition p3m.hpp:59
Utils::Vector3d pair_force(double q1q2, Utils::Vector3d const &d, double dist) const
Calculate real-space contribution of p3m Coulomb pair forces.
Definition p3m.hpp:146
void on_boxl_change()
Recalculate all box-length-dependent parameters.
Definition p3m.hpp:72
double pair_energy(double q1q2, double dist) const
Calculate real-space contribution of Coulomb pair energy.
Definition p3m.hpp:168
void sanity_checks_boxl() const
Checks for correctness of the k-space cutoff.
Definition p3m.cpp:793
virtual void prepare_fft_mesh(bool reset_weights)=0
void sanity_checks_cell_structure() const
Definition p3m.cpp:832
P3MParameters const & p3m_params
Definition p3m.hpp:56
virtual Utils::Vector9d long_range_pressure(ParticleRange const &)=0
Compute the k-space part of the pressure tensor.
virtual void on_activation()=0
void sanity_checks() const
Definition p3m.hpp:79
void on_periodicity_change() const
Definition p3m.hpp:74
virtual void adapt_epsilon_elc()=0
virtual void calc_influence_function_force()=0
void on_node_grid_change() const
Definition p3m.hpp:73
virtual bool is_double_precision() const noexcept=0
virtual bool is_tuned() const noexcept=0
virtual void count_charged_particles()=0
Count the number of charged particles and calculate the sum of the squared charges.
virtual double long_range_energy(ParticleRange const &)=0
Compute the k-space part of energies.
virtual void charge_assign(ParticleRange const &particles)=0
Assign the physical charges using the tabulated charge assignment function.
void sanity_checks_node_grid() const
Definition p3m.cpp:847
virtual void assign_charge(double q, Utils::Vector3d const &real_pos, bool skip_cache)=0
Assign a single charge into the current charge grid.
virtual void init()=0
Recalculate all derived parameters.
virtual void scaleby_box_l()=0
virtual void calc_influence_function_energy()=0
Structure to hold P3M parameters and some dependent variables.
double alpha
unscaled alpha_L for use with fast inline functions only
double r_cut
unscaled r_cut_iL for use with fast inline functions only