Show install with tag

Reviewed-on: #6

README.md

@@ -26,9 +26,6 @@ time, msd = md.correlation.shifted_correlation(
 
 ## Installation
 
-The package requires the Python package [pygmx](https://github.com/mdevaluate/pygmx),
-which handles reading of Gromacs file formats.
-Installation of pygmx is described in its own repository.
 
 The mdevaluate package itself is plain Python code and, hence, can be imported from its directory directly, 
 or may be installed via setuptools to the local Python environment by running
@@ -36,6 +33,14 @@ or may be installed via setuptools to the local Python environment by running
 	python setup.py install
 
 
+When you are using `uv` you can install it with:
+
+    uv venv some_folder
+    source some_folder/bin/activate
+    uv pip install git+https://gitea.pkm.physik.tu-darmstadt.de/IPKM/mdevaluate.git@2026.09
+
+Note: you can append a tag to get a specific release.
+
 ## Running the tests
 
 Mdevaluate includes a test suite that can be used to check if the installation was succesful.

create_mdevaluate_module.sh

@@ -0,0 +1,71 @@
+#!/bin/bash
+
+CONDA_VERSION=2024.10
+PYTHON_VERSION=3.12
+
+if [ -z "$1" ]; then
+    echo "No argument supplied, version to create expected"
+    exit 1
+fi
+
+
+if [ ! -w "/nfsopt/mdevaluate"]; then
+	echo "Please remount /nfsopt writable"
+	exit 2
+fi
+
+MD_VERSION=$1
+
+# purge evtl. loaded modules
+module purge
+
+
+echo "Create mdevaluate Python environemnt using conda"
+echo "Using conda version: $CONDA_VERSION"
+echo "Using Python version: $PYTHON_VERSION"
+
+module load anaconda3/$CONDA_VERSION
+conda create -y --prefix /nfsopt/mdevaluate/mdevaluate-${MD_VERSION} \
+             python=$PYTHON_VERSION
+module purge
+
+echo "Create modulefile for mdevaluate/$MD_VERSION"
+cat > /nfsopt/modulefiles/mdevaluate/$MD_VERSION <<EOF
+#%Module1.0#####################################################################
+##
+## dot modulefile
+##
+## modulefiles/dot.  Generated from dot.in by configure.
+##
+
+module-whatis	"Enables the mdevaluate Python environment."
+
+set version	${MD_VERSION}
+set module_path /nfsopt/mdevaluate/mdevaluate-\$version/bin
+
+prepend-path PATH \$module_path
+
+EOF
+
+echo "Loading mdevaluate environment and install packages"
+module load mdevaluate/${MD_VERSION}
+pip install jupyter \
+spyder \
+mdanalysis \
+pathos \
+pandas \
+dask \
+sqlalchemy \
+psycopg2-binary \
+trimesh \
+pyvista \
+seaborn \
+black \
+black[jupyter] \
+tables \
+pyedr \
+pytest
+
+pip install git+https://gitea.pkm.physik.tu-darmstadt.de/IPKM/mdevaluate.git
+pip install git+https://gitea.pkm.physik.tu-darmstadt.de/IPKM/python-store.git
+pip install git+https://gitea.pkm.physik.tu-darmstadt.de/IPKM/python-tudplot.git

pyproject.toml

@@ -4,12 +4,15 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "mdevaluate"
-version = "24.02"
+version = "2026.09"
+requires-python = ">=3.10"
 dependencies = [
     "mdanalysis",
     "pandas",
     "dask",
     "pathos",
     "tables",
-    "pyedr"
+    "pyedr",
+    "numpy<2.3",
+    "scipy<1.17",
 ]

src/mdevaluate/__init__.py

@@ -16,7 +16,7 @@ from . import reader
 from . import system
 from . import utils
 from . import extra
-from .logging_util import logger
+from .logging 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_util import logger
+from .logging import logger
 
 autosave_directory: Optional[str] = None
 load_autosave_data = False
@@ -166,10 +166,8 @@ 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_dir is not None:
+            if autosave_directory 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,14 +1,9 @@
 import functools
 import hashlib
-from .logging_util import logger
+from .logging 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
 
@@ -33,46 +28,19 @@ def version(version_nr: int, calls: Iterable = ()):
     return decorator
 
 
-def strip_comments(source: str) -> str:
+def strip_comments(s: str):
-    """Removes docstrings, comments, and irrelevant whitespace from Python source code."""
+    """Strips comment lines and docstring from Python source string."""
-
+    o = ""
-    # Step 1: Remove docstrings using AST
+    in_docstring = False
-    def remove_docstrings(node):
+    for l in s.split("\n"):
-        if isinstance(node, (ast.FunctionDef, ast.AsyncFunctionDef, ast.ClassDef, ast.Module)):
+        if l.strip().startswith(("#", '"', "'")) or in_docstring:
-            if (doc := ast.get_docstring(node, clean=False)):
+            in_docstring = l.strip().startswith(('"""', "'''")) + in_docstring == 1
-                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
-        if srow > last_lineno:
+        o += l + "\n"
-            last_col = 0
+    return o
-        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, _seen=None):
+def checksum(*args, csum=None):
     """
     Calculate a checksum of any object, by sha1 hash.
 
@@ -92,30 +60,10 @@ def checksum(*args, csum=None, _seen=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__"):
-            method = getattr(arg, "__checksum__")
+            logger.debug("Checksum via __checksum__: %s", str(arg))
-            if callable(method) and not isinstance(arg, type):
+            csum.update(str(arg.__checksum__()).encode())
-                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):
@@ -129,15 +77,15 @@ def checksum(*args, csum=None, _seen=None):
             for key in sorted(merged):  # deterministic ordering
                 v = merged[key]
                 if v is not arg:
-                    checksum(v, csum=csum, _seen=_seen)
+                    checksum(v, csum=csum)
         elif isinstance(arg, functools.partial):
             logger.debug("Checksum via partial for %s", str(arg))
-            checksum(arg.func, csum=csum, _seen=_seen)
+            checksum(arg.func, csum=csum)
             for x in arg.args:
-                checksum(x, csum=csum, _seen=_seen)
+                checksum(x, csum=csum)
             for k in sorted(arg.keywords.keys()):
                 csum.update(k.encode())
-                checksum(arg.keywords[k], csum=csum, _seen=_seen)
+                checksum(arg.keywords[k], csum=csum)
         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_util import logger
+from .logging import logger
 from typing import Optional, Callable, List, Tuple
 
 import numpy as np
@@ -436,34 +436,6 @@ def center_of_masses(
     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
 def pore_coordinates(
     frame: CoordinateFrame, origin: ArrayLike, sym_axis: str = "z"

src/mdevaluate/correlation.py

@@ -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(nargs=2, kwargs_keys=('selector', 'segments', 'skip', 'window', 'average', 'points',), version=1.0)
+@autosave_data(2)
 def shifted_correlation(
     function: Callable,
     frames: Coordinates,
@@ -147,8 +147,7 @@ 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)
-    dt = round(frames[1].time - frames[0].time, 6) # round to avoid bad floats
+    t = np.array([frames[i].time for i in idx]) - frames[0].time
-    t = idx * dt
 
     result = np.array(
         [
@@ -200,7 +199,7 @@ def msd(
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
 
 
-def isf_raw(
+def isf(
     start_frame: CoordinateFrame,
     end_frame: CoordinateFrame,
     q: float = 22.7,
@@ -217,59 +216,29 @@ def isf_raw(
         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)
+        return np.sinc(distance * q / np.pi).mean()
     elif axis == "xy" or axis == "yx":
         distance = (displacements[:, [0, 1]]**2).sum(axis=1) ** 0.5
-        return np.real(jn(0, distance * q))
+        return np.real(jn(0, distance * q)).mean()
     elif axis == "xz" or axis == "zx":
         distance = (displacements[:, [0, 2]]**2).sum(axis=1) ** 0.5
-        return np.real(jn(0, distance * q))
+        return np.real(jn(0, distance * q)).mean()
     elif axis == "yz" or axis == "zy":
         distance = (displacements[:, [1, 2]]**2).sum(axis=1) ** 0.5
-        return np.real(jn(0, distance * q))
+        return np.real(jn(0, distance * q)).mean()
     elif axis == "x":
         distance = np.abs(displacements[:, 0])
-        return np.cos(np.abs(q * distance))
+        return np.mean(np.cos(np.abs(q * distance)))
     elif axis == "y":
         distance = np.abs(displacements[:, 1])
-        return np.cos(np.abs(q * distance))
+        return np.mean(np.cos(np.abs(q * distance)))
     elif axis == "z":
         distance = np.abs(displacements[:, 2])
-        return np.cos(np.abs(q * distance))
+        return np.mean(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:
@@ -461,11 +430,10 @@ 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
     """
@@ -474,41 +442,27 @@ def non_gaussian_parameter(
     else:
         vectors = displacements_without_drift(start_frame, end_frame, trajectory)
     if axis == "all":
-        r2 = (vectors**2).sum(axis=1)
+        r = (vectors**2).sum(axis=1)
         dimensions = 3
     elif axis == "xy" or axis == "yx":
-        r2 = (vectors[:, [0, 1]]**2).sum(axis=1)
+        r = (vectors[:, [0, 1]]**2).sum(axis=1)
         dimensions = 2
     elif axis == "xz" or axis == "zx":
-        r2 = (vectors[:, [0, 2]]**2).sum(axis=1)
+        r = (vectors[:, [0, 2]]**2).sum(axis=1)
         dimensions = 2
     elif axis == "yz" or axis == "zy":
-        r2 = (vectors[:, [1, 2]]**2).sum(axis=1)
+        r = (vectors[:, [1, 2]]**2).sum(axis=1)
         dimensions = 2
     elif axis == "x":
-        r2 = vectors[:, 0] ** 2
+        r = vectors[:, 0] ** 2
         dimensions = 1
     elif axis == "y":
-        r2 = vectors[:, 1] ** 2
+        r = vectors[:, 1] ** 2
         dimensions = 1
     elif axis == "z":
-        r2 = vectors[:, 2] ** 2
+        r = vectors[:, 2] ** 2
         dimensions = 1
     else:
         raise ValueError('Parameter axis has to be ether "all", "x", "y", or "z"!')
 
-    m2 = np.mean(r2)
+    return (np.mean(r**2) / ((1 + 2 / dimensions) * (np.mean(r) ** 2))) - 1
-    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
-
-

src/mdevaluate/pbc.py

@@ -7,7 +7,7 @@ from numpy.typing import ArrayLike, NDArray
 
 from itertools import product
 
-from .logging_util import logger
+from .logging import logger
 
 if TYPE_CHECKING:
     from mdevaluate.coordinates import CoordinateFrame

src/mdevaluate/reader.py

@@ -19,15 +19,13 @@ import MDAnalysis
 from scipy import sparse
 
 from .checksum import checksum
-from .logging_util import logger
+from .logging 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]+) \]")
+Group_RE = re.compile("\[ ([-+\w]+) \]")
 
 
 class NojumpError(Exception):
@@ -187,52 +185,35 @@ def nojump_save_filename(reader: BaseReader):
         return full_path_fallback
 
 
-def parse_jumps(trajectory: Coordinates, whole: bool=True, fractional_inverted: bool=True):
+def parse_jumps(trajectory: Coordinates):
-    if whole:
+    prev = trajectory[0].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)
-        if not fractional_inverted:
+        r3 = np.subtract(curr, prev)
-            r3 = np.subtract(curr, prev)
+        delta_z = np.array(np.rint(np.divide(r3[:, 2], box[2][2])), dtype=np.int8)
-            delta_z = np.array(np.rint(np.divide(r3[:, 2], box[2][2])), dtype=np.int8)
+        r2 = np.subtract(
-            r2 = np.subtract(
+            r3,
-                r3,
+            (np.rint(np.divide(r3[:, 2], box[2][2])))[:, np.newaxis]
-                (np.rint(np.divide(r3[:, 2], box[2][2])))[:, np.newaxis]
+            * box[2][np.newaxis, :],
-                * box[2][np.newaxis, :],
+            )
-                )
+        delta_y = np.array(np.rint(np.divide(r2[:, 1], box[1][1])), dtype=np.int8)
-            delta_y = np.array(np.rint(np.divide(r2[:, 1], box[1][1])), dtype=np.int8)
+        r1 = np.subtract(
-            r1 = np.subtract(
+            r2,
-                r2,
+            (np.rint(np.divide(r2[:, 1], box[1][1])))[:, np.newaxis]
-                (np.rint(np.divide(r2[:, 1], box[1][1])))[:, np.newaxis]
+            * box[1][np.newaxis, :],
-                * box[1][np.newaxis, :],
+            )
-                )
+        delta_x = np.array(np.rint(np.divide(r1[:, 0], box[0][0])), dtype=np.int8)
-            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
-            delta = np.array([delta_x, delta_y, delta_z]).T
+        prev = curr
-            prev = curr
+        box = prev.box
-            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)
@@ -259,18 +240,7 @@ 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))}
@@ -283,7 +253,6 @@ 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)
@@ -306,7 +275,7 @@ def load_nojump_matrices(reader: BaseReader):
                 "Loaded Nojump matrices: {}".format(nojump_load_filename(reader))
             )
         else:
-            logger.info("Invalid Nojump Data: {}".format(nojump_load_filename(reader)))
+            logger.info("Invlaid 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_util import logger
+from .logging import logger
 from .functions import kww, kww_1e
 
 
@@ -334,11 +334,6 @@ 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
@@ -362,37 +357,6 @@ 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]: