GCC Code Coverage Report

Directory: ./
File: src/integrator.cpp
Date: 2025-07-07 11:13:23
Exec Total Coverage
Lines: 374 427 87.6%
Functions: 47 53 88.7%
Branches: 221 397 55.7%
Line Branch Exec Source
1 #include <algorithm>
2 #include <iomanip>
3 #include <iostream>
4 #ifndef NO_OMP
5 #include <omp.h>
6 #endif
7
8 #include "integrator.hpp"
9 #include "interaction.hpp"
10 #include "particle.hpp"
11 #include "rng.hpp"
12
13 #define RegisterIntegrator(A) \
14 { \
15 #A, [](System *s, const YAML::Node &c) -> Integrator * { return new A(s, c); } \
16 }
17
18 std::map<std::string, Integrator *(*)(System *s, const YAML::Node &c)>
19 IntegratorFactory::integratorMap = {
20 ✗ RegisterIntegrator(Static),
21
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 RegisterIntegrator(BrownianEuler),
22
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 RegisterIntegrator(BrownianHeunEulerRSwM2),
23
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 RegisterIntegrator(BrownianHeunEulerRSwM3),
24
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 RegisterIntegrator(MolecularVelocityVerlet),
25
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 RegisterIntegrator(MolecularNoseHooverVelocityVerlet),
26
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 RegisterIntegrator(LangevinSymplecticEuler),
27
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 RegisterIntegrator(MonteCarloCanonical),
28
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 2 RegisterIntegrator(MonteCarloGrandCanonical),
29 };
30
31 #undef RegisterIntegrator
32
33 19 Integrator *IntegratorFactory::create(System *system, const YAML::Node &config) {
34 19 std::string integratorName = "Static"; // default to Static
35
1/2
✓ Branch 0 taken 19 times.
✗ Branch 1 not taken.
 19 if (config) {
36
4/8
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 19 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 19 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 19 times.
✗ Branch 13 not taken.
 57 integratorName = config.begin()->first.as<std::string>();
37 }
38
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 const YAML::Node &c = config[integratorName];
39
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 if (auto it = integratorMap.find(integratorName); it != integratorMap.end()) {
40
1/2
✓ Branch 1 taken 19 times.
✗ Branch 2 not taken.
 19 return it->second(system, c);
41 } else {
42 ✗ throw std::invalid_argument("Integrator type " + integratorName + " is not implemented!");
43 }
44 38 }
45
46 19 Integrator::Integrator(System *system) : system(system) { std::clog << "\nIntegrator:"; }
47
48 ✗ double Integrator::weight() { return 1.; }
49
50 ✗ Static::Static(System *s, const YAML::Node &c) : Integrator(s) { std::clog << "\n\ttype: Static"; }
51
52 ✗ void Static::step() {}
53
54 RNG Brownian::rng = RNG();
55
56 9 Brownian::Brownian(System *s, const YAML::Node &c)
57
3/6
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 9 times.
✗ Branch 9 not taken.
 18 : Integrator(s), dt(YAML::parse<double>(c, "dt", 1e-4)),
58
4/8
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 9 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 9 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 9 times.
 27 randomSeed(YAML::parse<bool>(c, "randomSeed", true)) {
59
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if (!system->interaction->allDoF) {
60 ✗ throw std::invalid_argument(
61 "Interactions must implement the force calculation for use with a Brownian Dynamics "
62 ✗ "integrator. This is not the case for one of the chosen interactions.");
63 }
64
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 9 times.
 110 for (const auto &p : system->particles) {
65
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 101 times.
 101 if (!std::isfinite(p.gamma)) {
66 ✗ throw std::invalid_argument("You must specify the friction coefficient gamma of all "
67 ✗ "particles when using a Brownian Dynamics integrator.");
68 }
69 }
70 #ifndef NO_OMP
71 9 omp_set_dynamic(0); // this has to be off for consistent threadprivate behavior
72 #endif
73 // Workaround for IntelLLVM, which fails to copy the complete state of RNG into the threadprivate
74 // variable (the distributions remain uninitialized such that all threads except master produce a
75 // constant 0 as the random output). We therefore have to do the copying of the RNG state
76 // ourselves.
77
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 RNG rng_{};
78
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9 times.
 9 if (!randomSeed) {
79 #ifndef NO_OMP
80 ✗ if (omp_get_max_threads() > 1) {
81 ✗ std::clog << "\nWARNING: You have specified randomSeed = " << randomSeed
82 ✗ << " while using multithreading (number of OpenMP threads: "
83 << omp_get_max_threads()
84 << "). Consider setting the thread count to 1 (e.g. via the environment variable "
85 ✗ "OMP_NUM_THREADS) to ensure reproducible runs.";
86 }
87 #endif
88 ✗ rng_.engine.seed(0);
89 }
90 9 #pragma omp parallel
91 {
92 // Copy the complete state of the reference RNG
93 rng = rng_;
94 // Set individual streams (which ensure independent, non-overlapping random sequences) for each
95 // thread (feature of PCG)
96 #ifndef NO_OMP
97 rng.engine.set_stream(omp_get_thread_num());
98 #endif
99 }
100
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 9 std::clog << "\n\ttype: Brownian Dynamics";
101
2/4
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 9 times.
✗ Branch 5 not taken.
 9 std::clog << "\n\trandomSeed: " << this->randomSeed;
102 9 }
103
104 ✗ double Brownian::weight() { return dt; }
105
106 3 BrownianEuler::BrownianEuler(System *s, const YAML::Node &c) : Brownian(s, c) {
107
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 std::clog << "\n\tsubtype: Brownian Euler";
108
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
 3 std::clog << "\n\tdt: " << this->dt;
109 3 }
110
111 2 void BrownianEuler::step() {
112 2 #pragma omp parallel for
113 for (auto &p : system->particles) {
114 p.rRand = sqrt(2 * system->T * dt / p.gamma) * rng.vecNormal();
115 p.r += 1. / p.gamma * p.F() * dt + p.rRand;
116 }
117 2 system->applyPeriodic();
118 2 system->updateState();
119 2 system->t += dt;
120 2 }
121
122 6 BrownianHeunEulerRSwM::BrownianHeunEulerRSwM(System *s, const YAML::Node &c)
123
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
 12 : Brownian(s, c), dtTrial(dt), relTol(YAML::parse<double>(c, "relTol", 0.05)),
124
5/10
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 6 times.
✗ Branch 16 not taken.
 18 absTol(YAML::parse<double>(c, "absTol", 0.01)), qmin(YAML::parse<double>(c, "qmin", 1e-3)),
125
5/10
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
✓ Branch 15 taken 6 times.
✗ Branch 16 not taken.
 18 qmax(YAML::parse<double>(c, "qmax", 1.125)), alpha(YAML::parse<double>(c, "alpha", 2.0)),
126
4/8
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 6 times.
✗ Branch 12 not taken.
 18 beta(YAML::parse<double>(c, "beta", 0.9)), lNorm(YAML::parse<int>(c, "lNorm", -1)),
127
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 R(VecVec::Zero(3, system->particles.size())),
128
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 Rbridge(VecVec::Zero(3, system->particles.size())),
129
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 Rdiff(VecVec::Zero(3, system->particles.size())),
130
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 Rtmp(VecVec::Zero(3, system->particles.size())),
131
1/2
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
 6 drs0(VecVec::Zero(3, system->particles.size())),
132
2/4
✓ Branch 2 taken 6 times.
✗ Branch 3 not taken.
✓ Branch 8 taken 6 times.
✗ Branch 9 not taken.
 24 drs1(VecVec::Zero(3, system->particles.size())) {
133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 if (lNorm != 2 && lNorm != -1) {
134 ✗ throw std::invalid_argument("lNorm must be 2 (for l2 norm) or -1 (for maximum norm)");
135 }
136
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 generateRandVecs(R, dtTrial);
137
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 std::clog << "\n\tsubtype: Brownian Heun-Euler RSwM";
138
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\tdtInit: " << this->dt;
139
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\trelTol: " << this->relTol;
140
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\tabsTol: " << this->absTol;
141
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\tqmax: " << this->qmax;
142
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\tqmin: " << this->qmin;
143
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\talpha: " << this->alpha;
144
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
 6 std::clog << "\n\tbeta: " << this->beta;
145 6 }
146
147 12 void BrownianHeunEulerRSwM::generateRandVecs(VecVec &buffer, Real stddev) {
148 12 #pragma omp parallel for
149 for (int i = 0; i < buffer.cols(); ++i) {
150 buffer.col(i) = rng.vecNormal(stddev);
151 }
152 12 }
153
154 12 void BrownianHeunEulerRSwM::generateRandVecsBridge(VecVec &buffer, const VecVec &increment,
155 double q, double dt) {
156 12 #pragma omp parallel for
157 for (int i = 0; i < buffer.cols(); ++i) {
158 buffer.col(i) = rng.vecNormal(q * increment.col(i), (1 - q) * q * dt);
159 }
160 12 }
161
162 12 void BrownianHeunEulerRSwM::maybePush(RandomIncrementStack &S, const RandomIncrement &RI) {
163
1/2
✓ Branch 0 taken 12 times.
✗ Branch 1 not taken.
 12 if (RI.dt > dtCutoff) {
164 12 S.push(RI);
165 }
166 12 }
167
168 62 double BrownianHeunEulerRSwM::estimateParticleError(int i) {
169
1/2
✓ Branch 3 taken 62 times.
✗ Branch 4 not taken.
 62 return (drs1.col(i) - drs0.col(i)).norm() /
170
1/2
✓ Branch 3 taken 62 times.
✗ Branch 4 not taken.
 124 (absTol + 0.5 * (drs0.col(i) + drs1.col(i)).norm() * relTol);
171 }
172
173 8 double BrownianHeunEulerRSwM::estimateTotalError() {
174 8 double err = 0.;
175
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 size_t N = system->particles.size();
176
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (lNorm == 2) {
177 ✗ #pragma omp parallel for reduction(+ : err)
178 for (size_t i = 0; i < N; ++i) {
179 err += estimateParticleError(i);
180 }
181 ✗ err = sqrt(err / N);
182 } else {
183 8 #pragma omp parallel for reduction(max : err)
184 for (size_t i = 0; i < N; ++i) {
185 err = fmax(estimateParticleError(i), err);
186 }
187 }
188 8 return err;
189 }
190
191 8 bool BrownianHeunEulerRSwM::adaptStepsize() {
192 8 bool needSmallerStep = false;
193 8 error = estimateTotalError();
194 8 q = 1. / (alpha * alpha * error * error);
195
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 if (q < 1) {
196 4 needSmallerStep = true;
197 4 rejections++;
198 4 rejectStep();
199 } else {
200 4 needSmallerStep = false;
201 4 dt = dtTrial;
202 4 prepareNextStep();
203 }
204 8 return needSmallerStep;
205 }
206
207 4 void BrownianHeunEulerRSwM::step() {
208 4 rejections = 0;
209 4 particlesBefore = system->particles;
210 4 size_t N = system->particles.size();
211 8 do {
212
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 8 if (rejections > 0) {
213 4 system->particles = particlesBefore;
214 }
215 8 #pragma omp parallel for
216 for (size_t i = 0; i < N; ++i) {
217 Particle &p = system->particles[i];
218 p.rRand = sqrt(2 * system->T / p.gamma) * R.col(i);
219 drs0.col(i) = 1. / p.gamma * p.F() * dtTrial + p.rRand;
220 // Move particles for evaluation on second positions
221 p.r += drs0.col(i);
222 }
223 8 system->applyPeriodic();
224 8 system->updateState();
225 8 #pragma omp parallel for
226 for (size_t i = 0; i < N; ++i) {
227 Particle &p = system->particles[i];
228 drs1.col(i) = 1. / p.gamma * p.F() * dtTrial + p.rRand;
229 }
230
2/2
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 4 times.
 8 } while (adaptStepsize());
231 // Do the actual move based on the Heun algorithm
232 // Notice that we have already moved particles by drs0!
233 4 #pragma omp parallel for
234 for (size_t i = 0; i < N; ++i) {
235 system->particles[i].r += 0.5 * (drs1.col(i) - drs0.col(i));
236 }
237 4 system->applyPeriodic();
238 4 system->updateState();
239 4 system->t += dt;
240 4 std::ios coutStateBefore(nullptr);
241
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 coutStateBefore.copyfmt(std::cout);
242
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
 4 std::cout << ", dt: " << std::setprecision(3) << std::scientific << dt;
243
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::cout.copyfmt(coutStateBefore);
244 4 }
245
246 3 BrownianHeunEulerRSwM2::BrownianHeunEulerRSwM2(System *s, const YAML::Node &c)
247
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 : BrownianHeunEulerRSwM(s, c) {
248
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 std::clog << "\n\tRSwM type: RSwM2";
249 3 }
250
251 4 void BrownianHeunEulerRSwM2::rejectStep() {
252 4 q = fmax(qmin, q);
253 4 generateRandVecsBridge(Rbridge, R, q, dtTrial);
254 4 Rdiff = R - Rbridge;
255 12 maybePush(Sf, {(1 - q) * dtTrial, Rdiff});
256 4 dtTrial *= q;
257 4 R = Rbridge;
258
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 }
259
260 5 void BrownianHeunEulerRSwM2::prepareNextStep() {
261 5 dtTrial *= fmin(qmax, beta * q);
262 // dtTrial = fmin(0.0005, dtTrial);
263 5 double dts = 0;
264 5 R.setZero();
265
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
 7 while (!Sf.empty()) {
266
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 RandomIncrement &L = Sf.top();
267 4 double dtRest = L.dt;
268
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 if (dts + dtRest <= dtTrial) {
269 2 dts += dtRest;
270 2 R += L.R;
271 2 Sf.pop();
272 } else {
273 2 double qTmp = (dtTrial - dts) / dtRest;
274 2 generateRandVecsBridge(Rbridge, L.R, qTmp, dtRest);
275 2 Rdiff = L.R - Rbridge;
276 2 Sf.pop();
277 6 maybePush(Sf, {(1 - qTmp) * dtRest, Rdiff});
278 2 R += Rbridge;
279 2 dts += qTmp * dtRest;
280 2 break;
281 }
282 }
283
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
 5 if (dtTrial - dts > dtCutoff) {
284 3 generateRandVecs(R, dtTrial - dts);
285 }
286
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 7 }
287
288 3 BrownianHeunEulerRSwM3::BrownianHeunEulerRSwM3(System *s, const YAML::Node &c)
289
1/2
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
 3 : BrownianHeunEulerRSwM(s, c) {
290
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 Su.push({dtTrial, R});
291
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 std::clog << "\n\tRSwM type: RSwM3";
292
2/4
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
 9 }
293
294 4 void BrownianHeunEulerRSwM3::rejectStep() {
295 4 q = fmax(qmin, q);
296 4 double dts = 0;
297 4 Rtmp.setZero();
298
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 while (!Su.empty()) {
299 4 RandomIncrement &L = Su.top();
300 4 dts += L.dt;
301 4 Rtmp += L.R;
302
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (dts < (1 - q) * dtTrial) {
303 ✗ Sf.push(L);
304 ✗ Su.pop();
305 } else {
306 4 double dtu = L.dt;
307 4 double dtM = dts - (1 - q) * dtTrial;
308 4 double qM = dtM / dtu;
309 4 generateRandVecsBridge(Rbridge, L.R, qM, dtu);
310 4 Rdiff = L.R - Rbridge;
311 4 R += Rbridge - Rtmp;
312 ✗ maybePush(Sf, {(1 - qM) * dtu, Rdiff});
313 4 Su.pop();
314 4 Su.push({qM * dtu, Rbridge});
315 4 break;
316 }
317 }
318 4 dtTrial *= q;
319
2/4
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 4 times.
✗ Branch 7 not taken.
 8 }
320
321 5 void BrownianHeunEulerRSwM3::prepareNextStep() {
322 5 dtTrial *= fmin(qmax, beta * q);
323 // dtTrial = fmin(0.0005, dtTrial);
324
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 5 times.
 12 while (!Su.empty()) {
325 7 Su.pop();
326 }
327 5 double dts = 0;
328 5 R.setZero();
329
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 3 times.
 7 while (!Sf.empty()) {
330
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 RandomIncrement &L = Sf.top();
331 4 double dtf = L.dt;
332
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
 4 if (dts + dtf <= dtTrial) {
333 2 dts += dtf;
334 2 R += L.R;
335 2 Su.push(L);
336 2 Sf.pop();
337 } else {
338 2 double qM = (dtTrial - dts) / dtf;
339 2 generateRandVecsBridge(Rbridge, L.R, qM, dtf);
340 2 Rdiff = L.R - Rbridge;
341 2 R += Rbridge;
342 2 Sf.pop();
343 ✗ maybePush(Sf, {(1 - qM) * dtf, Rdiff});
344 2 Su.push({qM * dtf, Rbridge});
345 2 dts += qM * dtf;
346 2 break;
347 }
348 }
349
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 2 times.
 5 if (dtTrial - dts > dtCutoff) {
350 3 generateRandVecs(Rtmp, dtTrial - dts);
351 3 R += Rtmp;
352 3 Su.push({dtTrial - dts, Rtmp});
353 }
354
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 10 taken 3 times.
✗ Branch 11 not taken.
 10 }
355
356 4 Molecular::Molecular(System *s, const YAML::Node &c)
357
3/6
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 4 times.
 8 : Integrator(s), dt(YAML::parse<double>(c, "dt", 1e-4)) {
358
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (!system->interaction->allDoF) {
359 ✗ throw std::invalid_argument(
360 "Interactions must implement the force calculation for use with a Molecular Dynamics "
361 ✗ "integrator. This is not the case for one of the chosen interactions.");
362 }
363
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 4 times.
 26 for (const auto &p : system->particles) {
364
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 22 times.
 22 if (!std::isfinite(p.m)) {
365 ✗ throw std::invalid_argument("You must specify the mass m of all particles when using a "
366 ✗ "Molecular Dynamics integrator.");
367 }
368 }
369
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::clog << "\n\ttype: Molecular Dynamics";
370
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::clog << "\n\tdt: " << this->dt;
371 4 }
372
373 2 MolecularVelocityVerlet::MolecularVelocityVerlet(System *s, const YAML::Node &c) : Molecular(s, c) {
374
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::clog << "\n\tsubtype: Molecular Velocity-Verlet";
375 2 }
376
377 2 void MolecularVelocityVerlet::step() {
378 2 #pragma omp parallel for
379 for (auto &p : system->particles) {
380 p.r += p.v * dt + 0.5 * p.F() / p.m * dt * dt; // r(t + dt)
381 p.v += 0.5 * p.F() / p.m * dt; // v(t + dt/2)
382 }
383 2 system->applyPeriodic();
384 2 system->updateState(); // f(t + dt)
385 2 #pragma omp parallel for
386 for (auto &p : system->particles) {
387 p.v += 0.5 * p.F() / p.m * dt; // v(t + dt)
388 }
389 2 system->t += dt;
390 2 }
391
392 2 MolecularNoseHooverVelocityVerlet::MolecularNoseHooverVelocityVerlet(System *s, const YAML::Node &c)
393 2 : Molecular(s, c), dt2(dt / 2.), dt4(dt / 4.), dt8(dt / 8.), Ekin(0.),
394
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
 4 Q(YAML::parse<double>(c, "Q", 100)), xi1(0.), vxi1(0.), xi2(0.), vxi2(0.) {
395
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::clog << "\n\tsubtype: Molecular Nose-Hoover Velocity-Verlet";
396
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\ttDamp: " << 1. / this->Q;
397 2 }
398
399 4 void MolecularNoseHooverVelocityVerlet::chain() {
400 4 double L = DIM * system->particles.size();
401 4 double G2 = (Q * vxi1 * vxi1 - system->T);
402 4 vxi2 += G2 * dt4;
403 4 vxi1 *= std::exp(-vxi2 * dt8);
404 4 double G1 = (2. * Ekin - L * system->T) / Q;
405 4 vxi1 += G1 * dt4;
406 4 vxi1 *= std::exp(-vxi2 * dt8);
407 4 xi1 += vxi1 * dt2;
408 4 xi2 += vxi2 * dt2;
409 4 double s = std::exp(-vxi1 * dt2);
410
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 4 times.
 26 for (Particle &p : system->particles) {
411 22 p.v *= s;
412 }
413 4 Ekin *= s * s;
414 4 vxi1 *= std::exp(-vxi2 * dt8);
415 4 G1 = (2. * Ekin - L * system->T) / Q;
416 4 vxi1 += G1 * dt4;
417 4 vxi1 *= std::exp(-vxi2 * dt8);
418 4 G2 = (Q * vxi1 * vxi1 - system->T) / Q;
419 4 vxi2 += G2 * dt4;
420 4 }
421
422 2 void MolecularNoseHooverVelocityVerlet::posvel() {
423 2 Ekin = 0.;
424 2 #pragma omp parallel for
425 for (Particle &p : system->particles) {
426 p.r += p.v * dt2;
427 }
428 2 system->applyPeriodic();
429 2 system->updateState();
430 2 #pragma omp parallel for reduction(+ : Ekin)
431 for (Particle &p : system->particles) {
432 p.v += p.F() / p.m * dt;
433 p.r += p.v * dt2;
434 Ekin += 0.5 * p.m * p.v.squaredNorm();
435 }
436 2 system->applyPeriodic();
437 2 system->updateState();
438 2 }
439
440 2 void MolecularNoseHooverVelocityVerlet::step() {
441 2 chain();
442 2 posvel();
443 2 chain();
444 2 system->t += dt;
445 2 }
446
447 RNG Langevin::rng = RNG();
448
449 2 Langevin::Langevin(System *s, const YAML::Node &c)
450
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 4 : Integrator(s), dt(YAML::parse<double>(c, "dt", 1e-4)),
451
4/8
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 2 times.
 6 randomSeed(YAML::parse<bool>(c, "randomSeed", true)) {
452
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (!system->interaction->allDoF) {
453 ✗ throw std::invalid_argument(
454 "Interactions must implement the force calculation for use with a Langevin Dynamics "
455 ✗ "integrator. This is not the case for one of the chosen interactions.");
456 }
457
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 2 times.
 13 for (const auto &p : system->particles) {
458
2/4
✓ Branch 0 taken 11 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 11 times.
 11 if (!std::isfinite(p.gamma) || !std::isfinite(p.m)) {
459 ✗ throw std::invalid_argument("You must specify the friction coefficient gamma and the mass m "
460 ✗ "of all particles when using a Langevin Dynamics integrator.");
461 }
462 }
463 #ifndef NO_OMP
464 2 omp_set_dynamic(0); // this has to be off for consistent threadprivate behavior
465 #endif
466 // Workaround for IntelLLVM, which fails to copy the complete state of RNG into the threadprivate
467 // variable (the distributions remain uninitialized such that all threads except master produce a
468 // constant 0 as the random output). We therefore have to do the copying of the RNG state
469 // ourselves.
470
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 RNG rng_{};
471
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (!randomSeed) {
472 #ifndef NO_OMP
473 ✗ if (omp_get_max_threads() > 1) {
474 ✗ std::clog << "\nWARNING: You have specified randomSeed = " << randomSeed
475 ✗ << " while using multithreading (number of OpenMP threads: "
476 << omp_get_max_threads()
477 << "). Consider setting the thread count to 1 (e.g. via the environment variable "
478 ✗ "OMP_NUM_THREADS) to ensure reproducible runs.";
479 }
480 #endif
481 ✗ rng_.engine.seed(0);
482 }
483 2 #pragma omp parallel
484 {
485 // Copy the complete state of the reference RNG
486 rng = rng_;
487 // Set individual streams (which ensure independent, non-overlapping random sequences) for each
488 // thread (feature of PCG)
489 #ifndef NO_OMP
490 rng.engine.set_stream(omp_get_thread_num());
491 #endif
492 }
493
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::clog << "\n\ttype: Langevin Dynamics";
494
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\trandomSeed: " << this->randomSeed;
495 2 }
496
497 2 LangevinSymplecticEuler::LangevinSymplecticEuler(System *s, const YAML::Node &c) : Langevin(s, c) {
498
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::clog << "\n\tsubtype: Langevin Euler";
499 2 }
500
501 2 void LangevinSymplecticEuler::step() {
502 2 #pragma omp parallel for
503 for (auto &p : system->particles) {
504 p.v += p.F() * dt - p.gamma * p.m * p.v * dt +
505 sqrt(2. * system->T * p.gamma * p.m * dt) * rng.vecNormal();
506 p.r += p.v * dt;
507 }
508 2 system->applyPeriodic();
509 2 system->updateState();
510 2 system->t += dt;
511 2 }
512
513 4 MonteCarlo::MonteCarlo(System *s, const YAML::Node &c)
514
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
 8 : Integrator(s), randomSeed(YAML::parse<bool>(c, "randomSeed", true)),
515
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
 8 dx(YAML::parse<double>(c, "dx", 0.1)), nSweep(0), acceptCount(0), acceptRatio(0.),
516
4/8
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
 12 globalUpdateAfterSweep(YAML::parse<bool>(c, "globalUpdateAfterSweep", true)) {
517
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (!system->interaction->allDoE) {
518 ✗ throw std::invalid_argument(
519 "Interactions must implement the energy calculation for use with a Monte Carlo integrator. "
520 ✗ "This is not the case for one of the chosen interactions.");
521 }
522
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if (!randomSeed) {
523 ✗ rng.engine.seed(0);
524 }
525
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 std::clog << "\n\ttype: Monte Carlo";
526
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::clog << "\n\tdx: " << this->dx;
527
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::clog << "\n\tglobalUpdateAfterSweep: " << this->globalUpdateAfterSweep;
528
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 std::clog << "\n\trandomSeed: " << this->randomSeed;
529 4 }
530
531 4 void MonteCarlo::step() {
532 4 acceptCount = 0;
533 4 sweep(); // do one sweep (must also set nSweep and acceptCount)
534 4 acceptRatio = (double)acceptCount / nSweep;
535
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 4 if (globalUpdateAfterSweep) {
536 4 system->updateState();
537 }
538 4 system->t += 1.;
539 4 std::cout << ", acceptance: " << acceptRatio;
540 4 std::cout << ", E: " << system->E;
541 4 }
542
543 2 MonteCarloCanonical::MonteCarloCanonical(System *s, const YAML::Node &c)
544
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
 4 : MonteCarlo(s, c), du(YAML::parse<double>(c, "du", 0)) {
545
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tdu: " << du;
546
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::clog << "\n\tensemble: Canonical";
547 2 }
548
549 2 void MonteCarloCanonical::sweep() {
550 2 nSweep = system->particles.size();
551 // Sweep over all particles
552
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 2 times.
 13 for (auto &p : system->particles) {
553 11 Particle pBefore = p;
554 11 double E0 = system->updateParticleState(p);
555 // Trial move
556 11 p.r += dx * rng.vecNormal();
557
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
 11 if (du > 0) {
558 ✗ p.u += du * rng.vecNormal();
559 ✗ p.u.normalize();
560 }
561 11 system->applyPeriodic(p);
562 11 double E1 = system->updateParticleState(p);
563 11 double dE = E1 - E0;
564
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 9 times.
 11 if (rng.realUniform() > exp(-dE / system->T)) {
565 // Rejected: reset particle
566 2 p = pBefore;
567 2 system->updateParticleMeta(p);
568 } else {
569 // Accepted: update energy
570 9 system->E += dE;
571 9 acceptCount++;
572 }
573 }
574 2 }
575
576 2 MonteCarloGrandCanonical::MonteCarloGrandCanonical(System *s, const YAML::Node &c)
577
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 4 : MonteCarlo(s, c), nSweepMin(YAML::parse<int>(c, "nSweepMin", 10)),
578
4/8
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
 6 inoutProb(YAML::parse<double>(c, "inoutProb", 0.1)), V(system->L.prod()) {
579
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if (!std::isfinite(system->mu)) {
580 ✗ throw std::invalid_argument("You must specify the chemical potential mu in the system "
581 ✗ "properties to use Grand Canonical Monte Carlo.");
582 }
583
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 std::clog << "\n\tensemble: Grand Canonical";
584
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tnSweepMin: " << this->nSweepMin;
585
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 std::clog << "\n\tinoutProb: " << this->inoutProb;
586 2 }
587
588 ✗ double MonteCarloGrandCanonical::weight() { return nSweep / V; }
589
590 86 Particle &MonteCarloGrandCanonical::getRandomParticle() {
591 86 return system->particles[rng.intUniform(0, system->particles.size() - 1)];
592 }
593
594 2 void MonteCarloGrandCanonical::sweep() {
595
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 size_t N = system->particles.size();
596
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 nSweep = std::max(N, nSweepMin);
597
2/2
✓ Branch 0 taken 110 times.
✓ Branch 1 taken 2 times.
 112 for (int i = 0; i < nSweep; i++) {
598
2/2
✓ Branch 1 taken 70 times.
✓ Branch 2 taken 40 times.
 110 if (rng.realUniform() > inoutProb) {
599
1/2
✓ Branch 0 taken 70 times.
✗ Branch 1 not taken.
 70 if (!system->particles.empty()) {
600 70 trialMove();
601 }
602 } else {
603
2/2
✓ Branch 1 taken 24 times.
✓ Branch 2 taken 16 times.
 40 if (rng.realUniform() > 0.5) {
604 24 trialInsert();
605 } else {
606
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if (!system->particles.empty()) {
607 16 trialRemove();
608 }
609 }
610 }
611 }
612 2 std::cout << ", N: " << N;
613 2 }
614
615 70 void MonteCarloGrandCanonical::trialMove() {
616 70 Particle &p = getRandomParticle();
617 70 Particle pBefore = p;
618 70 double E0 = system->updateParticleState(p);
619 // Trial move
620 70 p.r += dx * rng.vecNormal();
621 70 system->applyPeriodic(p);
622 70 double E1 = system->updateParticleState(p);
623 70 double dE = E1 - E0;
624
2/2
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 68 times.
 70 if (rng.realUniform() > exp(-dE / system->T)) {
625 // Rejected: reset particle
626 2 p = pBefore;
627 2 system->updateParticleMeta(p);
628 } else {
629 // Accepted: update energy
630 68 system->E += dE;
631 68 acceptCount++;
632 }
633 70 }
634
635 24 void MonteCarloGrandCanonical::trialInsert() {
636 24 system->particles.emplace_back(rng.vecUniform(system->L));
637 24 Particle &p = system->particles.back();
638 24 double dE = system->updateParticleState(p, ParticleUpdateType::INSERT);
639
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 24 times.
 24 if (rng.realUniform() > V / system->particles.size() * exp((-dE + system->mu) / system->T)) {
640 // Rejected: particle must be removed again, so update meta info and really pop the element
641 ✗ system->updateParticleMeta(p, ParticleUpdateType::DELETE);
642 ✗ system->particles.pop_back();
643 } else {
644 // Accepted: change energy
645 24 system->E += dE;
646 24 acceptCount++;
647 }
648 24 }
649
650 16 void MonteCarloGrandCanonical::trialRemove() {
651 16 Particle &p = getRandomParticle(); // this is particle A
652 16 double dE = system->updateParticleState(p);
653 // Note the flipped sign in the Bolzmann factor because we remove a particle
654
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if (rng.realUniform() > system->particles.size() / V * exp((dE - system->mu) / system->T)) {
655 // Rejected: nothing to do
656 } else {
657 // Accepted
658 // Remove meta info of last particle
659 ✗ system->updateParticleMeta(system->particles.back(), ParticleUpdateType::DELETE);
660 // If A is the last particle, we can just pop the last particle element and update the energy
661 // (see below).
662 // If A is not the last particle ...
663 ✗ if (&p != &system->particles.back()) {
664 // ... also remove meta info of A,
665 ✗ system->updateParticleMeta(p, ParticleUpdateType::DELETE);
666 // swap with B,
667 ✗ std::swap(p, system->particles.back());
668 // and restore the meta info of B, which now sits at the previous position of A
669 ✗ system->updateParticleMeta(p, ParticleUpdateType::INSERT);
670 }
671 // Finally pop the last particle (which at this moment is guaranteed to be A) and update energy
672 ✗ system->particles.pop_back();
673 ✗ system->E -= dE;
674 ✗ acceptCount++;
675 }
676 16 }
677