ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
rattle.cpp
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#include "rattle.hpp"
23
24#ifdef BOND_CONSTRAINT
25
26#include "BoxGeometry.hpp"
27#include "Particle.hpp"
28#include "ParticleRange.hpp"
32#include "communication.hpp"
33#include "errorhandling.hpp"
34
35#include <boost/mpi/collectives/all_reduce.hpp>
36#include <boost/range/algorithm.hpp>
37
38#include <cmath>
39#include <functional>
40#include <span>
41
42/**
43 * @brief copy current position
44 *
45 * @param particles particle range
46 * @param ghost_particles ghost particle range
47 */
48void save_old_position(const ParticleRange &particles,
49 const ParticleRange &ghost_particles) {
50 auto save_pos = [](Particle &p) { p.pos_last_time_step() = p.pos(); };
51
52 boost::for_each(particles, save_pos);
53 boost::for_each(ghost_particles, save_pos);
54}
55
56/**
57 * @brief reset correction vectors to zero
58 *
59 * @param particles particle range
60 * @param ghost_particles ghost particle range
61 */
62static void init_correction_vector(const ParticleRange &particles,
63 const ParticleRange &ghost_particles) {
64 auto reset_force = [](Particle &p) { p.rattle_params().correction.fill(0); };
65
66 boost::for_each(particles, reset_force);
67 boost::for_each(ghost_particles, reset_force);
68}
69
70/**
71 * @brief Calculate the positional correction for the particles.
72 *
73 * @param ia_params Parameters
74 * @param box_geo Box geometry.
75 * @param p1 First particle.
76 * @param p2 Second particle.
77 * @return True if there was a correction.
78 */
79static bool calculate_positional_correction(RigidBond const &ia_params,
80 BoxGeometry const &box_geo,
81 Particle &p1, Particle &p2) {
82 auto const r_ij = box_geo.get_mi_vector(p1.pos(), p2.pos());
83 auto const r_ij2 = r_ij.norm2();
84
85 if (std::abs(1.0 - r_ij2 / ia_params.d2) > ia_params.p_tol) {
86 auto const r_ij_t =
88 auto const r_ij_dot = r_ij_t * r_ij;
89 auto const G =
90 0.50 * (ia_params.d2 - r_ij2) / r_ij_dot / (p1.mass() + p2.mass());
91
92 auto const pos_corr = G * r_ij_t;
93 p1.rattle_params().correction += pos_corr * p2.mass();
94 p2.rattle_params().correction -= pos_corr * p1.mass();
95
96 return true;
97 }
98
99 return false;
100}
101
102/**
103 * @brief Compute the correction vectors using given kernel.
104 *
105 * @param cs cell structure
106 * @param box_geo Box geometry
107 * @param bonded_ias Bonded interactions
108 * @param kernel kernel function
109 * @return True if correction is necessary
110 */
111template <typename Kernel>
113 BoxGeometry const &box_geo,
114 BondedInteractionsMap const &bonded_ias,
115 Kernel kernel) {
116 bool correction = false;
117 cs.bond_loop([&correction, &kernel, &box_geo, &bonded_ias](
118 Particle &p1, int bond_id, std::span<Particle *> partners) {
119 auto const &iaparams = *bonded_ias.at(bond_id);
120
121 if (auto const *bond = boost::get<RigidBond>(&iaparams)) {
122 auto const corrected = kernel(*bond, box_geo, p1, *partners[0]);
123 if (corrected)
124 correction = true;
125 }
126
127 /* Rigid bonds cannot break */
128 return false;
129 });
130
131 return correction;
132}
133
134/**
135 * @brief Apply positional corrections
136 *
137 * @param particles particle range
138 */
139static void apply_positional_correction(const ParticleRange &particles) {
140 boost::for_each(particles, [](Particle &p) {
141 p.pos() += p.rattle_params().correction;
142 p.v() += p.rattle_params().correction;
143 });
144}
145
147 BondedInteractionsMap const &bonded_ias) {
150
151 auto particles = cs.local_particles();
152 auto ghost_particles = cs.ghost_particles();
153
154 int cnt;
155 for (cnt = 0; cnt < SHAKE_MAX_ITERATIONS; ++cnt) {
156 init_correction_vector(particles, ghost_particles);
157 bool const repeat_ = compute_correction_vector(
158 cs, box_geo, bonded_ias, calculate_positional_correction);
159 bool const repeat =
160 boost::mpi::all_reduce(comm_cart, repeat_, std::logical_or<bool>());
161
162 // no correction is necessary, skip communication and bail out
163 if (!repeat)
164 break;
165
167
170 }
171 if (cnt >= SHAKE_MAX_ITERATIONS) {
172 runtimeErrorMsg() << "RATTLE failed to converge after " << cnt
173 << " iterations";
174 }
175
176 auto const resort_level =
178 cs.set_resort_particles(resort_level);
179}
180
181/**
182 * @brief Calculate the velocity correction for the particles.
183 *
184 * The position correction is accumulated in the forces
185 * of the particles so that it can be reduced over the ghosts.
186 *
187 * @param ia_params Parameters
188 * @param box_geo Box geometry.
189 * @param p1 First particle.
190 * @param p2 Second particle.
191 * @return True if there was a correction.
192 */
193static bool calculate_velocity_correction(RigidBond const &ia_params,
194 BoxGeometry const &box_geo,
195 Particle &p1, Particle &p2) {
196 auto const v_ij = p1.v() - p2.v();
197 auto const r_ij = box_geo.get_mi_vector(p1.pos(), p2.pos());
198
199 auto const v_proj = v_ij * r_ij;
200 if (std::abs(v_proj) > ia_params.v_tol) {
201 auto const K = v_proj / ia_params.d2 / (p1.mass() + p2.mass());
202
203 auto const vel_corr = K * r_ij;
204
205 p1.rattle_params().correction -= vel_corr * p2.mass();
206 p2.rattle_params().correction += vel_corr * p1.mass();
207
208 return true;
209 }
210
211 return false;
212}
213
214/**
215 * @brief Apply velocity corrections
216 *
217 * @param particles particle range
218 */
219static void apply_velocity_correction(ParticleRange const &particles) {
220 boost::for_each(particles,
221 [](Particle &p) { p.v() += p.rattle_params().correction; });
222}
223
225 BondedInteractionsMap const &bonded_ias) {
227
228 auto particles = cs.local_particles();
229 auto ghost_particles = cs.ghost_particles();
230
231 int cnt;
232 for (cnt = 0; cnt < SHAKE_MAX_ITERATIONS; ++cnt) {
233 init_correction_vector(particles, ghost_particles);
234 bool const repeat_ = compute_correction_vector(
235 cs, box_geo, bonded_ias, calculate_velocity_correction);
236 bool const repeat =
237 boost::mpi::all_reduce(comm_cart, repeat_, std::logical_or<bool>());
238
239 // no correction is necessary, skip communication and bail out
240 if (!repeat)
241 break;
242
244
245 apply_velocity_correction(particles);
247 }
248
249 if (cnt >= SHAKE_MAX_ITERATIONS) {
250 runtimeErrorMsg() << "VEL RATTLE failed to converge after " << cnt
251 << " iterations";
252 }
253}
254
255#endif
Data structures for bonded interactions.
container for bonded interactions.
mapped_type at(key_type const &key) const
Utils::Vector< T, 3 > get_mi_vector(const Utils::Vector< T, 3 > &a, const Utils::Vector< T, 3 > &b) const
Get the minimum-image vector between two coordinates.
A range of particles.
boost::mpi::communicator comm_cart
The communicator.
#define SHAKE_MAX_ITERATIONS
Maximal number of iterations in the RATTLE algorithm before it bails out.
Definition config.hpp:88
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
@ DATA_PART_MOMENTUM
Particle::m.
@ DATA_PART_PROPERTIES
Particle::p.
@ DATA_PART_POSITION
Particle::r.
void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Correction of current velocities using RATTLE algorithm.
Definition rattle.cpp:224
void save_old_position(const ParticleRange &particles, const ParticleRange &ghost_particles)
copy current position
Definition rattle.cpp:48
void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias)
Propagate velocity and position while using SHAKE algorithm for bond constraint.
Definition rattle.cpp:146
static void init_correction_vector(const ParticleRange &particles, const ParticleRange &ghost_particles)
reset correction vectors to zero
Definition rattle.cpp:62
static void apply_positional_correction(const ParticleRange &particles)
Apply positional corrections.
Definition rattle.cpp:139
static bool calculate_velocity_correction(RigidBond const &ia_params, BoxGeometry const &box_geo, Particle &p1, Particle &p2)
Calculate the velocity correction for the particles.
Definition rattle.cpp:193
static void apply_velocity_correction(ParticleRange const &particles)
Apply velocity corrections.
Definition rattle.cpp:219
static bool compute_correction_vector(CellStructure &cs, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias, Kernel kernel)
Compute the correction vectors using given kernel.
Definition rattle.cpp:112
static bool calculate_positional_correction(RigidBond const &ia_params, BoxGeometry const &box_geo, Particle &p1, Particle &p2)
Calculate the positional correction for the particles.
Definition rattle.cpp:79
RATTLE algorithm ().
Definition of the rigid bond data type for the Rattle algorithm.
Describes a cell structure / cell system.
ParticleRange ghost_particles() const
void update_ghosts_and_resort_particle(unsigned data_parts)
Update ghost particles, with particle resort if needed.
void ghosts_update(unsigned data_parts)
Update ghost particles.
bool check_resort_required(Utils::Vector3d const &additional_offset={}) const
Check whether a particle has moved further than half the skin since the last Verlet list update,...
void bond_loop(BondKernel const &bond_kernel)
Bonded pair loop.
void set_resort_particles(Cells::Resort level)
Increase the local resort level at least to level.
ParticleRange local_particles() const
void ghosts_reduce_rattle_correction()
Add rattle corrections from ghost particles to real particles.
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & rattle_params() const
Definition Particle.hpp:567
auto const & mass() const
Definition Particle.hpp:452
auto const & pos_last_time_step() const
Definition Particle.hpp:565
auto const & v() const
Definition Particle.hpp:433
auto const & pos() const
Definition Particle.hpp:431
Parameters for the rigid_bond/SHAKE/RATTLE ALGORITHM.
double d2
Square of the length of Constrained Bond.
double v_tol
Velocity Tolerance/Accuracy for termination of RATTLE/SHAKE iterations during velocity corrections.
double p_tol
Positional Tolerance/Accuracy value for termination of RATTLE/SHAKE iterations during position correc...