Update sheet5.py

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")