Update sheet6.py

Jaslo/sheet6.py

@@ -137,15 +137,19 @@ def LJ(r: ArrayLike, L: float) -> (ArrayLike, float):
 r, p, U = verlet_thermostat(LJ, r0, p0, m, dt, len(t), L, k_BT, int(0.1 / dt))
 T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2))
 
+# plt.plot(t, T)
 # plt.plot(t, U)
 # plt.show()
 
 t_eq = 30 * tau
-dT2 = (T[int(t_eq / dt):] - T[int(t_eq / dt)])**2
-dU2 = (U[int(t_eq / dt):] - U[int(t_eq / dt)])**2
+dT = T[int(t_eq / dt):] - np.mean(T[int(t_eq / dt):])
+c_V = 3/2 / (1 - 3/2 * N * np.mean((dT / T[int(t_eq / dt):])**2))
+dc_V = 3/2 / (1 - 3/2 * N * np.std((dT / T[int(t_eq / dt):])**2))
+dU = U[int(t_eq / dt):] - np.mean(U[int(t_eq / dt):])
 
-print(f"<dT²> = {np.average(dT2)} ± {np.std(dT2)} ϵ")
-print(f"<dU²> = {np.average(dU2)} ± {np.std(dU2)} ϵ")
+print(f"<dT²> = {np.average(dT**2)} ± {np.std(dT**2)} ϵ")
+print(f"c_V = {c_V} ± {dc_V} k_B")
+print(f"<dU²> = {np.average(dU**2)} ± {np.std(dU**2)} ϵ")
 
 # b)
 # i)
@@ -172,11 +176,11 @@ print(f"k_B T = {temp}")
 r, p, U = verlet(LJ, r[-1], p[-1], m, dt, int(500 / dt), L)
 T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2))
 
-dT2 = (T - T[0])**2
-dU2 = (U - U[0])**2
-
-print(f"<dT²> = {np.average(dT2)} ± {np.std(dT2)} ϵ")
-print(f"<dU²> = {np.average(dU2)} ± {np.std(dU2)} ϵ")
+H = U + T
+dH = H - np.mean(H)
+c_V = np.mean(dH**2) / (N * k_BT**2)
+dc_V = bp.std(dH**2) / (N * k_BT**2)
+print(f"c_V = {c_V} ± {dc_V} k_B")
 
 # c)
 k_BT = 1.35
@@ -203,8 +207,8 @@ print(f"k_B T = {temp}")
 r, p, U = verlet(LJ, r[-1], p[-1], m, dt, int(500 / dt), L)
 T = 1 / (2 * m0) * np.sum(p**2, axis=(1, 2))
 
-dT2 = (T - T[0])**2
-dU2 = (U - U[0])**2
-
-print(f"<dT²> = {np.average(dT2)} ± {np.std(dT2)} ϵ")
-print(f"<dU²> = {np.average(dU2)} ± {np.std(dU2)} ϵ")
+H = U + T
+dH = H - np.mean(H)
+c_V = np.mean(dH**2) / (N * k_BT**2)
+dc_V = bp.std(dH**2) / (N * k_BT**2)
+print(f"c_V = {c_V} ± {dc_V} k_B")