ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
mmm1d.hpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 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/**
23 * @file
24 * MMM1D algorithm for long-range Coulomb interactions on the CPU.
25 *
26 * Implementation of the MMM1D method for the calculation of the electrostatic
27 * interaction in one-dimensionally periodic systems. For details on the
28 * method see MMM in general. The MMM1D method works only with the N-squared
29 * cell system since neither the near nor far formula can be decomposed.
30 *
31 * MMM1D uses polygamma expansions for the near formula.
32 * The expansion of the polygamma functions is fairly easy and follows
33 * directly from @cite abramowitz65a. For details, see @cite arnold02a.
34 */
35
36#pragma once
37
38#include "config/config.hpp"
39
40#ifdef ELECTROSTATICS
41
43
44#include "Particle.hpp"
45
46#include <utils/Vector.hpp>
47
48#include <array>
49
50/** @brief Parameters for the MMM1D electrostatic interaction */
51struct CoulombMMM1D : public Coulomb::Actor<CoulombMMM1D> {
52 /**
53 * @brief Maximal allowed pairwise error for the potential and force.
54 * This error ignores prefactors, i.e. this is for a pure lattice 1/r-sum.
55 */
56 double maxPWerror;
57 /**
58 * @brief Far switch radius. Represents the xy-distance at which
59 * the calculation switches from the far to the near formula.
60 */
64
65 CoulombMMM1D(double prefactor, double maxPWerror, double switch_rad,
66 int tune_timings, bool tune_verbose);
67
68 /** Compute the pair force.
69 * @param[in] q1q2 Product of the charges on p1 and p2.
70 * @param[in] d Vector pointing from p1 to p2.
71 * @param[in] dist Distance between p1 and p2.
72 */
73 Utils::Vector3d pair_force(double q1q2, Utils::Vector3d const &d,
74 double dist) const;
75
76 /** Compute the pair energy.
77 * @param[in] q1q2 Product of the charges on p1 and p2.
78 * @param[in] d Vector pointing from p1 to p2.
79 * @param[in] dist Distance between p1 and p2.
80 */
81 double pair_energy(double q1q2, Utils::Vector3d const &d, double dist) const;
82
83 void tune();
84 bool is_tuned() const { return m_is_tuned; }
85
88 tune();
89 }
90 /** @brief Recalculate all box-length-dependent parameters. */
91 void on_boxl_change() { recalc_boxl_parameters(); }
92 void on_node_grid_change() const {}
93 void on_periodicity_change() const { sanity_checks_periodicity(); }
94 void on_cell_structure_change() { sanity_checks_cell_structure(); }
95 /** @brief Recalculate all derived parameters. */
96 void init() { recalc_boxl_parameters(); }
97
98 void sanity_checks() const {
99 sanity_checks_periodicity();
100 sanity_checks_cell_structure();
102 }
103
104private:
105 bool m_is_tuned;
106 /** @brief Square of the far switch radius. */
107 double far_switch_radius_sq;
108 /** @brief Squared inverse box length in z-direction. */
109 double uz2;
110 /** @brief Squared inverse box length in z-direction times prefactor. */
111 double prefuz2;
112 /** @brief Cubed inverse box length in z-direction times prefactor. */
113 double prefL3_i;
114 /**
115 * @brief Largest numerically stable cutoff for Bessel function.
116 * Don't change without improving the formulas.
117 */
118 static constexpr auto MAXIMAL_B_CUT = 30;
119 /** @brief From which distance a certain Bessel cutoff is valid. */
120 std::array<double, MAXIMAL_B_CUT> bessel_radii;
121 /** @brief Table of Taylor expansions of the modified polygamma functions. */
122 std::vector<std::vector<double>> modPsi;
123
124 /** @brief Create even and odd polygamma functions up to order `2 * new_n`. */
125 void create_mod_psi_up_to(int new_n);
126 void determine_bessel_radii();
127 void prepare_polygamma_series();
128 void recalc_boxl_parameters();
129 void sanity_checks_periodicity() const;
130 void sanity_checks_cell_structure() const;
131};
132
133#endif // ELECTROSTATICS
Vector implementation and trait types for boost qvm interoperability.
double prefactor
Electrostatics prefactor.
This file contains the defaults for ESPResSo.
Parameters for the MMM1D electrostatic interaction.
Definition mmm1d.hpp:51
bool tune_verbose
Definition mmm1d.hpp:63
void on_node_grid_change() const
Definition mmm1d.hpp:92
void sanity_checks() const
Definition mmm1d.hpp:98
double far_switch_radius
Far switch radius.
Definition mmm1d.hpp:61
Utils::Vector3d pair_force(double q1q2, Utils::Vector3d const &d, double dist) const
Compute the pair force.
Definition mmm1d.cpp:182
void on_periodicity_change() const
Definition mmm1d.hpp:93
double pair_energy(double q1q2, Utils::Vector3d const &d, double dist) const
Compute the pair energy.
Definition mmm1d.cpp:273
bool is_tuned() const
Definition mmm1d.hpp:84
void init()
Recalculate all derived parameters.
Definition mmm1d.hpp:96
double maxPWerror
Maximal allowed pairwise error for the potential and force.
Definition mmm1d.hpp:56
int tune_timings
Definition mmm1d.hpp:62
void tune()
Definition mmm1d.cpp:337
void on_cell_structure_change()
Definition mmm1d.hpp:94
void on_boxl_change()
Recalculate all box-length-dependent parameters.
Definition mmm1d.hpp:91
void on_activation()
Definition mmm1d.hpp:86