GCC Code Coverage Report

Directory: ./
File: include/integrator.hpp
Date: 2025-07-07 11:13:23
Exec Total Coverage
Lines: 5 8 62.5%
Functions: 2 2 100.0%
Branches: 2 4 50.0%
Line Branch Exec Source
1 #pragma once
2
3 #include <stack>
4 #include <yaml-cpp/yaml.h>
5
6 #include "rng.hpp"
7
8 class System;
9 struct Particle;
10
11 /**
12 * Abstract base class for all sorts of integrators. Those are time evolution or phase space
13 * sampling methods which are used to propagate the particles in a system.
14 * All integrators have to define the methods step() and weight().
15 */
16 class Integrator {
17 public:
18 System *system;
19
20 Integrator(System *s);
21 6 virtual ~Integrator() = default;
22
23 /**
24 * Do one step, i.e. propagate the particles of the system and their properties according to the
25 * respective method of the chosen integrator class.
26 */
27 virtual void step() = 0;
28
29 /**
30 * Return a weight that measures the relative importance of the current step to yield correct
31 * sampling of observables. The returned value is 1 for simple methods, but can be any real number
32 * e.g. in the case of variable timestepping algorithms.
33 * @return sampling weight of current step
34 */
35 virtual double weight();
36 };
37
38 /**
39 * An integrator that does nothing. Good for debugging other stuff and viewing a certain particle
40 * configuration in the graphics output.
41 */
42 class Static : public Integrator {
43 public:
44 Static(System *s, const YAML::Node &c);
45
46 void step();
47 };
48
49 /**
50 * Abstract base class for Brownian Dynamics simulations.
51 */
52 class Brownian : public Integrator {
53 public:
54 double dt;
55 bool randomSeed;
56 static RNG rng;
57 #pragma omp threadprivate(rng)
58
59 Brownian(System *s, const YAML::Node &c);
60 6 virtual ~Brownian() = default;
61
62 virtual void step() = 0;
63 double weight();
64 };
65
66 /**
67 * Brownian Dynamics via the simple Euler-Maruyama method and fixed timestepping.
68 */
69 class BrownianEuler : public Brownian {
70 public:
71 BrownianEuler(System *s, const YAML::Node &c);
72 4 virtual ~BrownianEuler() = default;
73
74 void step();
75 };
76
77 /**
78 * Abstract base class for Brownian Dynamics simulations with the embedded Heun-Euler method and a
79 * version of RSwM.
80 */
81 class BrownianHeunEulerRSwM : public Brownian {
82 public:
83
2/4
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 12 taken 3 times.
✗ Branch 13 not taken.
 72 using RandomIncrement = struct {
84 double dt;
85 VecVec R;
86 };
87 using RandomIncrementStack = std::stack<RandomIncrement, std::vector<RandomIncrement>>;
88
89 double dtTrial;
90 double relTol;
91 double absTol;
92 double qmin;
93 double qmax;
94 double alpha;
95 double beta;
96 int lNorm;
97
98 // Discard random increments beyond a cutoff threshold. This avoids storing very small
99 // random increments (and results of mere floating point rounding errors) on the
100 // stack(s). If chosen small enough, it should not alter the random process
101 double dtCutoff = 1e-13;
102
103 double error;
104 double q;
105 int rejections;
106
107 VecVec R;
108 VecVec Rbridge;
109 VecVec Rdiff;
110 VecVec Rtmp;
111
112 std::vector<Particle> particlesBefore;
113 VecVec drs0;
114 VecVec drs1;
115
116 BrownianHeunEulerRSwM(System *s, const YAML::Node &c);
117 8 virtual ~BrownianHeunEulerRSwM() = default;
118
119 void generateRandVecs(VecVec &buffer, Real stddev);
120 void generateRandVecsBridge(VecVec &buffer, const VecVec &increment, double q, double dt);
121 double estimateParticleError(int i);
122 double estimateTotalError();
123 void maybePush(RandomIncrementStack &S, const RandomIncrement &RI);
124
125 bool adaptStepsize();
126 void step();
127 virtual void rejectStep() = 0;
128 virtual void prepareNextStep() = 0;
129 };
130
131 /**
132 * Brownian Dynamics via the embedded Heun-Euler method and RSwM2 adaptive timestepping.
133 */
134 class BrownianHeunEulerRSwM2 : public BrownianHeunEulerRSwM {
135 public:
136 RandomIncrementStack Sf;
137
138 BrownianHeunEulerRSwM2(System *s, const YAML::Node &c);
139
140 void rejectStep();
141 void prepareNextStep();
142 };
143
144 /**
145 * Brownian Dynamics via the embedded Heun-Euler method and RSwM3 adaptive timestepping.
146 */
147 class BrownianHeunEulerRSwM3 : public BrownianHeunEulerRSwM {
148 public:
149 RandomIncrementStack Sf;
150 RandomIncrementStack Su;
151
152 BrownianHeunEulerRSwM3(System *s, const YAML::Node &c);
153
154 void rejectStep();
155 void prepareNextStep();
156 };
157
158 /**
159 * Abstract base class for Molecular Dynamics simulations.
160 */
161 class Molecular : public Integrator {
162 public:
163 double dt;
164
165 Molecular(System *s, const YAML::Node &c);
166 ✗ virtual ~Molecular() = default;
167
168 virtual void step() = 0;
169 };
170
171 /**
172 * Molecular Dynamics via the Velocity Verlet method.
173 */
174 class MolecularVelocityVerlet : public Molecular {
175 public:
176 MolecularVelocityVerlet(System *s, const YAML::Node &c);
177
178 void step();
179 };
180
181 /**
182 * Molecular Dynamics in the canonical ensemble via the Nose-Hoover thermostat and the Velocity
183 * Verlet method.
184 */
185 class MolecularNoseHooverVelocityVerlet : public Molecular {
186 void chain();
187 void posvel();
188
189 public:
190 double dt2;
191 double dt4;
192 double dt8;
193 double Ekin;
194 double Q;
195 double xi1;
196 double vxi1;
197 double xi2;
198 double vxi2;
199
200 MolecularNoseHooverVelocityVerlet(System *s, const YAML::Node &c);
201
202 void step();
203 };
204
205 /**
206 * Abstract base class for Langevin Dynamics simulations.
207 */
208 class Langevin : public Integrator {
209 public:
210 double dt;
211 bool randomSeed;
212 static RNG rng;
213 #pragma omp threadprivate(rng)
214
215 Langevin(System *s, const YAML::Node &c);
216 ✗ virtual ~Langevin() = default;
217
218 virtual void step() = 0;
219 };
220
221 /**
222 * Langevin Dynamics with the symplectic (or semi-implicit) Euler method
223 */
224 class LangevinSymplecticEuler : public Langevin {
225 public:
226 LangevinSymplecticEuler(System *s, const YAML::Node &c);
227
228 void step();
229 };
230
231 /**
232 * Monte Carlo base class.
233 * Note that MC simulations do not produce physical trajectories and that the simulation time has no
234 * further meaning (it is set equal to the number of steps for convenience).
235 */
236 class MonteCarlo : public Integrator {
237 public:
238 bool randomSeed;
239 Real dx;
240 int nSweep;
241 int acceptCount;
242 double acceptRatio;
243 bool globalUpdateAfterSweep;
244 RNG rng;
245
246 MonteCarlo(System *s, const YAML::Node &c);
247 ✗ virtual ~MonteCarlo() = default;
248
249 void step();
250 virtual void sweep() = 0;
251 };
252
253 /**
254 * Canonical (N, V, T) Monte Carlo via standard importance sampling.
255 */
256 class MonteCarloCanonical : public MonteCarlo {
257 public:
258 Real du;
259
260 MonteCarloCanonical(System *s, const YAML::Node &c);
261
262 void sweep();
263 };
264
265 /**
266 * Grand canonical (mu, V, T) Monte Carlo.
267 */
268 class MonteCarloGrandCanonical : public MonteCarlo {
269 public:
270 size_t nSweepMin;
271 double inoutProb;
272 double V;
273
274 MonteCarloGrandCanonical(System *s, const YAML::Node &c);
275 double weight();
276
277 Particle &getRandomParticle();
278 void sweep();
279 void trialMove();
280 void trialInsert();
281 void trialRemove();
282 };
283
284 class IntegratorFactory {
285 public:
286 static std::map<std::string, Integrator *(*)(System *s, const YAML::Node &c)> integratorMap;
287 static Integrator *create(System *system, const YAML::Node &config);
288 };
289