From 02fed343f098ec35856381e66176e56fa07f0ffa Mon Sep 17 00:00:00 2001
From: Sebastian Kloth <sebastian.kloth@pkm.tu-darmstadt.de>
Date: Tue, 23 Jan 2024 11:01:43 +0100
Subject: [PATCH] Fixed rdf

---
 src/mdevaluate/coordinates.py  | 8 ++++++--
 src/mdevaluate/distribution.py | 7 +++++--
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/src/mdevaluate/coordinates.py b/src/mdevaluate/coordinates.py
index 1e3b623..9841875 100755
--- a/src/mdevaluate/coordinates.py
+++ b/src/mdevaluate/coordinates.py
@@ -646,7 +646,7 @@ def next_neighbors(
         return distances_new, indices_new
     else:
         atoms_pbc, atoms_pbc_index = pbc_points(
-            query_atoms, box, thickness=distance_upper_bound + 0.1, index=True, **kwargs
+            atoms, box, thickness=distance_upper_bound + 0.1, index=True, **kwargs
         )
         tree = KDTree(atoms_pbc)
         distances, indices = tree.query(
@@ -692,12 +692,16 @@ def number_of_neighbors(
     elif not distinct:
         dnn = 1
 
+    print(atoms[:5])
+    print(query_atoms[:5])
+    print(" ")
+
     box = atoms.box
     if np.all(np.diag(np.diag(box)) == box):
         atoms = atoms % np.diag(box)
         tree = KDTree(atoms, boxsize=np.diag(box))
     else:
-        atoms_pbc = pbc_points(query_atoms, box, thickness=r_max + 0.1, **kwargs)
+        atoms_pbc = pbc_points(atoms, box, thickness=r_max + 0.1, **kwargs)
         tree = KDTree(atoms_pbc)
 
     num_of_neighbors = tree.query_ball_point(query_atoms, r_max, return_length=True)
diff --git a/src/mdevaluate/distribution.py b/src/mdevaluate/distribution.py
index cb20d21..6aa09fd 100644
--- a/src/mdevaluate/distribution.py
+++ b/src/mdevaluate/distribution.py
@@ -126,6 +126,9 @@ def rdf(
     particles_in_volume = int(
         np.max(number_of_neighbors(atoms_a, query_atoms=atoms_b, r_max=bins[-1])) * 1.1
     )
+    print(atoms_a[:5])
+    print(atoms_b[:5])
+    print(" ")
     distances, indices = next_neighbors(
         atoms_a,
         atoms_b,
@@ -147,9 +150,9 @@ def rdf(
         distances = [d for dist in distances for d in dist]
 
     hist, bins = np.histogram(distances, bins=bins, range=(0, bins[-1]), density=False)
-    hist = hist / len(atoms_a)
+    hist = hist / len(atoms_b)
     hist = hist / (4 / 3 * np.pi * bins[1:] ** 3 - 4 / 3 * np.pi * bins[:-1] ** 3)
-    n = len(atoms_b) / np.prod(np.diag(atoms_b.box))
+    n = len(atoms_a) / np.prod(np.diag(atoms_a.box))
     hist = hist / n
 
     return hist