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 << ", E: " << system->E; | |
540 | 4 | } | |
541 | |||
542 | 2 | MonteCarloCanonical::MonteCarloCanonical(System *s, const YAML::Node &c) | |
543 |
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)) { |
544 |
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; |
545 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | std::clog << "\n\tensemble: Canonical"; |
546 | 2 | } | |
547 | |||
548 | 2 | void MonteCarloCanonical::sweep() { | |
549 | 2 | nSweep = system->particles.size(); | |
550 | // Sweep over all particles | ||
551 |
2/2✓ Branch 0 taken 11 times.
✓ Branch 1 taken 2 times.
|
13 | for (auto &p : system->particles) { |
552 | 11 | Particle pBefore = p; | |
553 | 11 | double E0 = system->updateParticleState(p); | |
554 | // Trial move | ||
555 | 11 | p.r += dx * rng.vecNormal(); | |
556 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 11 times.
|
11 | if (du > 0) { |
557 | ✗ | p.u += du * rng.vecNormal(); | |
558 | ✗ | p.u.normalize(); | |
559 | } | ||
560 | 11 | system->applyPeriodic(p); | |
561 | 11 | double E1 = system->updateParticleState(p); | |
562 | 11 | double dE = E1 - E0; | |
563 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 9 times.
|
11 | if (rng.realUniform() > exp(-dE / system->T)) { |
564 | // Rejected: reset particle | ||
565 | 2 | p = pBefore; | |
566 | 2 | system->updateParticleMeta(p); | |
567 | } else { | ||
568 | // Accepted: update energy | ||
569 | 9 | system->E += dE; | |
570 | 9 | acceptCount++; | |
571 | } | ||
572 | } | ||
573 | 2 | } | |
574 | |||
575 | 2 | MonteCarloGrandCanonical::MonteCarloGrandCanonical(System *s, const YAML::Node &c) | |
576 |
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)), |
577 |
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()) { |
578 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | if (!std::isfinite(system->mu)) { |
579 | ✗ | throw std::invalid_argument("You must specify the chemical potential mu in the system " | |
580 | ✗ | "properties to use Grand Canonical Monte Carlo."); | |
581 | } | ||
582 |
1/2✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
|
2 | std::clog << "\n\tensemble: Grand Canonical"; |
583 |
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; |
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\tinoutProb: " << this->inoutProb; |
585 | 2 | } | |
586 | |||
587 | ✗ | double MonteCarloGrandCanonical::weight() { return nSweep / V; } | |
588 | |||
589 | 81 | Particle &MonteCarloGrandCanonical::getRandomParticle() { | |
590 | 81 | return system->particles[rng.intUniform(0, system->particles.size() - 1)]; | |
591 | } | ||
592 | |||
593 | 2 | void MonteCarloGrandCanonical::sweep() { | |
594 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | size_t N = system->particles.size(); |
595 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
|
2 | nSweep = std::max(N, nSweepMin); |
596 |
2/2✓ Branch 0 taken 110 times.
✓ Branch 1 taken 2 times.
|
112 | for (int i = 0; i < nSweep; i++) { |
597 |
2/2✓ Branch 1 taken 55 times.
✓ Branch 2 taken 55 times.
|
110 | if (rng.realUniform() > inoutProb) { |
598 |
1/2✓ Branch 0 taken 55 times.
✗ Branch 1 not taken.
|
55 | if (!system->particles.empty()) { |
599 | 55 | trialMove(); | |
600 | } | ||
601 | } else { | ||
602 |
2/2✓ Branch 1 taken 29 times.
✓ Branch 2 taken 26 times.
|
55 | if (rng.realUniform() > 0.5) { |
603 | 29 | trialInsert(); | |
604 | } else { | ||
605 |
1/2✓ Branch 0 taken 26 times.
✗ Branch 1 not taken.
|
26 | if (!system->particles.empty()) { |
606 | 26 | trialRemove(); | |
607 | } | ||
608 | } | ||
609 | } | ||
610 | } | ||
611 | 2 | std::cout << ", N: " << N; | |
612 | 2 | } | |
613 | |||
614 | 55 | void MonteCarloGrandCanonical::trialMove() { | |
615 | 55 | Particle &p = getRandomParticle(); | |
616 | 55 | Particle pBefore = p; | |
617 | 55 | double E0 = system->updateParticleState(p); | |
618 | // Trial move | ||
619 | 55 | p.r += dx * rng.vecNormal(); | |
620 | 55 | system->applyPeriodic(p); | |
621 | 55 | double E1 = system->updateParticleState(p); | |
622 | 55 | double dE = E1 - E0; | |
623 |
2/2✓ Branch 1 taken 2 times.
✓ Branch 2 taken 53 times.
|
55 | if (rng.realUniform() > exp(-dE / system->T)) { |
624 | // Rejected: reset particle | ||
625 | 2 | p = pBefore; | |
626 | 2 | system->updateParticleMeta(p); | |
627 | } else { | ||
628 | // Accepted: update energy | ||
629 | 53 | system->E += dE; | |
630 | 53 | acceptCount++; | |
631 | } | ||
632 | 55 | } | |
633 | |||
634 | 29 | void MonteCarloGrandCanonical::trialInsert() { | |
635 | 29 | system->particles.emplace_back(rng.vecUniform(system->L)); | |
636 | 29 | Particle &p = system->particles.back(); | |
637 | 29 | double dE = system->updateParticleState(p, ParticleUpdateType::INSERT); | |
638 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 29 times.
|
29 | if (rng.realUniform() > V / system->particles.size() * exp((-dE + system->mu) / system->T)) { |
639 | // Rejected: particle must be removed again, so update meta info and really pop the element | ||
640 | ✗ | system->updateParticleMeta(p, ParticleUpdateType::DELETE); | |
641 | ✗ | system->particles.pop_back(); | |
642 | } else { | ||
643 | // Accepted: change energy | ||
644 | 29 | system->E += dE; | |
645 | 29 | acceptCount++; | |
646 | } | ||
647 | 29 | } | |
648 | |||
649 | 26 | void MonteCarloGrandCanonical::trialRemove() { | |
650 | 26 | Particle &p = getRandomParticle(); // this is particle A | |
651 | 26 | double dE = system->updateParticleState(p); | |
652 | // Note the flipped sign in the Bolzmann factor because we remove a particle | ||
653 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 26 times.
|
26 | if (rng.realUniform() > system->particles.size() / V * exp((dE - system->mu) / system->T)) { |
654 | // Rejected: nothing to do | ||
655 | } else { | ||
656 | // Accepted | ||
657 | // Remove meta info of last particle | ||
658 | ✗ | system->updateParticleMeta(system->particles.back(), ParticleUpdateType::DELETE); | |
659 | // If A is the last particle, we can just pop the last particle element and update the energy | ||
660 | // (see below). | ||
661 | // If A is not the last particle ... | ||
662 | ✗ | if (&p != &system->particles.back()) { | |
663 | // ... also remove meta info of A, | ||
664 | ✗ | system->updateParticleMeta(p, ParticleUpdateType::DELETE); | |
665 | // swap with B, | ||
666 | ✗ | std::swap(p, system->particles.back()); | |
667 | // and restore the meta info of B, which now sits at the previous position of A | ||
668 | ✗ | system->updateParticleMeta(p, ParticleUpdateType::INSERT); | |
669 | } | ||
670 | // Finally pop the last particle (which at this moment is guaranteed to be A) and update energy | ||
671 | ✗ | system->particles.pop_back(); | |
672 | ✗ | system->E -= dE; | |
673 | ✗ | acceptCount++; | |
674 | } | ||
675 | 26 | } | |
676 |