ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
BoxGeometry.hpp
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#pragma once
21
24
25#include <utils/Vector.hpp>
26#include <utils/math/sgn.hpp>
27
28#include <bitset>
29#include <cassert>
30#include <cmath>
31#include <limits>
32#include <stdexcept>
33#include <utility>
34
35namespace detail {
36/**
37 * @brief Get the minimum-image distance between two coordinates.
38 * @param a Coordinate of the terminal point.
39 * @param b Coordinate of the initial point.
40 * @param box_length Box length.
41 * @param box_length_inv Inverse box length
42 * @param box_length_half Half box length
43 * @param periodic Box periodicity.
44 * @return Shortest distance from @p b to @p a across periodic images,
45 * i.e. <tt>a - b</tt>. Can be negative.
46 */
47template <typename T>
48T get_mi_coord(T a, T b, T box_length, T box_length_inv, T box_length_half,
49 bool periodic) {
50 auto const dx = a - b;
51
52 if (periodic && (std::abs(dx) > box_length_half)) {
53 return dx - std::round(dx * box_length_inv) * box_length;
54 }
55
56 return dx;
57}
58
59/**
60 * @brief Get the minimum-image distance between two coordinates.
61 * @param a Coordinate of the terminal point.
62 * @param b Coordinate of the initial point.
63 * @param box_length Box length.
64 * @param periodic Box periodicity.
65 * @return Shortest distance from @p b to @p a across periodic images,
66 * i.e. <tt>a - b</tt>. Can be negative.
67 */
68template <typename T> T get_mi_coord(T a, T b, T box_length, bool periodic) {
69 return get_mi_coord(a, b, box_length, 1. / box_length, 0.5 * box_length,
70 periodic);
71}
72
73/** @brief Calculate image box shift vector.
74 * @param image_box image box offset
75 * @param box box length
76 * @return Image box coordinates.
77 */
78inline auto image_shift(Utils::Vector3i const &image_box,
79 Utils::Vector3d const &box) {
80 return hadamard_product(image_box, box);
81}
82
83/** @brief Unfold particle coordinates to image box.
84 * @param pos coordinate to unfold
85 * @param image_box image box offset
86 * @param box box length
87 * @return Unfolded coordinates.
88 */
89inline auto unfolded_position(Utils::Vector3d const &pos,
90 Utils::Vector3i const &image_box,
91 Utils::Vector3d const &box) {
92 return pos + image_shift(image_box, box);
93}
94} // namespace detail
95
96enum class BoxType { CUBOID = 0, LEES_EDWARDS = 1 };
97
99public:
101 set_length(Utils::Vector3d{1., 1., 1.});
102 set_periodic(0u, true);
103 set_periodic(1u, true);
104 set_periodic(2u, true);
106 }
108 m_type = rhs.type();
109 set_length(rhs.length());
110 set_periodic(0u, rhs.periodic(0u));
111 set_periodic(1u, rhs.periodic(1u));
112 set_periodic(2u, rhs.periodic(2u));
113 m_lees_edwards_bc = rhs.m_lees_edwards_bc;
114 }
115
116private:
117 BoxType m_type = BoxType::CUBOID;
118 /** Flags for all three dimensions whether pbc are applied (default). */
119 std::bitset<3> m_periodic = 0b111;
120 /** Side lengths of the box */
121 Utils::Vector3d m_length = {1., 1., 1.};
122 /** Inverse side lengths of the box */
123 Utils::Vector3d m_length_inv = {1., 1., 1.};
124 /** Half side lengths of the box */
125 Utils::Vector3d m_length_half = {0.5, 0.5, 0.5};
126
127 /** Lees-Edwards boundary conditions */
128 LeesEdwardsBC m_lees_edwards_bc;
129
130public:
131 /**
132 * @brief Set periodicity for direction
133 *
134 * @param coord The coordinate to set the periodicity for.
135 * @param val True if this direction should be periodic.
136 */
137 void set_periodic(unsigned coord, bool val) { m_periodic.set(coord, val); }
138
139 /**
140 * @brief Check periodicity in direction.
141 *
142 * @param coord Direction to check
143 * @return true iff periodic in direction.
144 */
145 constexpr bool periodic(unsigned coord) const {
146 assert(coord <= 2u);
147 return m_periodic[coord];
148 }
149
150 /**
151 * @brief Box length
152 * @return Return vector of side-lengths of the box.
153 */
154 Utils::Vector3d const &length() const { return m_length; }
155
156 /**
157 * @brief Inverse box length
158 * @return Return vector of inverse side-lengths of the box.
159 */
160 Utils::Vector3d const &length_inv() const { return m_length_inv; }
161
162 /**
163 * @brief Half box length
164 * @return Return vector of half side-lengths of the box.
165 */
166 Utils::Vector3d const &length_half() const { return m_length_half; }
167
168 /**
169 * @brief Set box side lengths.
170 * @param box_l Length that should be set.
171 */
172 void set_length(Utils::Vector3d const &box_l) {
173 m_length = box_l;
174 m_length_inv = {1. / box_l[0], 1. / box_l[1], 1. / box_l[2]};
175 m_length_half = 0.5 * box_l;
176 }
177
178 /**
179 * @brief Box volume
180 * @return Return the volume of the box.
181 */
182 double volume() const { return Utils::product(m_length); }
183
184 /**
185 * @brief Get the minimum-image distance between two coordinates.
186 * @param a Coordinate of the terminal point.
187 * @param b Coordinate of the initial point.
188 * @param coord Direction
189 * @return Shortest distance from @p b to @p a across periodic images,
190 * i.e. <tt>a - b</tt>. Can be negative.
191 */
192 template <typename T> T inline get_mi_coord(T a, T b, unsigned coord) const {
193 assert(coord <= 2u);
194
195 return detail::get_mi_coord(a, b, m_length[coord], m_length_inv[coord],
196 m_length_half[coord], m_periodic[coord]);
197 }
198
199 /**
200 * @brief Get the minimum-image vector between two coordinates.
201 *
202 * @tparam T Floating point type.
203 *
204 * @param a Coordinate of the terminal point.
205 * @param b Coordinate of the initial point.
206 * @return Vector from @p b to @p a that minimizes the distance across
207 * periodic images, i.e. <tt>a - b</tt>.
208 */
209 template <typename T>
211 const Utils::Vector<T, 3> &b) const {
212 if (type() == BoxType::LEES_EDWARDS) {
213 auto const shear_plane_normal = lees_edwards_bc().shear_plane_normal;
214 auto a_tmp = a;
215 auto b_tmp = b;
216 a_tmp[shear_plane_normal] = Algorithm::periodic_fold(
217 a_tmp[shear_plane_normal], m_length[shear_plane_normal]);
218 b_tmp[shear_plane_normal] = Algorithm::periodic_fold(
219 b_tmp[shear_plane_normal], m_length[shear_plane_normal]);
220 return lees_edwards_bc().distance(a_tmp - b_tmp, m_length, m_length_half,
221 m_length_inv, m_periodic);
222 }
223 assert(type() == BoxType::CUBOID);
224 return {get_mi_coord(a[0], b[0], 0), get_mi_coord(a[1], b[1], 1),
225 get_mi_coord(a[2], b[2], 2)};
226 }
227
228 BoxType type() const { return m_type; }
229 void set_type(BoxType type) { m_type = type; }
230
231 LeesEdwardsBC const &lees_edwards_bc() const { return m_lees_edwards_bc; }
232 void set_lees_edwards_bc(LeesEdwardsBC bc) { m_lees_edwards_bc = bc; }
233
234 /**
235 * @brief Update the Lees-Edwards parameters of the box geometry
236 * for the current simulation time.
237 */
238 void lees_edwards_update(double pos_offset, double shear_velocity) {
239 assert(type() == BoxType::LEES_EDWARDS);
240 m_lees_edwards_bc.pos_offset = pos_offset;
241 m_lees_edwards_bc.shear_velocity = shear_velocity;
242 }
243
244 /** Calculate the velocity difference including the Lees-Edwards velocity */
246 Utils::Vector3d const &y,
247 Utils::Vector3d const &u,
248 Utils::Vector3d const &v) const {
249 auto ret = u - v;
250 if (type() == BoxType::LEES_EDWARDS) {
251 auto const &le = m_lees_edwards_bc;
252 auto const shear_plane_normal = le.shear_plane_normal;
253 auto const shear_direction = le.shear_direction;
254 auto const dy = x[shear_plane_normal] - y[shear_plane_normal];
255 if (fabs(dy) > 0.5 * length_half()[shear_plane_normal]) {
256 ret[shear_direction] -= Utils::sgn(dy) * le.shear_velocity;
257 }
258 }
259 return ret;
260 }
261
262 /** @brief Fold coordinates to primary simulation box in-place.
263 * Lees-Edwards offset is ignored.
264 * @param[in,out] pos coordinates to fold
265 * @param[in,out] image_box image box offset
266 */
267 void fold_position(Utils::Vector3d &pos, Utils::Vector3i &image_box) const {
268 for (auto i = 0u; i < 3u; i++) {
269 if (m_periodic[i]) {
270 auto const result =
271 Algorithm::periodic_fold(pos[i], image_box[i], m_length[i]);
272 if (result.second == std::numeric_limits<int>::min() or
273 result.second == std::numeric_limits<int>::max()) {
274 throw std::runtime_error(
275 "Overflow in the image box count while folding a particle "
276 "coordinate into the primary simulation box. Maybe a particle "
277 "experienced a huge force.");
278 }
279 std::tie(pos[i], image_box[i]) = result;
280 }
281 }
282 }
283
284 /**
285 * @brief Calculate coordinates folded to primary simulation box.
286 * @param[in] pos coordinates to fold
287 * @return Folded coordinates.
288 */
289 auto folded_position(Utils::Vector3d const &pos) const {
290 auto pos_folded = pos;
291 for (unsigned int i = 0u; i < 3u; i++) {
292 if (m_periodic[i]) {
293 pos_folded[i] = Algorithm::periodic_fold(pos[i], m_length[i]);
294 }
295 }
296
297 return pos_folded;
298 }
299
300 /**
301 * @brief Calculate image box of coordinates folded to primary simulation box.
302 * @param[in] pos coordinates
303 * @param[in] image_box image box to fold
304 * @return Folded image box.
305 */
307 Utils::Vector3i const &image_box) const {
308 auto image_box_folded = image_box;
309 for (auto i = 0u; i < 3u; i++) {
310 if (m_periodic[i]) {
311 image_box_folded[i] =
312 Algorithm::periodic_fold(pos[i], image_box[i], m_length[i]).second;
313 }
314 }
315
316 return image_box_folded;
317 }
318
319 /** @brief Calculate image box shift vector */
320 auto image_shift(Utils::Vector3i const &image_box) const {
321 return detail::image_shift(image_box, m_length);
322 }
323
324 /** @brief Unfold particle coordinates to image box. */
326 Utils::Vector3i const &image_box) const {
327 return detail::unfolded_position(pos, image_box, m_length);
328 }
329};
BoxType
@ LEES_EDWARDS
static int coord(std::string const &s)
Vector implementation and trait types for boost qvm interoperability.
auto unfolded_position(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Unfold particle coordinates to image box.
T get_mi_coord(T a, T b, unsigned coord) const
Get the minimum-image distance between two coordinates.
void lees_edwards_update(double pos_offset, double shear_velocity)
Update the Lees-Edwards parameters of the box geometry for the current simulation time.
auto folded_position(Utils::Vector3d const &pos) const
Calculate coordinates folded to primary simulation box.
Utils::Vector3d const & length() const
Box length.
LeesEdwardsBC const & lees_edwards_bc() const
BoxGeometry(BoxGeometry const &rhs)
constexpr bool periodic(unsigned coord) const
Check periodicity in direction.
double volume() const
Box volume.
BoxType type() const
auto image_shift(Utils::Vector3i const &image_box) const
Calculate image box shift vector.
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.
Utils::Vector3d const & length_half() const
Half box length.
Utils::Vector3d const & length_inv() const
Inverse box length.
void set_periodic(unsigned coord, bool val)
Set periodicity for direction.
void set_length(Utils::Vector3d const &box_l)
Set box side lengths.
auto folded_image_box(Utils::Vector3d const &pos, Utils::Vector3i const &image_box) const
Calculate image box of coordinates folded to primary simulation box.
void set_lees_edwards_bc(LeesEdwardsBC bc)
void fold_position(Utils::Vector3d &pos, Utils::Vector3i &image_box) const
Fold coordinates to primary simulation box in-place.
Utils::Vector3d velocity_difference(Utils::Vector3d const &x, Utils::Vector3d const &y, Utils::Vector3d const &u, Utils::Vector3d const &v) const
Calculate the velocity difference including the Lees-Edwards velocity.
void set_type(BoxType type)
std::pair< T, I > periodic_fold(T x, I i, T const l)
Fold value into primary interval.
T product(Vector< T, N > const &v)
Definition Vector.hpp:374
constexpr int sgn(T val)
Calculate signum of val.
Definition sgn.hpp:27
auto hadamard_product(Vector< T, N > const &a, Vector< U, N > const &b)
Definition Vector.hpp:379
unsigned int shear_plane_normal
Utils::Vector3d distance(Utils::Vector3d const &d, Utils::Vector3d const &l, Utils::Vector3d const &, Utils::Vector3d const &l_inv, std::bitset< 3 > const periodic) const