ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
dlc.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 "config/config.hpp"
23
24#ifdef DIPOLES
25
27
31#include "p3m/common.hpp"
32
33#include "BoxGeometry.hpp"
34#include "Particle.hpp"
36#include "communication.hpp"
37#include "errorhandling.hpp"
38#include "system/System.hpp"
39
40#include <utils/math/sqr.hpp>
41
42#include <boost/mpi/collectives/all_reduce.hpp>
43#include <boost/mpi/collectives/reduce.hpp>
44#include <boost/mpi/operations.hpp>
45
46#include <mpi.h>
47
48#include <algorithm>
49#include <cassert>
50#include <cmath>
51#include <cstdio>
52#include <functional>
53#include <numbers>
54#include <numeric>
55#include <stdexcept>
56#include <variant>
57#include <vector>
58
59void DipolarLayerCorrection::check_gap(Particle const &p) const {
60 if (p.dipm() != 0.) {
61 auto const z = p.pos()[2];
62 if (z < 0. or z > dlc.box_h) {
63 runtimeErrorMsg() << "Particle " << p.id() << " entered DLC gap region "
64 << "by " << z - ((z < 0.) ? 0. : dlc.box_h);
65 }
66 }
67}
68
69/** Calculate the maximal dipole moment in the system */
70static double calc_mu_max(System::System const &system) {
71 auto const local_particles = system.cell_structure->local_particles();
72 auto const mu_max_local = std::accumulate(
73 local_particles.begin(), local_particles.end(), 0.,
74 [](double mu, Particle const &p) { return std::max(mu, p.dipm()); });
75
76 return boost::mpi::all_reduce(comm_cart, mu_max_local,
77 boost::mpi::maximum<double>());
78}
79
80static double g1_DLC_dip(double g, double x) {
81 auto const c = g / x;
82 auto const cc2 = c * c;
83 auto const x3 = x * x * x;
84 auto const a = g * g * g / x + 1.5 * cc2 + 1.5 * g / x3 + 0.75 / (x3 * x);
85 return a;
86}
87
88static double g2_DLC_dip(double g, double x) {
89 auto const x2 = x * x;
90 auto const a = g * g / x + 2.0 * g / x2 + 2.0 / (x2 * x);
91 return a;
92}
93
94/** Compute the box magnetic dipole. */
96
97 Utils::Vector3d local_dip{};
98 for (auto const &p : particles) {
99 if (p.dipm() != 0.) {
100 local_dip += p.calc_dip();
101 }
102 }
103
104 return boost::mpi::all_reduce(comm_cart, local_dip, std::plus<>());
105}
106
107/**
108 * @brief Compute the dipolar force and torque corrections.
109 * Algorithm implemented accordingly to @cite brodka04a.
110 */
111static void dipolar_force_corrections(int kcut,
112 std::vector<Utils::Vector3d> &fs,
113 std::vector<Utils::Vector3d> &ts,
114 ParticleRange const &particles,
115 BoxGeometry const &box_geo) {
116 auto const facux = 2. * std::numbers::pi * box_geo.length_inv()[0];
117 auto const facuy = 2. * std::numbers::pi * box_geo.length_inv()[1];
118
119 auto const n_local_particles = particles.size();
120 std::vector<double> ReSjp(n_local_particles);
121 std::vector<double> ReSjm(n_local_particles);
122 std::vector<double> ImSjp(n_local_particles);
123 std::vector<double> ImSjm(n_local_particles);
124 std::vector<double> ReGrad_Mup(n_local_particles);
125 std::vector<double> ImGrad_Mup(n_local_particles);
126 std::vector<double> ReGrad_Mum(n_local_particles);
127 std::vector<double> ImGrad_Mum(n_local_particles);
128
129 for (int ix = -kcut; ix <= +kcut; ix++) {
130 for (int iy = -kcut; iy <= +kcut; iy++) {
131 if (ix == 0 and iy == 0) {
132 continue;
133 }
134 auto const gx = static_cast<double>(ix) * facux;
135 auto const gy = static_cast<double>(iy) * facuy;
136 auto const gr = sqrt(gx * gx + gy * gy);
137
138 // We assume short slab direction is in the z-direction
139 auto const fa1 = 1. / (gr * (exp(gr * box_geo.length()[2]) - 1.));
140
141 // ... Compute S+,(S+)*,S-,(S-)*, and Spj,Smj for the current g
142
143 std::size_t ip = 0;
144 double S[4] = {0., 0., 0., 0.}; // S of Brodka method, or is S[4] =
145 // {Re(S+), Im(S+), Re(S-), Im(S-)}
146 for (auto const &p : particles) {
147 if (p.dipm() != 0.) {
148 auto const &pos = p.pos();
149 auto const dip = p.calc_dip();
150
151 auto const a = gx * dip[0] + gy * dip[1];
152 auto const b = gr * dip[2];
153 auto const er = gx * pos[0] + gy * pos[1];
154 auto const ez = gr * pos[2];
155 auto const c = cos(er);
156 auto const d = sin(er);
157 auto const f = exp(ez);
158
159 ReSjp[ip] = (b * c - a * d) * f;
160 ImSjp[ip] = (c * a + b * d) * f;
161 ReSjm[ip] = (-b * c - a * d) / f;
162 ImSjm[ip] = (c * a - b * d) / f;
163 ReGrad_Mup[ip] = c * f;
164 ReGrad_Mum[ip] = c / f;
165 ImGrad_Mup[ip] = d * f;
166 ImGrad_Mum[ip] = d / f;
167
168 S[0] += ReSjp[ip];
169 S[1] += ImSjp[ip];
170 S[2] += ReSjm[ip];
171 S[3] += ImSjm[ip];
172 }
173 ip++;
174 }
175
176 MPI_Allreduce(MPI_IN_PLACE, S, 4, MPI_DOUBLE, MPI_SUM, comm_cart);
177
178 // ... Now we can compute the contributions to E,Fj,Ej for the current
179 // g-value
180 ip = 0;
181 for (auto &p : particles) {
182 if (p.dipm() != 0.) {
183 {
184 // compute contributions to the forces
185 auto const s1 = -(-ReSjp[ip] * S[3] + ImSjp[ip] * S[2]);
186 auto const s2 = +(ReSjm[ip] * S[1] - ImSjm[ip] * S[0]);
187 auto const s3 = -(-ReSjm[ip] * S[1] + ImSjm[ip] * S[0]);
188 auto const s4 = +(ReSjp[ip] * S[3] - ImSjp[ip] * S[2]);
189
190 auto const s1z = +(ReSjp[ip] * S[2] + ImSjp[ip] * S[3]);
191 auto const s2z = -(ReSjm[ip] * S[0] + ImSjm[ip] * S[1]);
192 auto const s3z = -(ReSjm[ip] * S[0] + ImSjm[ip] * S[1]);
193 auto const s4z = +(ReSjp[ip] * S[2] + ImSjp[ip] * S[3]);
194
195 auto const ss = s1 + s2 + s3 + s4;
196 fs[ip][0] += fa1 * gx * ss;
197 fs[ip][1] += fa1 * gy * ss;
198 fs[ip][2] += fa1 * gr * (s1z + s2z + s3z + s4z);
199 }
200 {
201 // compute contributions to the electrical field
202 auto const s1 = -(-ReGrad_Mup[ip] * S[3] + ImGrad_Mup[ip] * S[2]);
203 auto const s2 = +(ReGrad_Mum[ip] * S[1] - ImGrad_Mum[ip] * S[0]);
204 auto const s3 = -(-ReGrad_Mum[ip] * S[1] + ImGrad_Mum[ip] * S[0]);
205 auto const s4 = +(ReGrad_Mup[ip] * S[3] - ImGrad_Mup[ip] * S[2]);
206
207 auto const s1z = +(ReGrad_Mup[ip] * S[2] + ImGrad_Mup[ip] * S[3]);
208 auto const s2z = -(ReGrad_Mum[ip] * S[0] + ImGrad_Mum[ip] * S[1]);
209 auto const s3z = -(ReGrad_Mum[ip] * S[0] + ImGrad_Mum[ip] * S[1]);
210 auto const s4z = +(ReGrad_Mup[ip] * S[2] + ImGrad_Mup[ip] * S[3]);
211
212 auto const ss = s1 + s2 + s3 + s4;
213 ts[ip][0] += fa1 * gx * ss;
214 ts[ip][1] += fa1 * gy * ss;
215 ts[ip][2] += fa1 * gr * (s1z + s2z + s3z + s4z);
216 }
217 }
218 ++ip;
219 }
220 }
221 }
222
223 // Convert from the corrections to the electrical field to the corrections
224 // for the torques
225 std::size_t ip = 0;
226 for (auto const &p : particles) {
227 if (p.dipm() != 0.) {
228 ts[ip] = vector_product(p.calc_dip(), ts[ip]);
229 }
230 ++ip;
231 }
232
233 // Multiply by the factors we have left during the loops
234
235 auto const piarea =
236 std::numbers::pi * box_geo.length_inv()[0] * box_geo.length_inv()[1];
237 for (std::size_t j = 0; j < n_local_particles; ++j) {
238 fs[j] *= piarea;
239 ts[j] *= piarea;
240 }
241}
242
243/**
244 * @brief Compute the dipolar DLC energy correction.
245 * Algorithm implemented accordingly to @cite brodka04a.
246 */
247static double dipolar_energy_correction(int kcut,
248 ParticleRange const &particles,
249 BoxGeometry const &box_geo) {
250 auto const facux = 2. * std::numbers::pi * box_geo.length_inv()[0];
251 auto const facuy = 2. * std::numbers::pi * box_geo.length_inv()[1];
252
253 double energy = 0.;
254 double sum_S[4] = {0., 0., 0., 0.};
255 for (int ix = -kcut; ix <= +kcut; ix++) {
256 for (int iy = -kcut; iy <= +kcut; iy++) {
257 if ((ix == 0 && iy == 0)) {
258 continue;
259 }
260 auto const gx = static_cast<double>(ix) * facux;
261 auto const gy = static_cast<double>(iy) * facuy;
262 auto const gr = sqrt(gx * gx + gy * gy);
263
264 // We assume short slab direction is in the z-direction
265 auto const fa1 = 1. / (gr * (exp(gr * box_geo.length()[2]) - 1.));
266
267 // ... Compute S+,(S+)*,S-,(S-)*, and Spj,Smj for the current g
268
269 double S[4] = {0., 0., 0., 0.}; // S of Brodka method, or is S[4] =
270 // {Re(S+), Im(S+), Re(S-), Im(S-)}
271 for (auto const &p : particles) {
272 if (p.dipm() != 0.) {
273 auto const &pos = p.pos();
274 auto const dip = p.calc_dip();
275
276 auto const a = gx * dip[0] + gy * dip[1];
277 auto const b = gr * dip[2];
278 auto const er = gx * pos[0] + gy * pos[1];
279 auto const ez = gr * pos[2];
280 auto const c = cos(er);
281 auto const d = sin(er);
282 auto const f = exp(ez);
283
284 S[0] += (b * c - a * d) * f;
285 S[1] += (c * a + b * d) * f;
286 S[2] += (-b * c - a * d) / f;
287 S[3] += (c * a - b * d) / f;
288 }
289 }
290 boost::mpi::reduce(comm_cart, S, 4, sum_S, std::plus<>(), 0);
291
292 // compute contribution to the energy
293 // s2=(ReSm*ReSp+ImSm*ImSp); s2=s1!!!
294 energy += fa1 * 2. * (sum_S[0] * sum_S[2] + sum_S[1] * sum_S[3]);
295 }
296 }
297
298 auto const piarea =
299 std::numbers::pi * box_geo.length_inv()[0] * box_geo.length_inv()[1];
300 energy *= -piarea;
301 return (this_node == 0) ? energy : 0.;
302}
303
305 ParticleRange const &particles) const {
306 assert(dlc.far_cut > 0.);
307 auto const &box_geo = *get_system().box_geo;
308 auto const volume = box_geo.volume();
309 auto const correc = 4. * std::numbers::pi / volume;
310
311 // --- Create arrays that should contain the corrections to
312 // the forces and torques, and set them to zero.
313 std::vector<Utils::Vector3d> dip_DLC_f(particles.size());
314 std::vector<Utils::Vector3d> dip_DLC_t(particles.size());
315
316 //---- Compute the corrections ----------------------------------
317
318 // First the DLC correction
319 auto const kcut = static_cast<int>(std::round(dlc.far_cut));
320 dipolar_force_corrections(kcut, dip_DLC_f, dip_DLC_t, particles, box_geo);
321
322 // Now we compute the correction like Yeh and Klapp to take into account
323 // the fact that you are using a 3D PBC method which uses spherical
324 // summation instead of slab-wise summation.
325 // Slab-wise summation is the one required to apply DLC correction.
326 // This correction is often called SDC = Shape Dependent Correction.
327 // See @cite brodka04a.
328
329 auto const box_dip = calc_slab_dipole(particles);
330
331 // --- Transfer the computed corrections to the Forces, Energy and torques
332 // of the particles
333
334 std::size_t ip = 0;
335 for (auto &p : particles) {
336 check_gap(p);
337
338 if (p.dipm() != 0.) {
339 // SDC correction term is zero for the forces
340 p.force() += prefactor * dip_DLC_f[ip];
341
342 auto const dip = p.calc_dip();
343 // SDC correction for the torques
344 auto d = Utils::Vector3d{0., 0., -correc * box_dip[2]};
345#ifdef DP3M
347 d += correc * epsilon_correction * box_dip;
348 }
349#endif
350 p.torque() += prefactor * (dip_DLC_t[ip] + vector_product(dip, d));
351 }
352 ++ip;
353 }
354}
355
357 ParticleRange const &particles) const {
358 assert(dlc.far_cut > 0.);
359 auto const &box_geo = *get_system().box_geo;
360 auto const volume = box_geo.volume();
361 auto const pref = prefactor * 2. * std::numbers::pi / volume;
362
363 // Check if particles aren't in the forbidden gap region
364 // This loop is needed, because there is no other guaranteed
365 // single pass over all particles in this function.
366 for (auto const &p : particles) {
367 check_gap(p);
368 }
369
370 //---- Compute the corrections ----------------------------------
371
372 // First the DLC correction
373 auto const k_cut = static_cast<int>(std::round(dlc.far_cut));
374 auto dip_DLC_energy =
375 prefactor * dipolar_energy_correction(k_cut, particles, box_geo);
376
377 // Now we compute the correction like Yeh and Klapp to take into account
378 // the fact that you are using a 3D PBC method which uses spherical
379 // summation instead of slab-wise summation.
380 // Slab-wise summation is the one required to apply DLC correction.
381 // This correction is often called SDC = Shape Dependent Correction.
382 // See @cite brodka04a.
383
384 auto const box_dip = calc_slab_dipole(particles);
385
386 if (this_node == 0) {
387 dip_DLC_energy += pref * Utils::sqr(box_dip[2]);
388#ifdef DP3M
390 dip_DLC_energy -= pref * epsilon_correction * box_dip.norm2();
391 }
392#endif
393 return dip_DLC_energy;
394 }
395 return 0.;
396}
397
398static int count_magnetic_particles(System::System const &system) {
399 int local_n = 0;
400
401 for (auto const &p : system.cell_structure->local_particles()) {
402 if (p.dipm() != 0.) {
403 local_n++;
404 }
405 }
406
407 return boost::mpi::all_reduce(comm_cart, local_n, std::plus<>());
408}
409
410/** Compute the cut-off in the DLC dipolar part to get a certain accuracy.
411 * We assume particles to have all the same value of the dipolar momentum
412 * modulus, which is taken as the largest value of mu inside the system.
413 * If we assume the gap has a width @c gap_size (within
414 * which there is no particles): <tt>Lz = h + gap_size</tt>
415 */
416double DipolarLayerCorrection::tune_far_cut() const {
417 auto const &system = get_system();
418 /* we take the maximum dipole in the system, to be sure that the errors
419 * in the other case will be equal or less than for this one */
420 auto const mu_max_sq = Utils::sqr(calc_mu_max(system));
421 auto const &box_geo = *system.box_geo;
422 auto const lx = box_geo.length()[0];
423 auto const ly = box_geo.length()[1];
424 auto const lz = box_geo.length()[2];
425
426 if (std::abs(lx - ly) > 0.001) {
427 throw std::runtime_error("DLC tuning: box size in x direction is "
428 "different from y direction. The tuning "
429 "formula requires both to be equal.");
430 }
431
432 auto constexpr limitkc = 200;
433 auto const piarea = std::numbers::pi / (lx * ly);
434 auto const nmp = static_cast<double>(count_magnetic_particles(system));
435 auto const h = dlc.box_h;
436 auto far_cut = -1.;
437 for (int kc = 1; kc < limitkc; kc++) {
438 auto const gc = kc * 2. * std::numbers::pi / lx;
439 auto const fa0 = sqrt(9. * exp(+2. * gc * h) * g1_DLC_dip(gc, lz - h) +
440 9. * exp(-2. * gc * h) * g1_DLC_dip(gc, lz + h) +
441 22. * g1_DLC_dip(gc, lz));
442 auto const fa1 = sqrt(0.125 * piarea) * fa0;
443 auto const fa2 = g2_DLC_dip(gc, lz);
444 auto const de = nmp * mu_max_sq / (4. * (exp(gc * lz) - 1.)) * (fa1 + fa2);
445 if (de < dlc.maxPWerror) {
446 far_cut = static_cast<double>(kc);
447 break;
448 }
449 }
450 if (far_cut <= 0.) {
451 throw std::runtime_error("DLC tuning failed: maxPWerror too small");
452 }
453 return far_cut;
454}
455
456/** @brief Lock an actor and modify its parameters. */
458 explicit AdaptSolver(DipolarLayerCorrection *this_ptr) : m_actor{this_ptr} {}
459
460 void operator()(std::shared_ptr<DipolarDirectSum> const &solver) {
461 m_actor->prefactor = solver->prefactor;
462 m_actor->epsilon = P3M_EPSILON_METALLIC;
463 }
464
465#ifdef DP3M
466 void operator()(std::shared_ptr<DipolarP3M> const &solver) {
467 m_actor->prefactor = solver->prefactor;
468 m_actor->epsilon = solver->dp3m_params.epsilon;
469 }
470#endif
471
472private:
473 DipolarLayerCorrection *m_actor;
474};
475
477 auto visitor = AdaptSolver{this};
478 std::visit(visitor, base_solver);
480 (epsilon == P3M_EPSILON_METALLIC) ? 0. : 1. / (2. * epsilon + 1.);
481}
482
483void DipolarLayerCorrection::sanity_checks_node_grid() const {
484 auto const &box_geo = *get_system().box_geo;
485 if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) {
486 throw std::runtime_error("DLC: requires periodicity (True, True, True)");
487 }
488}
489
490dlc_data::dlc_data(double maxPWerror, double gap_size, double far_cut)
491 : maxPWerror{maxPWerror}, gap_size{gap_size}, box_h{-1.}, far_cut{far_cut},
492 far_calculated{far_cut == -1.} {
493 if (far_cut <= 0. and not far_calculated) {
494 throw std::domain_error("Parameter 'far_cut' must be > 0");
495 }
496 if (maxPWerror <= 0.) {
497 throw std::domain_error("Parameter 'maxPWerror' must be > 0");
498 }
499 if (gap_size <= 0.) {
500 throw std::domain_error("Parameter 'gap_size' must be > 0");
501 }
502}
503
505 auto const lz = get_system().box_geo->length()[2];
506 auto const new_box_h = lz - dlc.gap_size;
507 if (new_box_h < 0.) {
508 throw std::runtime_error("DLC gap size (" + std::to_string(dlc.gap_size) +
509 ") larger than box length in z-direction (" +
510 std::to_string(lz) + ")");
511 }
512 dlc.box_h = new_box_h;
513}
514
516 BaseSolver &&solver)
517 : dlc{parameters}, base_solver{solver} {
518 adapt_solver();
519}
520
521#endif // DIPOLES
Utils::Vector3d const & length() const
Box length.
Utils::Vector3d const & length_inv() const
Inverse box length.
A range of particles.
base_type::size_type size() const
Main system class.
std::shared_ptr< CellStructure > cell_structure
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
This file contains the defaults for ESPResSo.
__device__ void vector_product(float const *a, float const *b, float *out)
static void dipolar_force_corrections(int kcut, std::vector< Utils::Vector3d > &fs, std::vector< Utils::Vector3d > &ts, ParticleRange const &particles, BoxGeometry const &box_geo)
Compute the dipolar force and torque corrections.
Definition dlc.cpp:111
static double g1_DLC_dip(double g, double x)
Definition dlc.cpp:80
static double dipolar_energy_correction(int kcut, ParticleRange const &particles, BoxGeometry const &box_geo)
Compute the dipolar DLC energy correction.
Definition dlc.cpp:247
static double g2_DLC_dip(double g, double x)
Definition dlc.cpp:88
static Utils::Vector3d calc_slab_dipole(ParticleRange const &particles)
Compute the box magnetic dipole.
Definition dlc.cpp:95
static double calc_mu_max(System::System const &system)
Calculate the maximal dipole moment in the system.
Definition dlc.cpp:70
static int count_magnetic_particles(System::System const &system)
Definition dlc.cpp:398
P3M algorithm for long-range magnetic dipole-dipole interaction.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
Vector< T, N > sqrt(Vector< T, N > const &a)
Definition Vector.hpp:357
Common functions for dipolar and charge P3M.
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
Lock an actor and modify its parameters.
Definition dlc.cpp:457
AdaptSolver(DipolarLayerCorrection *this_ptr)
Definition dlc.cpp:458
void operator()(std::shared_ptr< DipolarDirectSum > const &solver)
Definition dlc.cpp:460
void operator()(std::shared_ptr< DipolarP3M > const &solver)
Definition dlc.cpp:466
Adapt a magnetostatics solver to remove contributions from the z-direction.
Definition dlc.hpp:78
double epsilon_correction
Definition dlc.hpp:89
double energy_correction(ParticleRange const &particles) const
Calculate the dipolar energy correction.
Definition dlc.cpp:356
BaseSolver base_solver
Magnetostatics solver that is adapted.
Definition dlc.hpp:94
void add_force_corrections(ParticleRange const &particles) const
Add the dipolar force and torque corrections.
Definition dlc.cpp:304
DipolarLayerCorrection(dlc_data &&parameters, BaseSolver &&solver)
Definition dlc.cpp:515
std::variant< std::shared_ptr< DipolarP3M >, std::shared_ptr< DipolarDirectSum > > BaseSolver
Definition dlc.hpp:83
Struct holding all information for one particle.
Definition Particle.hpp:395
auto calc_dip() const
Definition Particle.hpp:495
auto const & torque() const
Definition Particle.hpp:479
auto const & dipm() const
Definition Particle.hpp:493
auto const & pos() const
Definition Particle.hpp:431
auto const & force() const
Definition Particle.hpp:435
auto const & id() const
Definition Particle.hpp:414
Parameters for the DLC method.
Definition dlc.hpp:47
bool far_calculated
Flag whether far_cut was set by the user, or calculated by ESPResSo.
Definition dlc.hpp:71
dlc_data(double maxPWerror, double gap_size, double far_cut)
Definition dlc.cpp:490
double maxPWerror
maximal pairwise error of the potential and force
Definition dlc.hpp:51
double far_cut
Cutoff of the exponential sum.
Definition dlc.hpp:65
double gap_size
Size of the empty gap.
Definition dlc.hpp:56
double box_h
Up to where particles can be found.
Definition dlc.hpp:59