Update

PS11_Jung/akira_1.py

+import numpy as np
+import matplotlib.pyplot as plt
+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
+dt = 0.0005
+t_max = 10**6*dt
+
+t = np.arange(0,t_max,dt)
+M = N*m
+
+# 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))
+#     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])
+
+#     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 # 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
+
+#         f, U[i] = force(x[i])
+#         p[i] = p[i] + 1/2 * f * dt
+
+#     return x, p, U
+
+
+# def fene(r):
+#     dr = r[1:] - r[:-1]
+#     dr_abs = np.linalg.norm(dr, axis=1)
+
+#     fval = -K*dr / (1-(dr_abs / r_max)**2)[:,np.newaxis]
+
+#     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
+
+
+# r0= np.zeros((N,3))
+# r0[:, 0] = 0.5*np.arange(N)
+# p0 = np.random.normal(0, np.sqrt(m * k_BT), (N, d))
+
+# r, p, U = BAOAB(fene, r0,p0, m, dt, N)
+# E = 1 / (2 * m) * np.sum(p**2, axis=(1,2)) + U
+# plt.plot(t,E)
+# plt.show()
+
+# t_eq = 100
+# t_eq = int(t_eq/dt)
+# r_cm = 1/M * np.sum(m * r, axis=1)
+# r_cm = np.tile(r_cm[:, np.newaxis, :], (1, N, 1))
+# R_g = (1/M * np.sum(np.linalg.norm(r-r_cm, axis=2)**2, axis=1))**(1/2)
+# R_g_eq = np.mean(R_g[t_eq:])
+# d_R_g_eq = 0.5 / R_g_eq * np.std((1/M * np.sum(np.linalg.norm(r-r_cm, axis=2)**2, axis=1))**(1/2))
+
+# R_e_eq = np.mean(np.linalg.norm(r[t_eq:,-1] - r[t_eq:,0], axis=1)**2)**(1/2)
+# d_R_e_eq = 0.5 / R_e_eq * np.std(np.linalg.norm(r[t_eq:,-1] - r[t_eq:,0], axis=1)**2)
+
+# T = 1 / (2 * m) * np.sum(p**2, axis=(1, 2))
+
+# H = U[t_eq:] + T[t_eq:]
+# dH = H - np.mean(H)
+# c_V = np.mean(dH**2) / (N * k_BT**2)
+# dc_V = np.std(dH**2) / (N * k_BT**2)
+
+
+# print(f'R_g at equilibrium: {R_g_eq} +/- {d_R_g_eq*100/R_g_eq} %')
+# print(f'R_e at equilibrium: {R_e_eq} +/- {d_R_e_eq*100/R_e_eq} %')
+# print(f"c_V = {c_V} ± {dc_V*100/c_V} k_B")
+
+
+# b)
+L = N * r_max/3
+
+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
+
+def BAOAB_pbc(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
+        # enforce periodic boundaries
+        x[i] -= (x[i] > L / 2) * L
+        x[i] += (x[i] < -L / 2) * L
+
+        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
+
+r0_0= np.zeros((N,3))
+r0_0[:, 0] = 0.5*np.arange(N)
+p0_0 = np.random.normal(0, np.sqrt(m * k_BT), (N, d))
+
+r, p, U = BAOAB_pbc(fene_chain, r0_0,p0_0, m, dt, N)
+
+
+
+## mean square displacement
+msdv_0 = np.zeros(len(r))
+for i in tqdm(range(1, len(msdv_0))):
+    # msdv_0[i] = np.mean(np.sum((r[i:] - r[:-i])**2, axis=2))
+    msdv_0[i] = np.mean(np.linalg.norm(r[i:] - r[:-i], axis=2)**2)    
+
+
+N=96
+r0_1= np.zeros((N,3))
+r0_1[:, 0] = 0.5*np.arange(N)
+p0_1 = np.random.normal(0, np.sqrt(m * k_BT), (N, d))
+
+r1, p1, U1 = BAOAB_pbc(fene_chain, r0_1,p0_1, m, dt, N)
+msdv_1 = np.zeros(len(r1))
+for i in tqdm(range(1, len(msdv_1))):
+    # msdv_1[i] = np.mean(np.sum((r1[i:] - r1[:-i])**2, axis=2))
+    msdv_1[i] = np.mean(np.linalg.norm(r1[i:] - r1[:-i], axis=2)**2)    
+
+plt.figure()
+plt.loglog(t, msdv_0, label=r'$N=48$')
+plt.loglog(t, msdv_1, label=r'$N=96$')
+
+plt.xlabel(r'$\delta r^2(t) / r_{max}^{2}$')
+plt.ylabel(r'$t / \tau$')
+plt.legend()
+plt.tight_layout()
+plt.savefig('2b_N_48.pdf')
\ No newline at end of file

PS11_Jung/odyssey.py

@@ -160,17 +160,34 @@ def LJ_pair(r, r_new, L):
 t_max = 10**3
 t = np.arange(0,t_max*dt, dt)
 r_bef, p_bef, 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)
-# # 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')
+# E = 1 / (2 * m) * np.sum(p_bef**2, axis=(1,2)) + U
+
+# # root mean square deviation of energy
+# rmsq_dev_E = np.sqrt(np.mean((E-E[0])**2) / np.mean(E)**2)
+# # root mean square deviation of Impuls
+rmsq_dev_p = np.sqrt(np.mean(np.linalg.norm(p_bef - p_bef[0])**2))
+deta_p =np.abs(p_bef - p_bef[0])
+dev_p = np.einsum('...i,...i->...', deta_p, deta_p)
+rmsq_dev_p = np.sqrt(np.mean(dev_p))
+# print(f'{rmsq_dev_E=}')
+rmsq_dev_p_2 = (np.mean(np.linalg.norm(p_bef-p_bef[0], axis=(1,2))**2))**(1/2)
+rmsq_dev_p_2 = (np.mean(np.linalg.norm(np.sum(p_bef-p_bef[0],axis=1), axis=1)**2))**(1/2)
+
+print(f'{rmsq_dev_p=}')
+print(f'{rmsq_dev_p_2=}')
+
+# delta_p = np.sum(p_bef - p_bef[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
@@ -185,41 +202,51 @@ r_bef, p_bef, U, f = verlet_pot(LJ, x0, p0, m, dt, t_max, L)
 t_eq = 100
 t_max = 10**4
 
-
-no_of_simulations = 6
-p_ensemble = np.zeros(no_of_simulations)
-std_p_ensemble = np.zeros(no_of_simulations)
-U_ensemble = np.zeros(no_of_simulations)
-std_U_ensemble = np.zeros(no_of_simulations)
-mu_ensemble = np.zeros(no_of_simulations)
-std_mu_ensemble = np.zeros(no_of_simulations)
-
-## for different ensembles with different seeds !
-for i in range(no_of_simulations):
-    random.seed(i)
-    p0 = np.random.normal(0, np.sqrt(m * k_BT), (N, d))    
-    r,p,U,f = verlet_thermostat(LJ, r_bef[t_eq], p0, m, dt, t_max,L,int(1/dt))
-    V = np.sum(r*f, axis=(1,2))
-    p_samples = rho*k_BT + (1 / (3 * L**3))* V
-    U_by_N = U/N
-
-    p_ensemble[i] = np.mean(p_samples)
-    std_p_ensemble[i] = np.std(p_samples)
-    U_ensemble[i] = np.mean(U_by_N)
-    std_U_ensemble[i] = np.std(U_by_N)
-
-    #mu_ex
-    k = 7
-    r_test = r[int(1 / dt)::int(3 / dt)]
-
-    samples = np.zeros((len(r_test), 10**k))
-    for y, 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[y,j] = np.exp(-beta * LJ_pair(x, np.random.uniform(-L/2, L/2, d), L))
-
-    mean = np.mean(samples)
-    std = np.std(samples)
+r,p,U,f = verlet_thermostat(LJ, r_bef[t_eq], p0, m, dt, t_max,L,int(1/dt))
+V = np.sum(r*f, axis=(1,2))
+E_kin = 1 / (2 * m) * np.sum(p**2, axis=(1,2))
+# plt.figure()
+# plt.plot(np.arange(0,t_max,t_max/len(V)), V)
+# plt.xlabel(r'$t / \tau$')
+# plt.ylabel(r'$V / \epsilon \tau$')
+# plt.tight_layout()
+# plt.savefig('virial_function_1_protocol.pdf')
+
+
+# no_of_simulations = 6
+# p_ensemble = np.zeros(no_of_simulations)
+# std_p_ensemble = np.zeros(no_of_simulations)
+# U_ensemble = np.zeros(no_of_simulations)
+# std_U_ensemble = np.zeros(no_of_simulations)
+# mu_ensemble = np.zeros(no_of_simulations)
+# std_mu_ensemble = np.zeros(no_of_simulations)
+
+# ## for different ensembles with different seeds !
+# for i in range(no_of_simulations):
+#     random.seed(i)
+#     p0 = np.random.normal(0, np.sqrt(m * k_BT), (N, d))    
+#     r,p,U,f = verlet_thermostat(LJ, r_bef[t_eq], p0, m, dt, t_max,L,int(1/dt))
+#     V = np.sum(r*f, axis=(1,2))
+#     p_samples = rho*k_BT + (1 / (3 * L**3))* V
+#     U_by_N = U/N
+
+#     p_ensemble[i] = np.mean(p_samples)
+#     std_p_ensemble[i] = np.std(p_samples)
+#     U_ensemble[i] = np.mean(U_by_N)
+#     std_U_ensemble[i] = np.std(U_by_N)
+
+#     #mu_ex
+#     k = 5
+#     r_test = r[int(1 / dt)::int(3 / dt)]
+
+#     samples = np.zeros((len(r_test), 10**k))
+#     for y, 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[y,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
     # std_mu = np.abs(1 / (mean * beta)) * std
@@ -228,14 +255,14 @@ for i in range(no_of_simulations):
 
 # a = [-1.851747546190128, -1.878954082130441, -1.9081508511080745]
 # a_std = [5.412283418498762, 5.93628399967974, 5.53353222778447] 
-    mu_ensemble[i] = - np.log(mean) / beta
-    std_mu_ensemble[i] = np.abs(1 / (mean * beta)) * std
+#     mu_ensemble[i] = - np.log(mean) / beta
+#     std_mu_ensemble[i] = np.abs(1 / (mean * beta)) * std
 
 
-print(f'Pressure at equilibrium: {np.mean(p_ensemble)} +/- {np.std(std_p_ensemble)*100/np.mean(p_ensemble)} %')
-print(f'Internal energy per particle at equilibrium: {np.mean(U_ensemble)} +/- {np.std(std_U_ensemble)*100/np.mean(U_ensemble)} %')
-print(f'Excess chemical potential at equilibrium: {np.mean(mu_ensemble)} +/- {np.std(std_mu_ensemble)*100/np.mean(mu_ensemble)} %')
-print(f'Excess chemical potential at equilibrium_2: {np.mean(mu_ensemble)} +/- {np.std(mu_ensemble)*100/np.mean(mu_ensemble)} %')
+# print(f'Pressure at equilibrium: {np.mean(p_ensemble)} +/- {np.std(std_p_ensemble)*100/np.mean(p_ensemble)} %')
+# print(f'Internal energy per particle at equilibrium: {np.mean(U_ensemble)} +/- {np.std(std_U_ensemble)*100/np.mean(U_ensemble)} %')
+# print(f'Excess chemical potential at equilibrium: {np.mean(mu_ensemble)} +/- {np.std(std_mu_ensemble)*100/np.mean(mu_ensemble)} %')
+# print(f'Excess chemical potential at equilibrium_2: {np.mean(mu_ensemble)} +/- {np.std(mu_ensemble)*100/np.mean(mu_ensemble)} %')
 
 
 ## mu_ex