diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py
index abe659730..3526c56b4 100644
--- a/benchmarks/eko/benchmark_evol_to_unity.py
+++ b/benchmarks/eko/benchmark_evol_to_unity.py
@@ -1,4 +1,3 @@
-import pathlib
 from dataclasses import dataclass
 
 import numpy as np
@@ -47,7 +46,6 @@ def test_operator_grid(
         self,
         theory_card: TheoryCard,
         operator_card: OperatorCard,
-        tmp_path: pathlib.Path,
     ):
         """Test that eko_forward @ eko_backward gives ID matrix or zeros."""
         update_cards(theory_card, operator_card)
diff --git a/doc/source/code/Operators.rst b/doc/source/code/Operators.rst
index c58a64b37..b7c6a7250 100644
--- a/doc/source/code/Operators.rst
+++ b/doc/source/code/Operators.rst
@@ -19,10 +19,7 @@ The classes are nested as follows:
         PhysicalOperator [label="PhysicalOperator"];
         Operator [label="Operator" ];
         OME [label="OME" ];
-        OperatorGrid [label="OperatorGrid"];
 
-        OperatorGrid -> Operator;
-        OperatorGrid -> OME;
         Operator -> PhysicalOperator [weight=100,style=dashed];
         PhysicalOperator -> ndarray [style=dashed];
         OME -> MatchingCondition [weight=100,style=dashed];
@@ -31,17 +28,8 @@ The classes are nested as follows:
         OpMember -> PhysicalOperator [dir=back];
         OME -> OpMember;
         OpMember -> MatchingCondition [dir=back];
-
-        OperatorGrid -> OpMember -> ndarray [style=invis];
     }
 
-- :class:`~eko.evolution_operator.grid.OperatorGrid`
-
-    * is the master class which administrates all operator tasks
-    * is instantiated once for each run
-    * holds all necessary :doc:`configurations </code/IO>`
-    * holds all necessary instances of the :doc:`/code/Utilities`
-
 - :class:`~eko.evolution_operator.Operator`
 
     * represents a configuration for a fixed final scale :math:`Q_1^2`
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index bd1b19d65..366f934c9 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -7,6 +7,7 @@
 import logging
 import os
 import time
+from dataclasses import dataclass
 from multiprocessing import Pool
 from typing import Dict, Tuple
 
@@ -21,6 +22,8 @@
 from .. import basis_rotation as br
 from .. import interpolation, mellin
 from .. import scale_variations as sv
+from ..couplings import Couplings
+from ..interpolation import InterpolatorDispatcher
 from ..io.types import EvolutionMethod, OperatorLabel
 from ..kernels import ev_method
 from ..kernels import non_singlet as ns
@@ -28,7 +31,7 @@
 from ..kernels import singlet as s
 from ..kernels import singlet_qed as qed_s
 from ..kernels import valence_qed as qed_v
-from ..matchings import Segment, lepton_number
+from ..matchings import Atlas, Segment, lepton_number
 from ..member import OpMember
 from ..scale_variations import expanded as sv_expanded
 from ..scale_variations import exponentiated as sv_exponentiated
@@ -608,6 +611,15 @@ def quad_ker_qed(
 """Map of all operators."""
 
 
+@dataclass(frozen=True)
+class Managers:
+    """Set of steering objects."""
+
+    atlas: Atlas
+    couplings: Couplings
+    interpolator: InterpolatorDispatcher
+
+
 class Operator(sv.ScaleVariationModeMixin):
     """Internal representation of a single EKO.
 
@@ -637,7 +649,12 @@ class Operator(sv.ScaleVariationModeMixin):
     full_labels_qed: Tuple[OperatorLabel, ...] = br.full_unified_labels
 
     def __init__(
-        self, config, managers, segment: Segment, mellin_cut=5e-2, is_threshold=False
+        self,
+        config,
+        managers: Managers,
+        segment: Segment,
+        mellin_cut=5e-2,
+        is_threshold=False,
     ):
         self.config = config
         self.managers = managers
@@ -669,7 +686,7 @@ def xif2(self):
         return self.config["xif2"]
 
     @property
-    def int_disp(self):
+    def int_disp(self) -> InterpolatorDispatcher:
         """Return the interpolation dispatcher."""
         return self.managers.interpolator
 
diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py
deleted file mode 100644
index 9279e1f66..000000000
--- a/src/eko/evolution_operator/grid.py
+++ /dev/null
@@ -1,204 +0,0 @@
-"""Define operators container and computing workflow.
-
-The first is the driver class of eko as it is the one that collects all
-the previously instantiated information and does the actual computation
-of the Q2s.
-"""
-
-import logging
-from dataclasses import dataclass
-from typing import Any, Dict, List, Optional
-
-import numpy as np
-import numpy.typing as npt
-
-from .. import member
-from .. import scale_variations as sv
-from ..couplings import Couplings
-from ..interpolation import InterpolatorDispatcher
-from ..io.runcards import Configs, Debug
-from ..io.types import EvolutionPoint as EPoint
-from ..io.types import Order, SquaredScale
-from ..matchings import Atlas, Segment, flavor_shift, is_downward_path
-from . import Operator, OpMembers, flavors, matching_condition, physical
-from .operator_matrix_element import OperatorMatrixElement
-
-logger = logging.getLogger(__name__)
-
-OpDict = Dict[str, Optional[npt.NDArray]]
-"""In particular, only the ``operator`` and ``error`` fields are expected."""
-
-
-@dataclass(frozen=True)
-class Managers:
-    """Set of steering objects."""
-
-    atlas: Atlas
-    couplings: Couplings
-    interpolator: InterpolatorDispatcher
-
-
-class OperatorGrid(sv.ScaleVariationModeMixin):
-    """Collection of evolution operators for several scales.
-
-    The operator grid is the driver class of the evolution.
-
-    It receives as input a threshold holder and a generator of a_s.
-    From that point onwards it can compute any operator at any q2.
-
-    Attributes
-    ----------
-    config: dict
-    q2_grid: np.ndarray
-    managers: dict
-    """
-
-    def __init__(
-        self,
-        mu2grid: List[EPoint],
-        order: Order,
-        masses: List[float],
-        mass_scheme,
-        thresholds_ratios: List[float],
-        xif: float,
-        n3lo_ad_variation: tuple,
-        matching_order: Order,
-        configs: Configs,
-        debug: Debug,
-        atlas: Atlas,
-        couplings: Couplings,
-        interpol_dispatcher: InterpolatorDispatcher,
-        use_fhmruvv: bool,
-    ):
-        # check
-        config: Dict[str, Any] = {}
-        config["order"] = order
-        config["xif2"] = xif**2
-        config["HQ"] = mass_scheme
-        config["ModSV"] = configs.scvar_method
-        config["n3lo_ad_variation"] = n3lo_ad_variation
-        config["use_fhmruvv"] = use_fhmruvv
-
-        for i, q in enumerate("cbt"):
-            config[f"m{q}"] = masses[i]
-        config["thresholds_ratios"] = thresholds_ratios
-        method = config["method"] = configs.evolution_method.value
-        config["backward_inversion"] = configs.inversion_method
-        config["ev_op_max_order"] = configs.ev_op_max_order
-        config["ev_op_iterations"] = configs.ev_op_iterations
-        config["n_integration_cores"] = configs.n_integration_cores
-        config["debug_skip_singlet"] = debug.skip_singlet
-        config["debug_skip_non_singlet"] = debug.skip_non_singlet
-        config["polarized"] = configs.polarized
-        config["time_like"] = configs.time_like
-        config["matching_order"] = matching_order
-
-        if order == (1, 0) and method != "iterate-exact":
-            logger.warning("Evolution: In LO we use the exact solution always!")
-
-        logger.info(dict(polarized=configs.polarized))
-        logger.info(dict(time_like=configs.time_like))
-
-        self.config = config
-        self.q2_grid = mu2grid
-        self.managers = Managers(
-            atlas=atlas,
-            couplings=couplings,
-            interpolator=interpol_dispatcher,
-        )
-        self._threshold_operators: Dict[Segment, Operator] = {}
-        self._matching_operators: Dict[SquaredScale, OpMembers] = {}
-
-    def get_threshold_operators(self, path: List[Segment]) -> List[Operator]:
-        """Generate the threshold operators.
-
-        This method is called everytime the OperatorGrid is asked for a
-        grid on Q^2 with a list of the relevant areas. If new threshold
-        operators need to be computed, they will be cached in an
-        internal dictionary.
-
-        The internal dictionary is self._threshold_operators and its
-        structure is: (q2_from, q2_to) -> eko.operators.Operator
-
-        It computes and stores the necessary macthing operators.
-        """
-        # The base area is always that of the reference q
-        thr_ops = []
-        # is_downward point to smaller nf
-        is_downward = is_downward_path(path)
-        shift = flavor_shift(is_downward)
-        for seg in path[:-1]:
-            kthr = self.config["thresholds_ratios"][seg.nf - shift]
-            ome = OperatorMatrixElement(
-                self.config,
-                self.managers,
-                seg.nf - shift + 3,
-                seg.target,
-                is_downward,
-                np.log(kthr),
-                self.config["HQ"] == "MSBAR",
-            )
-            if seg not in self._threshold_operators:
-                # Compute the operator and store it
-                logger.info("Prepare threshold operator")
-                op_th = Operator(self.config, self.managers, seg, is_threshold=True)
-                op_th.compute()
-                self._threshold_operators[seg] = op_th
-            thr_ops.append(self._threshold_operators[seg])
-
-            # Compute the matching conditions and store it
-            if seg.target not in self._matching_operators:
-                ome.compute()
-                self._matching_operators[seg.target] = ome.op_members
-        return thr_ops
-
-    def compute(self) -> Dict[EPoint, dict]:
-        """Compute all ekos for the `q2grid`."""
-        return {q2: self.generate(q2) for q2 in self.q2_grid}
-
-    def generate(self, q2: EPoint) -> OpDict:
-        r"""Compute a single EKO.
-
-        eko :math:`\mathbf E(Q^2 \leftarrow Q_0^2)` in flavor basis as
-        numpy array.
-        """
-        # The lists of areas as produced by the thresholds
-        path = self.managers.atlas.path(q2)
-        # Prepare the path for the composition of the operator
-        thr_ops = self.get_threshold_operators(path)
-        # we start composing with the highest operator ...
-        operator = Operator(self.config, self.managers, path[-1])
-        operator.compute()
-
-        is_downward = is_downward_path(path)
-        qed = self.config["order"][1] > 0
-
-        final_op = physical.PhysicalOperator.ad_to_evol_map(
-            operator.op_members, operator.nf, operator.q2_to, qed
-        )
-        # and multiply the lower ones from the right
-        for op in reversed(list(thr_ops)):
-            phys_op = physical.PhysicalOperator.ad_to_evol_map(
-                op.op_members, op.nf, op.q2_to, qed
-            )
-
-            # join with the basis rotation, since matching requires c+ (or likewise)
-            nf_match = op.nf - 1 if is_downward else op.nf
-            matching = matching_condition.MatchingCondition.split_ad_to_evol_map(
-                self._matching_operators[op.q2_to],
-                nf_match,
-                op.q2_to,
-                qed=qed,
-            )
-            if is_downward:
-                invrot = member.ScalarOperator.promote_names(
-                    flavors.rotate_matching_inverse(op.nf, qed), op.q2_to
-                )
-                final_op = final_op @ matching @ invrot @ phys_op
-            else:
-                rot = member.ScalarOperator.promote_names(
-                    flavors.rotate_matching(op.nf + 1, qed), op.q2_to
-                )
-                final_op = final_op @ rot @ matching @ phys_op
-        values, errors = final_op.to_flavor_basis_tensor(qed)
-        return {"operator": values, "error": errors}
diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py
index 098b52db2..e7d51df9f 100644
--- a/src/eko/evolution_operator/operator_matrix_element.py
+++ b/src/eko/evolution_operator/operator_matrix_element.py
@@ -5,6 +5,7 @@
 import enum
 import functools
 import logging
+from typing import Optional
 
 import numba as nb
 import numpy as np
@@ -18,7 +19,7 @@
 from ..io.types import InversionMethod
 from ..matchings import Segment
 from ..scale_variations.exponentiated import gamma_variation
-from . import Operator, QuadKerBase
+from . import Managers, Operator, QuadKerBase
 
 logger = logging.getLogger(__name__)
 
@@ -31,7 +32,7 @@ class MatchingMethods(enum.IntEnum):
     BACKWARD_EXPANDED = enum.auto()
 
 
-def matching_method(s: InversionMethod) -> MatchingMethods:
+def matching_method(s: Optional[InversionMethod]) -> MatchingMethods:
     """Return the matching method.
 
     Parameters
@@ -241,7 +242,16 @@ class OperatorMatrixElement(Operator):
     # still valid in QED since Sdelta and Vdelta matchings are diagonal
     full_labels_qed = copy.deepcopy(full_labels)
 
-    def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar):
+    def __init__(
+        self,
+        config,
+        managers: Managers,
+        nf: int,
+        q2: float,
+        is_backward: bool,
+        L: float,
+        is_msbar: bool,
+    ):
         super().__init__(config, managers, Segment(q2, q2, nf))
         self.backward_method = matching_method(
             config["backward_inversion"] if is_backward else None
diff --git a/src/eko/runner/parts.py b/src/eko/runner/parts.py
index 002caddde..96c4d9ace 100644
--- a/src/eko/runner/parts.py
+++ b/src/eko/runner/parts.py
@@ -15,14 +15,13 @@
 from .. import evolution_operator as evop
 from ..evolution_operator import matching_condition, physical
 from ..evolution_operator import operator_matrix_element as ome
-from ..evolution_operator.grid import Managers
 from ..io import EKO
 from ..io.items import Evolution, Matching, Operator
 from ..quantities.heavy_quarks import QuarkMassScheme
 from . import commons
 
 
-def _managers(eko: EKO) -> Managers:
+def _managers(eko: EKO) -> evop.Managers:
     """Collect managers for operator computation.
 
     .. todo::
@@ -31,7 +30,7 @@ def _managers(eko: EKO) -> Managers:
     """
     tcard = eko.theory_card
     ocard = eko.operator_card
-    return Managers(
+    return evop.Managers(
         atlas=commons.atlas(tcard, ocard),
         couplings=commons.couplings(tcard, ocard),
         interpolator=commons.interpolator(ocard),
diff --git a/tests/eko/evolution_operator/test_grid.py b/tests/eko/evolution_operator/test_grid.py
deleted file mode 100644
index f06c152fc..000000000
--- a/tests/eko/evolution_operator/test_grid.py
+++ /dev/null
@@ -1,55 +0,0 @@
-"""Checks that the operator grid works as intended.
-
-These test can be slow as they require the computation of several values
-of Q But they should be fast as the grid is very small. It does *not*
-test whether the result is correct, it can just test that it is sane
-"""
-
-
-# from eko.runner import legacy
-
-# def test_compute_mu2grid(theory_ffns, operator_card, tmp_path):
-#     mugrid = [(10.0, 5), (100.0, 5)]
-#     operator_card.mugrid = mugrid
-#     opgrid = legacy.Runner(
-#         theory_ffns(3), operator_card, path=tmp_path / "eko.tar"
-#     ).op_grid
-#     opg = opgrid.compute()
-#     assert len(opg) == len(mugrid)
-#     assert all(k in op for k in ["operator", "error"] for op in opg.values())
-
-
-# def test_grid_computation_VFNS(theory_card, operator_card, tmp_path):
-#     """Checks that the grid can be computed."""
-#     mugrid = [(3, 4), (5, 5), (5, 4)]
-#     operator_card.mugrid = mugrid
-#     opgrid = legacy.Runner(
-#         theory_card, operator_card, path=tmp_path / "eko.tar"
-#     ).op_grid
-#     operators = opgrid.compute()
-#     assert len(operators) == len(mugrid)
-
-
-# def test_mod_expanded(theory_card, theory_ffns, operator_card, tmp_path: pathlib.Path):
-#     mugrid = [(3, 4)]
-#     operator_card.mugrid = mugrid
-#     operator_card.configs.scvar_method = eko.io.types.ScaleVariationsMethod.EXPANDED
-#     epsilon = 1e-1
-#     path = tmp_path / "eko.tar"
-#     for is_ffns, nf0 in zip([False, True], [5, 3]):
-#         if is_ffns:
-#             theory = theory_ffns(nf0)
-#         else:
-#             theory = theory_card
-#         theory.order = (1, 0)
-#         operator_card.init = (operator_card.init[0], nf0)
-#         path.unlink(missing_ok=True)
-#         opgrid = legacy.Runner(theory, operator_card, path=path).op_grid
-#         opg = opgrid.compute()
-#         theory.xif = 1.0 + epsilon
-#         path.unlink(missing_ok=True)
-#         sv_opgrid = legacy.Runner(theory, operator_card, path=path).op_grid
-#         sv_opg = sv_opgrid.compute()
-#         np.testing.assert_allclose(
-#             opg[(9, 4)]["operator"], sv_opg[(9, 4)]["operator"], atol=2.5 * epsilon
-#         )