From c760d373efb286058e9cc582e6224389de0ccc05 Mon Sep 17 00:00:00 2001
From: Jaslo Ziska <jaslo.ziska@fu-berlin.de>
Date: Mon, 29 May 2023 23:19:45 +0200
Subject: [PATCH] Update sheet5.py
---
Jaslo/sheet5.py | 79 +++++++++++++++++++++++++++++++++----------------
1 file changed, 54 insertions(+), 25 deletions(-)
diff --git a/Jaslo/sheet5.py b/Jaslo/sheet5.py
index b601f0b..5bcd785 100644
--- a/Jaslo/sheet5.py
+++ b/Jaslo/sheet5.py
@@ -258,31 +258,31 @@ plt.savefig("problem5.2a-E.png")
plt.close()
# b)
-k = np.array([1, 2, 3, 4])
-N = 2**(3 * k)
-T = np.zeros(N.shape)
-
-for i in range(len(N)):
- L = 2**k[i] * sig
- t_max = 10 * tau
- t = np.arange(0, t_max, dt)
- m = np.full((N[i], d), 1)
- n = int(np.round(N[i]**(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[0,0] * k_BT), (N[i], d))
-
- T1 = time.time()
- r, p, U = verlet_pot(LJ, r0, p0, m, dt, len(t), L)
- T2 = time.time()
-
- T[i] = T2 - T1
-
-plt.plot(N, T)
-plt.xlabel(r"$N$")
-plt.ylabel(r"$T / s$")
-#plt.savefig("problem5.2b.png")
-plt.close()
+# k = np.array([1, 2, 3, 4])
+# N = 2**(3 * k)
+# T = np.zeros(N.shape)
+
+# for i in range(len(N)):
+# L = 2**k[i] * sig
+# t_max = 10 * tau
+# t = np.arange(0, t_max, dt)
+# m = np.full((N[i], d), 1)
+# n = int(np.round(N[i]**(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[0,0] * k_BT), (N[i], d))
+
+# T1 = time.time()
+# r, p, U = verlet_pot(LJ, r0, p0, m, dt, len(t), L)
+# T2 = time.time()
+
+# T[i] = T2 - T1
+
+# plt.plot(N, T)
+# plt.xlabel(r"$N$")
+# plt.ylabel(r"$T / s$")
+# #plt.savefig("problem5.2b.png")
+# plt.close()
# c)
N = 64
@@ -301,3 +301,32 @@ t_max = 100 * tau
t = np.arange(0, t_max, dt)
r, p, U = verlet_pot(LJ, r0, p0, m, dt, len(t), L)
+
+V = np.sum(r * p, axis=(1, 2))
+
+plt.plot(t, U / N)
+plt.xlabel(r"$t / \tau$")
+plt.ylabel(r"$u / \epsilon$")
+plt.savefig("problem5.2c-U.png")
+plt.close()
+
+plt.plot(t, V / N)
+plt.xlabel(r"$t / \tau$")
+plt.ylabel(r"$v / \epsilon \tau$")
+plt.savefig("problem5.2c-V.png")
+plt.close()
+
+t_eq = 40
+t_max = 1000
+t = np.arange(0, t_max, dt)
+r, p, U = verlet_pot(LJ, r0, p0, m, dt, len(t), L)
+V = np.sum(r * p, axis=(1, 2))
+K = np.sum(p**2 / (2 * m[0][0]), axis=(1,2))
+
+U_a = np.average(U[int(t_eq / dt):]) / N
+V_a = np.average(V[int(t_eq / dt):]) / N
+temp = 2 * np.average(K[int(t_eq / dt):]) / N
+
+print(f"<U> = {U_a} eps")
+print(f"<V> = {V_a}")
+print(f"k_B T = {temp} k_B")
--
GitLab