ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
relative.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 *
4 * This file is part of ESPResSo.
5 *
6 * ESPResSo is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * ESPResSo is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20#include "config/config.hpp"
21#ifdef VIRTUAL_SITES_RELATIVE
22
23#include "BoxGeometry.hpp"
24#include "Particle.hpp"
25#include "PropagationMode.hpp"
27#include "cells.hpp"
28#include "errorhandling.hpp"
29#include "forces.hpp"
31#include "rotation.hpp"
32
33#include <utils/Vector.hpp>
36#include <utils/quaternion.hpp>
37
38#include <cassert>
39
40/**
41 * @brief Vector pointing from the real particle to the virtual site.
42 *
43 * @return Relative distance.
44 */
45static auto connection_vector(Particle const &p_ref, Particle const &p_vs) {
46 // Calculate the quaternion defining the orientation of the vector connecting
47 // the virtual site and the real particle
48 // This is obtained, by multiplying the quaternion representing the director
49 // of the real particle with the quaternion of the virtual particle, which
50 // specifies the relative orientation.
51 auto const director = Utils::convert_quaternion_to_director(
52 p_ref.quat() * p_vs.vs_relative().rel_orientation)
53 .normalize();
54
55 return p_vs.vs_relative().distance * director;
56}
57
58/**
59 * @brief Velocity of the virtual site
60 * @param p_ref Reference particle for the virtual site.
61 * @param p_vs Virtual site.
62 * @return Velocity of the virtual site.
63 */
64static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs) {
65 auto const d = connection_vector(p_ref, p_vs);
66
67 // Get omega of real particle in space-fixed frame
68 auto const omega_space_frame =
69 convert_vector_body_to_space(p_ref, p_ref.omega());
70 // Obtain velocity from v = v_real particle + omega_real_particle * director
71 return vector_product(omega_space_frame, d) + p_ref.v();
72}
73
74/**
75 * @brief Get real particle tracked by a virtual site.
76 *
77 * @param cell_structure Cell structure.
78 * @param p Virtual site.
79 * @return Pointer to real particle.
80 */
82 Particle const &p) {
83 auto const &vs_rel = p.vs_relative();
84 if (vs_rel.to_particle_id == -1) {
85 runtimeErrorMsg() << "Particle with id " << p.id()
86 << " is a dangling virtual site";
87 return nullptr;
88 }
89 auto p_ref_ptr = cell_structure.get_local_particle(vs_rel.to_particle_id);
90 if (!p_ref_ptr) {
91 runtimeErrorMsg() << "No real particle with id " << vs_rel.to_particle_id
92 << " for virtual site with id " << p.id();
93 }
94 return p_ref_ptr;
95}
96
97/**
98 * @brief Constraint force to hold the particle at its prescribed position.
99 *
100 * @param p_ref Reference particle.
101 * @param p_vs Virtual site.
102 * @return Constraint force.
103 */
104static auto constraint_stress(Particle const &p_ref, Particle const &p_vs) {
105 /* The constraint force is minus the force on the particle, make it force
106 * free. The counter force is translated by the connection vector to the
107 * real particle, hence the virial stress is */
108 return tensor_product(-p_vs.force(), connection_vector(p_ref, p_vs));
109}
110
114static bool is_vs_relative_rot(Particle const &p) {
116}
117
118static bool is_vs_relative(Particle const &p) {
120}
121
123 BoxGeometry const &box_geo) {
126
127 auto const particles = cell_structure.local_particles();
128 for (auto &p : particles) {
129 if (!is_vs_relative(p)) {
130 continue;
131 }
132
133 auto const *p_ref_ptr = get_reference_particle(cell_structure, p);
134 if (!p_ref_ptr)
135 continue;
136
137 auto const &p_ref = *p_ref_ptr;
138
139 // position update
140 if (is_vs_relative_trans(p)) {
141 p.image_box() = p_ref.image_box();
142 p.pos() = p_ref.pos() + connection_vector(p_ref, p);
143 p.v() = velocity(p_ref, p);
144
145 if (box_geo.type() == BoxType::LEES_EDWARDS) {
146 auto push = LeesEdwards::Push(box_geo);
147 push(p, -1); // includes a position fold
148 } else {
149 box_geo.fold_position(p.pos(), p.image_box());
150 }
151 }
152
153 // Orientation update
154 if (is_vs_relative_rot(p)) {
155 p.quat() = p_ref.quat() * p.vs_relative().quat;
156 }
157 }
158
159 if (cell_structure.check_resort_required()) {
161 }
162}
163
164// Distribute forces that have accumulated on virtual particles to the
165// associated real particles
167 CellStructure &cell_structure) {
168 cell_structure.ghosts_reduce_forces();
169
170 init_forces_ghosts(cell_structure.ghost_particles());
171
172 // Iterate over all the particles in the local cells
173 for (auto &p : cell_structure.local_particles()) {
174 if (!is_vs_relative(p))
175 continue;
176
177 auto *p_ref_ptr = get_reference_particle(cell_structure, p);
178 assert(p_ref_ptr != nullptr);
179
180 auto &p_ref = *p_ref_ptr;
181 if (is_vs_relative_trans(p)) {
182 p_ref.force() += p.force();
183 p_ref.torque() += vector_product(connection_vector(p_ref, p), p.force());
184 }
185
186 if (is_vs_relative_rot(p)) {
187 p_ref.torque() += p.torque();
188 }
189 }
190}
191
192// Rigid body contribution to scalar pressure and pressure tensor
196
197 for (auto const &p : cell_structure.local_particles()) {
198 if (is_vs_relative_trans(p)) {
199 if (auto const *p_ref_ptr = cell_structure.get_local_particle(
200 p.vs_relative().to_particle_id)) {
201 pressure_tensor += constraint_stress(*p_ref_ptr, p);
202 }
203 }
204 }
205
206 return pressure_tensor;
207}
208#endif // VIRTUAL_SITES_RELATIVE
@ LEES_EDWARDS
Vector implementation and trait types for boost qvm interoperability.
This file contains everything related to the global cell structure / cell system.
BoxType type() const
void fold_position(Utils::Vector3d &pos, Utils::Vector3i &image_box) const
Fold coordinates to primary simulation box in-place.
This file contains the defaults for ESPResSo.
__device__ void vector_product(float const *a, float const *b, float *out)
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
void init_forces_ghosts(ParticleRange const &particles)
Set forces of all ghosts to zero.
Definition forces.cpp:103
Force calculation.
Quaternion algebra.
@ DATA_PART_MOMENTUM
Particle::m.
@ DATA_PART_POSITION
Particle::r.
Vector< T, 3 > convert_quaternion_to_director(Quaternion< T > const &quat)
Convert quaternion to director.
Quaternion implementation and trait types for boost qvm interoperability.
Utils::Matrix< double, 3, 3 > vs_relative_pressure_tensor(CellStructure const &cell_structure)
Definition relative.cpp:194
static bool is_vs_relative_rot(Particle const &p)
Definition relative.cpp:114
void vs_relative_back_transfer_forces_and_torques(CellStructure &cell_structure)
Definition relative.cpp:166
static auto connection_vector(Particle const &p_ref, Particle const &p_vs)
Vector pointing from the real particle to the virtual site.
Definition relative.cpp:45
void vs_relative_update_particles(CellStructure &cell_structure, BoxGeometry const &box_geo)
Definition relative.cpp:122
static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs)
Velocity of the virtual site.
Definition relative.cpp:64
static Particle * get_reference_particle(CellStructure &cell_structure, Particle const &p)
Get real particle tracked by a virtual site.
Definition relative.cpp:81
static bool is_vs_relative_trans(Particle const &p)
Definition relative.cpp:111
static bool is_vs_relative(Particle const &p)
Definition relative.cpp:118
static auto constraint_stress(Particle const &p_ref, Particle const &p_vs)
Constraint force to hold the particle at its prescribed position.
Definition relative.cpp:104
This file contains all subroutines required to process rotational motion.
Utils::Vector3d convert_vector_body_to_space(const Particle &p, const Utils::Vector3d &vec)
Definition rotation.hpp:57
Describes a cell structure / cell system.
ParticleRange ghost_particles() const
Particle * get_local_particle(int id)
Get a local particle by id.
void ghosts_update(unsigned data_parts)
Update ghost particles.
void ghosts_reduce_forces()
Add forces from ghost particles to real 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 set_resort_particles(Cells::Resort level)
Increase the local resort level at least to level.
ParticleRange local_particles() const
Struct holding all information for one particle.
Definition Particle.hpp:395
auto const & propagation() const
Definition Particle.hpp:421
auto const & quat() const
Definition Particle.hpp:477
auto const & vs_relative() const
Definition Particle.hpp:527
auto const & v() const
Definition Particle.hpp:433
auto const & omega() const
Definition Particle.hpp:481
auto const & force() const
Definition Particle.hpp:435
auto const & id() const
Definition Particle.hpp:414
Matrix representation with static size.
Definition matrix.hpp:65