Add sheet 11

PS11_Jung/PS11.py

+import numpy as np
+import itertools
+import matplotlib.pyplot as plt
+from tqdm import tqdm
+
+N = 64
+d = 3
+sig = 1
+eps = 1
+L = 6 * sig
+m = 1
+tau = np.sqrt(m * sig**2 / eps)
+a = 1.5 * sig
+k_BT = 1.5 * eps
+dt = 0.002 * tau
+gamma = 1
+
+# create lattice
+n = int(np.round(N**(1/3)))
+x = np.linspace(-(n-1) / 2, (n-1) / 2, n)
+r0 = a * np.array(list(itertools.product(x, repeat=3)))
+
+p0 = np.random.normal(0, np.sqrt(m * k_BT), (N, d))
+
+t_max = 100 * tau
+t = np.arange(0, t_max, dt)
+
+def BAOAB(force, x0, p0, m, dt, L):
+    assert(x0.shape == p0.shape)
+
+    x = np.zeros((len(t), *x0.shape))
+    p = np.zeros((len(t), *p0.shape))
+    U = np.zeros(len(t))
+
+    x[0] = x0
+    p[0] = p0
+    xi_2 = np.sqrt(k_BT*(1 - np.exp(-2 * gamma * dt)))
+
+    f, U[0] = force(x[0], L)
+
+    for i in tqdm(range(1, len(t))):
+        r = np.random.normal(0,1,(N,d))
+        p[i] = p[i-1] + 1/2 * f * dt
+        x[i] = x[i-1] + p[i] / (2*m) * dt
+        p[i] = np.exp(-gamma*dt) * p[i] + xi_2 * r/np.sqrt(m)
+        x[i] = x[i] + (dt/2) * p[i]/m
+
+        # enforce periodic boundaries
+        x[i] -= (x[i] > L / 2) * L
+        x[i] += (x[i] < -L / 2) * L
+
+        f, U[i] = force(x[i], L)
+        p[i] = p[i] + 1/2 * f * dt
+
+    return x, p, U
+
+def LJ(r, L):
+    N = len(r)
+    f = np.zeros(r.shape)
+    U = 0
+
+    sig6 = sig**6
+
+    for i in range(1, N):
+        dr = r[i:] - r[:N-i]
+        dr -= (dr > L / 2) * L
+        dr += (dr < -L / 2) * L
+
+        r2 = 1 / np.sum(dr**2, axis=1) # scalar product per element
+        r6 = r2 * r2 * r2
+
+        U += 4 * eps * sig6 * np.sum(r6 * (sig6 * r6 - 1))
+        fval = np.transpose(-24 * eps * sig6 * r6 * r2 * (2 * sig6 * r6 - 1) * dr.T) # cursed scalar multiplication
+
+        f[:N-i] += fval
+        f[i:] -= fval
+
+    return f, U
+
+# E = 1 / (2 * m) * np.sum(p**2, axis=(1,2)) + U1
+# plt.plot(E)
+# plt.show()
+
+# a)
+r, p, U1 = BAOAB(LJ, r0, p0, m, dt, L)
+equil_t = int(20*tau)
+r = r[equil_t:]
+ms_displacement = np.average(np.linalg.norm(r - r0, axis=2)**2, axis=1)
+plt.plot(ms_displacement)
+plt.yscale('log')
+plt.xscale('log')
+plt.show()
\ No newline at end of file