Skip to content
Snippets Groups Projects
Commit e0e9fb5b authored by jung_42's avatar jung_42
Browse files

Update

parent 072df35b
Branches
No related tags found
No related merge requests found
...@@ -3,16 +3,19 @@ import itertools ...@@ -3,16 +3,19 @@ import itertools
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from tqdm import tqdm from tqdm import tqdm
plt.rcParams.update({ plt.rcParams['axes.grid'] = True
"text.usetex": True, # plt.rcParams.update({
"font.family": "sans-serif", # "text.usetex": True,
"font.sans-serif": "Helvetica", # "font.family": "sans-serif",
}) # "font.sans-serif": "Helvetica",
# })
# # Problem 1 (T=1.5, N=128, rho = 0.4)
# # Problem 1
rho = 0.4
sig = 1 sig = 1
eps = 1 eps = 1
k_BT = 1.5 * eps k_BT = 1.5 * eps
beta = 1 / k_BT
N = 125 N = 125
d=3 d=3
L = 6.786 * sig L = 6.786 * sig
...@@ -37,25 +40,26 @@ def verlet_pot(force, x0, p0, m, dt, t, L): ...@@ -37,25 +40,26 @@ def verlet_pot(force, x0, p0, m, dt, t, L):
x = np.zeros((t, *x0.shape)) x = np.zeros((t, *x0.shape))
p = np.zeros((t, *p0.shape)) p = np.zeros((t, *p0.shape))
f = np.zeros((t, *x0.shape))
U = np.zeros(t) U = np.zeros(t)
x[0] = x0 x[0] = x0
p[0] = p0 p[0] = p0
f, U[0] = force(x[0], L) f[0], U[0] = force(x[0], L)
for i in tqdm(range(1, t)): for i in tqdm(range(1, t)):
p[i] = p[i-1] + 1/2 * f * dt p[i] = p[i-1] + 1/2 * f[i-1] * dt
x[i] = x[i-1] + p[i] / m * dt x[i] = x[i-1] + p[i] / m * dt
# enforce periodic boundaries # enforce periodic boundaries
x[i] -= (x[i] > L / 2) * L x[i] -= (x[i] > L / 2) * L
x[i] += (x[i] < -L / 2) * L x[i] += (x[i] < -L / 2) * L
f, U[i] = force(x[i], L) f[i], U[i] = force(x[i], L)
p[i] = p[i] + 1/2 * f * dt p[i] = p[i] + 1/2 * f[i] * dt
return x, p, U return x, p, U, f
def LJ(r, L): def LJ(r, L):
N = len(r) N = len(r)
...@@ -97,17 +101,66 @@ def LJ(r, L): ...@@ -97,17 +101,66 @@ def LJ(r, L):
# plt.savefig(f'problem1a_{k}_order_2.png') # plt.savefig(f'problem1a_{k}_order_2.png')
# b) # b)
t = np.arange(0,10**3*dt, dt) t_max = 10**4
r, p, U = verlet_pot(LJ, x0, p0, m, dt, 10**3, L) t = np.arange(0,t_max*dt, dt)
delta_p = np.sum(p - p[0], axis=1) 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)
# 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.figure()
plt.plot(t,delta_p[:,0], label=r'$\alpha=x$') plt.plot(t, E)
plt.plot(t,delta_p[:,1], label=r'$\alpha=y$') plt.show()
plt.plot(t,delta_p[:,2], label=r'$\alpha=z$')
plt.xlabel(r'$t / \tau$') ## Pressure
plt.ylabel(r'$\delta p / \frac{\tau \epsilon}{\sigma}$') t_eq = 50
plt.legend() V_eq = (1 / (3 * L**3)) * np.mean(np.sum(r[t_eq:]*f[t_eq:], axis=(1,2)))
plt.tight_layout() print(f'Pressure at equilibrium: {rho*k_BT + V_eq =}')
plt.savefig(f'problem1b.png')
## U_pot/N
U_eq = np.mean(U[t_eq:])/N
print(f'Potential energy per particle at equilibrium: {U_eq}')
## mu_ex
def LJ_pair(r, r_new, L):
N = len(r)
U = 0
sig6 = sig**6
dr = r - r_new
# periodic boundary conditions
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))
return U
r_test = r[int(2 / dt)::int(3 / dt)]
k = 3
samples = np.zeros((len(r_test), 10**k))
for i, x in enumerate(r_test):
for j in range(10**k):
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
print(f'Excess chemical potential at equilibrium: {mu} +/- {dmu}')
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment