stop reading many frames to get times for shifted_correlation

Reviewed-on: #5

pyproject.toml

@@ -4,11 +4,12 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "mdevaluate"
-version = "24.01"
+version = "24.02"
 dependencies = [
     "mdanalysis",
     "pandas",
     "dask",
     "pathos",
-    "tables"
+    "tables",
+    "pyedr"
 ]

requirements.txt

@@ -3,4 +3,5 @@ pandas
 dask
 pathos
 tables
-pytest
+pytest
+pyedr

src/mdevaluate/__init__.py

@@ -16,7 +16,7 @@ from . import reader
 from . import system
 from . import utils
 from . import extra
-from .logging import logger
+from .logging_util import logger
 
 
 def open(

src/mdevaluate/autosave.py

@@ -5,7 +5,7 @@ from typing import Optional, Callable, Iterable
 
 import numpy as np
 from .checksum import checksum
-from .logging import logger
+from .logging_util import logger
 
 autosave_directory: Optional[str] = None
 load_autosave_data = False
@@ -166,8 +166,10 @@ def autosave_data(
         @functools.wraps(function)
         def autosave(*args, **kwargs):
             description = kwargs.pop("description", "")
+            autosave_dir_overwrite = kwargs.pop("autosave_dir_overwrite", None)
+            autosave_dir = autosave_dir_overwrite if autosave_dir_overwrite is not None else autosave_directory
             autoload = kwargs.pop("autoload", True) and load_autosave_data
-            if autosave_directory is not None:
+            if autosave_dir is not None:
                 relevant_args = list(args[:nargs])
                 if kwargs_keys is not None:
                     for key in [*posargs_keys, *kwargs_keys]:

src/mdevaluate/checksum.py

@@ -1,9 +1,14 @@
 import functools
 import hashlib
-from .logging import logger
+from .logging_util import logger
 from types import ModuleType, FunctionType
 import inspect
 from typing import Iterable
+import ast
+import io
+import tokenize
+import re
+import textwrap
 
 import numpy as np
 
@@ -28,19 +33,46 @@ def version(version_nr: int, calls: Iterable = ()):
     return decorator
 
 
-def strip_comments(s: str):
-    """Strips comment lines and docstring from Python source string."""
-    o = ""
-    in_docstring = False
-    for l in s.split("\n"):
-        if l.strip().startswith(("#", '"', "'")) or in_docstring:
-            in_docstring = l.strip().startswith(('"""', "'''")) + in_docstring == 1
+def strip_comments(source: str) -> str:
+    """Removes docstrings, comments, and irrelevant whitespace from Python source code."""
+
+    # Step 1: Remove docstrings using AST
+    def remove_docstrings(node):
+        if isinstance(node, (ast.FunctionDef, ast.AsyncFunctionDef, ast.ClassDef, ast.Module)):
+            if (doc := ast.get_docstring(node, clean=False)):
+                first_stmt = node.body[0]
+                if isinstance(first_stmt, ast.Expr) and isinstance(first_stmt.value, ast.Constant):
+                    node.body.pop(0)  # Remove the docstring entirely
+        for child in ast.iter_child_nodes(node):
+            remove_docstrings(child)
+
+    tree = ast.parse(textwrap.dedent(source))
+    remove_docstrings(tree)
+    code_without_docstrings = ast.unparse(tree)
+
+    # Step 2: Remove comments using tokenize
+    tokens = tokenize.generate_tokens(io.StringIO(code_without_docstrings).readline)
+    result = []
+    last_lineno = -1
+    last_col = 0
+
+    for toknum, tokval, (srow, scol), (erow, ecol), line in tokens:
+        if toknum == tokenize.COMMENT:
             continue
-        o += l + "\n"
-    return o
+        if srow > last_lineno:
+            last_col = 0
+        if scol > last_col:
+            result.append(" " * (scol - last_col))
+        result.append(tokval)
+        last_lineno, last_col = erow, ecol
+
+    code_no_comments = ''.join(result)
+
+    # Step 3: Remove empty lines (whitespace-only or truly blank)
+    return "\n".join([line for line in code_no_comments.splitlines() if line.strip() != ""])
 
 
-def checksum(*args, csum=None):
+def checksum(*args, csum=None, _seen=None):
     """
     Calculate a checksum of any object, by sha1 hash.
 
@@ -60,7 +92,15 @@ def checksum(*args, csum=None):
         csum = hashlib.sha1()
         csum.update(str(SALT).encode())
 
+    if _seen is None:
+        _seen = set()
+
     for arg in args:
+        obj_id = id(arg)
+        if obj_id in _seen:
+            continue
+        _seen.add(obj_id)
+
         if hasattr(arg, "__checksum__"):
             logger.debug("Checksum via __checksum__: %s", str(arg))
             csum.update(str(arg.__checksum__()).encode())
@@ -73,17 +113,19 @@ 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)
+                    checksum(v, csum=csum, _seen=_seen)
         elif isinstance(arg, functools.partial):
             logger.debug("Checksum via partial for %s", str(arg))
-            checksum(arg.func, csum=csum)
+            checksum(arg.func, csum=csum, _seen=_seen)
             for x in arg.args:
-                checksum(x, csum=csum)
+                checksum(x, csum=csum, _seen=_seen)
             for k in sorted(arg.keywords.keys()):
                 csum.update(k.encode())
-                checksum(arg.keywords[k], csum=csum)
+                checksum(arg.keywords[k], csum=csum, _seen=_seen)
         elif isinstance(arg, np.ndarray):
             csum.update(arg.tobytes())
         else:

src/mdevaluate/coordinates.py

@@ -1,6 +1,6 @@
 from functools import partial, wraps
 from copy import copy
-from .logging import logger
+from .logging_util import logger
 from typing import Optional, Callable, List, Tuple
 
 import numpy as np
@@ -218,7 +218,7 @@ class Coordinates:
             self.get_frame.clear_cache()
 
     def __iter__(self):
-        for i in range(len(self))[self._slice]:
+        for i in range(len(self.frames))[self._slice]:
             yield self[i]
 
     @singledispatchmethod
@@ -232,7 +232,7 @@ class Coordinates:
         return sliced
 
     def __len__(self):
-        return len(self.frames)
+        return len(self.frames[self._slice])
 
     def __checksum__(self):
         return checksum(self.frames, self.atom_filter, self._slice, self.mode)
@@ -434,6 +434,34 @@ def center_of_masses(
         ]
     ).T[mask]
     return np.array(positions)
+    
+    
+@map_coordinates
+def center_of_atoms(
+    frame: CoordinateFrame, atom_indices=None, shear: bool = False
+) -> NDArray:
+    if atom_indices is None:
+        atom_indices = list(range(len(frame)))
+    res_ids = frame.residue_ids[atom_indices]
+    if shear:
+        coords = frame[atom_indices]
+        box = frame.box
+        sort_ind = res_ids.argsort(kind="stable")
+        i = np.concatenate([[0], np.where(np.diff(res_ids[sort_ind]) > 0)[0] + 1])
+        coms = coords[sort_ind[i]][res_ids - min(res_ids)]
+        cor = pbc_diff(coords, coms, box)
+        coords = coms + cor
+    else:
+        coords = frame.whole[atom_indices]
+    mask = np.bincount(res_ids)[1:] != 0
+    positions = np.array(
+        [
+            np.bincount(res_ids, weights=c)[1:]
+            / np.bincount(res_ids)[1:]
+            for c in coords.T
+        ]
+    ).T[mask]
+    return np.array(positions)
 
 
 @map_coordinates
@@ -646,7 +674,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(
@@ -697,7 +725,7 @@ def number_of_neighbors(
         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)

src/mdevaluate/correlation.py

@@ -2,7 +2,7 @@ from typing import Callable, Optional
 
 import numpy as np
 from numpy.typing import ArrayLike
-from scipy.special import legendre
+from scipy.special import legendre, jn
 import dask.array as darray
 from functools import partial
 from scipy.spatial import KDTree
@@ -18,7 +18,7 @@ def log_indices(first: int, last: int, num: int = 100) -> np.ndarray:
     return np.unique(np.int_(ls) - 1 + first)
 
 
-@autosave_data(2)
+@autosave_data(nargs=2, kwargs_keys=('selector', 'segments', 'skip', 'window', 'average', 'points',), version=1.0)
 def shifted_correlation(
     function: Callable,
     frames: Coordinates,
@@ -147,7 +147,8 @@ def shifted_correlation(
     num_frames = int(len(frames) * window)
     ls = np.logspace(0, np.log10(num_frames + 1), num=points)
     idx = np.unique(np.int_(ls) - 1)
-    t = np.array([frames[i].time for i in idx]) - frames[0].time
+    dt = round(frames[1].time - frames[0].time, 6) # round to avoid bad floats
+    t = idx * dt
 
     result = np.array(
         [
@@ -183,6 +184,12 @@ def msd(
         displacements = displacements_without_drift(start_frame, end_frame, trajectory)
     if axis == "all":
         return (displacements**2).sum(axis=1).mean()
+    elif axis == "xy" or axis == "yx":
+        return (displacements[:, [0, 1]]**2).sum(axis=1).mean()
+    elif axis == "xz" or axis == "zx":
+        return (displacements[:, [0, 2]]**2).sum(axis=1).mean()
+    elif axis == "yz" or axis == "zy":
+        return (displacements[:, [1, 2]]**2).sum(axis=1).mean()
     elif axis == "x":
         return (displacements[:, 0] ** 2).mean()
     elif axis == "y":
@@ -193,7 +200,7 @@ def msd(
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
 
 
-def isf(
+def isf_raw(
     start_frame: CoordinateFrame,
     end_frame: CoordinateFrame,
     q: float = 22.7,
@@ -210,20 +217,59 @@ def isf(
         displacements = displacements_without_drift(start_frame, end_frame, trajectory)
     if axis == "all":
         distance = (displacements**2).sum(axis=1) ** 0.5
-        return np.sinc(distance * q / np.pi).mean()
+        return np.sinc(distance * q / np.pi)
+    elif axis == "xy" or axis == "yx":
+        distance = (displacements[:, [0, 1]]**2).sum(axis=1) ** 0.5
+        return np.real(jn(0, distance * q))
+    elif axis == "xz" or axis == "zx":
+        distance = (displacements[:, [0, 2]]**2).sum(axis=1) ** 0.5
+        return np.real(jn(0, distance * q))
+    elif axis == "yz" or axis == "zy":
+        distance = (displacements[:, [1, 2]]**2).sum(axis=1) ** 0.5
+        return np.real(jn(0, distance * q))
     elif axis == "x":
         distance = np.abs(displacements[:, 0])
-        return np.mean(np.cos(np.abs(q * distance)))
+        return np.cos(np.abs(q * distance))
     elif axis == "y":
         distance = np.abs(displacements[:, 1])
-        return np.mean(np.cos(np.abs(q * distance)))
+        return np.cos(np.abs(q * distance))
     elif axis == "z":
         distance = np.abs(displacements[:, 2])
-        return np.mean(np.cos(np.abs(q * distance)))
+        return np.cos(np.abs(q * distance))
     else:
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
 
 
+def isf(
+    start_frame: CoordinateFrame,
+    end_frame: CoordinateFrame,
+    q: float = 22.7,
+    trajectory: Coordinates = None,
+    axis: str = "all",
+) -> float:
+    """
+    Incoherent intermediate scattering function averaged over all particles.
+    See isf_raw for details.
+    """
+    return isf_raw(start_frame, end_frame, q=q, trajectory=trajectory, axis=axis).mean()
+
+
+def isf_mean_var(
+    start_frame: CoordinateFrame,
+    end_frame: CoordinateFrame,
+    q: float = 22.7,
+    trajectory: Coordinates = None,
+    axis: str = "all",
+) -> float:
+    """
+    Incoherent intermediate scattering function averaged over all particles and the
+    variance.
+    See isf_raw for details.
+    """
+    values = isf_raw(start_frame, end_frame, q=q, trajectory=trajectory, axis=axis)
+    return values.mean(), values.var()
+
+
 def rotational_autocorrelation(
     start_frame: CoordinateFrame, end_frame: CoordinateFrame, order: int = 2
 ) -> float:
@@ -262,6 +308,12 @@ def van_hove_self(
         vectors = displacements_without_drift(start_frame, end_frame, trajectory)
     if axis == "all":
         delta_r = (vectors**2).sum(axis=1) ** 0.5
+    elif axis == "xy" or axis == "yx":
+        delta_r = (vectors[:, [0, 1]]**2).sum(axis=1) ** 0.5
+    elif axis == "xz" or axis == "zx":
+        delta_r = (vectors[:, [0, 2]]**2).sum(axis=1) ** 0.5
+    elif axis == "yz" or axis == "zy":
+        delta_r = (vectors[:, [1, 2]]**2).sum(axis=1) ** 0.5
     elif axis == "x":
         delta_r = np.abs(vectors[:, 0])
     elif axis == "y":
@@ -409,10 +461,11 @@ def non_gaussian_parameter(
     end_frame: CoordinateFrame,
     trajectory: Coordinates = None,
     axis: str = "all",
+    full_output = False,
 ) -> float:
-    """
+    r"""
     Calculate the non-Gaussian parameter.
-    ..math:
+    .. math:
       \alpha_2 (t) =
         \frac{3}{5}\frac{\langle r_i^4(t)\rangle}{\langle r_i^2(t)\rangle^2} - 1
     """
@@ -421,18 +474,41 @@ def non_gaussian_parameter(
     else:
         vectors = displacements_without_drift(start_frame, end_frame, trajectory)
     if axis == "all":
-        r = (vectors**2).sum(axis=1)
+        r2 = (vectors**2).sum(axis=1)
         dimensions = 3
+    elif axis == "xy" or axis == "yx":
+        r2 = (vectors[:, [0, 1]]**2).sum(axis=1)
+        dimensions = 2
+    elif axis == "xz" or axis == "zx":
+        r2 = (vectors[:, [0, 2]]**2).sum(axis=1)
+        dimensions = 2
+    elif axis == "yz" or axis == "zy":
+        r2 = (vectors[:, [1, 2]]**2).sum(axis=1)
+        dimensions = 2
     elif axis == "x":
-        r = vectors[:, 0] ** 2
+        r2 = vectors[:, 0] ** 2
         dimensions = 1
     elif axis == "y":
-        r = vectors[:, 1] ** 2
+        r2 = vectors[:, 1] ** 2
         dimensions = 1
     elif axis == "z":
-        r = vectors[:, 2] ** 2
+        r2 = vectors[:, 2] ** 2
         dimensions = 1
     else:
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
+    
+    m2 = np.mean(r2)
+    m4 = np.mean(r2**2)
+    if m2 == 0.0:
+        if full_output:
+            return 0.0, 0.0, 0.0
+        else:
+            return 0.0
+
+    alpha_2 = (m4 / ((1 + 2 / dimensions) * m2**2)) - 1
+    if full_output:
+        return alpha_2, m2, m4
+    else:
+        return alpha_2
+
 
-    return (np.mean(r**2) / ((1 + 2 / dimensions) * (np.mean(r) ** 2))) - 1

src/mdevaluate/distribution.py

@@ -147,9 +147,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
@@ -182,10 +182,10 @@ def tetrahedral_order(
     )
 
     # Connection vectors
-    neighbors_1 -= atoms
-    neighbors_2 -= atoms
-    neighbors_3 -= atoms
-    neighbors_4 -= atoms
+    neighbors_1 = pbc_diff(neighbors_1, atoms, box=atoms.box)
+    neighbors_2 = pbc_diff(neighbors_2, atoms, box=atoms.box)
+    neighbors_3 = pbc_diff(neighbors_3, atoms, box=atoms.box)
+    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)
@@ -306,7 +306,7 @@ def next_neighbor_distribution(
     )[1]
     resname_nn = reference.residue_names[nn]
     count_nn = (resname_nn == atoms.residue_names.reshape(-1, 1)).sum(axis=1)
-    return np.histogram(count_nn, bins=bins, normed=normed)[0]
+    return np.histogram(count_nn, bins=bins, density=normed)[0]
 
 
 def hbonds(

src/mdevaluate/extra/chill.py

@@ -11,7 +11,7 @@ from mdevaluate.coordinates import CoordinateFrame, Coordinates
 from mdevaluate.pbc import pbc_points
 
 
-def a_ij(atoms: ArrayLike, N: int = 4, l: int = 3) -> tuple[NDArray, NDArray]:
+def calc_aij(atoms: ArrayLike, N: int = 4, l: int = 3) -> tuple[NDArray, NDArray]:
     tree = KDTree(atoms)
 
     dist, indices = tree.query(atoms, N + 1)
@@ -84,18 +84,18 @@ def count_ice_types(iceTypes: NDArray) -> NDArray:
 
 
 def selector_ice(
-    start_frame: CoordinateFrame,
-    traj: Coordinates,
+    oxygen_atoms_water: CoordinateFrame,
     chosen_ice_types: ArrayLike,
     combined: bool = True,
+    next_neighbor_distance: float = 0.35,
 ) -> NDArray:
-    oxygen = traj.subset(atom_name="OW")
-    atoms = oxygen[start_frame.step]
-    atoms = atoms % np.diag(atoms.box)
-    atoms_PBC = pbc_points(atoms, thickness=1)
-    aij, indices = a_ij(atoms_PBC)
+    atoms = oxygen_atoms_water
+    atoms_PBC = pbc_points(atoms, thickness=next_neighbor_distance * 2.2)
+    aij, indices = calc_aij(atoms_PBC)
     tree = KDTree(atoms_PBC)
-    neighbors = tree.query_ball_point(atoms_PBC, 0.35, return_length=True)
+    neighbors = tree.query_ball_point(
+        atoms_PBC, next_neighbor_distance, return_length=True
+    ) - 1
     index_SOL = atoms_PBC.tolist().index(atoms[0].tolist())
     index_SOL = np.arange(index_SOL, index_SOL + len(atoms))
     ice_Types = classify_ice(aij, indices, neighbors, index_SOL)
@@ -117,9 +117,9 @@ def selector_ice(
 def ice_types(trajectory: Coordinates, segments: int = 10000) -> pd.DataFrame:
     def ice_types_distribution(frame: CoordinateFrame, selector: Callable) -> NDArray:
         atoms_PBC = pbc_points(frame, thickness=1)
-        aij, indices = a_ij(atoms_PBC)
+        aij, indices = calc_aij(atoms_PBC)
         tree = KDTree(atoms_PBC)
-        neighbors = tree.query_ball_point(atoms_PBC, 0.35, return_length=True)
+        neighbors = tree.query_ball_point(atoms_PBC, 0.35, return_length=True) - 1
         index = selector(frame, atoms_PBC)
         ice_types_data = classify_ice(aij, indices, neighbors, index)
         ice_parts_data = count_ice_types(ice_types_data)
@@ -161,8 +161,8 @@ def ice_types(trajectory: Coordinates, segments: int = 10000) -> pd.DataFrame:
 def ice_clusters_traj(
     traj: Coordinates, segments: int = 10000, skip: float = 0.1
 ) -> list:
-    def ice_clusters(frame: CoordinateFrame, traj: Coordinates) -> Tuple[float, list]:
-        selection = selector_ice(frame, traj, [0, 1, 2])
+    def ice_clusters(frame: CoordinateFrame) -> Tuple[float, list]:
+        selection = selector_ice(frame, [0, 1, 2])
         if len(selection) == 0:
             return frame.time, []
         else:
@@ -194,6 +194,6 @@ def ice_clusters_traj(
         np.int_(np.linspace(len(traj) * skip, len(traj) - 1, num=segments))
     )
     all_clusters = [
-        ice_clusters(traj[frame_index], traj) for frame_index in frame_indices
+        ice_clusters(traj[frame_index]) for frame_index in frame_indices
     ]
     return all_clusters

src/mdevaluate/extra/free_energy_landscape.py

@@ -1,9 +1,9 @@
 from functools import partial
+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
@@ -11,6 +11,56 @@ import multiprocessing as mp
 from ..coordinates import Coordinates
 
 
+def _pbc_points_reduced(
+    coordinates: ArrayLike,
+    pore_geometry: str,
+    box: Optional[NDArray] = None,
+    thickness: Optional[float] = None,
+) -> tuple[NDArray, NDArray]:
+    if box is None:
+        box = coordinates.box
+    if pore_geometry == "cylindrical":
+        grid = np.array([[i, j, k] for k in [-1, 0, 1] for j in [0] for i in [0]])
+        indices = np.tile(np.arange(len(coordinates)), 3)
+    elif pore_geometry == "slit":
+        grid = np.array(
+            [[i, j, k] for k in [0] for j in [1, 0, -1] for i in [-1, 0, 1]]
+        )
+        indices = np.tile(np.arange(len(coordinates)), 9)
+    else:
+        raise ValueError(
+            f"pore_geometry is {pore_geometry}, should either be "
+            f"'cylindrical' or 'slit'"
+        )
+    coordinates_pbc = np.concatenate([coordinates + v @ box for v in grid], axis=0)
+    size = np.diag(box)
+
+    if thickness is not None:
+        mask = np.all(coordinates_pbc > -thickness, axis=1)
+        coordinates_pbc = coordinates_pbc[mask]
+        indices = indices[mask]
+        mask = np.all(coordinates_pbc < size + thickness, axis=1)
+        coordinates_pbc = coordinates_pbc[mask]
+        indices = indices[mask]
+
+    return coordinates_pbc, indices
+
+
+def _build_tree(points, box, r_max, pore_geometry):
+    if np.all(np.diag(np.diag(box)) == box):
+        tree = KDTree(points % box, boxsize=box)
+        points_pbc_index = None
+    else:
+        points_pbc, points_pbc_index = _pbc_points_reduced(
+            points,
+            pore_geometry,
+            box,
+            thickness=r_max + 0.01,
+        )
+        tree = KDTree(points_pbc)
+    return tree, points_pbc_index
+
+
 def occupation_matrix(
     trajectory: Coordinates,
     edge_length: float = 0.05,
@@ -28,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.arange(len(frame_indices))[i : i + size]
-        for i in range(0, len(frame_indices), size)
-    ]
+    indices = np.array_split(frame_indices, nodes)
     pool = mp.Pool(nodes)
     results = pool.map(
         partial(_calc_histogram, trajectory=trajectory, bins=bins), indices
@@ -72,19 +118,19 @@ def _calc_histogram(
 
 
 def find_maxima(
-    occupation_df: pd.DataFrame, box: ArrayLike, edge_length: float = 0.05
+    occupation_df: pd.DataFrame, box: ArrayLike, radius: float, pore_geometry: str
 ) -> pd.DataFrame:
     maxima_df = occupation_df.copy()
     maxima_df["maxima"] = None
     points = np.array(maxima_df[["x", "y", "z"]])
-    tree = KDTree(points, boxsize=box)
-    all_neighbors = tree.query_ball_point(
-        points, edge_length * 3 ** (1 / 2) + edge_length / 100
-    )
+    tree, points_pbc_index = _build_tree(points, box, radius, pore_geometry)
     for i in range(len(maxima_df)):
         if maxima_df.loc[i, "maxima"] is not None:
             continue
-        neighbors = np.array(all_neighbors[i])
+        maxima_pos = maxima_df.loc[i, ["x", "y", "z"]]
+        neighbors = np.array(tree.query_ball_point(maxima_pos, radius))
+        if points_pbc_index is not None:
+            neighbors = points_pbc_index[neighbors]
         neighbors = neighbors[neighbors != i]
         if len(neighbors) == 0:
             maxima_df.loc[i, "maxima"] = True
@@ -104,23 +150,39 @@ def _calc_energies(
     maxima_df: pd.DataFrame,
     bins: ArrayLike,
     box: NDArray,
+    pore_geometry: str,
     T: float,
+    nodes: int = 8,
 ) -> NDArray:
     points = np.array(maxima_df[["x", "y", "z"]])
-    tree = KDTree(points, boxsize=box)
+    tree, points_pbc_index = _build_tree(points, box, bins[-1], pore_geometry)
     maxima = maxima_df.loc[maxima_indices, ["x", "y", "z"]]
     maxima_occupations = np.array(maxima_df.loc[maxima_indices, "occupation"])
     num_of_neighbors = np.max(
         tree.query_ball_point(maxima, bins[-1], return_length=True)
     )
-    distances, indices = tree.query(
-        maxima, k=num_of_neighbors, distance_upper_bound=bins[-1]
-    )
+    split_maxima = []
+    for i in range(0, len(maxima), 1000):
+        split_maxima.append(maxima[i : i + 1000])
+
+    distances = []
+    indices = []
+    for maxima in split_maxima:
+        distances_step, indices_step = tree.query(
+            maxima, k=num_of_neighbors, distance_upper_bound=bins[-1], workers=nodes
+        )
+        distances.append(distances_step)
+        indices.append(indices_step)
+    distances = np.concatenate(distances)
+    indices = np.concatenate(indices)
     all_energy_hist = []
     all_occupied_bins_hist = []
     if distances.ndim == 1:
         current_distances = distances[1:][distances[1:] <= bins[-1]]
-        current_indices = indices[1:][distances[1:] <= bins[-1]]
+        if points_pbc_index is None:
+            current_indices = indices[1:][distances[1:] <= bins[-1]]
+        else:
+            current_indices = points_pbc_index[indices[1:][distances[1:] <= bins[-1]]]
         energy = (
             -np.log(maxima_df.loc[current_indices, "occupation"] / maxima_occupations)
             * T
@@ -131,8 +193,12 @@ def _calc_energies(
         return result
     for i, maxima_occupation in enumerate(maxima_occupations):
         current_distances = distances[i, 1:][distances[i, 1:] <= bins[-1]]
-        current_indices = indices[i, 1:][distances[i, 1:] <= bins[-1]]
-
+        if points_pbc_index is None:
+            current_indices = indices[i, 1:][distances[i, 1:] <= bins[-1]]
+        else:
+            current_indices = points_pbc_index[
+                indices[i, 1:][distances[i, 1:] <= bins[-1]]
+            ]
         energy = (
             -np.log(maxima_df.loc[current_indices, "occupation"] / maxima_occupation)
             * T
@@ -168,9 +234,12 @@ def distance_resolved_energies(
     distance_bins: ArrayLike,
     r_bins: ArrayLike,
     box: NDArray,
+    pore_geometry: str,
     temperature: float,
+    nodes: int = 8,
 ) -> pd.DataFrame:
     results = []
+    distances = []
     for i in range(len(distance_bins) - 1):
         maxima_indices = np.array(
             maxima_df.index[
@@ -179,11 +248,22 @@ def distance_resolved_energies(
                 * (maxima_df["maxima"] == True)
             ]
         )
-        results.append(
-            _calc_energies(maxima_indices, maxima_df, r_bins, box, temperature)
-        )
+        try:
+            results.append(
+                _calc_energies(
+                    maxima_indices,
+                    maxima_df,
+                    r_bins,
+                    box,
+                    pore_geometry,
+                    temperature,
+                    nodes,
+                )
+            )
+            distances.append((distance_bins[i] + distance_bins[i + 1]) / 2)
+        except ValueError:
+            pass
 
-    distances = (distance_bins[:-1] + distance_bins[1:]) / 2
     radii = (r_bins[:-1] + r_bins[1:]) / 2
     d = np.array([d for d in distances for r in radii])
     r = np.array([r for d in distances for r in radii])
@@ -192,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 = []
@@ -201,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)
-        energies.append(np.max(p3(np.linspace(r_min, r_max, 1000))))
+        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})

src/mdevaluate/extra/water.py

@@ -265,7 +265,7 @@ def LSI(
     )
     distributions = np.array(
         [
-            LSI_distribution(trajectory[frame_index], trajectory, bins, selector=None)
+            LSI_distribution(trajectory[frame_index], bins, selector=None)
             for frame_index in frame_indices
         ]
     )

src/mdevaluate/pbc.py

@@ -7,7 +7,7 @@ from numpy.typing import ArrayLike, NDArray
 
 from itertools import product
 
-from .logging import logger
+from .logging_util import logger
 
 if TYPE_CHECKING:
     from mdevaluate.coordinates import CoordinateFrame
@@ -149,32 +149,21 @@ def nojump(frame: CoordinateFrame, usecache: bool = True) -> CoordinateFrame:
             i0 = 0
             delta = 0
 
-        delta = (
-            delta
-            + np.array(
-                np.vstack(
-                    [m[i0 : abstep + 1].sum(axis=0) for m in reader.nojump_matrices]
-                ).T
-            )
-            * frame.box.diagonal()
-        )
+        delta = (delta
+            + np.vstack(
+                [m[i0 : abstep + 1].sum(axis=0) for m in reader.nojump_matrices]
+            ).T)
 
         reader._nojump_cache[abstep] = delta
         while len(reader._nojump_cache) > NOJUMP_CACHESIZE:
             reader._nojump_cache.popitem(last=False)
-        delta = delta[selection, :]
     else:
-        delta = (
-            np.array(
-                np.vstack(
-                    [
-                        m[: frame.step + 1, selection].sum(axis=0)
-                        for m in reader.nojump_matrices
-                    ]
+        delta = np.vstack(
+                [m[: frame.step + 1, selection].sum(axis=0) for m in reader.nojump_matrices]
                 ).T
-            )
-            * frame.box.diagonal()
-        )
+    
+    delta = delta[selection, :]
+    delta = np.array(delta @ frame.box)
     return frame - delta
 
 

src/mdevaluate/reader.py

@@ -19,13 +19,13 @@ import MDAnalysis
 from scipy import sparse
 
 from .checksum import checksum
-from .logging import logger
+from .logging_util import logger
 from . import atoms
 from .coordinates import Coordinates
 
 CSR_ATTRS = ("data", "indices", "indptr")
 NOJUMP_MAGIC = 2016
-Group_RE = re.compile("\[ ([-+\w]+) \]")
+Group_RE = re.compile(r"\[ ([-+\w]+) \]")
 
 
 class NojumpError(Exception):
@@ -152,7 +152,7 @@ def nojump_load_filename(reader: BaseReader):
         )
         if os.path.exists(full_path_fallback):
             return full_path_fallback
-    if os.path.exists(fname) or is_writeable(directory):
+    if os.path.exists(full_path) or is_writeable(directory):
         return full_path
     else:
         user_data_dir = os.path.join("/data/", os.environ["HOME"].split("/")[-1])
@@ -187,19 +187,33 @@ def nojump_save_filename(reader: BaseReader):
 
 def parse_jumps(trajectory: Coordinates):
     prev = trajectory[0].whole
-    box = prev.box.diagonal()
+    box = prev.box
     SparseData = namedtuple("SparseData", ["data", "row", "col"])
     jump_data = (
         SparseData(data=array("b"), row=array("l"), col=array("l")),
         SparseData(data=array("b"), row=array("l"), col=array("l")),
         SparseData(data=array("b"), row=array("l"), col=array("l")),
     )
-
     for i, curr in enumerate(trajectory):
         if i % 500 == 0:
             logger.debug("Parse jumps Step: %d", i)
-        delta = ((curr - prev) / box).round().astype(np.int8)
+        r3 = np.subtract(curr, prev)
+        delta_z = np.array(np.rint(np.divide(r3[:, 2], box[2][2])), dtype=np.int8)
+        r2 = np.subtract(
+            r3,
+            (np.rint(np.divide(r3[:, 2], box[2][2])))[:, np.newaxis]
+            * box[2][np.newaxis, :],
+            )
+        delta_y = np.array(np.rint(np.divide(r2[:, 1], box[1][1])), dtype=np.int8)
+        r1 = np.subtract(
+            r2,
+            (np.rint(np.divide(r2[:, 1], box[1][1])))[:, np.newaxis]
+            * box[1][np.newaxis, :],
+            )
+        delta_x = np.array(np.rint(np.divide(r1[:, 0], box[0][0])), dtype=np.int8)
+        delta = np.array([delta_x, delta_y, delta_z]).T
         prev = curr
+        box = prev.box
         for d in range(3):
             (col,) = np.where(delta[:, d] != 0)
             jump_data[d].col.extend(col)
@@ -261,7 +275,7 @@ def load_nojump_matrices(reader: BaseReader):
                 "Loaded Nojump matrices: {}".format(nojump_load_filename(reader))
             )
         else:
-            logger.info("Invlaid Nojump Data: {}".format(nojump_load_filename(reader)))
+            logger.info("Invalid Nojump Data: {}".format(nojump_load_filename(reader)))
     except KeyError:
         logger.info("Removing zip-File: %s", zipname)
         os.remove(nojump_load_filename(reader))

src/mdevaluate/utils.py

@@ -14,7 +14,7 @@ from scipy.ndimage import uniform_filter1d
 from scipy.interpolate import interp1d
 from scipy.optimize import curve_fit
 
-from .logging import logger
+from .logging_util import logger
 from .functions import kww, kww_1e
 
 
@@ -357,6 +357,37 @@ def quick1etau(t: ArrayLike, C: ArrayLike, n: int = 7) -> float:
     return tau_est
 
 
+def quicknongaussfit(t, C, width=2):
+    """
+    Estimates the time and height of the peak in the non-Gaussian function.
+    C is C(t) the correlation function
+    """
+    def ffunc(t,y0,A_main,log_tau_main,sig_main):
+        main_peak = A_main*np.exp(-(t - log_tau_main)**2 / (2 * sig_main**2))
+        return y0 + main_peak
+    
+    # first rough estimate, the closest time. This is returned if the interpolation fails!
+    tau_est = t[np.argmax(C)]
+    nG_max = np.amax(C)
+    try:
+        with np.errstate(invalid='ignore'):
+            corr = C[t > 0]
+            time = np.log10(t[t > 0])
+            tau = time[np.argmax(corr)]
+            mask = (time>tau-width/2) & (time<tau+width/2)
+            time = time[mask] ; corr = corr[mask]
+            nG_min = C[t > 0].min()
+            guess = [nG_min, nG_max-nG_min, tau, 0.6]
+            popt = curve_fit(ffunc, time, corr, p0=guess, maxfev=10000)[0]
+            tau_est = 10**popt[-2]
+            nG_max = popt[0] + popt[1]
+    except:
+        pass
+    if np.isnan(tau_est):
+        tau_est = np.inf
+    return tau_est, nG_max
+
+
 def susceptibility(
     time: NDArray, correlation: NDArray, **kwargs
 ) -> tuple[NDArray, NDArray]:

test/test_free_energy_landscape.py

@@ -13,42 +13,24 @@ def trajectory(request):
 
 
 def test_get_fel(trajectory):
-    test_array = np.array(
-        [
-            174.46253634,
-            174.60905476,
-            178.57658092,
-            182.43001192,
-            180.57916378,
-            176.49886217,
-            178.96018547,
-            181.13561782,
-            178.31026314,
-            176.08903996,
-            180.71215345,
-            181.59703135,
-            180.34329368,
-            187.02474488,
-            197.99167477,
-            214.05788031,
-            245.58571282,
-            287.52457507,
-            331.53492965,
-        ]
-    )
+    test_array = np.array([210., 214., 209., 192., 200., 193., 230., 218., 266.])
 
     OW = trajectory.subset(atom_name="OW")
-
-    box = np.diag(trajectory[0].box)
-    box_voxels = (box // [0.05, 0.05, 0.05] + [1, 1, 1]) * [0.05, 0.05, 0.05]
-    occupation_matrix = fel.occupation_matrix(OW, skip=0, segments=1000)
-    maxima_matrix = fel.find_maxima(occupation_matrix, box=box_voxels, edge_length=0.05)
-    maxima_matrix = fel.add_distances(maxima_matrix, "cylindrical", box / 2)
-    r_bins = np.arange(0, 2, 0.02)
-    distance_bins = np.arange(0.05, 2.05, 0.1)
+    box = trajectory[0].box
+    box_voxels = (np.diag(box) // [0.05, 0.05, 0.05] + [1, 1, 1]) * [0.05, 0.05, 0.05]
+    occupation_matrix = fel.occupation_matrix(OW, skip=0, segments=10)
+    radius_maxima = 0.05 * 3 ** (1 / 2) + 0.05 / 100
+    maxima_matrix = fel.find_maxima(
+        occupation_matrix,
+        box=box_voxels,
+        radius=radius_maxima,
+        pore_geometry="cylindrical",
+    )
+    maxima_matrix = fel.add_distances(maxima_matrix, "cylindrical", np.diag(box) / 2)
+    r_bins = np.arange(0, 0.5, 0.02)
+    distance_bins = np.arange(1.8, 1.9, 0.01)
     energy_df = fel.distance_resolved_energies(
-        maxima_matrix, distance_bins, r_bins, box, 225
+        maxima_matrix, distance_bins, r_bins, box, "cylindrical", 225
     )
     result = fel.find_energy_maxima(energy_df, r_min=0.05, r_max=0.15)
-
     assert (np.round(np.array(result["energy"])) == np.round(test_array)).all()