ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
forces.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/** \file
22 * Force calculation.
23 *
24 * The corresponding header file is forces.hpp.
25 */
26
27#include "BoxGeometry.hpp"
28#include "Particle.hpp"
29#include "ParticleRange.hpp"
30#include "PropagationMode.hpp"
33#include "cells.hpp"
34#include "collision_detection/CollisionDetection.hpp"
35#include "communication.hpp"
36#include "constraints/Constraints.hpp"
38#include "forces_inline.hpp"
39#include "galilei/ComFixed.hpp"
46#include "npt.hpp"
47#include "rotation.hpp"
48#include "short_range_loop.hpp"
49#include "system/System.hpp"
50#include "thermostat.hpp"
52
53#include <utils/Vector.hpp>
54#include <utils/math/sqr.hpp>
55
56#include <boost/variant.hpp>
57
58#ifdef CALIPER
59#include <caliper/cali.h>
60#endif
61
62#include <cassert>
63#include <cmath>
64#include <memory>
65#include <span>
66#include <variant>
67
68/** External particle forces */
70 ParticleForce f = {};
71
72#ifdef EXTERNAL_FORCES
73 f.f += p.ext_force();
74#ifdef ROTATION
75 f.torque += p.ext_torque();
76#endif
77#endif
78
79#ifdef ENGINE
80 // apply a swimming force in the direction of
81 // the particle's orientation axis
82 if (p.swimming().swimming and !p.swimming().is_engine_force_on_fluid) {
83 f.f += p.swimming().f_swim * p.calc_director();
84 }
85#endif
86
87 return f;
88}
89
90static void init_forces(ParticleRange const &particles,
91 ParticleRange const &ghost_particles) {
92#ifdef CALIPER
93 CALI_CXX_MARK_FUNCTION;
94#endif
95
96 for (auto &p : particles) {
97 p.force_and_torque() = external_force(p);
98 }
99
100 init_forces_ghosts(ghost_particles);
101}
102
103void init_forces_ghosts(ParticleRange const &particles) {
104 for (auto &p : particles) {
105 p.force_and_torque() = {};
106 }
107}
108
109static void force_capping(ParticleRange const &particles, double force_cap) {
110 if (force_cap > 0.) {
111 auto const force_cap_sq = Utils::sqr(force_cap);
112 for (auto &p : particles) {
113 auto const force_sq = p.force().norm2();
114 if (force_sq > force_cap_sq) {
115 p.force() *= force_cap / std::sqrt(force_sq);
116 }
117 }
118 }
119}
120
122#ifdef CALIPER
123 CALI_CXX_MARK_FUNCTION;
124#endif
125#ifdef CUDA
126#ifdef CALIPER
127 CALI_MARK_BEGIN("copy_particles_to_GPU");
128#endif
129 gpu.update();
130#ifdef CALIPER
131 CALI_MARK_END("copy_particles_to_GPU");
132#endif
133#endif // CUDA
134
135#ifdef COLLISION_DETECTION
136 collision_detection->clear_queue();
137#endif
138 bond_breakage->clear_queue();
139 auto particles = cell_structure->local_particles();
140 auto ghost_particles = cell_structure->ghost_particles();
141#ifdef ELECTROSTATICS
142 if (coulomb.impl->extension) {
143 if (auto icc = std::get_if<std::shared_ptr<ICCStar>>(
144 get_ptr(coulomb.impl->extension))) {
145 (**icc).iteration(*cell_structure, particles, ghost_particles);
146 }
147 }
148#endif // ELECTROSTATICS
149#ifdef NPT
151#endif
152 init_forces(particles, ghost_particles);
153 thermostat_force_init();
154
155 calc_long_range_forces(particles);
156
157 auto const elc_kernel = coulomb.pair_force_elc_kernel();
158 auto const coulomb_kernel = coulomb.pair_force_kernel();
159 auto const dipoles_kernel = dipoles.pair_force_kernel();
160
161#ifdef ELECTROSTATICS
162 auto const coulomb_cutoff = coulomb.cutoff();
163#else
164 auto const coulomb_cutoff = INACTIVE_CUTOFF;
165#endif
166
167#ifdef DIPOLES
168 auto const dipole_cutoff = dipoles.cutoff();
169#else
170 auto const dipole_cutoff = INACTIVE_CUTOFF;
171#endif
172#ifdef COLLISION_DETECTION
173 auto const collision_detection_cutoff = collision_detection->cutoff();
174#else
175 auto const collision_detection_cutoff = INACTIVE_CUTOFF;
176#endif
177
179 [coulomb_kernel_ptr = get_ptr(coulomb_kernel), &bonded_ias = *bonded_ias,
180 &bond_breakage = *bond_breakage, &box_geo = *box_geo](
181 Particle &p1, int bond_id, std::span<Particle *> partners) {
182 return add_bonded_force(p1, bond_id, partners, bonded_ias,
183 bond_breakage, box_geo, coulomb_kernel_ptr);
184 },
185 [coulomb_kernel_ptr = get_ptr(coulomb_kernel),
186 dipoles_kernel_ptr = get_ptr(dipoles_kernel),
187 elc_kernel_ptr = get_ptr(elc_kernel), &nonbonded_ias = *nonbonded_ias,
188 &thermostat = *thermostat, &bonded_ias = *bonded_ias,
189#ifdef COLLISION_DETECTION
190 &collision_detection = *collision_detection,
191#endif
192 &box_geo = *box_geo](Particle &p1, Particle &p2, Distance const &d) {
193 auto const &ia_params =
194 nonbonded_ias.get_ia_param(p1.type(), p2.type());
195 add_non_bonded_pair_force(p1, p2, d.vec21, sqrt(d.dist2), d.dist2,
196 ia_params, thermostat, box_geo, bonded_ias,
197 coulomb_kernel_ptr, dipoles_kernel_ptr,
198 elc_kernel_ptr);
199#ifdef COLLISION_DETECTION
200 if (not collision_detection.is_off()) {
201 collision_detection.detect_collision(p1, p2, d.dist2);
202 }
203#endif
204 },
205 *cell_structure, maximal_cutoff(), bonded_ias->maximal_cutoff(),
206 VerletCriterion<>{*this, cell_structure->get_verlet_skin(),
207 get_interaction_range(), coulomb_cutoff, dipole_cutoff,
208 collision_detection_cutoff});
209
210 constraints->add_forces(particles, get_sim_time());
211 oif_global->calculate_forces();
212
213 // Must be done here. Forces need to be ghost-communicated
214 immersed_boundaries->volume_conservation(*cell_structure);
215
216 if (thermostat->lb and (propagation->used_propagations &
218 lb_couple_particles();
219 }
220
221#ifdef CUDA
222#ifdef CALIPER
223 CALI_MARK_BEGIN("copy_forces_from_GPU");
224#endif
225 gpu.copy_forces_to_host(particles, this_node);
226#ifdef CALIPER
227 CALI_MARK_END("copy_forces_from_GPU");
228#endif
229#endif // CUDA
230
231#ifdef VIRTUAL_SITES_RELATIVE
232 if (propagation->used_propagations &
235 }
236#endif
237
238 // Communication step: ghost forces
239 cell_structure->ghosts_reduce_forces();
240
241 // should be pretty late, since it needs to zero out the total force
242 comfixed->apply(particles);
243
244 // Needs to be the last one to be effective
245 force_capping(particles, force_cap);
246
247 // mark that forces are now up-to-date
248 propagation->recalc_forces = false;
249}
250
251void calc_long_range_forces(const ParticleRange &particles) {
252#ifdef CALIPER
253 CALI_CXX_MARK_FUNCTION;
254#endif
255
256#ifdef ELECTROSTATICS
257 /* calculate k-space part of electrostatic interaction. */
259#endif // ELECTROSTATICS
260
261#ifdef DIPOLES
262 /* calculate k-space part of the magnetostatic interaction. */
264#endif // DIPOLES
265}
266
267#ifdef NPT
272#endif
Vector implementation and trait types for boost qvm interoperability.
This file contains everything related to the global cell structure / cell system.
A range of particles.
void calculate_forces()
Calculate all forces.
Definition forces.cpp:121
Returns true if the particles are to be considered for short range interactions.
int this_node
The number of this node.
const T * get_ptr(std::optional< T > const &opt)
static void init_forces(ParticleRange const &particles, ParticleRange const &ghost_particles)
Definition forces.cpp:90
static void force_capping(ParticleRange const &particles, double force_cap)
Definition forces.cpp:109
void init_forces_ghosts(ParticleRange const &particles)
Set forces of all ghosts to zero.
Definition forces.cpp:103
static ParticleForce external_force(Particle const &p)
External particle forces.
Definition forces.cpp:69
void npt_add_virial_force_contribution(const Utils::Vector3d &force, const Utils::Vector3d &d)
Update the NpT virial.
Definition forces.cpp:268
void calc_long_range_forces(const ParticleRange &particles)
Calculate long range forces (P3M, ...).
Definition forces.cpp:251
Force calculation.
void add_non_bonded_pair_force(Particle &p1, Particle &p2, Utils::Vector3d const &d, double dist, double dist2, IA_parameters const &ia_params, Thermostat::Thermostat const &thermostat, BoxGeometry const &box_geo, BondedInteractionsMap const &bonded_ias, Coulomb::ShortRangeForceKernel::kernel_type const *coulomb_kernel, Dipoles::ShortRangeForceKernel::kernel_type const *dipoles_kernel, Coulomb::ShortRangeForceCorrectionsKernel::kernel_type const *elc_kernel)
Calculate non-bonded forces between a pair of particles and update their forces and torques.
bool add_bonded_force(Particle &p1, int bond_id, std::span< Particle * > partners, BondedInteractionsMap const &bonded_ia_params, BondBreakage::BondBreakage &bond_breakage, BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel)
ICC is a method that allows to take into account the influence of arbitrarily shaped dielectric inter...
Solver const & get_coulomb()
Definition coulomb.cpp:66
Solver const & get_dipoles()
Definition dipoles.cpp:49
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Various procedures concerning interactions between particles.
constexpr double INACTIVE_CUTOFF
Cutoff for deactivated interactions.
void npt_add_virial_contribution(double energy)
Definition npt.cpp:141
void npt_reset_instantaneous_virials()
reset virial part of instantaneous pressure
Definition npt.cpp:134
Exports for the NpT code.
void vs_relative_back_transfer_forces_and_torques(CellStructure &cell_structure)
Definition relative.cpp:166
This file contains all subroutines required to process rotational motion.
void short_range_loop(BondKernel bond_kernel, PairKernel pair_kernel, CellStructure &cell_structure, double pair_cutoff, double bond_cutoff, VerletCriterion const &verlet_criterion={})
void calc_long_range_force(ParticleRange const &particles) const
Definition coulomb.cpp:256
void calc_long_range_force(ParticleRange const &particles) const
Definition dipoles.cpp:199
Distance vector and length handed to pair kernels.
Force information on a particle.
Definition Particle.hpp:290
Utils::Vector3d torque
torque.
Definition Particle.hpp:320
Utils::Vector3d f
force.
Definition Particle.hpp:316
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & swimming() const
Definition Particle.hpp:561
auto const & ext_force() const
Definition Particle.hpp:554
auto const & type() const
Definition Particle.hpp:418
auto const & ext_torque() const
Definition Particle.hpp:484
auto calc_director() const
Definition Particle.hpp:487