-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathsearch.html
More file actions
397 lines (270 loc) · 210 KB
/
search.html
File metadata and controls
397 lines (270 loc) · 210 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8">
<meta http-equiv="x-ua-compatible" content="ie=edge">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>Search results</title>
<link rel="stylesheet" href="/assets/css/readdy_documentation.css">
<link rel="canonical" href="https://readdy.github.io/search.html">
<link href="https://fonts.googleapis.com/css?family=Inconsolata|Roboto+Mono|Lora|Lato|Source+Sans+Pro|Roboto+Slab|Merriweather" rel="stylesheet">
<link rel="stylesheet" href="https://readdy.github.io/libraries/perfect-scrollbar/css/perfect-scrollbar.min.css"/>
<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.10.1/dist/katex.min.css" integrity="sha384-dbVIfZGuN1Yq7/1Ocstc1lUEm+AT+/rCkibIcC/OmWo5f0EA48Vf8CytHzGrSwbQ" crossorigin="anonymous">
<script defer src="https://cdn.jsdelivr.net/npm/katex@0.10.1/dist/katex.min.js" integrity="sha384-2BKqo+exmr9su6dir+qCw08N2ZKRucY4PrGQPPWU1A7FtlCGjmEGFqXCv5nyM5Ij" crossorigin="anonymous"></script>
<script defer src="https://cdn.jsdelivr.net/npm/katex@0.10.1/dist/contrib/auto-render.min.js"></script>
<script>
document.addEventListener("DOMContentLoaded", function() {
renderMathInElement(document.body, {
delimiters: [
{left: "$$", right: "$$", display: true},
{left: "$", right: "$", display: false},
{left: "\\[", right: "\\]", display: true},
{left: "\\(", right: "\\)", display: false},
]
});
});
</script>
<link rel="icon" type="image/png" href="/assets/icon_black_32px.png">
</head>
<body>
<div class="side-wrapper" id="unique-side-container">
<div class="side">
<div class="side-logo logo-readdy"><a href="https://readdy.github.io/index.html"></a></div>
<div style="margin-right: 1.5rem; text-align: center;">
<a href="https://readdy.github.io/index.html">ReaDDy - A particle-based<br>reaction-diffusion simulator</a>
</div>
<div class="side-searchbar-wrapper">
<form action="https://readdy.github.io/search.html" method="get">
<input type="text" id="search-box" name="query" placeholder="Search ...">
</form>
</div>
<nav class="side-nav">
<a class="side-nav-item" href="https://readdy.github.io/index.html">Home</a>
<a class="side-nav-item" href="https://readdy.github.io/installation.html">Install readdy</a>
<div class="nav-supergroup-delimiter">
<b>API</b>
</div>
<a class="side-nav-item" href="https://readdy.github.io/system.html">System configuration</a>
<a class="side-nav-item" href="https://readdy.github.io/simulation.html">Simulation</a>
<a class="side-nav-item" href="https://readdy.github.io/results.html">Post-processing</a>
<div class="nav-supergroup-delimiter">
<b>Examples</b>
</div>
<a class="side-nav-item" href="https://readdy.github.io/demonstration.html">Demonstration</a>
<a class="side-nav-item" href="https://readdy.github.io/validation.html">Validation</a>
<a class="side-nav-item" href="https://readdy.github.io/benchmark.html">Benchmark</a>
<div class="nav-supergroup-delimiter">
<b>Advanced topics</b>
</div>
<a class="side-nav-item" href="https://readdy.github.io/cookbook.html">Cookbook</a>
<a class="side-nav-item" href="https://readdy.github.io/development.html">Development</a>
<div class="nav-supergroup-delimiter">
<b>Legal notices</b>
</div>
<a class="side-nav-item" href="https://readdy.github.io/imprint.html">Imprint</a>
<a class="side-nav-item" href="https://readdy.github.io/software_license.html">Software license</a>
</nav>
<div class="github-wrapper">
<a href="https://github.com/readdy/readdy" style="width: 16rem; height:100%; position:absolute; top:0; left:0;">
<div class="github-text">ReaDDy on GitHub</div>
<div class="side-logo logo-github"></div>
</a>
</div>
<div class="side-logo logo-cmb"><a href="https://www.mi.fu-berlin.de/en/math/groups/comp-mol-bio/"></a></div>
</div>
</div>
<div class="main-container" id="unique-main-container">
<div class="main">
<article>
<h1 style="font-size: 2.7rem; padding-left: 0.5rem;">Search results</h1>
<form action="https://readdy.github.io/search.html" method="get">
<input type="text" id="search-box" name="query">
<input type="submit" value="search">
</form>
<ul id="search-results"></ul>
<script>
window.store = {
"benchmark-cytosolic-reactions": {
"title": "Benchmark > Cytosolic reactions",
"content": "We consider the system of particle species A, B and C that undergo reactions $A + B\\rightleftharpoons C$ and also have optional volume exclusion effects due to repulsive potentials between all particles. A and B Particles will be initially uniformly distributed in a periodic box and no C particles are present.We first demonstrate how the concentration of particles evolves in time and compares to a description by the law of mass action (LMA). This demonstrates that ReaDDy produces exact results when parameters are in the well-mixed regime. We also show that including a repulsive potential between all particles, but not changing the intrinsic reaction rates, will lead to different macroscopic reaction kinetics.The focus of this demonstration is the measurement of computation time. In particular we observe the computation time $\\tau$ per particle and per integration step. The system is simulated for different initial particle numbers, such that the density is constant.In [1]: import osimport numpy as npimport readdyut = readdy.units In [2]: readdy.__version__ Out[2]:'v2.0.3-27'Add species A, B and C to the system. Diffusion constants are chosen to represent the size of particles A, B and C according to the Einstein-Stokes relation$$D_i = \\frac{k_B T}{6\\pi\\eta r_i}$$where $D_i$ is the diffusion coefficient of species $i$ and $r_i$ is its radius. $\\eta$ is the viscosity of a cellular environment, i.e. mostly water. $k_B T$ is the termal energy, here we choose it to represent a system at room temperature.Note that in ReaDDy, the concept of a particle radius does not exist. The volume of a particle only becomes important when setting the reaction radius of a reaction or when setting the interaction distance of a interaction potential. However in modeling a certain system, we have particles with certain sizes in mind. To implement these, one can set the radii of particles outside of ReaDDy and use them to choose the reaction radii and interaction distances.In [3]: particle_radii = {"A": 1.5, "B": 3., "C": 3.12} # in nanometersdt = 1e-1 * ut.nanosecond In [4]: def get_system(with_repulsion): system = readdy.ReactionDiffusionSystem( box_size=[60., 60., 60.], temperature=293., unit_system={"length_unit": "nanometer", "time_unit": "nanosecond", "energy_unit": "kilojoule / mol"} ) system.add_species("A", diffusion_constant=143.1 * ut.micrometer ** 2 / ut.second) system.add_species("B", diffusion_constant=71.6 * ut.micrometer ** 2 / ut.second) system.add_species("C", diffusion_constant=68.82 * ut.micrometer ** 2 / ut.second) if with_repulsion: force_constant = 10. * ut.kilojoule / ut.mol / (ut.nanometer ** 2) else: force_constant = 0. * ut.kilojoule / ut.mol / (ut.nanometer ** 2) if force_constant.magnitude > 0.: all_pairs = [("A","A"), ("B","B"), ("C", "C"), ("A", "B"), ("A", "C"), ("B", "C")] for pair in all_pairs: distance = particle_radii[pair[0]] + particle_radii[pair[1]] system.potentials.add_harmonic_repulsion(pair[0], pair[1], force_constant, interaction_distance=distance) return system Particles can undergo the two reactions $A + B \\rightarrow C$ and $C \\rightarrow A + B$. Note that the particles radii are chosen, such that the total volume of particles is conserved in this reaction, i.e.$$ r_A^3 + r_B^3 = r_C^3 $$The reaction radius $R$ is chosen as $R=r_A + r_B$. This is used as the maximum distance of educts in the Fusion reaction as well as the distance at which products will be placed in the Fission reaction.In [5]: reaction_radius = particle_radii["A"] + particle_radii["B"]def configure_reactions(system): system.reactions.add_fusion("fusion", "A", "B", "C", rate=1e6 / ut.second, educt_distance=reaction_radius * ut.nanometer) system.reactions.add_fission("fission", "C", "A", "B", rate=5e4 / ut.second, product_distance=reaction_radius * ut.nanometer) In [6]: number_a = 509number_b = 509def simulate(system, output_file): simulation = system.simulation(kernel="SingleCPU") simulation.output_file = output_file simulation.reaction_handler = "UncontrolledApproximation" edge_length = system.box_size[0] initial_positions_a = np.random.random(size=(number_a, 3)) * edge_length - .5 * edge_length initial_positions_b = np.random.random(size=(number_b, 3)) * edge_length - .5 * edge_length simulation.add_particles("A", initial_positions_a) simulation.add_particles("B", initial_positions_b) simulation.observe.number_of_particles(stride=1, types=["A", "B", "C"]) if os.path.exists(simulation.output_file): os.remove(simulation.output_file) n_steps = int(10000. / dt.magnitude) # simulate to 10 microseconds simulation.run(n_steps=n_steps, timestep=dt) In [7]: system_with_repulsion = get_system(True)configure_reactions(system_with_repulsion)system_without_repulsion = get_system(False)configure_reactions(system_without_repulsion) In [8]: outfile_with_repulsion = "cytosolic_reactions_with_repulsion.h5"simulate(system_with_repulsion, outfile_with_repulsion) 0%| | 5/10000 [00:00<03:40, 45.38it/s] Configured kernel context with:-------------------------------- - kBT = 2.43614 - periodic b.c. = (true, true, true) - box size = (60.0, 60.0, 60.0) - particle types: * particle type "C" with D=0.06882 * particle type "A" with D=0.1431 * particle type "B" with D=0.0716 - potentials of order 2: * for types "B" and "C" * Harmonic repulsion with Force constant k=10.0 * for types "B" and "B" * Harmonic repulsion with Force constant k=10.0 * for types "A" and "C" * Harmonic repulsion with Force constant k=10.0 * for types "A" and "A" * Harmonic repulsion with Force constant k=10.0 * for types "C" and "C" * Harmonic repulsion with Force constant k=10.0 * for types "A" and "B" * Harmonic repulsion with Force constant k=10.0 - unimolecular reactions: * Fission C -> A + B with a rate of 5.0e-05, a product distance of 4.5, and weights 0.5 and 0.5 - bimolecular reactions: * Fusion A + B -> C with a rate of 0.00100000, an educt distance of 4.5, and weights 0.5 and 0.5Configured simulation loop with:-------------------------------- - timeStep = 0.1 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 10000/10000 [01:53<00:00, 87.99it/s]In [9]: outfile_without_repulsion = "cytosolic_reactions_without_repulsion.h5"simulate(system_without_repulsion, outfile_without_repulsion) 0%| | 9/10000 [00:00<01:56, 85.84it/s] Configured kernel context with:-------------------------------- - kBT = 2.43614 - periodic b.c. = (true, true, true) - box size = (60.0, 60.0, 60.0) - particle types: * particle type "C" with D=0.06882 * particle type "A" with D=0.1431 * particle type "B" with D=0.0716 - unimolecular reactions: * Fission C -> A + B with a rate of 5.0e-05, a product distance of 4.5, and weights 0.5 and 0.5 - bimolecular reactions: * Fusion A + B -> C with a rate of 0.00100000, an educt distance of 4.5, and weights 0.5 and 0.5Configured simulation loop with:-------------------------------- - timeStep = 0.1 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 10000/10000 [01:04<00:00, 156.14it/s]Obtain outputIn [10]: c = dict()traj_with_repulsion = readdy.Trajectory(outfile_with_repulsion)times, counts = traj_with_repulsion.read_observable_number_of_particles()concentrations = counts / system_with_repulsion.box_volume.magnitudec["w/ repulsion"] = {"A": concentrations[:, 0], "B": concentrations[:, 1], "C": concentrations[:, 2]}traj_without_repulsion = readdy.Trajectory(outfile_without_repulsion)times2, counts2 = traj_without_repulsion.read_observable_number_of_particles()concentrations2 = counts2 / system_without_repulsion.box_volume.magnitudec["w/o repulsion"] = {"A": concentrations2[:, 0], "B": concentrations2[:, 1], "C": concentrations2[:, 2]}times = times * dt.magnitude Under the assumption that the system is well mixed (particles encounter often, but rarely react), and there are no interaction potentials, that would obstruct particles in their free diffusion, one can approximate the macroscopic association rate $k_\\mathrm{on}$ as$$k_\\mathrm{on} = \\lambda_\\mathrm{on} V_R $$where $\\lambda_\\mathrm{on}$ is the intrinsic association rate and $V_R$ is the reaction volume available to a pair of particles A and B.In [11]: from scipy.integrate import odeintlambda_on = 1e6 * 1e-9 # per nanosecondreaction_volume = 4. / 3. * np.pi * reaction_radius**3k_on = lambda_on * reaction_volume # dilute (large volume) and well mixed (fast diffusion) approximationk_off = 5e4 * 1e-9 # per nanoseconddef f(x, _): result = np.zeros_like(x) delta = - k_on * x[0] * x[0] + k_off * x[1] result[0] = delta result[1] = -1. * delta return resultx_init = np.array([number_a, 0.]) / system_with_repulsion.box_volume.magnitudeode_times = np.arange(0., 10000. , 1.) # in nanosecondsode_solution = odeint(f, x_init, t=ode_times) In [12]: %matplotlib inlineimport matplotlib.pyplot as pltplt.plot(times, c["w/o repulsion"]["A"], label="[A]=[B]")plt.plot(times, c["w/o repulsion"]["C"], label="[C]")plt.plot(times, c["w/ repulsion"]["A"], label="[A]=[B], w/ repulsion")plt.plot(times, c["w/ repulsion"]["C"], label="[C], w/ repulsion")plt.plot(ode_times, ode_solution[:, 0], "k--", label="law of mass action")plt.plot(ode_times, ode_solution[:, 1], "k--")plt.legend(loc="best")plt.xlabel(r"Time in $\\mathrm{ns}$")plt.ylabel(r"Concentration in $\\mathrm{nm}^{-3}$")plt.show() PerformanceThe system shown above is simulated for constant densities but different number of particles. For the measurements we always have particle repulsion present. The computation time is measured for the different compute-kernels SingleCPU and CPU. The most important value is the computation time $\\tau$ per particle and per integration step. This value is measured at the highest number of particles as it usually saturates when load-balancing improves. An exception is the legacy Java version, which was measured at 1500 particles since it does not scale well. The measurements were performed on a machine equipped with a Intel(R) Xeon(R) CPU E5-1630 v4 @ 3.70GHz that has 8 logical CPUs, all of which were used for profiling the CPU compute-kernel.$\\tau$ - computation time per particle and per integration stepReaDDy version$\\tau$ SingleCPU$\\tau$ CPUjava legacy(@1500 particles)3.77806e-06n.a.v1.0.0-py36_7_g50d9ee52.11896e-06 s8.21963e-07 sv1.0.0-py36_59_ga47c00481.43479e-06 s6.02905e-07 sv1.0.11.23611e-06 s5.13914e-07 s ",
"url": "/benchmark.html#"
},
"benchmark-michaelis-menten": {
"title": "Benchmark > Michaelis-Menten - well mixed",
"content": "In [1]: %matplotlib inlineimport osimport numpy as npimport readdy In [2]: readdy.__version__ Out[2]:'v2.0.3-27'Setup ReaDDy systemThe intrinsic forward rate $\\lambda_\\mathrm{f}$ is obtained from the macroscopic forward rate $k_\\mathrm{f} = 5.88 \\,\\mathrm{\\mu M}^{-1} \\mathrm{s}^{-1} = 0.98 \\times 10^{-2} \\mathrm{\\mu m}^{3} \\mathrm{s}^{-1}$ by application of$$k_\\mathrm{f} = 4\\pi (D_E + D_S) R \\left[1 - \\frac{\\tanh(\\kappa R)}{\\kappa R}\\right]$$where$$\\kappa = \\sqrt{\\frac{\\lambda_\\mathrm{f}}{D_E + D_S}}$$as described in [1]. Note that this will only work in the well mixed regime. The equation cannot be solved explicitly for $\\kappa$. Hence, the intrinsic rate was numerically approximated by solving a minimization problem.[1] R. Erban and J. Chapman, “Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions.,” Phys. Biol., vol. 6, no. 4, p. 46001, Jan. 2009.In [3]: system = readdy.ReactionDiffusionSystem( box_size=[0.3, 0.3, 0.3], unit_system={"length_unit": "micrometer", "time_unit": "second"})system.add_species("E", diffusion_constant=10.)system.add_species("S", diffusion_constant=10.)system.add_species("ES", diffusion_constant=10.)system.add_species("P", diffusion_constant=10.)system.reactions.add("fwd: E +(0.03) S -> ES", rate=86.78638438)system.reactions.add("back: ES -> E +(0.03) S", rate=1.)system.reactions.add("prod: ES -> E +(0.03) P", rate=1.) Simulate the systemIn [4]: simulation = system.simulation(kernel="CPU")simulation.output_file = "out.h5"simulation.reaction_handler = "UncontrolledApproximation"n_particles_e = 909n_particles_s = 9091edge_length = system.box_size[0]initial_positions_e = np.random.random(size=(n_particles_e, 3)) * edge_length - .5*edge_lengthinitial_positions_s = np.random.random(size=(n_particles_s, 3)) * edge_length - .5*edge_lengthsimulation.add_particles("E", initial_positions_e)simulation.add_particles("S", initial_positions_s)simulation.observe.number_of_particles(stride=1, types=["E", "S", "ES", "P"]) In [5]: if os.path.exists(simulation.output_file): os.remove(simulation.output_file)dt = 8e-4n_steps = int(10. / dt)simulation.run(n_steps=n_steps, timestep=dt) 0%| | 2/1250 [00:00<01:20, 15.59it/s] Configured kernel context with:-------------------------------- - kBT = 2.43614 - periodic b.c. = (true, true, true) - box size = (0.3, 0.3, 0.3) - particle types: * particle type "P" with D=10.0 * particle type "ES" with D=10.0 * particle type "E" with D=10.0 * particle type "S" with D=10.0 - unimolecular reactions: * Fission ES -> E + S with a rate of 1.0, a product distance of 0.03, and weights 0.5 and 0.5 * Fission ES -> E + P with a rate of 1.0, a product distance of 0.03, and weights 0.5 and 0.5 - bimolecular reactions: * Fusion E + S -> ES with a rate of 86.7864, an educt distance of 0.03, and weights 0.5 and 0.5Configured simulation loop with:-------------------------------- - timeStep = 0.000800000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 1250/1250 [01:19<00:00, 15.64it/s]Obtain simulation outputIn [6]: traj = readdy.Trajectory(simulation.output_file)time, counts = traj.read_observable_number_of_particles()counts = counts / system.box_volume.magnitudetime = time * dt Analytical solution (integrated ODE)In [7]: from scipy.integrate import odeintdef f(x, t0, kf, kr, kcat): """ x: state vector with concentrations of E, S, ES, P """ return np.array([ -kf * x[0] * x[1] + (kr+kcat)*x[2], -kf * x[0] * x[1] + kr * x[2], kf * x[0] * x[1]- (kr+kcat)*x[2], kcat*x[2] ])init_state = np.array([n_particles_e, n_particles_s, 0., 0.]) / system.box_volume.magnitudeode_time = np.linspace(0.,10,100000)ode_result = odeint(f, y0=init_state, t=ode_time, args=(0.98e-2, 1., 1.)) In [8]: import matplotlib.pyplot as pltstride = 100plt.plot(time[::stride], counts[:,0][::stride], "-", label='E')plt.plot(time[::stride], counts[:,1][::stride], "-", label='S')plt.plot(time[::stride], counts[:,2][::stride], "-", label='ES')plt.plot(time[::stride], counts[:,3][::stride], "-", label='P')plt.plot(ode_time, ode_result[:,0], "--", color="grey", label="LMA solution")plt.plot(ode_time, ode_result[:,1], "--", color="grey")plt.plot(ode_time, ode_result[:,2], "--", color="grey")plt.plot(ode_time, ode_result[:,3], "--", color="grey")plt.legend()plt.xlabel(r"time in $\\mathrm{s}$")plt.ylabel(r"concentration in $\\mathrm{\\mu m}^{-3}$")plt.show() Performance dataIn [10]: perf = simulation._simulation.performance_root() In [11]: print(perf) simulation: time 104.921 s, count 1\tinitNeighborList: time 0.003026 s, count 1\t\tset_up: time 0.002545 s, count 1\t\t\tsetUpCellNeighbors: time 0.001843 s, count 1\t\t\tsetUpBins: time 0.000693 s, count 1\t\t\t\tallocate: time 1.3e-05 s, count 1\t\t\t\tfillBins parallel: time 0.000673 s, count 1\t\tupdate: time 0.000469 s, count 1\t\t\tsetUpBins: time 0.000465 s, count 1\t\t\t\tallocate: time 1.1e-05 s, count 1\t\t\t\tfillBins parallel: time 0.000446 s, count 1\tforces: time 0.000563 s, count 12501\tintegrator: time 5.49103 s, count 12500\tneighborList: time 3.85435 s, count 25000\t\tupdate: time 3.80827 s, count 25000\t\t\tsetUpBins: time 3.76443 s, count 25000\t\t\t\tallocate: time 0.120096 s, count 25000\t\t\t\tfillBins parallel: time 3.56029 s, count 25000\treactionScheduler: time 93.8283 s, count 12500\tevaluateTopologyReactions: time 7.8e-05 s, count 12500\tclearNeighborList: time 0 s, count 1 ",
"url": "/benchmark.html#"
},
"demonstration-api": {
"title": "Demonstration > Simple example",
"content": "In [1]: import readdyprint(readdy.__version__) v2.0.3-27Step 1: Set up a reaction diffusion systemIn [2]: system = readdy.ReactionDiffusionSystem([25.,25.,25.], temperature=300.*readdy.units.kelvin)system.add_species("A", 1.0) Define reactionsIn [3]: system.reactions.add("myfusion: A +(2) A -> A", rate=10.)system.reactions.add("myfission: A -> A +(2) A", rate=3.) Define potentialsIn [4]: system.potentials.add_harmonic_repulsion("A", "A", force_constant=10., interaction_distance=2.) Step 2: Create a corresponding simulationIn [5]: simulation = system.simulation(kernel="CPU")simulation.output_file = "out.h5"simulation.reaction_handler = "UncontrolledApproximation"simulation.add_particle("A", [0.,0.,0.])simulation.record_trajectory(stride=1)simulation.observe.number_of_particles(stride=100, callback=lambda n: print("#A:", n)) Step 3: Run the simulationIn [6]: simulation.run(n_steps=3000, timestep=1e-2) 0%| | 0/300 [00:00<?, ?it/s] Configured kernel context with:-------------------------------- - kBT = 2.49434 - periodic b.c. = (true, true, true) - box size = (25.0, 25.0, 25.0) - particle types: * particle type "A" with D=1.0 - potentials of order 2: * for types "A" and "A" * Harmonic repulsion with Force constant k=10.0 - unimolecular reactions: * Fission A -> A + A with a rate of 3.0, a product distance of 2.0, and weights 0.5 and 0.5 - bimolecular reactions: * Fusion A + A -> A with a rate of 10.0, an educt distance of 2.0, and weights 0.5 and 0.5Configured simulation loop with:-------------------------------- - timeStep = 0.01 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true#A: [1]#A: [7] 5%|▌ | 16/300 [00:00<00:01, 153.23it/s] #A: [16] 13%|█▎ | 39/300 [00:00<00:02, 116.80it/s] #A: [36]#A: [94] 21%|██ | 63/300 [00:00<00:02, 80.95it/s] #A: [175]#A: [242] 26%|██▌ | 78/300 [00:00<00:03, 68.90it/s] #A: [328]#A: [403] 33%|███▎ | 99/300 [00:01<00:03, 65.20it/s] #A: [450]#A: [512] 40%|████ | 120/300 [00:01<00:02, 60.11it/s] #A: [529]#A: [525] 47%|████▋ | 141/300 [00:01<00:02, 58.78it/s] #A: [502]#A: [507] 54%|█████▎ | 161/300 [00:02<00:02, 61.82it/s] #A: [475]#A: [504] 61%|██████ | 182/300 [00:02<00:01, 63.12it/s] #A: [476]#A: [517] 65%|██████▌ | 196/300 [00:02<00:01, 60.30it/s] #A: [485]#A: [469] 72%|███████▏ | 217/300 [00:03<00:01, 61.38it/s] #A: [516]#A: [503] 79%|███████▉ | 237/300 [00:03<00:01, 56.02it/s] #A: [487]#A: [497] 86%|████████▌ | 257/300 [00:03<00:00, 62.03it/s] #A: [483]#A: [531] 93%|█████████▎| 278/300 [00:04<00:00, 62.94it/s] #A: [482]#A: [489] 100%|██████████| 300/300 [00:04<00:00, 65.86it/s] #A: [496]#A: [533] Step 4: Look at resultsIn [7]: trajectory = readdy.Trajectory('out.h5')trajectory.convert_to_xyz(particle_radii={'A': 1.}) Visualize in VMDIn [8]: !vmd -e out.xyz.tcl /usr/lib/vmd/vmd_LINUXAMD64: /usr/lib/libGL.so.1: no version information available (required by /usr/lib/vmd/vmd_LINUXAMD64)Info) VMD for LINUXAMD64, version 1.9.4a38 (October 20, 2019)Info) http://www.ks.uiuc.edu/Research/vmd/ Info) Email questions and bug reports to vmd@ks.uiuc.edu Info) Please include this reference in published work using VMD: Info) Humphrey, W., Dalke, A. and Schulten, K., `VMD - Visual Info) Molecular Dynamics', J. Molec. Graphics 1996, 14.1, 33-38.Info) -------------------------------------------------------------Info) Multithreading available, 12 CPUs detected.Info) CPU features: SSE2 AVX AVX2 FMA Info) Free system memory: 25GB (78%)Info) Creating CUDA device pool and initializing hardware...Info) Detected 1 available CUDA accelerator::Info) [0] GeForce GTX 1080 20 SM_6.1 1.9 GHz, 7.9GB RAM SP32 KT AE2 ZCWarning) Detected X11 'Composite' extension: if incorrect display occursWarning) try disabling this X server option. Most OpenGL driversWarning) disable stereoscopic display when 'Composite' is enabled.Info) OpenGL renderer: GeForce GTX 1080/PCIe/SSE2Info) Features: STENCIL MSAA(4) MDE CVA MTX NPOT PP PS GLSL(OVFGS) Info) Full GLSL rendering mode is available.Info) Textures: 2-D (32768x32768), 3-D (16384x16384x16384), Multitexture (4)Info) Detected 1 available TachyonL/OptiX ray tracing acceleratorInfo) Compiling 1 OptiX shaders on 1 target GPU...Info) Dynamically loaded 3 plugins in directory:Info) /usr/lib/vmd/plugins/LINUXAMD64/molfileERROR) Error opening file out.xyz.tclERROR) 1.0couldn't open "out.xyz.tcl": no such file or directoryvmd > Info) VMD for LINUXAMD64, version 1.9.4a38 (October 20, 2019)Info) Exiting normally. ",
"url": "/demonstration.html#"
},
"demonstration-constant-concentration": {
"title": "Demonstration > Constant concentration",
"content": "This will demonstrate how to achieve a constant concentration $[S]_\\mathrm{eq}$ of S particles, where the value of $[S]_\\mathrm{eq}$ can be roughly adjusted under some approximations.Define reactions of substrate with factory particles F$$F \\rightleftharpoons F + S$$with forward propensity (intrinsic rate constant) $f_+$ and backward propensity $f_-$. The according macroscopic rates of these processes would be $k_+$ and $k_-$. In equilibrium under assumption of the law of mass action [1]:$$K = \\frac{[F]_\\mathrm{eq}[S]_\\mathrm{eq}}{[F]_\\mathrm{eq}} = [S]_\\mathrm{eq}$$The equilibrium constant also relates $k_+$ and $k_-$, which are functions of the propensities $f_+$ and $f_-$ respectively:$$K = \\frac{k_+(f_+)}{k_-(f_-)}$$Since the \"$+$\" reaction is unimolecular we have $k_+=f_+$. The \"$-$\" reaction depends on the choice of reaction model. Here we choose the Doi model. Furthermore in equilibrium the effective \"$-$\" rate does not depend on diffusion and we can formulate it for a dilute system(well described by an isolated pair, see [2]) and in the absence of potentials yielding $k_-= f_- V_R$. The reaction volume is spherical $V_R=\\frac{4}{3}\\pi R^3$ with reaction radius $R$.Combining the equations above gives the constant concentration$$[S]_\\mathrm{eq} = \\frac{f_+}{f_- V_R}$$[1] Atkins, Peter; de Paula, Julio: Atkins’ physical chemistry, Oxford University Press, 2006 [2] Fröhner, Christoph; Noé, Frank: Reversible interacting-particle reaction dynamics. The Journal of Physical Chemistry B, 122.49 (2018): 11240-11250.In [1]: import osimport numpy as npimport matplotlib.pyplot as pltimport readdyprint(readdy.__version__) v2.0.3-27In [2]: factory_minus_rate = 0.08factory_plus_rate = 0.1factory_reaction_radius = 0.5factory_reaction_volume = 4./3. * np.pi * factory_reaction_radius**3concentration = factory_plus_rate / factory_minus_rate / factory_reaction_volumeprint(concentration) 2.3873241463784303In [3]: system = readdy.ReactionDiffusionSystem( box_size=[10.,10.,10.], periodic_boundary_conditions=[True, True, True], unit_system=None)system.add_species("S", 1.)system.add_species("F", 1.)system.reactions.add_fusion( "factory-fusion", "F", "S", "F", factory_minus_rate, factory_reaction_radius, weight1=0., weight2=1.)system.reactions.add_fission( "factory-fission", "F", "F", "S", factory_plus_rate, factory_reaction_radius, weight1=0., weight2=1.) In [4]: simulation = system.simulation("SingleCPU", reaction_handler="UncontrolledApproximation")simulation.observe.number_of_particles(1, types=["S", "F"])simulation.output_file = "const_concentration.h5"if os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.add_particles("F", np.random.uniform(size=(1000,3)) * 10. - 5.)timestep = 5e-3simulation.run(50000, timestep) 0%| | 11/5000 [00:00<00:48, 103.50it/s] Configured kernel context with:-------------------------------- - kBT = 1.0 - periodic b.c. = (true, true, true) - box size = (10.0, 10.0, 10.0) - particle types: * particle type "F" with D=1.0 * particle type "S" with D=1.0 - unimolecular reactions: * Fission F -> F + S with a rate of 0.1, a product distance of 0.5, and weights 0.0 and 1.0 - bimolecular reactions: * Fusion F + S -> F with a rate of 0.08, an educt distance of 0.5, and weights 0.0 and 1.0Configured simulation loop with:-------------------------------- - timeStep = 0.00500000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 5000/5000 [02:39<00:00, 31.29it/s]In [6]: traj = readdy.Trajectory(simulation.output_file)time, counts = traj.read_observable_number_of_particles()volume = 10**3plt.plot(time * timestep, counts[:,0]/volume, label=r"$[S]$")plt.hlines(concentration, *plt.xlim(), label=r"theoretical value $[S]_\\mathrm{eq}$")plt.legend(loc="lower right")plt.xlabel(r"Time in $[\\mathrm{time}]$")plt.ylabel(r"Concentration in $[\\mathrm{length}]^{-3}$")plt.show() ",
"url": "/demonstration.html#"
},
"demonstration-growing-polymer": {
"title": "Demonstration > Growing polymer",
"content": "In [1]: import numpy as npimport readdyimport os In [2]: readdy.__version__ Out[2]:'v2.0.3-27'Set up reaction-diffusion system¶In [3]: system = readdy.ReactionDiffusionSystem(box_size=[30., 30., 30.]) In [4]: system.add_species("Substrate", 1.)system.topologies.add_type("Polymer")system.add_topology_species("Head", 1.)system.add_topology_species("Core", 1.) In [5]: system.topologies.configure_harmonic_bond("Head", "Core", force_constant=70, length=1.)system.topologies.configure_harmonic_bond("Core", "Core", force_constant=70, length=1.)system.topologies.configure_harmonic_angle("Core", "Core", "Core", force_constant=5., equilibrium_angle=np.pi)system.potentials.add_harmonic_repulsion("Core", "Core", force_constant=70., interaction_distance=1.) Define reaction that is called \"Bind\" and attaches a \"Substrate\" particle to a \"Polymer\" topologies \"Head\" particle with a fixed rate once they are close enough, changing types from \"Head\" to \"Core\" and from \"Substrate\" to \"Head\".In [6]: system.topologies.add_spatial_reaction( "Bind: Polymer(Head) + (Substrate) -> Polymer(Core--Head)", rate=10.0, radius=1.5) Set up simulation¶In [7]: simulation = system.simulation(kernel="SingleCPU") One initial topology with three particles (Head--Core--Head).In [8]: init_top_pos = np.array([ [ 1. ,0. ,0.], [ 0. ,0. ,0.], [-1. ,0. ,0.]])top = simulation.add_topology("Polymer", ["Head", "Core", "Head"], init_top_pos)top.get_graph().add_edge(0, 1)top.get_graph().add_edge(1, 2) Some substrate particles, randomly distributed in the simulation box.In [9]: n_particles = 500positions = np.random.uniform(size=(500,3)) * system.box_size.magnitude - system.box_size.magnitude * 0.5simulation.add_particles("Substrate", positions) Save the trajectory and topologies observable to an output file (which is removed if it already exists).In [10]: simulation.output_file = 'growing_polymer.h5'if os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.observe.topologies(10)simulation.record_trajectory(10)simulation.progress_output_stride = 10 Run the simulation.In [11]: simulation.run(200000, .001) 0%| | 20/20000 [00:00<01:42, 195.66it/s] Configured kernel context with:-------------------------------- - kBT = 2.43614 - periodic b.c. = (true, true, true) - box size = (30.0, 30.0, 30.0) - particle types: * Topology particle type "Core" with D=1.0 * particle type "Substrate" with D=1.0 * Topology particle type "Head" with D=1.0 - potentials of order 2: * for types "Core" and "Core" * Harmonic repulsion with Force constant k=70.0 - topology potential configuration: - bonds (2): - Bonds for particle types Core and Core: * Harmonic bond with force constant 70.0 and length 1.0 - Bonds for particle types Head and Core: * Harmonic bond with force constant 70.0 and length 1.0 - angles (1): * Harmonic angle with force constant 5.0 and equilibrium angle 3.14159 - topology types: * topology type "Polymer" with 0 structural reactions - spatial topology reactions: * Topology-particle fusion reaction "Bind: Polymer(Head) + (Substrate) -> Polymer(Core--Head)"Configured simulation loop with:-------------------------------- - timeStep = 0.00100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 20000/20000 [01:35<00:00, 208.84it/s]Read back the recorded data.In [12]: traj = readdy.Trajectory(simulation.output_file)time, topology_records = traj.read_observable_topologies() The topology_records object is a list of lists, each outer element representing a frame and each frame containing a list of topologies at that point of the simulation.In [13]: chain_length = [ len(tops[0].particles) for tops in topology_records ] In [14]: %matplotlib inlineimport matplotlib.pyplot as pltplt.plot(time, chain_length)plt.show() ",
"url": "/demonstration.html#"
},
"demonstration-living-polymers": {
"title": "Demonstration > Living polymers",
"content": "In [1]: import numpy as npimport readdyimport matplotlib.pyplot as pltimport osfrom collections import defaultdict In [2]: readdy.__version__ Out[2]:'v2.0.3-27'In [3]: system = readdy.ReactionDiffusionSystem(box_size=[100, 100, 100]) In [4]: system.topologies.add_type("Polymer")system.add_topology_species("Head", .002)system.add_topology_species("Tail", .002) In [5]: system.topologies.configure_harmonic_bond("Head", "Tail", force_constant=50, length=1.)system.topologies.configure_harmonic_bond("Tail", "Tail", force_constant=50, length=1.)system.topologies.configure_harmonic_bond("Head", "Head", force_constant=50, length=1.)system.topologies.configure_harmonic_angle("Head", "Tail", "Tail", force_constant=10, equilibrium_angle=np.pi)system.topologies.configure_harmonic_angle("Tail", "Tail", "Tail", force_constant=10, equilibrium_angle=np.pi) In [6]: # dissociationdef dissociation_rate_function(topology): edges = topology.get_graph().get_edges() return .000005 * float(len(edges)) if len(edges) > 2 else 0.def dissociation_reaction_function(topology): recipe = readdy.StructuralReactionRecipe(topology) edges = topology.get_graph().get_edges() vertices = topology.get_graph().get_vertices() # at least a structure like # v1 -- v2 -- v3 -- v4 if len(edges) > 2: # find the end particles counts = defaultdict(int) for e in edges: pix1 = e[0].get().particle_index pix2 = e[1].get().particle_index counts[pix1] += 1 counts[pix2] += 1 # the end particles are the ones that have exactly one edge leading from/to them endpoints = [] for pix in counts.keys(): if counts[pix] == 1: endpoints.append(pix) assert len(endpoints) == 2, "the number of end points should always be 2 \\ but was {} (counts: {})".format(len(endpoints), counts) # draw an edge excluding the edges leading to the both ends edge_index = np.random.randint(0, len(edges) - 2) removed_edge = None # for each edge in the graph for e in edges: pix1 = e[0].get().particle_index pix2 = e[1].get().particle_index # check if it belongs to one of the end vertices if pix1 not in endpoints and pix2 not in endpoints: # if not, reduce edge_index until 0 if edge_index == 0: # remove this edge recipe.remove_edge(e[0], e[1]) removed_edge = e break else: edge_index -= 1 assert removed_edge is not None pix1 = removed_edge[0].get().particle_index pix2 = removed_edge[1].get().particle_index assert pix1 not in endpoints assert pix2 not in endpoints # since the edge was removed we now have two topologies and need to set the correct particle types recipe.change_particle_type([vix for vix, v in enumerate(vertices) if v.particle_index == pix1][0], "Head") recipe.change_particle_type([vix for vix, v in enumerate(vertices) if v.particle_index == pix2][0], "Head") else: print("this should not have happened") return recipesystem.topologies.add_structural_reaction("dissociation", "Polymer", dissociation_reaction_function, dissociation_rate_function) In [7]: system.topologies.add_spatial_reaction("Association: Polymer(Head) + Polymer(Head) -> Polymer(Tail--Tail)", rate=1.0, radius=1.0) In [8]: simulation = system.simulation(kernel="CPU") In [9]: # randomly place some polymers of length 4n_polymers = 500for head in range(n_polymers): head_position = 80. * np.random.random((1, 3)) - 40. offset1 = 2.*np.random.random((1, 3))-1. offset1 /= np.linalg.norm(offset1) tail1 = head_position + offset1 offset2 = 2.*np.random.random((1, 3))-1. offset2 /= np.linalg.norm(offset2) tail2 = head_position + offset1 + offset2 offset3 = 2.*np.random.random((1, 3)) - 1. offset3 /= np.linalg.norm(offset3) head_position2 = head_position + offset1 + offset2 + offset3 top = simulation.add_topology("Polymer", ["Head", "Tail", "Tail", "Head"], np.array([head_position, tail1, tail2, head_position2]).squeeze()) top.get_graph().add_edge(0, 1) top.get_graph().add_edge(1, 2) top.get_graph().add_edge(2, 3) In [10]: simulation.output_file = 'living_polymers.h5'simulation.observe.topologies(300)simulation.record_trajectory(300)simulation.progress_output_stride = 10simulation.show_progress = True In [11]: timestep = 1.if os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.run(50000, timestep) 0%| | 0/5000 [00:00<?, ?it/s] Configured kernel context with:-------------------------------- - kBT = 2.43614 - periodic b.c. = (true, true, true) - box size = (100.0, 100.0, 100.0) - particle types: * Topology particle type "Tail" with D=0.00200000 * Topology particle type "Head" with D=0.00200000 - topology potential configuration: - bonds (3): - Bonds for particle types Head and Head: * Harmonic bond with force constant 50.0 and length 1.0 - Bonds for particle types Head and Tail: * Harmonic bond with force constant 50.0 and length 1.0 - Bonds for particle types Tail and Tail: * Harmonic bond with force constant 50.0 and length 1.0 - angles (2): * Harmonic angle with force constant 10.0 and equilibrium angle 3.14159 * Harmonic angle with force constant 10.0 and equilibrium angle 3.14159 - topology types: * topology type "Polymer" with 1 structural reactions - spatial topology reactions: * Topology-topology fusion reaction "Association: Polymer(Head) + Polymer(Head) -> Polymer(Tail--Tail)" - structural topology reactions: - for topology type "Polymer" with 1 structural reactions: * reaction with create child topologies = trueConfigured simulation loop with:-------------------------------- - timeStep = 1.0 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 5000/5000 [03:15<00:00, 25.63it/s]In [12]: # read back topologiest = readdy.Trajectory(simulation.output_file)time, topology_records = t.read_observable_topologies() In [13]: # gather average length of topologies for respective time stepavg_length = []# for each time stepfor topologies in topology_records: # gather average polymer length avg_length.append(0) for top in topologies: avg_length[-1] += len(top.edges) avg_length[-1] /= len(topologies)f, ax = plt.subplots(nrows=1, ncols=1)ax.plot(time, avg_length)ax.set_title('average length')ax.set_xlabel('time')ax.set_ylabel('# beads')plt.show() In [14]: from IPython.display import YouTubeVideoYouTubeVideo('1fZqbZRQnEQ') Out[14]: ",
"url": "/demonstration.html#"
},
"demonstration-topologies": {
"title": "Demonstration > Topologies",
"content": "In [1]: import numpy as npimport readdy In [2]: readdy.__version__ Out[2]:'v2.0.3-27'Step 1: Define the systemIn [3]: system = readdy.ReactionDiffusionSystem(box_size=[150, 150, 150])system.periodic_boundary_conditions = False, False, Falsesystem.add_species("Ligand", diffusion_constant=3.)system.add_species("Decay", diffusion_constant=1.)system.add_topology_species("T", diffusion_constant=1.)system.add_topology_species("unstable T", diffusion_constant=1.) Add a decay reaction for \"Decay\" praticles and harmonic repulsion between \"unstable T\" and \"Decay\"In [4]: system.reactions.add("decay: Decay ->", .1)system.potentials.add_box("Ligand", 10., [-70, -70, -70], [130, 130, 130])system.potentials.add_box("Decay", 10., [-70, -70, -70], [130, 130, 130])system.potentials.add_box("T", 10., [-70, -70, -70], [130, 130, 130])system.potentials.add_box("unstable T", 10., [-70, -70, -70], [130, 130, 130])system.potentials.add_harmonic_repulsion("Decay", "unstable T", force_constant=20., interaction_distance=2.) Addharmonic bonds and angles between pairs/triples of \"T\" and \"unstable T\"three topology types: stable (with no topology reactions), intermediate, and unstableIn [5]: system.topologies.configure_harmonic_bond("T", "T", force_constant=20., length=2.)system.topologies.configure_harmonic_bond("unstable T", "unstable T", force_constant=20., length=2.)system.topologies.add_type("stable")system.topologies.add_type("intermediate")system.topologies.add_type("unstable") Change the topology type from stable to intermediateIn [6]: system.topologies.add_spatial_reaction( "encounter: stable(T) + (Ligand) -> intermediate(T) + (Ligand)", rate=10.0, radius=2.0) Change the topology type from intermediate to unstable, change particle types from T to unstable TIn [8]: def intermediate_rate_function(topology): return 1e3def intermediate_reaction_function(topology): recipe = readdy.StructuralReactionRecipe(topology) for i in range(len(topology.get_graph().get_vertices())): recipe.change_particle_type(i, "unstable T") recipe.change_topology_type("unstable") return recipesystem.topologies.add_structural_reaction(name='to_intermediate', topology_type="intermediate", reaction_function=intermediate_reaction_function, rate_function=intermediate_rate_function) Randomly select a vertex and separate it from the graph, change it's particle type to DecayIn [9]: def unstable_rate_function(topology): return .1def unstable_reaction_function(topology): recipe = readdy.StructuralReactionRecipe(topology) index = np.random.randint(0, len(topology.particles)) recipe.separate_vertex(index) recipe.change_particle_type(index, "Decay") return recipesystem.topologies.add_structural_reaction(name="to_unstable", topology_type="unstable", reaction_function=unstable_reaction_function, rate_function=unstable_rate_function) Step 2: Create the simulation object out of the systemIn [10]: simulation = system.simulation(kernel="CPU") add topologyIn [11]: n_topology_particles = 70positions = [[0, 0, 0], np.random.normal(size=3)]for i in range(n_topology_particles-2): delta = positions[-1] - positions[-2] offset = np.random.normal(size=3) + delta offset = offset / np.linalg.norm(offset) positions.append(positions[-1] + 2.*offset)topology = simulation.add_topology(topology_type="stable", particle_types="T", positions=np.array(positions)) set up connectivity of topologyIn [12]: graph = topology.get_graph()for i in range(len(graph.get_vertices())-1): graph.add_edge(i, i+1) add ligandsIn [13]: simulation.add_particles("Ligand",-6 * np.ones((5, 3))) In [14]: simulation.output_file = "topology_simulation.h5"simulation.record_trajectory() In [15]: simulation.run(n_steps=10000, timestep=1e-2) 0%| | 0/1000 [00:00<?, ?it/s] Configured kernel context with:-------------------------------- - kBT = 2.43614 - periodic b.c. = (false, false, false) - box size = (150.0, 150.0, 150.0) - particle types: * Topology particle type "unstable T" with D=1.0 * Topology particle type "T" with D=1.0 * particle type "Ligand" with D=3.0 * particle type "Decay" with D=1.0 - potentials of order 1: * for type "unstable T" * Box potential with origin=(-70, -70, -70), extent=(130, 130, 130), and Force constant k=10.0 * for type "T" * Box potential with origin=(-70, -70, -70), extent=(130, 130, 130), and Force constant k=10.0 * for type "Ligand" * Box potential with origin=(-70, -70, -70), extent=(130, 130, 130), and Force constant k=10.0 * for type "Decay" * Box potential with origin=(-70, -70, -70), extent=(130, 130, 130), and Force constant k=10.0 - potentials of order 2: * for types "Decay" and "unstable T" * Harmonic repulsion with Force constant k=20.0 - unimolecular reactions: * Decay Decay -> ø with a rate of 0.1 - topology potential configuration: - bonds (2): - Bonds for particle types unstable T and unstable T: * Harmonic bond with force constant 20.0 and length 2.0 - Bonds for particle types T and T: * Harmonic bond with force constant 20.0 and length 2.0 - topology types: * topology type "stable" with 0 structural reactions * topology type "intermediate" with 1 structural reactions * topology type "unstable" with 1 structural reactions - spatial topology reactions: * Topology-particle enzymatic reaction "encounter: stable(T) + (Ligand) -> intermediate(T) + (Ligand)" - structural topology reactions: - for topology type "intermediate" with 1 structural reactions: * reaction with create child topologies = true - for topology type "unstable" with 1 structural reactions: * reaction with create child topologies = trueConfigured simulation loop with:-------------------------------- - timeStep = 0.01 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 1000/1000 [00:35<00:00, 27.98it/s]In [16]: traj = readdy.Trajectory(simulation.output_file) In [17]: traj.convert_to_xyz() In [18]: from IPython.display import YouTubeVideoYouTubeVideo('9cG2J1Nihnk') Out[18]: ",
"url": "/demonstration.html#"
},
"development-architecture-html": {
"title": "Development > Architecture",
"content": "Source treereaddy/├── README.md│ ...│├── cmake/│ ││ ├── Modules/│ │ ├── (cmake modules)│ └── sources/│ └── (readdy source file lists)│├── docs/│ └── (internal docs, doxygen configuration)│├── kernels/│ ││ ├── cpu/│ │ ├── include/│ │ │ └── (kernel includes)│ │ ├── src/│ │ │ └── (kernel sources)│ │ └── test/│ │ └── (kernel tests)│ ││ └── cuda/│ └── (yet to be added)│├── include/│ └── *.h (core library includes)│├── readdy/│ ││ ├── main/│ │ └── (core library sources)│ ├── test/│ │ └── (core library tests)│ └── examples/│ └── (cpp examples)│├── readdy_testing/│ └── (unit testing relevant tools)│├── libraries/│ └── (c-blosc, h5rd, pybind11, spdlog)│└── wrappers/ └── python/ └── src/ (code for python api) ├── cxx/ │ └── (c++ code for the python module) │ └── python/ └── (python api and tests)",
"url": "/development.html#architecture"
},
"development-build-html": {
"title": "Development > Building from source",
"content": "ReaDDy has the following dependencies (bold printed are expected to be present on the system): hdf5 cmake spdlog, included by git submodule c-blosc, included by git submodule h5rd, included by git submodule catch2, for testing, included by git submodule for python bindings pybind11, included by git submodule python 3 with numpy, h5py Build by using CMakeThis type of build is suggested if one is interested in development of the software. There are a number of ReaDDy specific CMake options that influence the type of build: CMake option = default Description READDY_CREATE_TEST_TARGET=ON Determining if the test targets should be generated. READDY_CREATE_MEMORY_CHECK_TEST_TARGET=OFF Determining if the test targets should be additionally called through valgrind. Requires the previous option to be ON and valgrind to be installed. READDY_INSTALL_UNIT_TEST_EXECUTABLE=OFF Determining if the unit test executables should be installed. This is option is mainly important for the conda recipe. READDY_BUILD_SHARED_COMBINED=OFF Determining if the core library should be built monolithically or as separated shared libraries. READDY_BUILD_PYTHON_WRAPPER=ON Determining if the python wrapper should be built. READDY_DEBUG_PYTHON_MODULES=OFF If this flag is set to ON, the generated python module will be placed in-source rather than in the output directory to enable faster development. READDY_DEBUG_CONDA_ROOT_DIR=\"\" This option is to be used in conjunction with the previous option and only has effect if it is set to ON. It should point to the conda environment which is used for development and then effects the output directory of the binary files such that they get compiled directly into the respective environment. READDY_GENERATE_DOCUMENTATION_TARGET=OFF Determines if the documentation target should be generated or not, which, if generated, can be called by “make doc”. READDY_GENERATE_DOCUMENTATION_TARGET_ONLY=OFF This option has the same effect as the previous option, just that it does not need any dependencies other than doxygen to be fulfilled and generates the documentation target exclusively. READDY_LOG_CMAKE_CONFIGURATION=OFF This option determines if the status of relevant cmake cache variables should be logged at configuration time or not. READDY_KERNELS_TO_TEST=\"SingleCPU,CPU\" Comma separated list of kernels against which the core library should be tested within the test targets. After having configured the cmake cache variables, one can invoke cmake and make and compile the project.Altogether, a shell script invoking cmake with modified parameters in an environment with multiple python versions could look like this.Build by using conda-buildconda install conda-buildconda-build PATH_TO_READDY/tools/conda-recipe",
"url": "/development.html#build"
},
"development-topologies-html": {
"title": "Development > Topologies",
"content": "A topology is a collection of particles, whose unique ids are stored in a std::vector. These particles are subject to bond-, angle- and dihedral-potentials. For a pair (or triple or quadruple) of particles to have a potential term, they have to be connected. The connection is defined by a graph, which is a linked list of vertices. Each vertex represents one particle in the topology and contains references to other vertices. The actual potential terms are yielded by a lookup table on pairs/triples/quadruples of particle types. All potentials holding the topology together are thus obtained from both the Graph and the PotentialConfiguration.Internal structure of topologiesTopology contains global particle indices and potential terms between them, referencing to local indices w.r.t. the global-particle-index-vectorGraphTopology is a derived class of Topology and has the following properties has a topology type which currently defines the possible structural reactions contains a graph consisting out of vertices that have a one-to-one correspondence to particle indices in the respective topology: the data structure for the vertices is a linked list, hence iterators can be used as vertex references and a tuple of vertex references denotes an edge in the graph a topology graph needs to be connected upon simulation start has PotentialConfiguration that contains particle-type pairs/triples/quadruples definitions for bonds/angles/torsion-potentials, respectively the graph’s connectivity together with the potential configuration gives the actual potential terms of the topology topology needs to be connected w.r.t. bonds as yielded by graph+potential config upon simulation start contains a vector of rates for structural reactions as registered for the topology typePotential termsPotential terms hold the topology together.They are configured in the readdy::api::PotentialConfiguration. This kernel-unique object holds three maps: pairPotentials: mapping (type1, type2) -> bonds anglePotentials: mapping (type1, type2, type3) -> angles torsionPotentials: mapping (type1, type2, type3, type4) -> torsion angles.These maps apply a hashing and equality operator that allows for asking for the reverse key, e.g., (type1, type2, type3) and (type3, type2, type1) should yield the same value.Topology reactionsThere are structurally-dependent and spatially-dependent reactions:Structural means that the reaction recipe as well as its rate only depend on the topology’s graph, no dependence on other particles or other topologies is possible. Spatial means that the reaction is occurring between two particles due to their spatial proximity. At least one of the particles must belong to a topology. structural reactions are defined on the topology-type. They require a rate-function and a reaction-function. With structural reactions one can realize Conversion of a topology, just a change of the graph Fission, change of the graph but ending up with disconnected components. Think of a linear polymer that breaks apart in the middle Special case: Fission of a topology that yields components which consist of only one particle with the flavor “NORMAL”. Here, the topology will be erased and the particle will be treated normally. spatial reactions are between two particle-types of which at least one has topology flavor with a certain topology-type TopologyParticleReaction TP, association of a normal particle and a topology particle, the normal particle becomes a topology flavor and joins the list of particles and an edge in the graph is established. TopologyTopologyReaction TT, association of two particles that already belong to a topology of a certain type, if these two are different topologies, their lists of particles and graphs are merged, and an edge is established between the given vertices. The type of the merged topology is determined by the reaction All topology reactions are stored in the TopologyRegistryStructural topology reactions detailA structural reaction may change the connectivity of the topologies’ graph and may also change the types of the vertices, i.e. change the particle-type of a particle.Rate functionA rate function has the signature scalar(const GraphTopology&) and is supposed to yield a rate depending on the current state of the topology. This rate is assumed to be constant as long as the topology keeps its type and its graph does not change.Reaction functionA reaction function has the signature reaction_recipe(GraphTopology&) and is supposed to yield a reaction recipe that is a set of operations which will then be applied to the topology by the selected computing kernel. Available operations are: addEdge(vertex1, vertex2): introduce an edge between two vertices removeEdge(vertex1, vertex2): remove an edge between two vertices separateVertex(vertex): remove all edges going from or to the given vertex changeParticleType(vertex, newType): change the type of the particle corresponding to the given vertex changeTopologyType(newType): change the topology’s typeSpatial topology reactions detailSpatial reactions are defined on both particle types and topology types, before and after the reaction.Additionally spatial reactions differ in behavior with respect to merging the topologies or not. Analogue to fusion-like or enzymatic-like reactions.Particle- and topology-types before and after the reactionTo simplify the definition of TP and TT reactions, one can use a describing string such as:$$\\text{label: } T_1 (P_1) + T_2 (P_2)\\to T_3(P_3) + T_4(P_4),$$where $T_i$ and $P_i$ denote topology- and particle-types, respectively. The above example would be for an enzymatic TT reaction. In case of an enzymatic TP reaction, one can omit $T_2$ and $T_4$. In case one would like to perform a fusion reaction, the notation is$$\\text{label: } T_1 (P_1) + T_2 (P_2)\\to T_3(P_3\\text{--}P_4),$$where again $T_2$ can be omitted if a TP reaction is wanted.In the special case of a TT-Fusion, one can additionally specify [self=true] at the end of the describing string to indicate that the topology may also form bonds with its own particles.Examples: TT-Fusion: T1 (p1) + T2 (p2) -> T3 (p3--p4) [self=true] TT-Fusion2: T1 (p1) + T2 (p2) -> T3 (p3--p4), where self=false is used as default TT-Enzymatic: T1 (p1) + T2 (p2) -> T3 (p3) + T4 (p4) TP-Fusion: T1 (p1) + (p2) -> T2 (p3--p4) TP-Enzymatic: T1 (p1) + (p2) -> T3 (p3) + (p4)Spatial reaction modeThe behavior of a reaction is summarized internally by an enum SpatialTopologyReactionMode, which is extracted from the user’s descriptor string: SpatialTopologyReactionMode meaning TT_ENZYMATIC reaction between two topologies, no edge is created TT_FUSION reaction between two topologies, edge can only be established between particles of different topology instances TT_FUSION_ALLOW_SELF reaction between two topologies, edge can be established between particles of different instances and even within the same instance of a topology TP_ENZYMATIC reaction between topology and particle changing types but not establishing a bond TP_FUSION reaction between topology and particle introducing a bond and possibly changing types Use cases: problem STRMode Linear polymers with two end-particles each. Polymers can fuse together at the ends, but not with themselves SpatialTopologyReaction with STRMode==TT_FUSION A complex/topology has an active site, that switches its type when a ligand/normal particle is close SpatialTopologyReaction with STRMode==TP_ENZYMATIC ",
"url": "/development.html#topologies"
},
"home-cite-html": {
"title": "Home > Citation",
"content": "When using ReaDDy in your work, please cite the following paper@article{hoffmann2019readdy, title={ReaDDy 2: Fast and flexible software framework for interacting-particle reaction dynamics}, author={Hoffmann, Moritz and Fr{\\\"o}hner, Christoph and No{\\'e}, Frank}, journal={PLoS computational biology}, volume={15}, number={2}, pages={e1006830}, year={2019}, publisher={Public Library of Science}}Also note these further publications related to ReaDDy: Schöneberg, Johannes, et al. Lipid-mediated PX-BAR domain recruitment couples local membrane constriction to endocytic vesicle fission. Nature communications 8 (2017): 15873. Schöneberg, Johannes, and Frank Noé. ReaDDy-a software for particle-based reaction-diffusion dynamics in crowded cellular environments. PloS one 8.9 (2013): e74261. Dibak, Manuel, et al. Diffusion-influenced reaction rates in the presence of pair interactions. The Journal of chemical physics 151.16 (2019): 164105. Fröhner, Christoph, and Frank Noé. Reversible interacting-particle reaction dynamics. The Journal of Physical Chemistry B 122.49 (2018): 11240-11250. Ullrich, Alexander, et al. Dynamical organization of syntaxin-1A at the presynaptic active zone. PLoS computational biology 11.9 (2015): e1004407. Gunkel, Monika, et al. Higher-order architecture of rhodopsin in intact photoreceptors and its implication for phototransduction kinetics. Structure 23.4 (2015): 628-638. Biedermann, Johann, et al. ReaDDyMM: Fast interacting particle reaction-diffusion simulations using graphical processing units. Biophysical journal 108.3 (2015): 457-461. Schöneberg, Johannes, et al. Explicit spatiotemporal simulation of receptor-G protein coupling in rod cell disk membranes. Biophysical journal 107.5 (2014): 1042-1053. Schöneberg, Johannes, Alexander Ullrich, and Frank Noé. Simulation tools for particle-based reaction-diffusion dynamics in continuous space. BMC biophysics 7.1 (2014): 11.",
"url": "/index.html#citation"
},
"home-get-started-html": {
"title": "Home > Get started",
"content": "import readdy# ----- Step 1: Set up reaction diffusion systemsystem = readdy.ReactionDiffusionSystem(box_size=(10, 10, 10))system.add_species(\"A\", diffusion_constant=2.0)system.reactions.add(\"mydecay: A ->\", rate=1.)system.reactions.add(\"myfission: A -> A +(1) A\", rate=3.)# ----- Step 2: Create simulation instance out of configured systemsimulation = system.simulation(kernel=\"CPU\")simulation.observe.number_of_particles(stride=5)simulation.output_file = \"out.h5\"simulation.add_particle(\"A\", [0.,0.,0.])# ------ Step 3: run the simulationsimulation.run(100, 0.01)The above snippet performs a ReaDDy simulation, which consists of three steps: Configure the system Setup and run the simulation Analyze resultsSee this ipython notebook for an example of the basic features",
"url": "/index.html#get_started"
},
"home-what-is-html": {
"title": "Home > What is ReaDDy?",
"content": " The logo simulation mimicks a predator prey system, i.e., a population growth process that frequently occurs in biology. Sometimes, this growth process is subjected to spatial constraints. There are three different particle types, referring to that biological model: Type 1, the red “logo particles“, serve as the spatial barriers. They have been given an attraction potential between them and start in a position that resembles the ReaDDy logo. Type 2, the purple “prey“. If there are no predators around, they will replicate. Type 3, the grey “predator” particles. They die out if there is no prey but replicate in their presence by consuming them.It is visible during the time course of the simulation, that the spatial distribution of the particles, their crowding inducing occurrence in masses as well as spatial constraints like barriers influence the growth of the populations dramatically. What is true for this simplified example is ubiquitous not only in molecular and cellular biology but in multiple other fields.ReaDDy has been designed to fit the modeling requirements of such processes: Particle (or agent) based reaction diffusion systems in which particle-particle interactions play an important role and where the systems are subjected to crowding or spatial constraints.",
"url": "/index.html#what_is"
},
"simulation-configuration-html": {
"title": "Simulation > Configuration",
"content": "In the following it will be explained how to add particles, add topologies,configure specifics of the selected kernel, how to record a trajectory, and how to perform checkpointing.Adding particlesFor adding particles to a system there are two separate methods. One is can be to place a single particle, oneis for bulk insertion.Adding a single particle of type A to a simulation box amounts tosimulation.add_particle(type=\"A\", position=pos)where pos can be a list [x, y, z], tuple (x, y, z), or a numpy.ndarray: np.array([x, y, z]) with three entriesrepresenting the x,y,z components.When one wants several particles of a certain type to the simulation, one can can exchange multiple calls tosimulation.add_particle by the better performing variantX = np.random.random((100, 3))simulation.add_particles(type=\"A\", positions=X)taking a (N, 3)-shaped numpy array as position argument, resulting in N particles with their respective positionsbeing added to the simulation. In this example, 100 particles of type A would be placed uniformly at randomin $[0,1)^3$.Adding topologiesA topology can be added by invokingmy_topology = simulation.add_topology(topology_type=\"My topology type\", particle_types=\"T\", positions=np.random.random((100, 3)))which requires a “My topology type” topology type and a topology species “T” to be registered in the ReactionDiffusionSystem. In this example the topology will contain 100 randomly placed topology particles of type “T” that are for now disconnected.If the topology should contain several different particle types one can pass a list of particle types to the particle_types argumentthat contains types for all the positions:my_topology = simulation.add_topology( topology_type=\"My topology type\", particle_types=[\"T1\", \"T2\", \"T3\", \"T2\", \"T1\"], positions=np.random.random((5, 3)))Unless the topology consists out of only one particle, one still needs to set up the connectivity graph before running the simulation. The returned object my_topology is a topology object as the ones described in topology reactions. Edges in the graph can be introduced likemy_graph = my_topology.get_graph()for i in range(len(graph.get_vertices())-1): my_graph.add_edge(i, i+1)where the indices that go into add_edge correspond to the particle positions that were entered in add_topology.The simulation can only be executed if the graph of each topology is connected, i.e., there are no independentcomponents (between each pair of vertices there is at least one path in the graph that connects them),and for each edge there is a bond, i.e., all particle type pairs that are contained in the graph have at least one entry in the topology potential configuration.Should one of these two conditions be not fulfilled, starting the simulation will raise an error.Kernel configurationIn case of selecting the CPU kernel with a parallelized implementation of the ReaDDy model, one can change certainaspects of its behavior: The number of threads to be used can be selected by simulation.kernel_configuration.n_threads = 4 Mainly due to pairwise interactions and bimolecular reactions there is a neighbor list to reduce the time needed for evaluating these. The neighbor list imposes a spatial discretization on the simulation box into cuboids. In thedefault case, each of these cuboids has an edge length of at least the maximal cutoff radius / reaction radius.This means that instead of naively looping over all particle pairs ($\\mathcal{O}(N^2)$), one can assign each particleto its cuboid and then loop over all particles in a cuboid and its 26 neighboring cuboids to find particle pairs. When collecting particle pairs in this fashion one effectively approximates a sphere with cuboids. The number ofpotential interaction or reaction partners can be further reduced by using only a fraction of the edge length butincreasing the search radius of the neighboring boxes so that one still covers at least the cutoff radius in eachspatial dimension. Reducing the edge length usually comes with a price, at some point the bookkeeping of neighboring boxes dominatesthe runtime. The edge length and therefore search radius can be controlled by simulation.kernel_configuration.cell_linked_list_radius = 4 which would yield cuboids with $\\frac{1}{4}$ the edge lengths of the default case. Recording a trajectoryReaDDy records trajectories and observable data in HDF5 files. For doing so one needs to set an output filesimulation.output_file = \"my_trajectory.h5\"and instruct the simulation to record a trajectory:simulation.record_trajectory(stride=1, name=\"\", chunk_size=1000)The stride arguments causes the trajectory to be recorded every stride time steps. If a name (other thanthe default one) is given, the trajectory data will be stored in a different data set. The chunk_size is mainlya performance parameter that has an effect on how large every chunk of data in the binary file is going to be,influencing the time needed for IO during the simulation and the resulting file size.For reading back the trajectory data, please refer to post-processing. NOTE: When running long simulations on a cluster it can happen that the process runs into a timeout causing the already recorded data to be corrupted. This can possibly be mitigated by configuring the job manager to send SIGINT before KILLing the process. That way the file can still be properly closed (see issue #220, thanks @jansteinkuehler). Please make sure this works in your environment before running long simulations.CheckpointingCheckpoints in ReaDDy consist out of the particles’ and topologies’ configurations at specific points in simulation time. They can be enabled by callingsimulation.make_checkpoints(stride=1000, output_directory=\"checkpoints/\", max_n_saves=10)which causes checkpoints to be made every 1000 steps. Each checkpoint is a separate file and all checkpoint files will besaved to the directory specified by output_directory. The option max_n_saves decides how many checkpoint filesare allowed to be saved to the directory, e.g. if max_n_saves=10 then only the last 10 most recent checkpointsare kept.Once the simulation has run its course and checkpoint files have been written, they can be listed bysimulation.list_checkpoint_files('checkpoints/')A particular checkpoint file can in principle also contain multiple checkpoints. They can be inspected bysimulation.list_checkpoints('checkpoints/checkpoint_10000.h5')and a system’s state can be restored by a call tosimulation.load_particles_from_checkpoint('checkpoints/checkpoint_10000.h5')amounting to restoring the latest checkpoint of that particular file. If the file contains multiple checkpoints, let’s say 5, you canselect the 5th checkpoint by supplying the optional argument n=4 (enumeration starts at n=0 per file).Oftentimes you just need the last checkpoint of all checkpoint files in a certain directory. This can be achieved bysimulation.load_particles_from_latest_checkpoint('checkpoints/')It should be noted that if the simulation should be continued and the output_directory for the new checkpoints is the same as of the original simulation, the old checkpoint files will be overwritten. If you want to keep the checkpointsof the original simulation, specify another output_directory.",
"url": "/simulation.html#simulation_configuration"
},
"simulation-observables-html": {
"title": "Simulation > Observables",
"content": "The currently available observables are: Radial distribution function Particles Particle positions Number of particles Energy Forces Reactions Reaction counts Virial PressureThere are three things that all observables have in common: The evaluation can be strided, they can have a callback function and they can be saved to the simulation.output_file.The callback is usually a function that takes one argument being the current value of the observable. During the course of the simulation this callback function will be evaluated whenever the particular observable is evaluated.Per default, whenever an output_file is given, the registered observables’ outputs are saved to that file. Eachobservable has a certain place in the group hierarchy of the HDF5 file, however this place can be modified so that, e.g., multiple observables of the same type can be recorded into the same file.To this end, the save argument of the respective observable can be modified. By providing None or False writing the results to file can be switched off, providing a dictionary with keys 'name' and 'chunk_size' can modify the name under which the observable data is storedin the group hierarchy and the hdf5 chunk size.The chunk_size is always to be considered into the “temporal direction”, i.e., if an observable yields data in the form of3x3 matrices each time it is evaluated, a chunk would be of shape (3, 3, chunk_size).Radial distribution functionThe radial distribution function for certain particle types can be observed bydef rdf_callback(current_value): print(current_value)simulation.observe.rdf( stride=1, bin_borders=np.linspace(0., 2., 10), types_count_from=[\"A\"], types_count_to=[\"B\"], particle_to_density=1./system.box_volume, callback=rdf_callback)which causes the observable to be evaluated in each time step (stride=1) and print the value (callback=rdf_callback).The RDF is determined by calculating the distance from all particles of a type contained in types_count_from toall particles of a type contained in types_count_to and then binning the distance into a histogram as given by bin_borders.The histogram is normalized with respect to $g(r) = 4\\pi r^2\\rho dr$, where $\\rho$ is the number density of particles with types contained in types_count_to, reflected by particle_to_density.ParticlesThis observable records all particles in the system, as in: Each particle’s type, (unique) id, and position.It can be registered bydef particles_callback(particles): types, ids, positions = particles print(\"Particle 5 has type {}, id {}, and position {}.\" .format(types[5], ids[5], positions[5])simulation.observe.particles( stride=5, callback=particles_callback, save=False)where the argument of the callback function is a 3-tuple containing a list of types, unique ids, and positionscorresponding to each particle in the system. In this example the callback function prints these properties of thefifth particle every fifth time step, the output of the observable is not saved into the trajectory file (save=False).Particle positionsThe particles’ positions can be recorded bysimulation.observe.particle_positions( stride=200, types=None, callback=lambda x: print(x))which makes this observable very similar to the particles one, however one can select specific types of particlesthat are recorded. In case of types=None, all particle positions will be recorded, in case of types=[\"A\", \"B\"]only positions of particles with type A or B are returned.In this case the callback will simply print x every 200 steps, where x is a list of three-dimensional vectors.Since save is not explicitly set to False or None the observed data will be recorded into the trajectory fileif n simulation.output_file was configured.Number of particlesWhen one is only interested in the sheer number of particles then one can use this observable. Depending on the input,it will either observe the total number of particles or the number of particles per selected type:simulation.observe.number_of_particles( stride=1, types=[\"A\", \"B\", \"C\"], callback=lambda x: print(x), save=False)This example records the numbers of particles with types A, B, C in each time step. The callback takesa list with three elements as argument where each element corresponds to a particle type as given in types and contains the respective counts. If types=None was given, the observable would record the total number of particles,regardless of their types.EnergyThe system’s current potential energy can be observed bysimulation.observe.energy( stride=123, callback=lambda x: print(\"Potential energy is {}\".format(x)), save=False)where stride=123 indicates that the observable will be evaluated every 123rd time step. The argument of the callbackfunction is a scalar value and the observable’s output is not saved to a potentially configured trajectory file.ForcesThe forces acting on particles can be observed bysimulation.observe.forces( stride=1, types=[\"A\"], callback=lambda x: print(sum(x)))yielding an observable that is evaluated every time step (stride=1) and collects the forces acting on all particlesof type A. If types=None, all types are considered. The callback function takes a list of vectors as argument.Since save is not further specified, this observable would be recorded into the trajectory file.ReactionsThis observable records all occurred reactions in the system for a particular time step. It can be registered by invokingdef reactions_callback(reactions): for r in reactions: print(\"{} reaction {} occurred: {} -> {}, position {}\" .format(r.type, r.reaction_label, r.educts, r.products, r.position)) print(\"----\")simulation.observe.reactions( stride=5, callback=reactions_callback) where stride=5 indicates that the observable is evaluated every fifth time step. The callback takes a list ofreaction records as argument, where each reaction record stores information about the type of reaction (r.type), i.e., one of conversion, fission, fusion, enzymatic, decay, reaction name (r.reaction_label), i.e., the name under which the reaction was registered in the system, educt unique particle ids (r.educts) as in the particles observable, product unique particle ids (r.products), and position (r.position) of the reaction event which is evaluated to the midpoint between educts in case of a bimolecular reaction.Since the save argument was left out, it is defaulted and given a simulation.output_file, the observed reactions are recorded.Reaction countsInstead of recording all reaction events one can also record the number of occurred reactions per registered reaction per time step. This can be achieved bysimulation.observe.reaction_counts( stride=2, callback=lambda x: print(x), save=False)where stride=2 causes the observable to be evaluated every second time step. The callback function takesa dictionary as argument where the keys are the reaction names as given in the system configuration and the valuesare the occurrences of the corresponding reaction in that time step.VirialThis observable evaluates the pressure virial according to pair-wise forces by evaluating$$\\mathbf{W}_t = \\sum_{i > j} \\mathbf{r}_{ij} \\otimes\\mathbf{f}_{ij},$$where $\\mathbf{r}_{ij}$ is the vector difference between the positions of particle i and j, $\\mathbf{f}$ the pair-wise force.It can be registered bysimulation.observe.virial( stride=5, callback=lambda x: print(x))which causes it to be evaluated every fifth time step (stride=5). The callback function takes a numpy array of shape (3,3) as argument.PressureThe pressure of a system can be understood in terms of the trace of a pressure tensor$$P_t = \\frac{1}{3}\\,\\mathrm{tr}(\\mathbf{P}_t)$$which is defined via the virial tensor$$\\mathbf{P}_tV = N_tk_BT + \\mathbf{W}_t,$$where $V$ is the system’s volume, $N_t$ the number of (physical) particles that somehow partake in the dynamics of the system,and $k_BT$ the thermal energy.Observing the pressure in a simulation amounts tosimulation.observe.pressure( stride=1, physical_particles=[\"A\", \"B\", \"C\"], callback=lambda p: print(\"current pressure: {}\".format(p)))where stride=1 causes the pressure to be evaluated in every time step, physical_particles are set to be particles oftype A, B, or C (physical_particles=None causes all particles to be considered physical), and the callbackfunction takes a scalar value as argument.Internally, the pressure observable builds up on the virial observableand the number_of_particles observable and when writing it to filenot the actual pressure is recorded but the output of these other two observables.In the HDF5 group hierarchy the observable’s group is per default postfixed by _pressure, causing the virialto be stored under virial_pressure and the number of particles under n_particles_pressure.Changing the postfix amounts to providing a dictionary to the save argument withpressure_save_options = { 'name': \"empty_or_other_postfix\", 'chunk_size': 500}simulation.observe.pressure(1, [\"A\", \"B\", \"C\"], save=pressure_save_options)where the chunk_size refers to the HDF5 chunk size of data sets.",
"url": "/simulation.html#observables"
},
"simulation-run-html": {
"title": "Simulation > Running the simulation",
"content": "Per default, the simulation loop looks like belowThis means that all observables are evaluated at the initial state regardless of their stride, after whichthe actual loop is performed.A simulation is started by invokingsimulation.run(n_steps=1000, timestep=1e-5)where n_steps is the number of time steps to perform and the timestep is the time step. Per default an overviewof the system configuration is printed upon simulation start, this can be disabled by providing the argument show_system=False.One can influence portions of the loop through the simulation object: Per default a progress bar is shown when the simulation is executed, however it can be hidden print(simulation.show_progress)simulation.show_progress = False If one does not want to evaluate topology reactions at all, one can invoke simulation.evaluate_topology_reactions = False Evaluation of forces can be deactivated by invoking simulation.evaluate_forces=False Evaluation of observables can be switched off altogether by simulation.evaluate_observables=False Note that setting this to False also causes the trajectory to not be recorded. In case of a large simulation box but small interaction radii one can sometimes boost performance by artifically increasing the cuboid size of the neighbor list’s spatial discretization by setting simulation.skin = s where $s\\geq 0$ is a scalar that will be added to the maximal cutoff. Furthermore, one can select the Integrator Reaction handlerIntegratorCurrently the only available integrator is the EulerBDIntegrator which is selected by default and can be selected bya call tosimulation.integrator = \"EulerBDIntegrator\"It integrates the isotropic brownian dynamics$$\\frac{d\\mathbf{x}(t)}{dt} = -D\\frac{\\nabla V(\\mathbf{x}(t))}{k_BT} + \\xi(t)$$with an Euler-Maruyama discretization$$\\mathbf{x}_{t+\\tau} = \\mathbf{x}_t - \\tau D\\frac{\\nabla V(\\mathbf{x}(t))}{k_BT} + \\sqrt{2D\\tau}\\eta_t,$$where$$\\eta_t \\sim \\begin{pmatrix}\\mathcal{N}(0, 1) & \\cdots & \\mathcal{N}(0, 1) \\end{pmatrix}^T.$$Reaction handlerReactions in ReaDDy are evaluated per time step, meaning that each particle might have multiple possible reactionpartners. In order to solve this, one can chose between two different reaction handlers: The UncontrolledApproximation reaction handler is the less rigorous of the two. It performs as follows: A list of all possible reaction events is gathered. For each element of this list a decision is made whether it is to be evaluated based on the time step and its rate as described in the section about reactions. The filtered list is shuffled. For each event in the list evaluate it unless its educts have already been consumed by another reaction event. A problem of this reaction handler is that it does not give preference to reaction events based on their rate in caseof a conflict, i.e., two events that share educts. However in the limit of small time steps this problem disappears. This reaction handler can be selected by invoking simulation.reaction_handler = \"UncontrolledApproximation\" The Gillespie reaction handler is the default choice and is statistically exact in the sense that it imposes a Gillespie reaction order on the events of a particular time step: A list of all possible reaction events is gathered. Then Each reaction event is assigned its normalized reaction rate $\\lambda_i/\\sum_i\\lambda_i$ A random event is picked from the list based on its normalized reaction rate, i.e., a uniform random number between 0 and 1 is drawn and then used together with the cumulative normalized reaction rates to determine the event Depending on its rate the reaction described by the event might be performed: if not the event is simply removed from the list if it was performed it is also removed and any other event that shared educts with this particular event if there are still events in the list go back to 1., otherwise stop An example of conflicting reaction events with expected outcome might be $$\\left\\{ \\begin{array}{rcc} A + B &\\xrightarrow{\\lambda_1}& C\\\\ A &\\xrightarrow{\\lambda_2}& D \\end{array} \\right. \\xrightarrow{\\lambda_1 \\gg \\lambda_2} \\left\\{ \\begin{array}{rcc} A + B &\\xrightarrow{\\lambda_1}& C\\\\ &\\mathrm{ignored} \\end{array} \\right.$$ This reaction handler can be selected by invoking simulation.reaction_handler = \"Gillespie\" ",
"url": "/simulation.html#simulation_run"
},
"system-external-potentials-html": {
"title": "System configuration > External potentials",
"content": "External potentials or first-order potentials are potentials that solely depend on the absolute position of each particle, i.e., the relative positioning of particles towards one another has no influence.They are registered with respect to a certain particle type. The potential willthen exert a force on each particle of that type individually.Currently available external potentials are: Box potential Spherical potential Spherical exclusion Spherical inclusion Spherical barrier Box potentialA box potential is a potential acting with a harmonic force on particles of the given type once they leave the area spanned by the cuboid that has origin as its front lower left and origin+extent as its back upper right vertex, respectively. Therefore, the cuboid is spanned by three intervals $C=\\prod_{i=1}^dC_i$ with $C_i := [\\text{origin}_i, \\text{origin}_i+\\text{extent}_i]$. Its energy term is given by$$V(\\mathbf{x}) = \\sum_{i=1}^d \\begin{cases} 0&\\text{, if } x_i \\in C_i\\\\ \\frac{1}{2}k \\,d(x_i, C_i)^2 &\\text{, otherwise,} \\end{cases}$$where $d(x_i, C_i)$ denotes the shortest distance between the set $C_i$ and the point $x_i$.Adding a box potential to the system amounts to:system.box_size=[3, 3, 3]system.potentials.add_box( particle_type=\"A\", force_constant=10., origin=[-1, -1, -1], extent=[2, 2, 2])Note that the simulation box and the box potential are completely independent.In the above example the simulation box is chosen larger than the full extent of the box potential. This is becauseparticles should never leave the simulation box if it is non-periodic. The box potential however is a soft potential,i.e., particles may penetrate the boundaries of it for a short time and then be pushed back inside. To make sure thatparticles do not penetrate the simulation box, it has a slightly larger extent.In particular there is a check upon simulation start that if the simulation box is not completely periodic, there must be a box potential for each particle type to keep it contained in the non-periodic directions, i.e., if there is no box potential such thatbox_lower_left[dim] < potential_lower_left[dim] and box_upper_right[dim] > potential_upper_right[dim]where dim is a non-periodic direction, an error is raised.Spherical potentialIn ReaDDy one can find three different types of spherical potentials: exclusion potentials to keep particles out of a spherical region, inclusion potentials to keep particles inside a spherical region, barriers that can function as a spatial separator or, if initialized with negative height, as a sticky spherical surface.All these potentials are harmonic, i.e., particles can potentially penetrate.Spherical exclusionAdds a spherical potential that keeps particles of a certain type excluded from the inside of the specified sphere. Its energy term is given by$$V(\\mathbf{x}) = \\begin{cases}\\frac{1}{2} k (\\|\\mathbf{x} - \\mathbf{c}\\|_2-r)^2 &\\text{, if } \\|\\mathbf{x} - \\mathbf{c}\\|_2 where $\\mathbf{c}\\in\\mathbb{R}^3$ denotes the center of the sphere and $r\\in\\mathbb{R}_{>0}$ the radius of the sphere.Adding such a potential to a reaction diffusion system amounts tosystem.box_size = [3, 3, 3]system.potentials.add_sphere_out( particle_type=\"A\", force_constant=10., origin=[0, 0, 0], radius=1.)yielding a spherical region of radius 1 in the center of the simulation box which keeps particles of type A from entering that region with a harmonic repulsion potential. Due to the harmonic nature of the potential and dependent on the force constant, particles are allowed to penetrate the sphere for a short period of time.Spherical inclusionAdds a spherical potential that keeps particles of a certain type restrained to the inside of the specified sphere. Its energy term is given by$$V(\\mathbf{x}) = \\begin{cases}\\frac{1}{2} k (\\|\\mathbf{x} - \\mathbf{c}\\|_2-r)^2 &\\text{, if } \\|\\mathbf{x} - \\mathbf{c}\\|_2 \\geq r,\\\\0 &\\text{, otherwise,} \\end{cases}$$where $\\mathbf{c}\\in\\mathbb{R}^3$ denotes the center of the sphere and $r\\in\\mathbb{R}_{>0}$ the radius of the sphere.Adding such a potential to a system amounts tosystem.box_size = [3, 3, 3]system.potentials.add_sphere_in( particle_type=\"A\", force_constant=10., origin=[0, 0, 0], radius=1.)which will cause all particles of type A to be contained in a sphere of radius 1 centered in the origin with a harmonic repulsion potential. Due to the harmonic nature of the potential and dependent on the force constant, particles are allowed to penetrate the sphere for a short period of time.Spherical barrierA potential that forms a concentric barrier at a certain radius around a given origin. It is given a height (in terms of energy) and a width. Note that the height can also be negative, then this potential acts as a ‘sticky’ sphere. The potential consists of harmonic snippets, such that the energy landscape is continuous and differentiable, the force is only continuous and not differentiable. Its energy term is given by$$V(\\mathbf{x}) = \\begin{cases}0 & \\text{, if } \\| \\mathbf{x} - \\mathbf{c}\\|_2 where $\\mathbf{c}\\in\\mathbb{R}^3$ is the center of the sphere, $r\\in\\mathbb{R}$ the sphere’s radius, $w\\in\\mathbb{R}$ the width of the barrier, and $h\\in\\mathbb{R}$ the height of the barrier.Adding such a potential to a system amounts tosystem.box_size = [3, 3, 3]# as a barriersystem.potentials.add_spherical_barrier( particle_type=\"A\", height=1.0, width=0.1, origin=[0, 0, 0], radius=1.)# stickysystem.potentials.add_spherical_barrier( particle_type=\"A\", height=-1.0, width=0.1, origin=[0, 0, 0], radius=1.)which will cause particles to be trapped inside or outside of the spherical barrier in the first case or will make them stick onto the surface of the sphere in the second case.",
"url": "/system.html#external_potentials"
},
"system-pair-potentials-html": {
"title": "System configuration > Pair potentials",
"content": "Pair potentials or second-order potentials are potentials that depend on the relative positioning of a pair of particles to one another. They are registered with respect to two not necessarily distinct particle types and then exert a force on each particle of these types individually.Currently available pair potentials are: Harmonic repulsion Weak interaction piecewise harmonic Lennard-Jones Screened electrostaticsHarmonic repulsionA harmonic repulsion potential causes each two particles of a certain type to repulse each other once they enter a certain radius. It can be used to emulate a radius of a particle type while still allowing for a relatively large time step. The potential term is given by$$V(\\mathbf{x}_1, \\mathbf{x}_2) = \\begin{cases}\\frac{1}{2}k(\\|\\mathbf{x_1} - \\mathbf{x_2}\\|_2 - r)^2 &\\text{, if } \\|\\mathbf{x_1} - \\mathbf{x_2}\\|_2 where $r$ is the distance at which particles begin to interact with respect to this potential.Adding such a potential to a system amounts to, e.g.,system.potentials.add_harmonic_repulsion( \"A\", \"A\", force_constant=10., interaction_distance=5.)system.potentials.add_harmonic_repulsion( \"B\", \"B\", force_constant=10., interaction_distance=6.)system.potentials.add_harmonic_repulsion( \"A\", \"B\", force_constant=10., interaction_distance=2.5+3.)which would cause all particles of type A or B to repulse from all other particles of type A or B. This set of potentials can be understood as inducing a “soft” radius of $r_A = 2.5$ and $r_B=3$ for particle types A and B, respectively. Soft meaning that the particles’ spheres can be penetrated for a short period of time by another particle before it is pushed out again.Weak interaction piecewise harmonicA weak interaction piecewise harmonic potential causes particles to crowd together once they are within a certain distance of one another. It is defined by three harmonic terms and described by a force constant $k$ which is responsible for the strength of the repulsive part of the potential, a ‘desired distance’ $d$, i.e., a distance at which the potential energy is lowest inside the interaction radius, a ‘depth’ $h$, denoting the depth of the potential well and therefore the stickiness of the potential, and a ‘cutoff’ $r_c$, denoting the distance at which particles begin to interact. The potential term reads$$V(\\|\\mathbf{x}_1- \\mathbf{x}_2\\|_2) = V(r) = \\begin{cases}\\frac{1}{2} k (r - d)^2 - h &\\text{, if } r It typically has the following shape:Adding such a potential to a system can be achieved by callingsystem.potentials.add_weak_interaction_piecewise_harmonic( \"A\", \"B\", force_constant=10., desired_distance=5., depth=2., cutoff=7.)yielding in this example a potential that lets all particle type pairings interact with one another given they are of type A and B and closer than the cutoff $r_c=7$. Once they are within range they would either by drawn into the potential well from the outside (third case in the potential term), kept in the potential well (second case in the potential term), or pushed back into the potential well if they come too close to one another (first case in the potential term).Lennard-JonesSimilarly to a weak interaction potential, the Lennard-Jones potential models the interaction between a pair of particles. However it is not a soft potential and therefore requires a relatively small time step in order to function correctly. The potential term reads$$V_\\text{LJ}(\\|\\mathbf{x_1}-\\mathbf{x_2}\\|_2) = V_\\text{LJ}(r) = k(\\varepsilon, n, m)\\left[ \\left(\\frac{\\sigma}{r}\\right)^m - \\left(\\frac{\\sigma}{r}\\right)^n \\right],$$where $k(\\varepsilon, n, m)\\in\\mathbb{R}$ is the force constant, $\\sigma\\in\\mathbb{R}$ the distance at which the inter-particle potential is zero, $\\varepsilon\\in\\mathbb{R}$ the depth of the potential well, i.e., $V_\\text{LJ}(r_\\text{min})=-\\varepsilon$, and $m,n\\in\\mathbb{N}, m>n$ exponents which determine the stiffness and range of the potential. For the classical Lennard-Jones potential the exponents are given by $m=12$ and $n=6$.The potential itself approaches but never reaches zero beyond the interaction well. Therefore, a cutoff is introduced, usually at $r_c=2.5\\sigma$ (for the 12-6 LJ potential), which is the point at which the potential as roughly $1/60$th of its minimal value $-\\varepsilon$. This however leads to a jump discontinuity at $r_c$ in the energy landscape, which can be avoided by shifting the potential by the value at the discontinuity:$$V_{\\text{LJ}_\\text{trunc}}(r) = \\begin{cases} V_\\text{LJ}(r) - V_\\text{LJ}(r_c) &\\text{, if } r \\leq r_c,\\\\ 0&\\text{, otherwise.} \\end{cases}$$The force constant $k$ is chosen such that the depth at the well is is $V(r_\\mathrm{min}) = -\\varepsilon$, i.e.,$$k = -\\frac{\\varepsilon}{\\left( \\frac{\\sigma}{r_\\mathrm{min}} \\right)^m - \\left( \\frac{\\sigma}{r_\\mathrm{min}} \\right)^n}.$$Different choices of exponents that can be found in the literature are, e.g., $9-3$, $9-6$, or $8-6$.Adding such a potential to a system can be achieved by callingsystem.potentials.add_lennard_jones( \"A\", \"B\", m=12, n=6, cutoff=2.5, shift=True, epsilon=1.0, sigma=1.0))yielding a truncated 12-6 Lennard-Jones potential between particles of type A and B with a zero inter-particle potential at $\\sigma=2.5$, a well depth of $\\varepsilon=1.0$, and a cutoff radius of $r_c=2.5 = 2.5\\cdot\\sigma$. If the shift in the energy landscape to avoid the jump discontinuity is not desired, it can be switched off by setting shift=False.Screened electrostaticsThe “screened electrostatics” (also: Yukawa or Debye-Hückel) potential is a potential that represents electrostatic interaction (both repulsive or attractive), which is screened with a certain screening depth. This kind of potential becomes important when dealing with particles representing proteins that have a net-charge. However, it is usually more expensive to evaluate than, e.g., harmonic repulsion, as it requires a larger cutoff and smaller time step. If the electrostatic term is attractive, a core repulsion term is added. The potential term reads$$V(\\|\\mathbf{x_1}-\\mathbf{x_2}\\|_2) = V(r) = \\begin{cases}C\\frac{\\exp(-\\kappa r)}{r} + D\\left( \\frac{\\sigma}{r} \\right)^n &\\text{, if } r \\leq r_c,\\\\0 &\\text{, otherwise,}\\end{cases}$$where $C\\in\\mathbb{R}$ is the electrostatic repulsion strength in units of energy times distance, $\\kappa\\in\\mathbb{R}$ is the inverse screening depth in units of 1/length, $D\\in\\mathbb{R}$ is the repulsion strength in units of energy, $\\sigma\\in\\mathbb{R}$ is the core repulsion radius or zero-interaction radius in units of length, $n\\in\\mathbb{N}$ is the core repulsion exponent (dimensionless), and $r_c\\in\\mathbb{R}$ the cutoff radius in units of length. It typically has the following shape:One can observe that the first term in the potential’s definition diverges towards $-\\infty$ for $r\\searrow 0$ which is an effect that gets countered by the second term, diverging towards $+\\infty$ for $r\\searrow 0$. The depth of the potential well increases with an increasing exponent $n$.Adding such a potential to a system amounts tosystem.potentials.add_screened_electrostatics( \"A\", \"B\", electrostatic_strength=-1., inverse_screening_depth=1., repulsion_strength=1., repulsion_distance=1., exponent=6, cutoff=6.)which would introduce an electrostatic interaction between particles of type A and B that resembles the above plot.",
"url": "/system.html#pair_potentials"
},
"system-potentials-html": {
"title": "System configuration > Potentials",
"content": "Potentials create an energy landscape in which particles diffuse in, subject to the corresponding forces.They can be used to build traps, obstacles or compartments for particles.One could also utilize them for clustering and crowding effects that are typically observed in biological fluid-like media.Potentials in ReaDDy are divided into first-order potentials/external potentials, i.e., those that depend only on the position of one particle, andsecond-order potentials/pair potentials, i.e., those that depend on the relative position of two particles.The topology functionality also provides higher order potentials like angles and dihedrals.All potentials are part of the ReactionDiffusionSystem and can be registered for certain particle types likesystem = readdy.ReactionDiffusionSystem()system.potentials.add(...)",
"url": "/system.html#potentials"
},
"system-reactions-html": {
"title": "System configuration > Reactions",
"content": "Reactions remove particles from, and add particles to the system. They typically have a microscopic/intrinsic rate $\\lambda$.This rate has units of inverse time and can be understood as the probability per unit time of the reaction occurring. Given an integrationstep $\\tau$ the probability of a reaction event is evaluated as $p = 1 - e^{-\\lambda \\tau}$.Additionally Fusion and Enzymatic reactions can only occur when particles are closer than a certain distance $R$.All reactions are added to the reaction registry, which is part of the ReactionDiffusionSystemsystem = readdy.ReactionDiffusionSystem()system.reactions.add(...)Each of the below listed reaction types can be registered in two ways: Either with the generic reactions.add(...) method which accepts a descriptor string, or by calling reactions.add_xxx(...), where xxx is to be replaced with one of conversion, decay, fusion, fission, or enzymatic.ConversionIn a conversion reaction, a particle of type $A$ will convert into type $B$ with the rate constant $\\lambda$$$A \\overset{\\lambda}{\\rightarrow} B$$Adding a conversion reaction to the system amounts tosystem.reactions.add_conversion(name=\"conv\", type_from=\"A\", type_to=\"B\", rate=0.1)which is equivalent tosystem.reactions.add(\"conv: A -> B\", rate=0.1)DecayIn a decay reaction, a particle of type $A$ will disappear with the rate constant $\\lambda$$$A \\overset{\\lambda}{\\rightarrow} \\emptyset$$Example of adding a decay reaction to the system:system.reactions.add_decay(name=\"decay of A\", particle_type=\"A\", rate=0.1)which is equivalent tosystem.reactions.add(\"decay of A: A ->\", rate=0.1)FusionIn a fusion reaction, a particle of type $A$ will associate with type $B$ to form a particle of type $C$.This will occur with the rate constant $\\lambda$, if the particles $A$ and $B$ are closer than the reaction radius$R$ (educt_distance).$$A \\overset{R}{+} B \\overset{\\lambda}{\\rightarrow} C$$Example of adding a fusion reaction to the system:system.reactions.add(\"fus: A +(2) B-> C\", rate=0.1)which is equivalent tosystem.reactions.add_fusion( name=\"fus\", type_from1=\"A\", type_from2=\"B\", type_to=\"C\", rate=0.1, educt_distance=2.)FissionIn a fission reaction, a particle of type $C$ will dissociate into two particles of type $B$ and $A$.This will occur with the rate constant $\\lambda$. The two products will be placed at a distance smaller thanthe reaction radius $R$ (product_distance).$$C \\overset{\\lambda}{\\rightarrow} A \\overset{R}{+} B$$Add a fission reaction to the systemsystem.reactions.add(\"fis: C -> A +(2) B\", rate=0.1)which is equivalent tosystem.reactions.add_fission( name=\"fis\", type_from=\"C\", type_to1=\"A\", type_to2=\"B\", rate=0.1, product_distance=2.)EnzymaticIn an enzymatic reaction, a particle of type $A$ convert into type $B$ in the presence of a particle of type $C$.This will occur with the rate constant $\\lambda$, if the particles $A$ and $C$ are closer than the reaction radius$R$ (educt_distance).$$A \\overset{R}{+} C \\overset{\\lambda}{\\rightarrow} B + C$$Add an enzymatic reaction to the systemsystem.reactions.add(\"enz: A +(2) C -> B + C\", rate=0.1)which is equivalent tosystem.reactions.add_enzymatic( name=\"enz\", type_catalyst=\"C\", type_from=\"A\", type_to=\"B\", rate=0.1, educt_distance=2.)",
"url": "/system.html#reactions"
},
"system-topologies-html": {
"title": "System configuration > Topologies",
"content": "Topologies are a way to include molecular structure in a reaction-diffusion simulation. More specifically, a topology is a group of particles with fixed potential terms like bonds and angles between certain particles. Topologies in ReaDDy consist of two ingredients: A connectivity graph, where the vertices of the graph correspond to the particles in the topology, a lookup table for potential terms between certain combinations of particle types.How to set up the actual connectivity graph can be found in the section about setting up and running the simulation, as it requires particles being added to the simulation.Since particles that are part of a topology are internally treated differently than “normal” particles, their species have to be configured by a call tosystem.add_topology_species(\"T\", diffusion_constant=2.0)Furthermore, for operations that function on the topology level, topologies have a “topology type” which can be seen as the generalization of a “particle type”. To add such a type, one can invokesystem.topologies.add_type(\"My topology type\")For a topology to be recognized as “valid”, both of the following conditions need to be fulfilled: The connectivity graph needs to be connected, i.e., there must not be independent components. Each edge in the connectivity graph needs to have a corresponding bond configured based on the respective particle types.If one of these conditions is not fulfilled, an exception is raised and the simulation will not start.",
"url": "/system.html#topologies"
},
"system-topologies-potentials-html": {
"title": "System configuration > Potentials",
"content": "Topologies are defined by a set of particles that are connected with a graph and a lookup table that defines what connectivities between what particle types yield which potentials.This section deals with the latter, i.e., with the lookup table. The lookup table is independent of the topology type, so all potentials that are defined here will be applied to pairs/triples/quadruples of particles which are connected in the respective topologies connectivity graph.In this picture the dashed lines denote the connectivity graph between the particles, the blue lines bond potentials, the green lines angle potentials, and the orange lines dihedral potentials. One can see that bonds are defined on pairs of particles, angles on triples, dihedrals on quadruples. In this particular case one has Bonds Angles Dihedrals $A\\leftrightarrow A$ $C\\leftrightarrow B\\leftrightarrow A$ $A\\leftrightarrow A\\leftrightarrow B \\leftrightarrow A$ $A\\leftrightarrow B$ $A\\leftrightarrow C$ In an actual instance of a topology one would also have to define a bond between particles of type $C\\leftrightarrow C$ or remove that edge from the graph, otherwise it would not be considered valid.ReaDDy supports three types of potentials within topologies: Harmonic bonds Harmonic angles Cosine dihedralsHarmonic bondsHarmonic bonds model, e.g., covalent bonds in a molecular structure. The potential term yields forces that push pairs of particles away from one another if they become closer than a certain distance and attracts them if they are further apart than that distance. It reads$$V(\\|\\mathbf{x}_1-\\mathbf{x}_2\\|_2) = V(r) = k(r-r_0)^2,$$where $r_0$ is the preferred distance and $k$ the force constant.Adding such a potential term to a system amounts to, e.g.,system.add_topology_species(\"T1\", diffusion_constant=2.)system.add_topology_species(\"T2\", diffusion_constant=4.)system.topologies.configure_harmonic_bond( \"T1\", \"T2\", force_constant=10., length=2.)which would have the effect of introducing for each topology a harmonic bond with force constant 10 and preferred distance 2 between each adjacent pair of particles with types “T1” and “T2”, respectively.Harmonic anglesHarmonic angles are potential terms that yield a preferred configuration for a triple of particles in terms of the spanned angle $\\theta$.Should the spanned angle be different from the preferred angle $\\theta_0$, a harmonic force acts on each of the particles toward the preferred angle. The potential energy term reads$$V(\\theta_{i,j,k}) = k(\\theta_{i,j,k} - \\theta_0)^2,$$where $\\theta_{i,j,k}$ corresponds to the angle spanned by three particle’s positions $\\mathbf{x}_i, \\mathbf{x}_j, \\mathbf{x}_k$ and $k$ is the force constant.Configuring such a potential for a system amounts to, e.g.,system.add_topology_species(\"T1\", diffusion_constant=2.)system.add_topology_species(\"T2\", diffusion_constant=4.)system.add_topology_species(\"T3\", diffusion_constant=4.)system.topologies.configure_harmonic_angle( \"T1\", \"T2\", \"T3\", force_constant=1., equilibrium_angle=3.141)yielding harmonic angle potential terms for each triple of particles with types (T1, T2, T3) (or equivalently types (T3, T2, T1)) that are contained in a topology and have edges between the particles corresponding to types (T1, T2) and (T2, T3), respectively.Cosine dihedralsThe sketch above depicts the definition of the proper dihedral angle $\\phi$ spanned by four particles with corrdinates $\\mathbf{x}_i$, $\\mathbf{x}_j$, $\\mathbf{x}_k$, and $\\mathbf{x}_l$, respectively. The potential energy is given by$$V(\\phi) = k(1+\\cos (n\\phi - \\phi_0)),$$where $\\phi_0\\in [ -\\pi,\\pi ]$ represents the offset angle, $k$ the force constant in units of energy/angle**2, and $n\\in\\mathbb{N}_{>0}$ the multiplicity, indicating the number of minima encountered when rotating the bond through $360^\\circ$.The $i$-th preferred particle configuration, i.e. the $i$-th minimum of $V$ are $\\phi_i = \\frac{1}{n}\\left( \\frac{\\pi}{2} - \\phi_0 + i\\pi \\right)$Configuring such a potential for a system amounts to, e.g.,system.add_topology_species(\"T1\", diffusion_constant=2.)system.add_topology_species(\"T2\", diffusion_constant=4.)system.add_topology_species(\"T3\", diffusion_constant=4.)system.add_topology_species(\"T4\", diffusion_constant=4.)system.topologies.configure_cosine_dihedral( \"T1\", \"T2\", \"T3\", \"T4\", force_constant=10, multiplicity=1., phi0=0.)yielding cosine dihedral potentials for each path of length 4 with vertices corresponding to particles of types (T1, T2, T3, T4) in the connectivity graph of a topology instance. The sequence in which the types are given can be reversed, i.e., system.topologies.configure_cosine_dihedral( \"T4\", \"T3\", \"T2\", \"T1\", force_constant=10, multiplicity=1., phi0=0.)results in the same potential terms.Angles are internally always expressed in radians.",
"url": "/system.html#topology_potentials"
},
"system-topologies-reactions-html": {
"title": "System configuration > Topology reactions",
"content": "Topology reactions provide means to change the structure of a topology during the course of a simulation. Changing the structure can involve: Changing particle types of particles inside topologies and therefore changing the force field, breaking and forming bonds inside a topology resulting in different connectivity or the separation of a topology in separated instances, attaching particles to topologies, connecting two topologies by introducing an edge between their graphs.These changes are divided into two types: Structural reactions and spatial reactions, where, as the name suggests, structural reactions change the internal structure of a topology and are independent of the actual spatial configuration of the system and spatial reactions represent local changes of the graph triggered by spatial events, i.e., attaching particles or forming bonds between two topology instances. The following sections are ordered accordingly: Structural reactions The reaction function The rate function Adding a structural reaction Spatial reactions Predefined reactions Topology dissociation Structural reactionsStructural topology reactions are defined on the topology type. They basically consist out of two functions: the reaction function, taking a topology object as input and returning a reaction recipe describing what the structural changes are to be applied to the topology the rate function, which takes a topology object as input and returns a corresponding fixed rate.The rate function is evaluated initially and then only when the topology has changed due to other reactions. The reaction function is only evaluated when a topology reaction is performed. It should be noted that these function evaluations can have negative performance impacts when occurring frequently.The above figure shows an example of a structural topology reaction. In the upper row there are particles $i,j,k,l$ from left to right with a graph that connects the pairs $(i,j)$, $(j,k)$, and $(k,l)$. Due to this adjacency, there are bonds $b_{01}$, $b_{12}$, and $b_{23}$ as well as angles $a_{012}$ and $a_{123}$. The lower row represents the configuration after a topology reaction that removed the edge $(k,l)$. In its absence the bond $b_{12}$ and the angle $a_{123}$ are removed and the topology originally consisting out of four particles decays into two topologies - one with three particles and one trivial topology containing just one particle.The reaction functionIn order to configure such a reaction for a reaction diffusion system, one first needs to set up a function that given a topology returns an instance of StructuralReactionRecipe, essentially describing the desired changes in structure:def no_op_reaction_function(topology): recipe = readdy.StructuralReactionRecipe(topology) return recipeOne can base the behavior of the reaction on the current state of the topology instance. It offers information about the contained particles configuration: topology.get_graph() or topology.graph yields the connectivity graph of the topology: graph.get_vertices() or graph.vertices yields a list of vertices that has a 1-1 correspondence to what is yielded by topology.particles. Each vertex itself can be iterated, yielding its adjacent vertices: # for every vertexfor v in graph.vertices: print(\"vertex {}\".format(v.particle_index)) # obtain number of its neighbors' neighbors n_neighbors_neighbors = 0 for neighbor in v: for neighbor_neighbor in neighbor.get(): n_neighbors_neighbors += 1 graph.get_edges() or graph.edges yields a list of edges contained in the graph, where each edge is represented by a tuple of vertices: for e in graph.get_edges(): v1, v2 = e[0], e[1] print(\"edge {} -- {}\".format(v1.particle_index, v2.particle_index)) topology.position_of_vertex(v) yields the position of the particle corresponding to the provided vertex, topology.particle_type_of_vertex(v) yields the type of the particle corresponding to the provided vertex, topology.particle_id_of_vertex(v) yields the unique id of the particle corresponding to the provided vertex.With these information, there are several operations that can be added to a recipe: recipe.change_particle_type(vertex_index, type): Changes the particle type of the to vertex_index associated particle to the given type, where the vertex index corresponds to the particle’s index. Can also be called with vertex instance directly. recipe.add_edge(v_index1, v_index2): Introduces an edge in the graph between vertices corresponding to indices v_index1 and v_index2. Can also be called with vertex instances directly. recipe.remove_edge(v_index1, v_index2): Attempts to remove an edge between vertices corresponding to the indices. Depending on the configuration of the topology reaction, this can lead to a failed state or multiple sub-topologies. Can also be called with vertex instance directly. recipe.remove_edge(edge): Same as with indices just that it takes an edge instance as contained in graph.get_edges(). recipe.separate_vertex(index): Removes all edges from the topology’s graph that contain the vertex corresponding to the provided index (or vertex instance). If no new edge is formed between the given vertex this call, depending on the configuration of the reaction, can lead to a failed state or to formation of a topology consisting out of only one particle. In the latter case, this call can be followed by a call to recipe.change_particle_type, where the target type is no topology type. Then, no one-particle topology will be formed but the particle will simply be emitted and treated as normal particle. recipe.change_topology_type(type): Changes the type of the topology to the given type, potentially changing its structural and spatial topology reactions. recipe.append_particle(list_of_neighbor_vertices, particle_type, position): Adds a new particle to the topology. It will have the specified particle_type and position and will be connected to the vertices specified in list_of_neighbor_vertices.The rate functionSame as ordinary reactions, also structural topology reactions have a rate with which they occur. This rate is microscopic, has units of inverse time and can be understood as the probability per unit time of the reaction taking place. Same as for normal reactions, the probability is evaluated as $p=1-e^{-\\lambda\\tau}$, where $\\lambda\\geq0$ is the rate and $\\tau$ the integration time step.In order to define a rate for a certain structural reaction, one needs to provide a rate function:def my_rate_function(topology): n = len(topology.get_graph().get_vertices()) if n > 3: return .5 * n else: return 20.The function takes a topology instance as argument and returns a floating point value representing the rate in terms of the magnitude w.r.t. the default units. In this example it returns half the number of vertices if the number of vertices is larger than three, otherwise a constant value of 20.For performance reasons it is only evaluated if the topology changes structurally, therefore the rate should optimally not depend on anything that can change a lot in the simulation time between evaluating the rate function and performing the reaction, e.g., the particles’ spatial configuration.Adding a structural reactionGiven these two functions, the reaction and the rate function, all that is left to do is to add them to a certain topology type in the system:system.topologies.add_structural_reaction( name=\"my_structural_reaction\", topology_type=\"TType\", reaction_function=no_op_reaction_function, rate_function=my_rate_function, raise_if_invalid=True, expect_connected=False)where name is a user-provided unique identifier for this reaction, which is used in observing reaction counts. The above snippet adds the structural reaction my_structural_reaction to all topologies of type TType with the provided reaction function and rate function.The option raise_if_invalid raises, if set to True, an error if the outcome of the reaction function is invalid, otherwise it will just roll back to the state of before the reaction and print a warning into the log.The other option expect_connected can trigger depending on the value of raise_if_invalid a raise if set to True and the topology’s connectivity graph decayed into two or more independent components.Spatial reactionsSpatial reactions are locally triggered by proximity of particles, therefore they are not only defined on topology types but also on particle types. In principle there are two kinds of spatial reactions: The kind that causes forming a bond between two particles and the kind that just changes a particle/topology type, corresponding to particle fusion and enzymatic reactions, respectively. Analogously spatial topology reactions also possess a rate determining how likely the reaction occurs per time step as well as a radius determining the volume which is scanned for potential reaction partners.To simplify the definition of these reactions a descriptor language is used, deciding about the nature of the reaction. It consists out of a label and a topology-particle-type reaction equation:$$\\begin{aligned} \\mathrm{label\\_enzymatic: }& T_1(P_1)+T_2(P_2)\\rightarrow T_3(P_3) + T_4(P_4)\\\\ \\mathrm{label\\_fusion: } & T_1(P_1)+T_2(P_2)\\rightarrow T_3(P_3\\mathrm{--}P_4)\\end{aligned}$$where $T_i$ denote topology types, $P_i$ denote particle types. The first reaction is of “enzymatic” type, changing the types of particles corresponding to $P_1$ to $P_3$ and $P_2$ to $P_4$ if they are contained in topologies of type $T_1$ and $T_2$ which are also changed to $T_3$ and $T_4$, respectively. The second reaction is of “fusion” type, merging two topologies of types $T_1$ and $T_2$ by forming a bond between a particle pair with types $P_1$ and $P_2$ in their respective topologies. The result is a topology of type $T_3$ and the particles between which the bond was formed are of types $P_3$ and $P_4$.Some of these reactions can also be performed with a topology and a free particle. In particular, the following types of reactions are possible: TT-Fusion: T1(p1)+T2(p2) -> T3(p3--p4): a fusion of a topology of type T1 with a topology of type T2 by forming a bond between a pair of particles with types p1 and p2, where the product is a topology of type T3 and the newly connected particles changed their types to p3 and p4, respectively. TT-Fusion-self: T1(p1)+T1(p2) -> T3(p3--p4) [self=true]: a fusion of two topologies of type T1 similarly to the first type with the difference that now also particles within the same topology can be reaction partners. TP-Fusion: T1(p1)+(p2) -> T2(p3--p4): attaching a free particle of type p2 to a topology of type T1 if it is close to a particle of type p1 in that topology, yielding a topology of type T2 in which the newly connected particles are now of type p3 and p4, respectively. TT-Enzymatic: T1(p1)+T2(p2) -> T3(p3)+T4(p4): not changing the structure of the graph of the reaction partners but changing particle types possibly locally influencing the force field and changing topology types possibly leading to different topology reactions. TP-Enzymatic: T1(p1)+(p2) -> T2(p3)+(p4): same as the TT-Enzymatic reaction just that here it is performed with one topology and one free particle.Adding such a reaction to a system amounts to, e.g.,system.topologies.add_spatial_reaction( 'TT-Fusion: T1(p1)+T2(p2) -> T3(p3--p4)', rate=1., radius=1.)where the rate is in units of 1/time and the radius is a length. The descriptor is always the first argument and can be of any of the above discussed types.It should be noted that while usually “normal” particle-particle reactions are not possible with topology-typed particles, one can define enzymatic reactions where the catalyst is a topology type as this leaves the topology untouched and therefore can be evaluated in the normal reaction procedure.system.add_species(\"A\", diffusion_constant=1.)system.add_species(\"B\", diffusion_constant=1.)system.add_topology_species(\"P\", diffusion_constant=.5)system.topologies.add_type(\"T\")# OK, this attaches the particle A to the topologysystem.topologies.add_spatial_reaction('label1: T1(P)+(A)->T1(P--P)')# Fails, A is not a topology species typesystem.topologies.add_spatial_reaction('label1: T1(P)+(A)->T1(P--A)')# Fails, P is not a normal particle typesystem.topologies.add_spatial_reaction('label1: T1(P)+(P)->T1(P--P)')# OK, this is a normal fusion reactionsystem.reactions.add('A +(2.) A -> A', rate=.1)# Fails, P is not a normal particle type but a topology particle typesystem.reactions.add('A +(2.) P -> A', rate=.2)# OK, this is the special case where P is the catalystsystem.reactions.add('A +(2.) P -> B + P', rate=.3)Predefined reactionsFor convenience there are predefined topology reactions that can be added to a certain topology type.Topology dissociationThis reaction causes a topology to break bonds with a rate of n_edges*bond_breaking_rate, causing it to dissociate. In this process it may decompose into multiple independent components of the same topology type. Consequently, each of these independent components again has a topology dissociation reaction.Adding such a reaction to a system amounts tosystem.topologies.add_type(\"T\")system.topologies.add_topology_dissociation(\"T\", bond_breaking_rate=10.)",
"url": "/system.html#topology_reactions"
},
"validation-bimolecular-reaction": {
"title": "Validation > Bimolecular reaction - well mixed",
"content": "In [1]: import osimport numpy as npimport matplotlib.pyplot as pltimport readdy In [2]: readdy.__version__ Out[2]:'v2.0.3-27'Setup ReaDDy systemIn [3]: system = readdy.ReactionDiffusionSystem([20.,20.,20.], temperature=300.*readdy.units.kelvin)system.add_species("A", diffusion_constant=1.0)system.add_species("B", diffusion_constant=1.0)system.add_species("C", diffusion_constant=1.0)lambda_on = 1.system.reactions.add("myfusion: A +(1) B -> C", rate=lambda_on/readdy.units.nanosecond) Simulate the systemIn [4]: simulation = system.simulation(kernel="CPU")simulation.output_file = "out.h5"simulation.reaction_handler = "Gillespie"n_particles = 2000initial_positions_a = np.random.random(size=(n_particles, 3)) * 20. - 10.initial_positions_b = np.random.random(size=(n_particles, 3)) * 20. - 10.simulation.add_particles("A", initial_positions_a)simulation.add_particles("B", initial_positions_b)simulation.observe.number_of_particles(stride=1, types=["A"]) In [5]: if os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.run(n_steps=5000, timestep=1e-3*readdy.units.nanosecond) 0%| | 2/500 [00:00<00:39, 12.64it/s] Configured kernel context with:-------------------------------- - kBT = 2.49434 - periodic b.c. = (true, true, true) - box size = (20.0, 20.0, 20.0) - particle types: * particle type "C" with D=1.0 * particle type "A" with D=1.0 * particle type "B" with D=1.0 - bimolecular reactions: * Fusion A + B -> C with a rate of 1.0, an educt distance of 1.0, and weights 0.5 and 0.5Configured simulation loop with:-------------------------------- - timeStep = 0.00100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 500/500 [00:17<00:00, 28.56it/s]In [6]: traj = readdy.Trajectory(simulation.output_file)time, counts = traj.read_observable_number_of_particles() Analytical solutionIn ReaDDy, one defines the intrinsic rate constant. In a well-mixed setting, one can use the followingequation (see [1]) to obtain the corresponding macroscopic rate (ODE rate equation/ law of mass action)$$k_\\mathrm{on} = 4\\pi (D_A + D_B) R \\left[1 - \\frac{\\tanh(\\kappa R)}{\\kappa R}\\right]$$where$$\\kappa = \\sqrt{\\frac{\\lambda_\\mathrm{on}}{D_A + D_B}}$$Parameters:intrinsic rate constant $\\lambda_\\mathrm{on} = 1\\,\\mathrm{ns}^{-1}$diffusion coefficient $D_A=D_B=1\\, \\mathrm{nm}^2 \\mathrm{ns}^{-1}$reaction radius $R = 1\\,\\mathrm{nm}$Law of mass action ODE for the concentration of A particles $[A](t) = a(t)$$$\\frac{\\mathrm{d}a}{\\mathrm{d}t} = -k_\\mathrm{on}\\,a^2 \\quad\\text{, with } a(0) = a_0$$which yields$$a(t) = \\frac{1}{a_0^{-1} + k_\\mathrm{on}t}$$[1]: R. Erban and J. Chapman, “Stochastic modelling of reaction-diffusion processes: algorithms for bimolecular reactions.,” Phys. Biol., vol. 6, no. 4, p. 46001, Jan. 2009.In [7]: kappa = np.sqrt(lambda_on / 2.)k_on = 4. * np.pi * 2. * 1. * (1. - np.tanh(kappa * 1.) / (kappa * 1.) )def a(t): return 1. / ((system.box_volume.magnitude / n_particles) + k_on * t)t_range = np.linspace(0., 5000 * 1e-3, 10000) In [8]: plt.plot(time[::200]*1e-3, counts[::200] / system.box_volume.magnitude, "x", label="ReaDDy")plt.plot(t_range, a(t_range), label=r"analytical $a(t)$")plt.legend(loc="best")plt.xlabel("time in nanoseconds")plt.ylabel(r"concentration in nm$^{-3}$")plt.show() Note that in the diffusion-limited regime, the law of mass action solution generally does not reflect what can be observed from the Doi model. ",
"url": "/validation.html#"
},
"validation-boltzmann-distribution-harmonic-repulsion": {
"title": "Validation > Harmonic repulsion - Boltzmann distribution",
"content": "In [1]: import osimport numpy as npimport matplotlib.pyplot as pltimport readdy In [2]: readdy.__version__ Out[2]:'v2.0.3-27'Prepare the systemIn [3]: system = readdy.ReactionDiffusionSystem( box_size=(4, 4, 4), periodic_boundary_conditions=[True, True, True], unit_system=None)system.add_species("A", 1.0)system.add_species("B", 1.0)system.potentials.add_harmonic_repulsion("A", "B", 1., 1.) Prepare the simulationIn [4]: simulation = system.simulation(kernel="SingleCPU")simulation.output_file = "out_rdf.h5"simulation.observe.rdf(1000, np.linspace(0., 2., 10), ["A"], ["B"], 1. / system.box_volume)simulation.add_particle("A", [0., 0., 0.])simulation.add_particle("B", [0., 0., 1.]) Run the simulationIn [5]: if os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.progress_output_stride = 1000simulation.run(n_steps=100000000, timestep=1e-4, show_summary=False) 100%|██████████| 100000/100000 [02:16<00:00, 734.66it/s]In [6]: traj = readdy.Trajectory(simulation.output_file)rdf_times, bin_centers, rdf_values = traj.read_observable_rdf() In [7]: def average_across_first_axis(values): n_values = len(values) mean = np.sum(values, axis=0) / n_values # shape = n_bins difference = values - mean # broadcasting starts with last axis std_dev = np.sqrt(np.sum(difference * difference, axis=0) / n_values) std_err = np.sqrt(np.sum(difference * difference, axis=0) / n_values ** 2) return mean, std_dev, std_errdef plot_boltzmann(force_const, interaction_radius): def potential(r, force_const, interaction_radius): if r < interaction_radius: return 0.5 * force_const * np.power(r - interaction_radius, 2.) else: return 0. boltz = lambda r: np.exp(-1. * potential(r, force_const, interaction_radius)) r_range = np.linspace(0.1, 2., 100) b_range = np.fromiter(map(boltz, r_range), dtype=float) plt.plot(r_range, b_range, label=r"Boltzmann correlation $e^{-\\beta U(r)}$")mean, std_dev, std_err = average_across_first_axis(rdf_values)plt.errorbar(bin_centers, mean, yerr=std_err, fmt=".", label="ReaDDy")plot_boltzmann(1., 1.)plt.legend()plt.xlabel(r"Distance $r$ of A and B")plt.ylabel(r"Radial distribution function $g(r)$")plt.show() ",
"url": "/validation.html#"
},
"validation-lennard-jones-fluid": {
"title": "Validation > Lennard-Jones fluid - thermodynamic state variables",
"content": "We consider a simple fluid of particles that interact due to the Lennard Jones potential$$U(r) = 4\\varepsilon \\left[ \\left(\\frac{\\sigma}{r}\\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^6 \\right]$$To validate the integration of the diffusion process and correct implementation of energies and forces, we compare observables to results from the literature. The observables are: the mean potential energy of the system per particle $U$, the pressure $P$ and the radial distribution function $g(r)$.The thermodynamic quantities of a Lennard Jones fluid are typically expressed in terms of the rescaled density$\\rho^* = (N/V)\\sigma^3$ and the rescaled temperature $T^* = k_BT/\\varepsilon$, where $N$ is the number of particles in the system, which is constant over time, and $V$ is the available volume.In simulation practice we set $\\sigma=1$ and $\\varepsilon=1$ to achieve the reduced units. The quantity $\\sigma^2 / 4D$ gives rise to a typical time scale, where $D$ is the self diffusion coefficient of the particles. In practice we set the diffusion coefficient to 1.We use an Euler-Maruyama scheme to integrate the positions of particles according to the overdamped Langevin equation of motion, in contrast to an integration of positions and momenta in the underdamped limit.The pressure can be measured from the acting forces according to [4]$$PV = Nk_BT + \\langle \\mathscr{W} \\rangle$$where$$\\mathscr{W} = \\frac{1}{3} \\sum_i \\sum_{j>i} \\mathbf{r}_{ij} \\mathbf{f}_{ij},$$is the virial that describes the deviation from ideal-gas-behavior. It is measured in terms of pairwise acting forces $\\mathbf{f}_{ij}$ between particles $i$ and $j$ which are separated by $\\mathbf{r}_{ij}$. This is implemented by ReaDDy's pressure observable.Resultsorigincutoff radius $r_c$density $\\rho$temperature $T$pressure $P$potential energy per particle $U$[1]40.331.023(2)-1.673(2)[2]40.331.0245-1.6717HALMD [3]40.331.0234(3)-1.6731(4)ReaDDy40.331.035(7)-1.679(3)[1]40.633.69(1)-3.212(3)[2]40.633.7165-3.2065HALMD [3]40.633.6976(8)-3.2121(2)ReaDDy40.633.70(2)-3.208(7)[1] Molecular dynamics simulations, J. K. Johnson, J. A. Zollweg, and K. E. Gubbins,The Lennard-Jones equation of state revisited, Mol. Phys. 78, 591 (1993)[2] Integral equations theory, A. Ayadim, M. Oettel, and S Amokrane,Optimum free energy in the reference functional approach for the integral equations theory, J. Phys.: Condens. Matter 21, 115103 (2009).[3] HAL's MD package, http://halmd.org/validation.html[4] Allen, M. P., & Tildesley, D. J. (1987). Computer Simulation of Liquids. New York: Oxford University Press.In [1]: import osimport numpy as npimport readdyprint(readdy.__version__) v2.0.3-27Utility methodsIn [2]: def average_across_first_axis(values): n_values = len(values) mean = np.sum(values, axis=0) / n_values # shape = n_bins difference = values - mean # broadcasting starts with last axis std_dev = np.sqrt(np.sum(difference * difference, axis=0) / n_values) std_err = np.sqrt(np.sum(difference * difference, axis=0) / n_values ** 2) return mean, std_dev, std_errdef lj_system(edge_length, temperature=1.): system = readdy.ReactionDiffusionSystem( box_size=[edge_length, edge_length, edge_length], unit_system=None ) system.kbt = temperature system.add_species("A", diffusion_constant=1.) system.potentials.add_lennard_jones("A", "A", m=12, n=6, epsilon=1., sigma=1., cutoff=4., shift=True) return system Wrap the whole simulation and analysis in a function and perform it for the two densities 0.3 and 0.6.In [3]: def equilibrate_and_measure(density=0.3): n_per_axis = 12 n_particles = n_per_axis ** 3 edge_length = (n_particles / density) ** (1. / 3.) pos_x = np.linspace(-edge_length/2., edge_length/2.-1., n_per_axis) pos = [] for x in pos_x: for y in pos_x: for z in pos_x: pos.append([x,y,z]) pos = np.array(pos) print("n_particles", len(pos), "box edge_length", edge_length) assert len(pos)==n_particles def pos_callback(x): nonlocal pos n = len(x) pos = np.zeros((n,3)) for i in range(n): pos[i][0] = x[i][0] pos[i][1] = x[i][1] pos[i][2] = x[i][2] print("saved positions") # create system system = lj_system(edge_length, temperature=3.) # equilibrate sim = system.simulation(kernel="CPU") sim.add_particles("A", pos) sim.observe.particle_positions(2000, callback=pos_callback, save=None) sim.observe.energy(500, callback=lambda x: print(x), save=None) sim.record_trajectory(stride=1) sim.output_file = "lj_eq.h5" if os.path.exists(sim.output_file): os.remove(sim.output_file) sim.run(n_steps=10000, timestep=1e-4) traj = readdy.Trajectory(sim.output_file) traj.convert_to_xyz(particle_radii={"A": 0.5}) # measure sim = system.simulation(kernel="CPU") sim.add_particles("A", pos) sim.observe.energy(200) sim.observe.pressure(200) sim.observe.rdf( 200, bin_borders=np.linspace(0.5, 4., 50), types_count_from="A", types_count_to="A", particle_to_density=density) sim.output_file = "lj_measure.h5" if os.path.exists(sim.output_file): os.remove(sim.output_file) sim.run(n_steps=10000, timestep=1e-4) # obtain results traj = readdy.Trajectory(sim.output_file) _, energy = traj.read_observable_energy() _, bin_centers, rdf = traj.read_observable_rdf() _, pressure = traj.read_observable_pressure() energy_mean, _, energy_err = average_across_first_axis(energy) # time average energy_mean /= n_particles energy_err /= n_particles pressure_mean, _, pressure_err = average_across_first_axis(pressure) # time average rdf_mean, _, rdf_err = average_across_first_axis(rdf) # time average return { "energy_mean": energy_mean, "energy_err": energy_err, "pressure_mean": pressure_mean, "pressure_err": pressure_err, "rdf_mean": rdf_mean, "rdf_err": rdf_err, "rdf_bin_centers": bin_centers } In [4]: result_low_dens = equilibrate_and_measure(density=0.3) 0%| | 0/1000 [00:00<?, ?it/s] n_particles 1728 box edge_length 17.925618986228656Configured kernel context with:-------------------------------- - kBT = 3.0 - periodic b.c. = (true, true, true) - box size = (17.9256, 17.9256, 17.9256) - particle types: * particle type "A" with D=1.0 - potentials of order 2: * for types "A" and "A" * 12.0-6.0-Lennard-Jones potential with cutoff=4.0, epsilon=1.0, k=4.0, and with energy shiftConfigured simulation loop with:-------------------------------- - timeStep = 0.000100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? truesaved positions-2006.092638276083 5%|▌ | 50/1000 [00:16<05:38, 2.81it/s] -2663.1424920832724 10%|█ | 100/1000 [00:32<05:05, 2.95it/s] -2759.8158335119415 15%|█▌ | 150/1000 [00:48<04:22, 3.24it/s] -2769.480985654246 20%|██ | 200/1000 [01:02<03:31, 3.78it/s] saved positions-2758.3520150659533 25%|██▌ | 250/1000 [01:15<03:08, 3.97it/s] -2838.185778357589 30%|███ | 300/1000 [01:29<03:00, 3.89it/s] -2841.100358293537 35%|███▌ | 350/1000 [01:43<02:51, 3.79it/s] -2801.9811061164173 40%|████ | 400/1000 [01:57<02:44, 3.66it/s] saved positions-2865.1272257328483 45%|████▌ | 450/1000 [02:11<02:37, 3.49it/s] -2891.940879474006 50%|█████ | 500/1000 [02:26<02:25, 3.43it/s] -2824.1108673381364 55%|█████▌ | 550/1000 [02:41<02:13, 3.38it/s] -2930.0205689454847 60%|██████ | 600/1000 [02:55<01:55, 3.46it/s] saved positions-2901.9845693204093 65%|██████▌ | 650/1000 [03:10<01:46, 3.28it/s] -2861.9925088341834 70%|███████ | 700/1000 [03:25<01:33, 3.21it/s] -2829.326983694836 75%|███████▌ | 750/1000 [03:41<01:16, 3.26it/s] -2874.095379626244 80%|████████ | 800/1000 [03:57<01:04, 3.09it/s] saved positions-2918.8710282053125 85%|████████▌ | 850/1000 [04:13<00:47, 3.14it/s] -2809.1883820136936 90%|█████████ | 900/1000 [04:28<00:30, 3.25it/s] -2862.796625209495 95%|█████████▌| 950/1000 [04:45<00:15, 3.15it/s] -2996.5870985783413 100%|██████████| 1000/1000 [05:01<00:00, 3.32it/s] saved positions-2969.2603788072392 0%| | 0/1000 [00:00<?, ?it/s] Configured kernel context with:-------------------------------- - kBT = 3.0 - periodic b.c. = (true, true, true) - box size = (17.9256, 17.9256, 17.9256) - particle types: * particle type "A" with D=1.0 - potentials of order 2: * for types "A" and "A" * 12.0-6.0-Lennard-Jones potential with cutoff=4.0, epsilon=1.0, k=4.0, and with energy shiftConfigured simulation loop with:-------------------------------- - timeStep = 0.000100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 1000/1000 [05:45<00:00, 2.89it/s]In [5]: result_low_dens Out[5]:{'energy_mean': -1.6663065999844973, 'energy_err': 0.004068511442129084, 'pressure_mean': 1.0241356362812084, 'pressure_err': 0.007471245837546726, 'rdf_mean': array([0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 8.73771465e-04, 1.42643258e-01, 7.90076755e-01, 1.37058562e+00, 1.48700339e+00, 1.40774728e+00, 1.26934966e+00, 1.17417650e+00, 1.09760049e+00, 1.03248774e+00, 1.00018151e+00, 9.77339310e-01, 9.52078308e-01, 9.54373365e-01, 9.58670740e-01, 9.78288430e-01, 9.85144126e-01, 9.96332472e-01, 9.96971138e-01, 1.01546690e+00, 1.00858485e+00, 1.01263937e+00, 1.00667808e+00, 1.00287432e+00, 1.00332166e+00, 1.00407493e+00, 1.00003326e+00, 9.97703466e-01, 1.00256106e+00, 1.00126730e+00, 1.00149326e+00, 9.99601950e-01, 9.99207132e-01, 1.00320129e+00, 1.00385635e+00, 1.00234372e+00, 1.00496500e+00, 9.96759322e-01, 1.00211504e+00, 1.00032322e+00, 9.97838840e-01, 1.00435081e+00, 1.00073824e+00, 9.97157011e-01, 1.00660962e+00]), 'rdf_err': array([0. , 0. , 0. , 0. , 0.00035392, 0.00397702, 0.00752209, 0.00932349, 0.00951939, 0.00999279, 0.00882471, 0.00700311, 0.00657965, 0.00635551, 0.00550131, 0.00465666, 0.00450317, 0.00426148, 0.00446077, 0.00411222, 0.00375464, 0.00418982, 0.00456293, 0.00339979, 0.00380776, 0.00383472, 0.00351115, 0.00340765, 0.0029781 , 0.00335204, 0.00309728, 0.00392425, 0.003419 , 0.00321209, 0.0030782 , 0.00278284, 0.00280004, 0.00268631, 0.00244541, 0.0022105 , 0.00284042, 0.00229337, 0.00223074, 0.00268411, 0.00236706, 0.00272272, 0.00256921, 0.00230784, 0.00235968]), 'rdf_bin_centers': array([0.53571429, 0.60714286, 0.67857143, 0.75 , 0.82142857, 0.89285714, 0.96428571, 1.03571429, 1.10714286, 1.17857143, 1.25 , 1.32142857, 1.39285714, 1.46428571, 1.53571429, 1.60714286, 1.67857143, 1.75 , 1.82142857, 1.89285714, 1.96428571, 2.03571429, 2.10714286, 2.17857143, 2.25 , 2.32142857, 2.39285714, 2.46428571, 2.53571429, 2.60714286, 2.67857143, 2.75 , 2.82142857, 2.89285714, 2.96428571, 3.03571429, 3.10714286, 3.17857143, 3.25 , 3.32142857, 3.39285714, 3.46428571, 3.53571429, 3.60714286, 3.67857143, 3.75 , 3.82142857, 3.89285714, 3.96428571])}In [6]: result_hi_dens = equilibrate_and_measure(density=0.6) 0%| | 0/1000 [00:00<?, ?it/s] n_particles 1728 box edge_length 14.227573217960249Configured kernel context with:-------------------------------- - kBT = 3.0 - periodic b.c. = (true, true, true) - box size = (14.2276, 14.2276, 14.2276) - particle types: * particle type "A" with D=1.0 - potentials of order 2: * for types "A" and "A" * 12.0-6.0-Lennard-Jones potential with cutoff=4.0, epsilon=1.0, k=4.0, and with energy shiftConfigured simulation loop with:-------------------------------- - timeStep = 0.000100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? truesaved positions-6881.589189751791 5%|▌ | 50/1000 [01:03<19:58, 1.26s/it] -5551.336656562339 10%|█ | 100/1000 [02:07<18:26, 1.23s/it] -5517.792907415611 15%|█▌ | 150/1000 [03:10<18:16, 1.29s/it] -5578.622215635525 20%|██ | 200/1000 [04:15<17:41, 1.33s/it] saved positions-5565.460833824768 25%|██▌ | 250/1000 [05:19<16:23, 1.31s/it] -5530.96678538631 30%|███ | 300/1000 [06:22<14:44, 1.26s/it] -5748.91818758591 35%|███▌ | 350/1000 [07:27<14:08, 1.31s/it] -5522.570473647219 40%|████ | 400/1000 [08:31<12:17, 1.23s/it] saved positions-5427.976867616478 45%|████▌ | 450/1000 [09:40<13:40, 1.49s/it] -5404.535270232334 50%|█████ | 500/1000 [10:53<11:03, 1.33s/it] -5425.220376798604 55%|█████▌ | 550/1000 [11:57<10:02, 1.34s/it] -5598.445519678548 60%|██████ | 600/1000 [13:03<08:40, 1.30s/it] saved positions-5534.771794722609 65%|██████▌ | 650/1000 [14:08<07:34, 1.30s/it] -5598.408378384099 70%|███████ | 700/1000 [15:12<06:14, 1.25s/it] -5391.389721014838 75%|███████▌ | 750/1000 [16:18<05:27, 1.31s/it] -5525.528134909822 80%|████████ | 800/1000 [17:25<04:24, 1.32s/it] saved positions-5454.222681817357 85%|████████▌ | 850/1000 [18:30<03:13, 1.29s/it] -5578.676448117959 90%|█████████ | 900/1000 [19:37<02:11, 1.31s/it] -5409.8892046793335 95%|█████████▌| 950/1000 [20:43<01:06, 1.33s/it] -5757.165020514167 100%|██████████| 1000/1000 [21:49<00:00, 1.31s/it] saved positions-5447.436977709065 0%| | 0/1000 [00:00<?, ?it/s] Configured kernel context with:-------------------------------- - kBT = 3.0 - periodic b.c. = (true, true, true) - box size = (14.2276, 14.2276, 14.2276) - particle types: * particle type "A" with D=1.0 - potentials of order 2: * for types "A" and "A" * 12.0-6.0-Lennard-Jones potential with cutoff=4.0, epsilon=1.0, k=4.0, and with energy shiftConfigured simulation loop with:-------------------------------- - timeStep = 0.000100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 1000/1000 [26:54<00:00, 1.61s/it]In [7]: result_hi_dens Out[7]:{'energy_mean': -3.193799119175668, 'energy_err': 0.005697042916386074, 'pressure_mean': 3.736313759014431, 'pressure_err': 0.0193156366839875, 'rdf_mean': array([0. , 0. , 0. , 0. , 0.00318302, 0.21544415, 1.09583401, 1.71095452, 1.71660056, 1.5030973 , 1.28105097, 1.11817947, 1.00079171, 0.9281504 , 0.88249208, 0.87185074, 0.87061258, 0.89063383, 0.92622229, 0.97184412, 1.01878889, 1.05056475, 1.06129247, 1.04824337, 1.03340413, 1.019903 , 1.00605256, 0.99656021, 0.98725976, 0.98180169, 0.98345413, 0.98252596, 0.99277386, 0.99484756, 0.99802621, 1.00303123, 1.00645232, 1.00655447, 1.00674462, 1.00682027, 1.00066395, 0.99939614, 0.99703547, 0.99774215, 0.9951048 , 0.99884033, 0.99708233, 0.99665095, 0.99983404]), 'rdf_err': array([0. , 0. , 0. , 0. , 0.000414 , 0.00362726, 0.00552918, 0.0078111 , 0.00731 , 0.00614851, 0.00501308, 0.0051015 , 0.00379222, 0.00358177, 0.00398928, 0.00347434, 0.00301331, 0.00324984, 0.00282947, 0.00364745, 0.00313543, 0.00295688, 0.00287935, 0.00295896, 0.00219868, 0.00230414, 0.00289815, 0.00289699, 0.00200497, 0.00210439, 0.00214856, 0.00240868, 0.00188535, 0.00211364, 0.00206456, 0.00200478, 0.00194404, 0.0020108 , 0.00211606, 0.0020575 , 0.00168416, 0.001548 , 0.00144649, 0.00169588, 0.00170068, 0.00159604, 0.00160979, 0.00145555, 0.0016807 ]), 'rdf_bin_centers': array([0.53571429, 0.60714286, 0.67857143, 0.75 , 0.82142857, 0.89285714, 0.96428571, 1.03571429, 1.10714286, 1.17857143, 1.25 , 1.32142857, 1.39285714, 1.46428571, 1.53571429, 1.60714286, 1.67857143, 1.75 , 1.82142857, 1.89285714, 1.96428571, 2.03571429, 2.10714286, 2.17857143, 2.25 , 2.32142857, 2.39285714, 2.46428571, 2.53571429, 2.60714286, 2.67857143, 2.75 , 2.82142857, 2.89285714, 2.96428571, 3.03571429, 3.10714286, 3.17857143, 3.25 , 3.32142857, 3.39285714, 3.46428571, 3.53571429, 3.60714286, 3.67857143, 3.75 , 3.82142857, 3.89285714, 3.96428571])}In [9]: %matplotlib inlineimport matplotlib.pyplot as pltprint("density 0.3:")print("mean energy per particle {}\\nerr energy per particle {}".format( result_low_dens["energy_mean"], result_low_dens["energy_err"]))print("pressure {}\\nerr pressure {}".format( result_low_dens["pressure_mean"], result_low_dens["pressure_err"]))print("density 0.6:")print("mean energy per particle {}\\nerr energy per particle {}".format( result_hi_dens["energy_mean"], result_hi_dens["energy_err"]))print("pressure {}\\nerr pressure {}".format( result_hi_dens["pressure_mean"], result_hi_dens["pressure_err"]))plt.plot(result_low_dens["rdf_bin_centers"], result_low_dens["rdf_mean"], label=r"density $\\rho=0.3$")plt.plot(result_hi_dens["rdf_bin_centers"], result_hi_dens["rdf_mean"], label=r"density $\\rho=0.6$")plt.xlabel(r"distance $r/\\sigma$")plt.ylabel(r"radial distribution $g(r)$")plt.legend(loc="best")plt.show() density 0.3:mean energy per particle -1.6663065999844973err energy per particle 0.004068511442129084pressure 1.0241356362812084err pressure 0.007471245837546726density 0.6:mean energy per particle -3.193799119175668err energy per particle 0.005697042916386074pressure 3.736313759014431err pressure 0.0193156366839875 ",
"url": "/validation.html#"
},
"validation-mean-squared-displacement": {
"title": "Validation > Diffusion - Mean squared displacement",
"content": "In [1]: import osimport numpy as npimport matplotlib.pyplot as pltimport readdyprint(readdy.__version__) v2.0.3-27Set up a cubic periodic box and register one species A with diffusion constant $D=1$. The simulation adds 800 particles in the origin. We observe the particle positions at every timestep and run the simulation for a given amount of steps.In [2]: system = readdy.ReactionDiffusionSystem( box_size=[10, 10, 10], periodic_boundary_conditions=[True, True, True], unit_system=None)system.add_species("A", 1.0)simulation = system.simulation(kernel="SingleCPU")simulation.output_file = "out_msd.h5"simulation.observe.particle_positions(stride=1)init_pos = np.zeros((800, 3))simulation.add_particles("A", init_pos)if os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.run(n_steps=5000, timestep=1e-3) 6%|▋ | 32/500 [00:00<00:01, 316.09it/s] Configured kernel context with:-------------------------------- - kBT = 1.0 - periodic b.c. = (true, true, true) - box size = (10.0, 10.0, 10.0) - particle types: * particle type "A" with D=1.0Configured simulation loop with:-------------------------------- - timeStep = 0.00100000 - evaluateObservables = true - progressOutputStride = 100 - context written to file = true - Performing actions: * Initialize neighbor list? true * Update neighbor list? true * Clear neighbor list? true * Integrate diffusion? true * Calculate forces? true * Handle reactions? true * Handle topology reactions? true 100%|██████████| 500/500 [00:01<00:00, 319.21it/s]Load the trajectory containing the observables. Since there were no reactions in the system we can assume that the particle positions is of fixed shape (T, N, 3), where N does not change over time. We may want to cast the positions into a numpy array.In [3]: traj = readdy.Trajectory(simulation.output_file)times, positions = traj.read_observable_particle_positions()times = np.array(times) * 1e-3# convert to pure numpy array to make use of fancy operationsT = len(positions)N = len(positions[0])pos = np.zeros(shape=(T, N, 3))for t in range(T): for n in range(N): pos[t, n, 0] = positions[t][n][0] pos[t, n, 1] = positions[t][n][1] pos[t, n, 2] = positions[t][n][2] Wrap the trajectories back to account for periodic boundariesfind the box index of each point in time for each particle and each coordinatewrap the trajectory back to absolute positions for each particle using the box indicesIn [4]: # step 1.box_size = 10.box_indices = np.zeros(shape=(T, N, 3), dtype=int)for t in range(1, T): for n in range(N): for coord in [0, 1, 2]: delta = pos[t, n, coord] - pos[t-1, n, coord] if delta > 0.5 * box_size: box_indices[t, n, coord] = box_indices[t-1, n, coord] - 1 elif delta < - 0.5 * box_size: box_indices[t, n, coord] = box_indices[t-1, n, coord] + 1 else: box_indices[t, n, coord] = box_indices[t-1, n, coord]# step 2.absolute_pos = np.zeros_like(pos)for t in range(T): for n in range(N): absolute_pos[t, n] = pos[t, n] + box_indices[t, n].astype(float) * box_size Using the absolute positions we can calculate the displacement from the initial positions and obtain the quadratic variation (the mean squared displacement) and its standard deviation. Since the particles do not interact their trajectories are uncorrelated and we obtain the standard error of the mean from the same experiment.In [5]: difference = absolute_pos - init_pos# sum over coordinates, per particle per timestepsquared_displacements = np.sum(difference * difference, axis=2)# squared_displacements has shape (T,N)mean = np.mean(squared_displacements, axis=1)std_dev = np.std(squared_displacements, axis=1)std_err = np.std(squared_displacements, axis=1) / np.sqrt(squared_displacements.shape[1])# plotstride = 100plt.errorbar(times[::stride], mean[::stride], yerr=std_err[::stride], fmt=".", label="Particle diffusion")plt.plot(times[::stride], 6. * times[::stride], label=r"$6 D t$")plt.legend(loc="best")plt.xlabel(r"Time $t$")plt.ylabel(r"Mean squared displacement $\\langle (x(t) - x_0)^2 \\rangle_N$")plt.show() ",
"url": "/validation.html#"
},
"workshop-sessions-cookbook-html": {
"title": "Practical sessions > Cookbook",
"content": "This section has some solutions to common problems you might run into. The order in which you create and manipulate the system and simulation matters!Remember that the workflow should always besystem = readdy.ReactionDiffusionSystem(...)# ... system configurationsimulation = system.simulation(...)# ... simulation setupsimulation.run(...)# ... analysis of resultsIf you made a mistake while registering species, potentials, etc., just create a new system and run the configuration again. The same goes for the simulation, just create a new one. In the jupyter notebooks it is sufficient, to just run all the cells again. Simulation box and box potential for confining particlesIf you want to confine particles to some cube without periodic boundaries, the box potential must be fully inside the simulation box. Remember that the box is centered around (0,0,0). For example the following will confine A particles to a cube of edge length 10system = readdy.ReactionDiffusionSystem( box_size=[12., 12., 12.], unit_system=None, periodic_boundary_conditions=[False, False, False])system.add_species(\"A\", 1.)origin = np.array([-5., -5., -5.])extent = np.array([10., 10., 10.])system.potentials.add_box( \"A\", force_constant=10., origin=origin, extent=extent)The two vectors origin and extent span a cube in 3D space, the picture you should have in mind is the following Initial placement of particles inside a certain volumeIf you have already defined origin and extent of your confining box potential, it is easy to generate random positions uniformly in this volumen_particles = 30uniform = np.random.uniform(size=(n_particles,3))init_pos = uniform * extent + originHere uniform is a matrix Nx3 where each row is a vector in the unit cube $\\in{[0,1]\\times[0,1]\\times[0,1]}$ . This is multiplied with extent, yielding a uniform cube ${[0,\\mathrm{extent}_0]\\times[0,\\mathrm{extent}_1]\\times[0,\\mathrm{extent}_2]}$. If you add the origin to this you get this cube at the right position (with respect to our box coordinates centered around (0,0,0)), i.e.$$\\begin{aligned}\\{&[\\mathrm{origin}_0,\\mathrm{extent}_0+\\mathrm{origin}_0]\\\\\\times&[\\mathrm{origin}_1,\\mathrm{extent}_1+\\mathrm{origin}_1]\\\\\\times&[\\mathrm{origin}_2,\\mathrm{extent}_2+\\mathrm{origin}_2]\\}\\end{aligned}$$ 2D planeIf you want to confine particles to a 2D plane, just use the box potential but make the extent in one dimension very small, i.e.system = readdy.ReactionDiffusionSystem( box_size=[12., 12., 3.], unit_system=None, periodic_boundary_conditions=[False, False, False])system.add_species(\"A\", 1.)origin = np.array([-5., -5., -0.01])extent = np.array([10., 10., 0.02])system.potentials.add_box( \"A\", force_constant=10., origin=origin, extent=extent)Having defined origin and extent it is now easy to add particles to this 2D plane. Note that I also made the box_size in z direction smaller, however it should be large enough. Output file sizePlease make use of your /storage/mi/user directories, your home will fill up quicker.data_dir = \"/storage/mi/user\" # replace user with your usernamesimulation.output_file = os.path.join(data_dir, \"myfile.h5\")Additionally, use a stride on your observables, e.g.simulation.record_trajectory(stride=100)# orsimulation.observe.particle_positions(stride=100) Reaction descriptor languageIn expressions likesystem.reactions.add(\"fus: A +(2) B-> C\", rate=0.1)the value in the parentheses +(2) is the reaction distance. look at VMD output if something is fishyThis might reveal some obvious mistakes. Therefore you must have registered the according observablesimulation.record_trajectory(stride=100) # use appropriate stride# ... run simulationtraj = readdy.Trajectory(simulation.output_file)traj.convert_to_xyz(particle_radii={\"A\": 0.5, \"B\": 1.})# particle_radii is optional hereThen in a bash shell dovmd -e myfile.h5.xyz.tclor prefix with ! in the jupyter notebook.",
"url": "/workshop_sessions.html#cookbook"
},
"workshop-sessions-friday-html": {
"title": "Practical sessions > Friday",
"content": "The idea for the friday session is to either Continue to create spatial patterns with the diffusion-limited LV model Try to simulate a system that you are working on / you know well Think of any other scientifically related/unrelated systemAt ca. 15:30h everyone will present their system: a short vmd movie and an explanation of what your system does, how it works.",
"url": "/workshop_sessions.html#friday"
},
"workshop-sessions-monday-html": {
"title": "Practical sessions > Monday - ReaDDy intro",
"content": "Task 0) installationGet miniconda and install itwget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.shbash Miniconda3-latest-Linux-x86_64.shCheck your quota -s. If your home directory is limited in space, you might want to install it under storage, i.e. when prompted for install location, choose either /home/mi/user/miniconda3 or /storage/mi/user/miniconda3 (replace user with your username). Your ~/.bashrc should contain the line. /storage/mi/user/miniconda3/etc/profile.d/conda.shto enable the conda command. If you try out which conda, you should see its function definition, and you’re good to go.Create a conda environment workshopconda create -n workshopand activate it.conda activate workshopAdd conda-forge as default channelconda config --add channels conda-forge --envand install the latest readdy in it(workshop) $ conda install -c readdy/label/dev readdyCheck if it worked, start a python interpreter and doimport readdyprint(readdy.__version__)If this does not return an error, you are ready to readdy (almost).Make sure you also have vmd installed. It should be installed on the universities machines.Task 1) Python, ipython notebook, numpy, matplotlib crash courseFollow along the ipython-intro. You should be able to reproduce the presented usage in your own ipython notebook.The notebook covers the following: usage of an ipython notebook filesystem operations with os save and load data objects with pickle numerical operations with numpy plotting with matplotlib.pyplotTask 2) ReaDDy intro I: Particles, diffusion and potentialsFollow along the readdy-intro. You should be able to reproduce the presented usage in your own ipython notebook.The notebook covers the following: principle workflow of readdy: system and simulation adding particle species to system adding potentials to system spatial layout of simulation box and box potentials adding particle instances at given positions to simulation convert trajectory output to VMD viewable formatTask 3) ReaDDy intro II: Reactions, observables and checkpointsFollow along the readdy-observablesThe notebook covers the following: adding reaction to system adding observable to simulation reading back the observable from trajectory file using checkpoints to continue a simulationTask 4) Crowded mixture, MSDThe time dependent mean squared displacement (MSD) is defined as$$\\langle(\\mathbf{x}_t-\\mathbf{x}_0)^2\\rangle_N$$where $\\mathbf{x}_t$ is the position of a particle at time $t$, and $\\mathbf{x}_0$ is its initial position. The difference of these is then squared, yielding a squared travelled distance. This quantity is then averaged over all particles.The task is to set up a crowded mixture of particles and measure the MSD. There is one species A with diffusion coefficient 0.1. The simulation box shall be $20\\times20\\times20$ and non-periodic, i.e. define a box potential that keeps particles inside, e.g. with an extent $19\\times19\\times19$ and appropriate origin.Define a harmonic repulsion between A particles with force constant 100 and interaction distance 2.Add 1000 A particles to the simulation, uniformly distributed within the box potential.Observe the positions of particles with a stride of 1.Run the simulation with a time step size of 0.01 for 20000 steps.In the post-processing you have to calculate the MSD from the observed positions. Implement the following steps: Convert positions list into numpy array $T\\times N\\times3$ From every particles position subtract the initial position Square this Average over particles np.mean(...)Since the positions observable returns a list of list, it might come in handy to convert this to a numpy array of shape $T\\times N\\times3$ (step 1). You may use the following brute force method to do thisT = len(positions)N = len(positions[0])pos = np.zeros(shape=(T, N, 3))for t in range(T): for n in range(N): pos[t, n, 0] = positions[t][n][0] pos[t, n, 1] = positions[t][n][1] pos[t, n, 2] = positions[t][n][2]You shall implement steps 2.-4. by yourself.4a) Finally plot the MSD as a function of time in a log-log plot, let’s call this the crowded result. (Hint: You may find that there is an initial jump in the squared displacements. Equilibrating the particle positions before starting the measurement may help you there. Make use of checkpointing to use an already equilibrated state)4b) Repeat the whole procedure, but do not register the harmonic repulsion between particles, this shall be the free result. Compare the MSD to the previous result, possibly in the same plot.4c) Additionally plot the analytical solution for freely diffusing particles$$\\langle(\\mathbf{x}_t-\\mathbf{x}_0)^2\\rangle_N = 6 D t$$From your resulting plot identify “finite size saturation”, “subdiffusion”, “reduced normal diffusion”, “ballistic/normal diffusion”Task 5) Crowded mixture, RDFThe radial distribution function (RDF) $g(r)$ describes how likely it is to find two particles at distance $r$. Compare the RDF of harmonically repelling particles and the RDF of non-repelling particles.Therefore set up the same system as in Task 4) but this time the system shall be periodic and there is no box potential.You may want to equilibrate the initial particle positions, use checkpointing.Instead of observing the particle positions, observe the radial distribution functionsimulation.observe.rdf(stride=1000, bin_borders=np.arange(0.,10.,0.2), types_count_from=\"A\", types_count_to=\"A\", particle_to_density=n_particles/system.box_volume)In the post-processing, you shall use_, bin_centers, distribution = traj.read_observable_rdf()to obtain the observable. The distribution contains multiple $g(r)$, one for each time it was recorded. Average them over time.5a) Plot the RDF as a function of the distance (i.e. meandistribution as function of bin_centers)5b) Perform the same procedure but for non-interacting particles, compare with the previous result, preferably in the same plot. What are your observations?",
"url": "/workshop_sessions.html#monday"
},
"workshop-sessions-thursday-html": {
"title": "Practical sessions > Thursday",
"content": "This session will deal with the self assembly of macromolecular structures due to reactions.Task 1) linear filament (e.g. actin) assemblyIn this task we will look at a linear polymer chain, but instead of placing all beads initially, we will let it self-assemble from monomers in solution. The polymer that we will build will have the following structurehead--core--(...)--core--tailwhere (...) means that there can be many core particles, but the structure is always linear.The simulation box is periodic with box_size=[20., 20., 20.].We define three topology particle species and one normal particle speciessystem.add_species(\"substrate\", 0.1)system.add_topology_species(\"head\", 0.1)system.add_topology_species(\"core\", 0.1)system.add_topology_species(\"tail\", 0.1)We also define one topology typesystem.topologies.add_type(\"filament\")There should be the following potentials for topology particlessystem.topologies.configure_harmonic_bond( \"head\", \"core\", force_constant=100, length=1.)system.topologies.configure_harmonic_bond( \"core\", \"core\", force_constant=100, length=1.)system.topologies.configure_harmonic_bond( \"core\", \"tail\", force_constant=100, length=1.)The polymer should be rather stiff, which is the case for Actin filaments in biological cells. We can compactly write this for all triplets of particle types:triplets = [ (\"head\", \"core\", \"core\"), (\"core\", \"core\", \"core\"), (\"core\", \"core\", \"tail\"), (\"head\", \"core\", \"tail\")]for (t1, t2, t3) in triplets: system.topologies.configure_harmonic_angle( t1, t2, t3, force_constant=50., equilibrium_angle=np.pi)We now introduce a topology reaction. They allow changes to the graph of a topology in form of a reaction. Here we will use the following definition.system.topologies.add_spatial_reaction( \"attach: filament(head) + (substrate) -> filament(core--head)\", rate=5.0, radius=1.5)Using the documentation, familiarize yourself, what this means. Do not hesitate to ask, since topology reactions can become quite a tricky concept!Next create a simulation object. We want to observe the followingsimulation.record_trajectory(stride=100)simulation.observe.topologies(stride=100)Add one filament topology to the simulationinit_top_pos = np.array([ [ 1. ,0. ,0.], [ 0. ,0. ,0.], [-1. ,0. ,0.]])top = simulation.add_topology( \"filament\", [\"head\", \"core\", \"tail\"], init_top_pos)top.get_graph().add_edge(0, 1)top.get_graph().add_edge(1, 2)Additionally we need substrate particles, that can attach themselves to the filamentn_substrate = 300origin = np.array([-10., -10., -10.])extent = np.array([20., 20., 20.])init_pos = np.random.uniform(size=(n_substrate,3)) * extent + originsimulation.add_particles(\"substrate\", positions=init_pos)Then, run the simulationif os.path.exists(simulation.output_file): os.remove(simulation.output_file)dt = 5e-3simulation.run(400000, dt)One important observable will be the length of the filament as a function of time. It can be obtained from the trajectory as follows:times, topology_records = traj.read_observable_topologies()chain_length = [ len(tops[0].particles) for tops in topology_records ]The last line is a list comprehension. tops is a list of topologies for a given time step. Hence, tops[0] is the first (and in this case, the only) topology in the system. tops[0].particles is a list of particles belonging to this topology. Thus, its length yields the length of the filament.1a) Have a look at the VMD output. Describe what happens? Additionally plot the length of the filament as a function of time. Note that you shall now plot the simulation time and not the time step indices, i.e. do the followingtimes = np.array(times) * dtwhere dt is the time step size.1b) Using your data of the filament-length, fit a function of the form$$f(t)=a(1-e^{-bt})+3$$to your data. You should use scipy.optimize.curve_fit to do soimport scipy.optimize as sodef func(t, a, b): return a*(1. - np.exp(-b * t)) + 3.popt, pcov = so.curve_fit(func, times, chain_length)print(\"popt\", popt)print(\"pcov\", pcov)f = lambda t: func(t, popt[0], popt[1])plt.plot(times, chain_length, label=\"data\")plt.plot(times, f(times), label=r\"fit $f(t)=a(1-e^{-bt})+3$\")plt.xlabel(\"Time\")plt.ylabel(\"Filament length\")plt.legend()plt.show()Question: Given the result of the fitting parameters How large is the equilibration rate? What will be the length of the filament for $t\\to\\infty$?1c) We now introduce a disassembly reaction for the tail particle. This is done by adding the following to your system configuration.def rate_function(topology): \"\"\" if the topology has at least (head, core, tail) the tail shall be removed with a fixed probability per time \"\"\" vertices = topology.get_graph().get_vertices() if len(vertices) > 3: return 0.05 else: return 0.def reaction_function(topology): \"\"\" find the tail and remove it, and make the adjacent core particle the new tail \"\"\" recipe = readdy.StructuralReactionRecipe(topology) vertices = topology.get_graph().get_vertices() tail_idx = None adjacent_core_idx = None for v in vertices: if topology.particle_type_of_vertex(v) == \"tail\": adjacent_core_idx = v.neighbors()[0].get().particle_index tail_idx = v.particle_index recipe.separate_vertex(tail_idx) recipe.change_particle_type(tail_idx, \"substrate\") recipe.change_particle_type(adjacent_core_idx, \"tail\") return recipesystem.topologies.add_structural_reaction( \"detach\", \"filament\", reaction_function=reaction_function, rate_function=rate_function)Familiarize yourself with this kind of structural topology reactionRepeat the same analysis as before, and also observe your VMD output. How large is the equilibration rate? What will be the length of the filament for $t\\to\\infty$?Task 2) Assembly of virus capsidsThis task will suggest a model for the assembly of monomers into hexamers, and in the bonus task: the assembly of these hexamers into even larger superstructures.First we need a model for one monomer, which should look as followsIt is essentially one bigger core particle with two sites attached. This is sometimes called a patchy particle, i.e. a particle with reaction patches. In the language of ReaDDy, this group of particles is a topology, let’s give it the name CAsystem.topologies.add_type(\"CA\")system.add_topology_species(\"core\", 0.1)system.add_topology_species(\"site\", 0.1)The angle between the triplet site--core--site, shall be fixed to 120°. The sites shall react with the sites of other monomers to form a dimer that looks likeFor this dimer, the four particles shall be confined to a single 2D plane, which we do using topology potentials. Summarizing all potentials that you need to configure for the particles: harmonic bond between core--site with force constant 100 and length 1 harmonic bond between core--core with force constant 100 and length 2 harmonic bond between site--site with force constant 100 and length 0.1 harmonic angle between site--core--site with force constant 200 and equilibrium angle 120 degrees harmonic angle between site--core--core with force constant 200 and equilibrium angle 120 degrees harmonic angle between core--core--core with force constant 200 and equilibrium angle 120 degrees dihedral between core--core--core--core with force constant 200, mutliplicity of 1 and equilibrium angle of 0 dihedral between site--core--core--core with force constant 200, mutliplicity of 1 and equilibrium angle of 0 dihedral between site--core--core--site with force constant 200, mutliplicity of 1 and equilibrium angle of 0 (normal) harmonic repulsion between core and core with force constant 80 and interaction distance 2The formation of the dimer will be done using a spatial topology reactionof the formsystem.topologies.add_spatial_reaction( \"attach: CA(site)+CA(site)->CA(site--site) [self=true]\", rate=10., radius=0.4)After such a reaction the connectivity looks like (...)--core--site--site--core--(...). In a second step after the reaction we want to get rid of the two site particles in the middle. This will be done using a structural topology reaction of the following kind (make sure that you understand what happens in these code snippets. Do ask, if you have problems understanding!). The first ingredient is the rate function for the structural reaction, i.e. given a toplogy this function shall return a very high rate, if there is a site--site connection, and shall return 0 otherwisedef clean_sites_rate_function(topology): edges = topology.get_graph().get_edges() vertices = topology.get_graph().get_vertices() if len(vertices) > 3: for e in edges: v1_ref, v2_ref = e[0], e[1] v1 = v1_ref.get() v2 = v2_ref.get() v1_type = topology.particle_type_of_vertex(v1) v2_type = topology.particle_type_of_vertex(v2) if v1_type == \"site\" and v2_type == \"site\": return 1e12 else: return 0. return 0.The second ingredient is the reaction function, that performs the removing of the two site particles, after the rate function has returned a very high ratedef clean_sites_reaction_function(topology): recipe = readdy.StructuralReactionRecipe(topology) vertices = topology.get_graph().get_vertices() def search_configuration(): # dfs for finding configuration core-site-site-core for v1 in vertices: if topology.particle_type_of_vertex(v1) == \"core\": for v2_ref in v1.neighbors(): v2 = v2_ref.get() if topology.particle_type_of_vertex(v2) == \"site\": for v3_ref in v2.neighbors(): v3 = v3_ref.get() if v3.particle_index != v1.particle_index: if topology.particle_type_of_vertex(v3) == \"site\": for v4_ref in v3.neighbors(): v4 = v4_ref.get() if v4.particle_index != v2.particle_index: if topology.particle_type_of_vertex(v4) == \"core\": return v1.particle_index, v2.particle_index, v3.particle_index, v4.particle_index core1_p_idx, site1_p_idx, site2_p_idx, core2_p_idx = search_configuration() # find corresponding vertex indices from particle indices core1_v_idx = None site1_v_idx = None site2_v_idx = None core2_v_idx = None for i, v in enumerate(vertices): if v.particle_index == core1_p_idx and core1_v_idx is None: core1_v_idx = i elif v.particle_index == site1_p_idx and site1_v_idx is None: site1_v_idx = i elif v.particle_index == site2_p_idx and site2_v_idx is None: site2_v_idx = i elif v.particle_index == core2_p_idx and core2_v_idx is None: core2_v_idx = i else: pass if (core1_v_idx is not None) and (core2_v_idx is not None) and (site1_v_idx is not None) and ( site2_v_idx is not None): recipe.add_edge(core1_v_idx, core2_v_idx) recipe.separate_vertex(site1_v_idx) recipe.separate_vertex(site2_v_idx) recipe.change_particle_type(site1_v_idx, \"dummy\") recipe.change_particle_type(site2_v_idx, \"dummy\") else: raise RuntimeError(\"core-site-site-core wasn't found\") return recipeFinally add the structural reaction to the systemsystem.topologies.add_structural_reaction( \"clean_sites\", topology_type=\"CA\", reaction_function=clean_sites_reaction_function, rate_function=clean_sites_rate_function, raise_if_invalid=True, expect_connected=False)2a) Simulate the system described above in a periodic box of size [25, 25, 25] for n_steps=50000 steps with a timestep of 0.005. Initially place 150 CA patchy particles uniformly distributed in the box. While simulating, observe the trajectory and topologies with the same stride.sim.record_trajectory(n_steps//2000)sim.observe.topologies(n_steps//2000)What do you observe in the VMD output? Do particles assemble in the way you expected?Hints: You should place the particles such that the site particles are already at their prescribed 120 degree angle and a distance of 1 away from the core particle. Adding one such particle can be done in the following way core = np.array([0., 0., 0.])site1 = np.array([0., 0., 1.])site2 = np.array([np.sin(np.pi * 60. / 180.), 0., - 1. * np.cos(np.pi * 60. / 180.)]) top = sim.add_topology(\"CA\", [\"site\", \"core\", \"site\"], np.array([site1, core, site2]))top.get_graph().add_edge(0, 1)top.get_graph().add_edge(1, 2) To distribute the particles uniformly you should add a random translation vector to all positions core, site1 and site2. If you want to be super cool, you can rotate the patchy particle by a random amount before translating it. Ask google how to generate a random rotation matrix and how to apply it to your vectors core, site1 and site22b) From your output file and using the topologies and trajectory observable, calculate the time dependent distribution of molecular mass. This means: Given an instance of a topology, the degree of polymerization is the number of connected core particles in this topology. For one polymer the molecular mass is equal to the degree of polymerization. Obtain such a value for all topologies in a given timestep and make a histogram of that. Now that histogram only counts the occurrence of how many times a topology with a certain molecular mass shows up. To convert that into a distribution of molecular mass itself, you have to multiply the number of occurrence for each degree of polymerization by the degree of polymerization itself. Repeat this for all observed times to obtain as many histograms as there are timesteps.Hints The actual trajectory can be obtained from the trajectory file like so traj_file = readdy.Trajectory(out_file)traj = traj_file.read() The particle type (string) of a particle with index v at time t is traj[t][v].type Construct the histogram for each time using np.histogram(current_sizes, bins=bin_edges), where current_sizes is the list of the molecular masses you have obtained, and bin_edges=np.arange(1,10,1) For plotting it might come in handy to convert the bin_edges to bin_centers, by calculating the midpoints for each bin Plot the histograms using the following snippet xs = np.array(times) * dtys = bin_centersX, Y = np.meshgrid(xs, ys-1)Z = all_histograms.transpose()plt.pcolor(X, Y, Z, cmap=plt.cm.viridis_r)plt.xlabel(\"Time\")plt.ylabel(\"Degree of polymerization\")plt.title(\"Distribution of molecular mass\") 2c) Calculate a similar distribution of molecular mass, but now only for completely assembled topologies, i.e. topologies with no open site particles left. What is the percentage of “misfolded” topologies?2d) From looking at your distribution of pentamers, hexamers, and heptamers. Can you form a full capsid out of the patchy particles we have used? Have a look at the introductory image with the viruses, what do you notice about the capsomers?2e) Bonus task: Introduce a third reactive patch for each patchy particle called offsite, which allows binding to other offsite particles. In this way try to assemble a larger super structure out of the hexamers.",
"url": "/workshop_sessions.html#thursday"
},
"workshop-sessions-tuesday-html": {
"title": "Practical sessions > Tuesday",
"content": "This session will cover a well studied reaction diffusion system, the Lotka Volterra system, that describes a predator prey dynamics. We will investigate how parameters in these system affect spatiotemporal correlations.Task 1) Diffusion influenced reactionConsider the reaction$$A + A\\to B\\quad\\text{with rate: }\\lambda=0.1$$The reaction distance should be $R=1$.Set up a periodic box with dimensions $20\\times20\\times20$Add 1000 A particles with $D=1$. The diffusion constant of B should be $D=1$. They should be uniformly distributed in the box, usen_particles = 1000init_pos = np.random.uniform(size=(n_particles,3)) * extent + originsimulation.add_particles(\"A\", init_pos)Observe the number of particles with a stride of 1. Additionally you can print the current number of particles using the callback functionalitysimulation.observe.number_of_particles( 1, [\"A\",\"B\"], callback=lambda x: print(\"A {}, B {}\".format(x[0],x[1])))Simulate the system for 10000 steps with time step size 0.01.1a) Obtain the number of particles and plot them as a function of time.The analytic solution of the concentration of particles subject to the reaction above is obtained by solving the ODE$$\\frac{\\mathrm{d}a}{\\mathrm{d}t} = -k a^2\\quad a(0)=a_0$$which yields$$a(t) = \\frac{1}{a_0^{-1}+kt}$$1b) Fit the function $a(t)$ to your counts data, to obtain the parameter $k$. For $a_0$ use n_particles. Make use of scipy.optimize.curve_fit. Can the function $a(t)$ describe the data well?1c) Repeat the procedures a) and b) but change the diffusion coefficient to $D=0.01$ and change the microscopic reaction rate to $\\lambda=1$. What happened to your fit?Task 2) well mixed predatory-preySimulate a Lotka-Volterra/predator-prey system for the given parameters and determine the period of oscillation.There are two particle species A (prey) and B (predator) with the same diffusion coefficient $D=0.7$. Both particle species shall be confined to a quadratic 2D plane with an edge length of roughly 50, and non periodic boundaries, i.e. the 2D plane must be fully contained in the simulation box. Choose a very small thickness for the plane, e.g. 0.02, and a force constant of 50.Particles of the same species interact via harmonic repulsion, with a force constant of 50, and an interaction distance of 2.There are three reactions for the Lotka Volterra system: birth, eat, and decay.$$\\mathrm{birth}: \\mathrm{A}\\to \\mathrm{A} +\\mathrm{A}\\quad\\mathrm{with }~ \\lambda=4\\times 10^{-2}, R=2$$$$\\mathrm{eat}: \\mathrm{A}+\\mathrm{B}\\to \\mathrm{B} +\\mathrm{B}\\quad\\mathrm{with }~ \\lambda=10^{-2}, R=2$$$$\\mathrm{decay}: \\mathrm{B}\\to \\emptyset\\quad\\mathrm{with }~ \\lambda=3\\times 10^{-2}$$Initialize a system by randomly positioning 250 particles of each species on the 2D plane. Run the simulation for 100,000 integration steps with a step size of 0.01. Observe the number of particles and the trajectory, with a stride of 100.From the observable, plot the number of particles for the two species as a function of time in the same plot. Make sure to label the lines accordingly. Additionally plot the phase space trajectory, i.e. plot the number of predator particles against the number of prey particles.Question: What is the period of oscillation?Task 3) diffusion influenced predator-preyWe simulate the same system as in Task 2) but now introduce a scaling parameter $\\ell$. This scaling parameter is used to control the ratio of reaction rates $\\lambda$ and the rate of diffusion $D$.For a given parameter $\\ell$ change all reaction rate constants $\\tilde\\lambda=\\lambda\\times\\sqrt{\\ell}$. change all diffusion coefficients $\\tilde D= D/\\sqrt{\\ell}$.where the tilde ~ denotes the actually used value in the simulation and the value without tilde is the one from Task 1).This means that the ratio$$\\frac{\\tilde\\lambda}{\\tilde D} = \\ell\\,\\frac{\\lambda}{D}=\\ell\\times\\mathrm{const}$$is controlled by $\\ell$.Perform the same simulation as in Task 1) but for different scaling parameters $\\ell$, each time saving the resulting plots (time series and phase plot) to a file, use plt.savefig(...). Run each simulation for $n$ number of integration steps, where$$n(\\ell) = 10^4 + \\frac{10^5}{\\sqrt{\\ell}}$$Vary the scaling parameter from 1 to 400.Question: For which $\\ell$ do you see a qualitative change of behavior in the system, both from looking at the plots and the VMD visualization? Which cases are reaction-limited and which cases are diffusion-limited?Task 4)Starting from the parameters of a diffusion influenced system from Task 3), set up a simulation, where prey particles form a traveling wavefront, closely followed by predators. Therefore you might want to change the simulation box and box potential parameters to get a 2D rectangular plane.Feel free to adjust all parameters. You can experiment with other spatial patterns as well, have a look at this paper.",
"url": "/workshop_sessions.html#tuesday"
},
"workshop-sessions-wednesday-html": {
"title": "Practical sessions > Wednesday",
"content": "This session will cover macromolecules and their dynamics. In particular we want to model short RNA chains, represented by linear chains of beads.Assume $N$ beads of particles at positions $\\mathbf{x}_i$, the $i$th bead is connected with the $i+1$th bead by a spring which has a fixed length $l$. Thus the whole chain of particles is linearly connected. The vector $\\mathbf{r}_i=\\mathbf{x}_{i+1}-\\mathbf{x}_i$ is the segment that connnects adjacent beads.One can measure how strongly the $i$th and the $j$th segment correlate by considering the scalar product $\\mathbf{r}_i\\cdot\\mathbf{r}_j$.Since this value alone is quite meaningless, one can consider its average over a whole ensemble of segments. This average can be taken over all segments, i.e. $\\forall i,j$ in the linear chain, and also over many times if the linear chain evolves over time. For realistic polymers one typically observes that the $i$th segment strongly correlates with the $j=i+1$th segment. However for $j\\gg i$ the correlation vanishes. Phenomenologically this is understood by an exponential decay of the correlation$$\\langle \\mathbf{r}_i\\cdot\\mathbf{r}_j \\rangle=l^2\\exp\\left(-|j-i|\\,l/l_p\\right)$$where we have defined the persistence length $l_p$. The value of $l_p$ is determined by the interaction of the beads. For example structured RNA molecules typically show a persistence length of 72nm. Today we will instead model unstructed RNA molecules, which show a persistence length of roughly 2nm.We will consider two models for the polymer, namely the freely jointed chain (FJC), and the freely rotating chain (FRC).In the FJC, the beads are connected by segments of fixed length $l=0.48$. Other than that there is no interaction.In the FRC, the beads are also connected by segments of fixed length $l=0.48$. Additionally the angle between neighbouring segments is fixed to $\\theta=35^\\circ$.Task 1) Equilibrate polymersIn the first task you will a) equilibrate a freely jointed chain (FJC) of $N=50$ beads. Equilibration is ensured by measuring the radius of gyration of a polymer over time. b) Once the polymer is equilibrated, you will measure the persistence length $l_p$. c) Repeat a) and b) for the freely rotating chain (FRC)We will only need one particle species monomer with diffusion constant $D=0.1$. Note that this will not be a normal species but will be a topology species:system.add_topology_species(\"monomer\", 0.1)The simulation box shall have box_size= [102.,102.,102.], and be non-periodic. This means that there has to be a box potential registered for the type monomer with force constant 50, origin=np.array([-50,-50,-50]) and extent=np.array([100,100,100]).In order to build a polymer we use topologies. At first we need to register a type of toplogysystem.topologies.add_type(\"polymer\")The monomers in a polymer must be held together by harmonic bonds, defined by pairs of particle typessystem.topologies.configure_harmonic_bond( \"monomer\", \"monomer\", force_constant=50., length=0.48)Only in the case of a FRC we also specify an angular potential, that is defined for a triplet of particle typessystem.topologies.configure_harmonic_angle( \"monomer\", \"monomer\", \"monomer\", force_constant=20., equilibrium_angle=(180.-35.)/180.*np.pi)where an equilibrium_angle is given in radians (Note that the equilibrium angle here is not the same as the $\\theta$ as defined above, thus the conversion by 180 degrees).The next step is creating the simulation. Then we specify the observables, the trajectory and the particle positionssimulation.record_trajectory(stride=10000)simulation.observe.particle_positions(stride=1000)We also want to make use of checkpointing to continue simulation for an already equilibrated polymer. If there are no checkpoints, we want to create new positions for the polymer. The new positions represent a random walk in three dimensions with fixed step length bond_length.if os.path.exists(checkpoint_dir): # load checkpoint simulation.load_particles_from_latest_checkpoint(checkpoint_dir)else: # new positions init_pos = [np.zeros(3)] for i in range(1, chain_length): displacement = np.random.normal(size=(3)) displacement /= np.sqrt(np.sum(displacement * displacement)) displacement *= bond_length init_pos.append(init_pos[i - 1] + displacement) init_pos = np.array(init_pos) # subtract center of mass init_pos -= np.mean(init_pos, axis=0) # add all particles for the topology at once top = simulation.add_topology(\"polymer\", len(init_pos) * [\"monomer\"], init_pos) # set up the linear connectivity for i in range(1, chain_length): top.get_graph().add_edge(i - 1, i)# this also creates the directorysimulation.make_checkpoints(n_steps // 100, output_directory=checkpoint_dir, max_n_saves=10)Tip: Keep two separate checkpoint directories output files and for the FJC and the FRC model. This means you may want to have the following defined in the beginning of your notebookchain_type = \"fjc\" # fjc or frcout_dir = \"/some/place/on/your/drive\"out_file = os.path.join(out_dir, f\"polymer_{chain_type}.h5\")checkpoint_dir = os.path.join(out_dir, f\"ckpts_{chain_type}\")Now that we have defined the simulation object we can run the simulationif os.path.exists(simulation.output_file): os.remove(simulation.output_file)simulation.run(n_steps, dt)and also observe the outputtraj = readdy.Trajectory(out_file)traj.convert_to_xyz( particle_radii={\"monomer\": bond_length / 2.}, draw_box=True)1a) The radius of gyration is a measure of how ‘extended’ in space a polymer is. To calculate it, we must have observed the particle positions. As a first step we convert the readdy output to a numpy arraytimes, positions = traj.read_observable_particle_positions()# convert to numpy arrayT = len(positions)N = len(positions[0])pos = np.zeros(shape=(T, N, 3))for t in range(T): for n in range(N): pos[t, n, 0] = positions[t][n][0] pos[t, n, 1] = positions[t][n][1] pos[t, n, 2] = positions[t][n][2]Then from the pos array you may use the following to calculate the radius of gyration (the assertion statements may help you understand how the arrays are shaped)# calculate radius of gyrationmean_pos = np.mean(pos, axis=1)difference = np.zeros_like(pos)for i in range(n_particles): difference[:,i] = pos[:,i] - mean_posassert difference.shape == (T,N,3)# square and sum over coordinates (axis=2)squared_radius_g = np.sum(difference * difference, axis=2)assert squared_radius_g.shape == (T,N)# average over particles (axis=1)squared_radius_g = np.mean(squared_radius_g, axis=1)radius_g = np.sqrt(squared_radius_g)assert radius_g.shape == times.shape == (T,)Plot the radius of gyration as a function of time, is it equilibrated? If not, simulate for a longer time.1b) The mean correlation of segments $\\langle \\mathbf{r}_i\\cdot\\mathbf{r}_j \\rangle$ shall be calculated from the pos array. You will average over all pairs of the linear chain and also over all times.Use the following snippet to calculate the segments vectorassert pos.shape == (T, N, 3)# calculate segmentssegments = pos[:, 1:, :] - pos[:, :-1, :]The correlation between $i$ and $j$ shall be measured as a function of their separation $s=\\mid j-i\\mid$, which is a value between 0 and $N-1$, e.g.n_beads = pos.shape[1]separations = np.arange(0, n_beads - 1, 1)corrs = None # your taskNow for every separation $s$, calculate the average correlation, averaged over all pairs $(i,j)$ that lead to $s=\\mid j-i\\mid$.Hints: The calculation may involve a double loop over all segments for i in range(n_beads-1): for j in range(i, n_beads-1): # something The scalar product of the $i$th and the $j$th bead for all times is np.sum(segments[:, i, :] * segments[:, j, :], axis=1) where the summation over axis=1 is over the x,y,z coordinates Then plot the mean correlation as a function of the separation. What do you observe?To determine the persistence length $l_p$, fit a function of the form$$f(s)=c_1e^{-sc_2}$$to the data using scipy.optimize.curve_fit. From the constant $c_2$ you should be able to obtain the persistence length. What is its value?1c) Repeat a) and b) for the FRC. Note the value for $l_p$ for both models.Given the values in the introduction text, which of the models (FJC, FRC) is better suited to model unstructured RNA?Task 2) Identify given structuresObtain the two data files a.npy and b.npy and save them to a directory of your liking. Each of them contains 100 positions of beads, i.e. a and b are two polymer configurations. You can load them as followsa = np.load(\"a.npy\")b = np.load(\"b.npy\")assert a.shape == (100, 3)assert b.shape == (100, 3)Your task is to identify which of them is the FJC model and which is the FRC model, from what you’ve learned in task 1.Task 3) First-passage times of finding targetYou shall now use the configuration x.npy (where x is either a or b)from task 2 that corresponds to the FRC, to set up another simulation, in which one bead of the polymer is of target type. Freely diffusing A particles have to find the target particle. The application you should have in mind is proteins docking to a certain part of nucleid acids, which is crucial for the function of each biological cell. To determine when an A particle has found the target we implement the following kind of reaction$$\\mathrm{A} + \\mathrm{target} \\to \\mathrm{B} + \\mathrm{target}$$The time when the first B particle is created, then corresponds to the first-passage time of that reaction.3a) Perform a simulation that initializes the polymer from x.npy where the 10th bead is of type target, and places 50 A particles normally distributed (with variance $\\sigma=0.5$) around the origin.Therefore use the following system configurationsystem = readdy.ReactionDiffusionSystem( box_size=[16., 16., 16.], periodic_boundary_conditions=[False, False, False], unit_system=None)system.add_topology_species(\"monomer\", 0.1)system.add_topology_species(\"target\", 0.1)system.add_species(\"A\", 0.5)system.add_species(\"B\", 0.5)system.topologies.add_type(\"polymer\")origin = np.array([-7.5, -7.5, -7.5])extent = np.array([15., 15., 15.])system.potentials.add_box(\"monomer\", force_constant=50., origin=origin, extent=extent)system.potentials.add_box(\"target\", force_constant=50., origin=origin, extent=extent)system.potentials.add_box(\"A\", force_constant=50., origin=origin, extent=extent)system.potentials.add_box(\"B\", force_constant=50., origin=origin, extent=extent)with the following topology potentialssystem.topologies.configure_harmonic_bond( \"monomer\", \"monomer\", force_constant=50., length=bond_length)system.topologies.configure_harmonic_bond( \"monomer\", \"target\", force_constant=50., length=bond_length)system.topologies.configure_harmonic_angle( \"monomer\", \"monomer\", \"monomer\", force_constant=20., equilibrium_angle=2.530727415391778)system.topologies.configure_harmonic_angle( \"monomer\", \"monomer\", \"target\", force_constant=20., equilibrium_angle=2.530727415391778)system.topologies.configure_harmonic_angle( \"monomer\", \"target\", \"monomer\", force_constant=20., equilibrium_angle=2.530727415391778)Define a boolean flag interaction = True.If the bool interaction is True then there should be a weak interaction between A and monomer particles with a force_constant of 50, desired_distance=bond_length, a depth of 1.4, and a cutoff of 2.2*bond_length.What does such an interaction result in?Additionally there will be repulsion potentials between monomers and between A particles.system.potentials.add_harmonic_repulsion( \"monomer\", \"monomer\", force_constant=50., interaction_distance=1.1*bond_length)system.potentials.add_harmonic_repulsion( \"A\", \"A\", force_constant=50., interaction_distance=1.5*bond_length)Finally the system needs the reactionsystem.reactions.add(\"found: A +(0.48) target -> B + target\", rate=10000.)where we use a very high rate, such that the reaction will happen directly on contact, where 0.48 is the contact distance.Observe the number of B particles with a stride of 1.Load the polymer configuration and turn the 10th monomer into a target.init_pos = np.load(\"x.npy\")types = len(init_pos) * [\"monomer\"]types[10] = \"target\" # define the target to be the 10-th monomertop = sim.add_topology(\"polymer\", types, init_pos)for i in range(1, len(init_pos)): top.get_graph().add_edge(i - 1, i)Place 50 A particles normally distributed, with variance $\\sigma=0.5$, around the origin.Simulate the system with a timestep of 0.001 for 30000 steps. Have a look at the VMD output. How do the A particles interact with the polymer? Calculate the first passage time from the observed number of particles, i.e. the time when the first B was created.3b) Combine the simulation procedure above into a function of the signaturedef find_target(interaction=False): ... # Since we will run many simulations # you may want to supress the textual output by # setting the two options show_progress # and show_summary sim.show_progress = False sim.run(..., show_summary=False) ... return passage_timeHint: One such simulation should not take much longer than 10 seconds.Gather passage times in a listts_int = []Repeat the simulation many times (50-100 should suffice) and append the result to the list.from tqdm import tqdm_notebook as tqdmn=50for _ in tqdm(range(n)): ts_int.append(find_target(interaction=True))As this might take a while, you will want to observe how long the whole process takes, which is done here using tqdm.Do the same for the case of no interaction, i.e. set interaction=False and gather the results in another list ts_noint.ProTip: To save computation time, run this second case in a copy of your notebook (i.e. at the same time) and save the resulting list ts_noint into a pickle file, which you can read in in the original notebook.For both cases interaction=True and interaction=False, calculate the distribution of first passage times, i.e. plot a histogram of the lists you constructed using plt.hist(). Use bins=np.logspace(0,2,20) and density=True.When assuming a memory less (Poisson) process, the only relevant parameter is the mean rate $\\lambda=N/\\sum_{i=1}^N\\tau_i$. Plot the distribution of first-passage-times with mean rate $\\lambda$, i.e. the Poisson probability density (not the cumulative) of 1 event occurring before time t. Compare against your measured distribution of waiting times?Is the process of finding the target with or without interaction well suited to be modeled as a memory-less process?Now additionally mark the mean first passage time for each case using plt.vlines(). Is the difference of the two cases well described by the mean first passage time?",
"url": "/workshop_sessions.html#wednesday"
},
"blub": {
"title": "",
"content": "",
"url": ""
}
};
</script>
<script src="https://readdy.github.io/libraries/lunr/js/lunr.min.js"></script>
<script src="https://readdy.github.io/assets/js/search.js"></script>
</article>
<div class="foot">
© Copyright 2020 <a href="https://www.mi.fu-berlin.de/en/math/groups/comp-mol-bio/">AI4Science (former CMB) Group</a>
</div>
</div>
</div>
<script src="https://readdy.github.io/libraries/perfect-scrollbar/js/perfect-scrollbar.min.js"></script>
<script src="https://readdy.github.io/assets/js/scrollbar.js"></script>
</body>
</html>