using pbc_diff now in tetrahedral_order parameter calculation, since reference positions will not be periodic images in the default use case

Reviewed-on: #2

src/mdevaluate/checksum.py

@@ -73,7 +73,9 @@ def checksum(*args, csum=None):
         elif isinstance(arg, FunctionType):
             csum.update(strip_comments(inspect.getsource(arg)).encode())
             c = inspect.getclosurevars(arg)
-            for v in {**c.nonlocals, **c.globals}.values():
+            merged = {**c.nonlocals, **c.globals}
+            for key in sorted(merged):  # deterministic ordering
+                v = merged[key]
                 if v is not arg:
                     checksum(v, csum=csum)
         elif isinstance(arg, functools.partial):

src/mdevaluate/distribution.py

@@ -182,10 +182,10 @@ def tetrahedral_order(
     )
 
     # Connection vectors
-    neighbors_1 -= atoms
+    neighbors_1 = pbc_diff(neighbors_1, atoms, box=atoms.box)
-    neighbors_2 -= atoms
+    neighbors_2 = pbc_diff(neighbors_2, atoms, box=atoms.box)
-    neighbors_3 -= atoms
+    neighbors_3 = pbc_diff(neighbors_3, atoms, box=atoms.box)
-    neighbors_4 -= atoms
+    neighbors_4 = pbc_diff(neighbors_4, atoms, box=atoms.box)
 
     # Normed Connection vectors
     neighbors_1 /= np.linalg.norm(neighbors_1, axis=-1).reshape(-1, 1)

src/mdevaluate/extra/free_energy_landscape.py

@@ -4,7 +4,6 @@ from typing import Optional
 import numpy as np
 from numpy.typing import ArrayLike, NDArray
 from numpy.polynomial.polynomial import Polynomial as Poly
-import math
 from scipy.spatial import KDTree
 import pandas as pd
 import multiprocessing as mp
@@ -49,7 +48,7 @@ def _pbc_points_reduced(
 
 def _build_tree(points, box, r_max, pore_geometry):
     if np.all(np.diag(np.diag(box)) == box):
-        tree = KDTree(points, boxsize=box)
+        tree = KDTree(points % box, boxsize=box)
         points_pbc_index = None
     else:
         points_pbc, points_pbc_index = _pbc_points_reduced(
@@ -79,11 +78,7 @@ def occupation_matrix(
     z_bins = np.arange(0, box[2][2] + edge_length, edge_length)
     bins = [x_bins, y_bins, z_bins]
     # Trajectory is split for parallel computing
-    size = math.ceil(len(frame_indices) / nodes)
+    indices = np.array_split(frame_indices, nodes)
-    indices = [
-        np.arange(len(frame_indices))[i : i + size]
-        for i in range(0, len(frame_indices), size)
-    ]
     pool = mp.Pool(nodes)
     results = pool.map(
         partial(_calc_histogram, trajectory=trajectory, bins=bins), indices
@@ -277,7 +272,11 @@ def distance_resolved_energies(
 
 
 def find_energy_maxima(
-    energy_df: pd.DataFrame, r_min: float, r_max: float
+    energy_df: pd.DataFrame,
+    r_min: float,
+    r_max: float,
+    r_eval: float = None,
+    degree: int = 2,
 ) -> pd.DataFrame:
     distances = []
     energies = []
@@ -286,6 +285,9 @@ def find_energy_maxima(
         x = np.array(data_d["r"])
         y = np.array(data_d["energy"])
         mask = (x >= r_min) * (x <= r_max)
-        p3 = Poly.fit(x[mask], y[mask], deg=2)
+        p3 = Poly.fit(x[mask], y[mask], deg=degree)
+        if r_eval is None:
             energies.append(np.max(p3(np.linspace(r_min, r_max, 1000))))
+        else:
+            energies.append(p3(r_eval))
     return pd.DataFrame({"d": distances, "energy": energies})

test/test_free_energy_landscape.py

@@ -13,19 +13,7 @@ def trajectory(request):
 
 
 def test_get_fel(trajectory):
-    test_array = np.array(
+    test_array = np.array([210., 214., 209., 192., 200., 193., 230., 218., 266.])
-        [
-            184.4909136,
-            233.79320471,
-            223.12003988,
-            228.49746397,
-            200.5626769,
-            212.82484221,
-            165.10818396,
-            170.74123681,
-            175.86672931,
-        ]
-    )
 
     OW = trajectory.subset(atom_name="OW")
     box = trajectory[0].box