diff --git a/PS11_Jung/akira.py b/PS11_Jung/akira.py new file mode 100644 index 0000000000000000000000000000000000000000..7a36d8408f3a0feba64382d33289729be6283871 --- /dev/null +++ b/PS11_Jung/akira.py @@ -0,0 +1,121 @@ +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 diff --git a/PS11_Jung/jung.py b/PS11_Jung/odyssey.py similarity index 56% rename from PS11_Jung/jung.py rename to PS11_Jung/odyssey.py index 366c6592198108a0eab168542a867d26168fbafd..3b3b151f7dfb9b20ad6a3e927be6ce34bbb9e7af 100644 --- a/PS11_Jung/jung.py +++ b/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') - -# c) -## Determine equilibration time -## TODO: Really strange behaviour !! No equilibration time ?? -plt.figure() -plt.plot(t, E) -plt.show() - -## 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 =}') +# 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 + + + +# ## Pressure +t_eq = 50 +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}')