Figure 5.13:

The equilibrium reaction extent's probability density for Reactions \ref {rxn:rxn2} at system volume \Omega =20 (top) and \Omega = 200 (bottom). Notice the decrease in variance in the reaction extent as system volume increases.

Code for Figure 5.13

Text of the GNU GPL.

main.py

 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
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import expm

def simulate_system(Vol):
    x = np.array([1, 5, 0]) * Vol
    nmol = x[0]
    k = np.zeros(2)
    k[0] = 1 / Vol
    k[1] = 3
    extent = np.linspace(1, 0, num=int(nmol+1))

    A = np.zeros((int(nmol+1), int(nmol+1)))
    for i in range(int(nmol+1)):
        n = i
        nu = nmol - n
        a = n
        b = x[1] - nu
        c = nu
        A[i, i] = -(k[0] * a * b + k[1] * c)
        if i > 0:
            A[i, i-1] = k[1] * (c+1)
        if i < nmol:
            A[i, i+1] = k[0] * (a+1) * (b+1)

    p = np.zeros(int(nmol+1))
    p[-1] = nmol

    deltat = 0.01
    tfinal = 1
    time = np.arange(0, tfinal + deltat, deltat)
    Ad = expm(A * deltat)
    for t in time:
        p = Ad @ p
        p = p / p.sum() * nmol

    peq = p / nmol
    table = np.column_stack((extent, peq))

    return table

table_small = simulate_system(20)
table_large = simulate_system(200)

plt.figure()
plt.plot(table_large[:,0], table_large[:,1],'^')
plt.figure()
plt.plot(table_small[:,0], table_small[:,1],'^')
plt.show(block=False)

with open('prob_dis_scale.dat', 'w') as f:
    # First block: Mean square displacement
    np.savetxt(f, table_small, fmt='%f', header="table_small")
    f.write("\n\n")
    np.savetxt(f, table_large, fmt='%f', header="table_large")