parse jumps now correctly uses fractional coordinates and should work for triclinic boxes. this should lead to correct behavior in all cases except for maybe changing angles in triclinic boxes. old behavior can still be used

src/mdevaluate/autosave.py

@@ -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

@@ -102,8 +102,20 @@ def checksum(*args, csum=None, _seen=None):
         _seen.add(obj_id)
 
         if hasattr(arg, "__checksum__"):
-            logger.debug("Checksum via __checksum__: %s", str(arg))
-            csum.update(str(arg.__checksum__()).encode())
+            method = getattr(arg, "__checksum__")
+            if callable(method) and not isinstance(arg, type):
+                logger.debug("Checksum via __checksum__: %s", str(arg))
+                csum.update(str(method()).encode())
+            elif isinstance(arg, type):
+                try:
+                    src = inspect.getsource(arg)
+                    csum.update(strip_comments(src).encode())
+                    logger.debug("Checksum via class source for %s", arg.__name__)
+                except (OSError, TypeError):
+                    csum.update(arg.__name__.encode())
+                    logger.debug("Checksum via class name for %s", arg.__name__)
+            else:
+                logger.debug("Skipping unbound __checksum__ on %s", type(arg))
         elif isinstance(arg, bytes):
             csum.update(arg)
         elif isinstance(arg, str):

src/mdevaluate/coordinates.py

@@ -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

src/mdevaluate/correlation.py

@@ -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(
         [

src/mdevaluate/pbc.py

@@ -163,7 +163,7 @@ def nojump(frame: CoordinateFrame, usecache: bool = True) -> CoordinateFrame:
                 ).T
     
     delta = delta[selection, :]
-    delta = delta @ frame.box
+    delta = np.array(delta @ frame.box)
     return frame - delta
 
 

src/mdevaluate/reader.py

@@ -23,6 +23,8 @@ from .logging_util import logger
 from . import atoms
 from .coordinates import Coordinates
 
+from unittest.mock import MagicMock
+
 CSR_ATTRS = ("data", "indices", "indptr")
 NOJUMP_MAGIC = 2016
 Group_RE = re.compile(r"\[ ([-+\w]+) \]")
@@ -185,35 +187,52 @@ def nojump_save_filename(reader: BaseReader):
         return full_path_fallback
 
 
-def parse_jumps(trajectory: Coordinates):
-    prev = trajectory[0].whole
+def parse_jumps(trajectory: Coordinates, whole: bool=True, fractional_inverted: bool=True):
+    if whole:
+        prev = trajectory[0].whole
+    else:
+        prev = trajectory[0]
     box = prev.box
+    if fractional_inverted:
+        s_prev = prev @ np.linalg.inv(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)
-        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
+        if not fractional_inverted:
+            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
+        else:
+            s_curr = curr @ np.linalg.inv(curr.box)
+
+            ds = s_curr - s_prev
+            delta = np.array(np.rint(ds), dtype=np.int8)
+
+            s_prev = s_curr
+
+
         for d in range(3):
             (col,) = np.where(delta[:, d] != 0)
             jump_data[d].col.extend(col)
@@ -240,7 +259,18 @@ def generate_nojump_matrices(trajectory: Coordinates):
     save_nojump_matrices(trajectory.frames)
 
 
+def _ensure_xdr(reader: BaseReader):
+    """Patch missing _xdr attribute for non-XDR readers (e.g. LAMMPS DumpReader)
+    with a stable mock so checksums are consistent across runs."""
+    if not hasattr(reader.rd, '_xdr'):
+        mock_xdr = MagicMock()
+        mock_xdr.offsets = np.arange(len(reader))
+        print(f"Adding mock _xdr attribute for to reader of length {len(reader)}.")
+        reader.rd._xdr = mock_xdr
+
+
 def save_nojump_matrices(reader: BaseReader, matrices: npt.ArrayLike = None):
+    _ensure_xdr(reader)
     if matrices is None:
         matrices = reader.nojump_matrices
     data = {"checksum": checksum(NOJUMP_MAGIC, checksum(reader))}
@@ -253,6 +283,7 @@ def save_nojump_matrices(reader: BaseReader, matrices: npt.ArrayLike = None):
 
 
 def load_nojump_matrices(reader: BaseReader):
+    _ensure_xdr(reader)
     zipname = nojump_load_filename(reader)
     try:
         data = np.load(zipname, allow_pickle=True)

src/mdevaluate/utils.py

@@ -334,6 +334,11 @@ def quick1etau(t: ArrayLike, C: ArrayLike, n: int = 7) -> float:
     C is C(t) the correlation function
     n is the minimum number of points around 1/e required
     """
+    # norm, if t=0 provided
+    if t[0] == 0:
+        C /= C[0]
+    C, t = C[t>0], t[t>0] # make sure t=0 is dropped
+
     # first rough estimate, the closest time. This is returned if the interpolation fails!
     tau_est = t[np.argmin(np.fabs(C - np.exp(-1)))]
     # reduce the data to points around 1/e