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 |