Update sheet 9

Jaslo/sheet9.py

@@ -82,15 +82,14 @@ def LJ_pair(r: ArrayLike, r_new: ArrayLike, L: float) -> float:
 
     sig6 = sig**6
 
-    for i in range(N):
-        dr = r[i:] - 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
+    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))
+    U = 4 * eps * sig6 * np.sum(r6 * (sig6 * r6 - 1))
 
     return U
 
@@ -130,7 +129,7 @@ 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, ), L))
+        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)
@@ -138,4 +137,5 @@ std = np.std(samples)
 mu = - np.log(mean) / beta
 dmu = np.abs(1 / (mean * beta)) * std
 
+print(mu)
 print(np.abs(dmu / mu))