Update

PS11_Jung/akira.py

+import numpy as np
+from tqdm import tqdm
+
+
+N = 48
+k = 1 # spring constant
+m = 1
+gamma = 0.2 
+k_BT = 0.1
+d = 3
+r_max = 1
+eps = 1
+
+## a)
+def BAOAB(force, x0, p0, m, dt, N):
+    assert(x0.shape == p0.shape)
+
+    x = np.zeros((len(t), *x0.shape))
+    p = np.zeros((len(t), *p0.shape))
+    x[0] = x0
+    p[0] = p0
+    xi_2 = np.sqrt(k_BT*(1 - np.exp(-2 * gamma * dt)))
+
+    for i in range(1, len(t)):
+        r = np.random.normal(0,1,(N,d))
+        p[i] = p[i-1] + 1/2 * force(x[i-1]) * dt # same
+        x[i] = x[i-1] + p[i] / (2*m) * dt # diff
+        p[i] = np.exp(-gamma*dt) * p[i] + xi_2 * r * np.sqrt(m)
+        x[i] = x[i] + (dt/2) * p[i]/m
+        p[i] = p[i] + 1/2 * force(x[i]) * dt
+
+    return x, p
+
+# 1D
+def FENE(x):
+    r1 = x[:3]
+    r2 = x[3:]
+
+    f1 = -epsilon * (r1 - r2) / (a**2 - np.linalg.norm(r1 - r2)**2)
+    f2 = -f1
+
+    return np.array([
+        f1,
+        f2,
+    ]).flatten()
+
+
+## b) 
+
+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
+
+
+
+
+##### noch mal anschauen
+def fene(r, L):
+    dr = r[1-] - r[:-1]
+    dr -= (dr > L/2) * L
+    dr += (dr < -L/2) * L
+    assert np.all(np.abs(dr) <= L/2)
+
+    dr_abs = np.linalg.norm(dr, axis=1)
+    
+    fval = -K*dr / (1-(dr_abs / r_max)**2)[:,None]
+
+    f = np.zeros(r.shape)
+    f[:-1] = -fval
+    f[1:] += fval
+
+    U = -1/2 * K * r_max**2 *np.sum(np.log(1 - (dr_abs / r_max)**2))
+
+    return f, U
+
+def fene_chain(r, L):
+    dr = np.roll(r, -1, axis=0) -r
+    dr -= (dr > L/2) * L
+    dr += (dr < -L/2) * L
+    assert np.all(np.abs(dr) <= L/2)
+
+    dr_abs = np.linalg.norm(dr, axis=1)
+    
+    fval = -K*dr / (1-(dr_abs / r_max)**2)[:,None]
+
+    f = np.zeros(r.shape)
+    f = -fval.copy()
+    f += np.roll(fval, 1, axis=0)
+
+    U = -1/2 * K * r_max**2 *np.sum(np.log(1 - (dr_abs / r_max)**2))
+
+    return f, U
+
+# R_g = 1.3
+# R_l = 3.3
+# c_v = 2.5
\ No newline at end of file

PS11_Jung/jung.py → PS11_Jung/odyssey.py

@@ -84,16 +84,15 @@ def LJ(r, L):
 
     return f, U
 
-# a)
+# a) p=2 stimmt übereinander für alle delta_t
 # for k in range(4):
 #     dt_prime = 2**k * dt
-#     t_prime = 10**3
-#     r, p, U = verlet_pot(LJ, x0, p0, m, dt_prime, t_prime, L)
+#     t = np.arange(0,10**3*dt, dt_prime)
+#     r, p, U, f = verlet_pot(LJ, x0, p0, m, dt_prime, len(t), L)
 #     E = 1 / (2 * m) * np.sum(p**2, axis=(1,2)) + U
 #     delta_E= E - E[0]
 #     plt.figure()
-#     plt.plot(np.arange(0,10**3*dt_prime, dt_prime), delta_E / (dt_prime)**2 )
-#     # plt.plot(np.arange(0,10**3*dt_prime, dt_prime), delta_E  )
+#     plt.plot(t, delta_E / (dt_prime)**2 )
 
 #     plt.ylabel(r'$\frac{\Delta E(t)}{(\Delta t)} / \frac{\epsilon}{\tau}$')
 #     plt.xlabel(r'$t / \tau$')
@@ -105,38 +104,85 @@ t_max = 10**4
 t = np.arange(0,t_max*dt, dt)
 r, p, U, f = verlet_pot(LJ, x0, p0, m, dt, t_max, L)
 E = 1 / (2 * m) * np.sum(p**2, axis=(1,2)) + U
-# delta_p = np.sum(p - p[0], axis=1)
+# # delta_p = np.sum(p - p[0], axis=1)
+# # plt.figure()
+# # plt.plot(t,delta_p[:,0], label=r'$\alpha=x$')
+# # plt.plot(t,delta_p[:,1], label=r'$\alpha=y$')
+# # plt.plot(t,delta_p[:,2], label=r'$\alpha=z$')
+# # plt.xlabel(r'$t / \tau$')
+# # plt.ylabel(r'$\delta p / \frac{\tau \epsilon}{\sigma}$')
+# # plt.legend()
+# # plt.tight_layout()
+# # plt.savefig(f'problem1b.png')
+
+# # c) 
+# ## Determine equilibration time
 # plt.figure()
-# plt.plot(t,delta_p[:,0], label=r'$\alpha=x$')
-# plt.plot(t,delta_p[:,1], label=r'$\alpha=y$')
-# plt.plot(t,delta_p[:,2], label=r'$\alpha=z$')
-# plt.xlabel(r'$t / \tau$')
-# plt.ylabel(r'$\delta p / \frac{\tau \epsilon}{\sigma}$')
-# plt.legend()
-# plt.tight_layout()
-# plt.savefig(f'problem1b.png')
+# plt.plot(t, E)
+# plt.show()
+
+def verlet_thermostat(
+    force,
+    r0,
+    p0,
+    m,
+    dt,
+    t,
+    L,
+    N_T
+):
+    assert(r0.shape == p0.shape)
+
+    r = np.zeros((t, *r0.shape))
+    p = np.zeros((t, *p0.shape))
+    f = np.zeros((t, *r0.shape))
+    U = np.zeros(t)
+
+    r[0] = r0
+    p[0] = p0
+
+    f[0], U[0] = force(r[0], L)
+
+    for i in tqdm(range(1, t)):
+        # first half step
+        p[i] = p[i-1] + 1/2 * f[i-1] * dt
+        r[i] = r[i-1] + p[i] / m * dt
+
+        # enforce periodic boundaries
+        r[i] -= (r[i] > L / 2) * L
+        r[i] += (r[i] < -L / 2) * L
+
+        # recompute force
+        f[i], U[i] = force(r[i], L)
+        # second half step
+        p[i] = p[i] + 1/2 * f[i] * dt
+
+        # thermostat (every N_T steps):
+        if i % N_T == 0:
+            p[i] = np.random.normal(0, np.sqrt(m * k_BT), p0.shape)
+
+    return r, p, U, f
+
 
-# c) 
-## Determine equilibration time
-## TODO: Really strange behaviour !! No equilibration time ??
-plt.figure()
-plt.plot(t, E)
-plt.show()
 
-## Pressure 
+# ## Pressure 
 t_eq = 50
-V_eq = (1 / (3 * L**3)) * np.mean(np.sum(r[t_eq:]*f[t_eq:], axis=(1,2)))
-print(f'Pressure at equilibrium: {rho*k_BT + V_eq =}')
+t_max = 10**5
+r,p,U,f = verlet_thermostat(LJ, r[50], p[50], m, dt, t_max,L,int(1/dt))
+# V = np.sum(r[t_eq:]*f[t_eq:], axis=(1,2))
+# p_int = (1 / (3 * L**3))* np.mean(V)
+V = np.sum(r*f, axis=(1,2))
+p_int = (1 / (3 * L**3))* np.mean(V)
+print(f'Pressure at equilibrium: {rho*k_BT + p_int} +/- {np.std(V)/(3*L**3)}')
 
 ## U_pot/N
-U_eq = np.mean(U[t_eq:])/N
-print(f'Potential energy per particle at equilibrium: {U_eq}')
+# U_eq = np.mean(U[t_eq:])/N
+U_eq = np.mean(U)/N
+
+print(f'Potential energy per particle at equilibrium: {U_eq} +/- {np.std(U[t_eq:]) / N}')
 
 ## mu_ex
 def LJ_pair(r, r_new, L):
-    N = len(r)
-    U = 0
-
     sig6 = sig**6
 
     dr = r - r_new
@@ -150,17 +196,19 @@ def LJ_pair(r, r_new, L):
 
     return U
 
-r_test = r[int(2 / dt)::int(3 / dt)]
+r_test = r[int(1 / dt)::int(3 / dt)]
 
-k = 3
+k = 4
 samples = np.zeros((len(r_test), 10**k))
 for i, x in enumerate(r_test):
     for j in range(10**k):
+        # while np.exp(-beta * LJ_pair(x, np.random.uniform(-L/2, L/2, d), L)) == 0:
         samples[i,j] = np.exp(-beta * LJ_pair(x, np.random.uniform(-L/2, L/2, d), L))
 
-mean = np.mean(samples)
-std = np.std(samples)
 
-mu = - np.log(mean) / beta
-dmu = np.abs(1 / (mean * beta)) * std
+mean = np.mean(samples, axis=1)
+std = np.std(mean)
+
+mu = - np.log(np.mean(mean)) / beta
+dmu = np.abs(1 / (np.mean(mean) * beta)) * std
 print(f'Excess chemical potential at equilibrium: {mu} +/- {dmu}')