diff --git a/benchmarks/bench_checks.py b/benchmarks/bench_checks.py
index c75759a0..6d88d002 100644
--- a/benchmarks/bench_checks.py
+++ b/benchmarks/bench_checks.py
@@ -1,7 +1,7 @@
 import eko
-import numpy as np
 import pineappl
 import pytest
+from eko.io import manipulate
 
 import pineko.check
 
@@ -16,12 +16,27 @@ def benchmark_check_grid_and_eko_compatible(test_files, tmp_path):
     with eko.EKO.edit(
         test_files / "data/ekos/400/HERA_NC_225GEV_EP_SIGMARED.tar"
     ) as ekoop:
-        with pytest.raises(ValueError):
-            pineko.check.check_grid_and_eko_compatible(wrong_grid, ekoop, 1.0, 3, 3)
-        pineko.check.check_grid_and_eko_compatible(grid, ekoop, 1.0, 3, 3)
-        eko.io.manipulate.xgrid_reshape(
-            ekoop, targetgrid=eko.interpolation.XGrid([0.0001, 0.001, 0.1, 0.5, 1.0])
-        )
-        with pytest.raises(ValueError):
-            pineko.check.check_grid_and_eko_compatible(grid, ekoop, 1.0, 10, 10)
-        eko.io.manipulate.xgrid_reshape(ekoop, targetgrid=ekoop.xgrid)
+        eko_xgrid = ekoop.xgrid.tolist()
+        for (eko_mu2, _), _ in ekoop.items():
+            with pytest.raises(ValueError):
+                pineko.check.check_grid_and_eko_compatible(
+                    wrong_grid, eko_xgrid, eko_mu2, 1.0, 3, 3
+                )
+            pineko.check.check_grid_and_eko_compatible(
+                grid, eko_xgrid, eko_mu2, 1.0, 3, 3
+            )
+
+        # test a wrong rotation
+        wrong_xgrid = eko.interpolation.XGrid([0.0001, 0.001, 0.1, 0.5, 1.0])
+        for (eko_mu2, _), op in ekoop.items():
+            # TODO: here we can only check inputgrid as this eko has dimension (14,40,14,50)
+            # and ekoop.xgrid has 50
+            op = manipulate.xgrid_reshape(op, ekoop.xgrid, 4, inputgrid=wrong_xgrid)
+            assert op.operator.shape[-1] == len(wrong_xgrid)
+            with pytest.raises(ValueError):
+                pineko.check.check_grid_and_eko_compatible(
+                    grid, wrong_xgrid.tolist(), eko_mu2, 1.0, 10, 10
+                )
+            # restore xgrid
+            op = manipulate.xgrid_reshape(op, wrong_xgrid, 4, inputgrid=ekoop.xgrid)
+            assert op.operator.shape[-1] == len(ekoop.xgrid)
diff --git a/benchmarks/bench_cli.py b/benchmarks/bench_cli.py
index 54b3ddec..3944fe76 100644
--- a/benchmarks/bench_cli.py
+++ b/benchmarks/bench_cli.py
@@ -1,7 +1,6 @@
 import pathlib
 import shutil
 
-import lhapdf
 from click.testing import CliRunner
 
 from pineko.cli._base import command
@@ -19,7 +18,7 @@ def benchmark_check_cli(test_files):
     result = runner.invoke(
         command, ["check", "compatibility", str(grid_path), str(eko_path)]
     )
-    assert "Success: grids are compatible" in result.output
+    assert "Success: grids and eko are compatible" in result.output
     wrong_result = runner.invoke(
         command, ["check", "compatibility", str(wrong_grid_path), str(eko_path)]
     )
diff --git a/benchmarks/bench_evolve.py b/benchmarks/bench_evolve.py
index 79bdfee0..6404fc5f 100644
--- a/benchmarks/bench_evolve.py
+++ b/benchmarks/bench_evolve.py
@@ -2,7 +2,6 @@
 import pathlib
 
 import eko
-import eko.io.legacy
 import numpy as np
 import pineappl
 import pytest
diff --git a/benchmarks/bench_regression.py b/benchmarks/bench_regression.py
index b6e2bf7f..57d8dc9c 100644
--- a/benchmarks/bench_regression.py
+++ b/benchmarks/bench_regression.py
@@ -72,9 +72,10 @@ def xfxQ2(self, pid, x, q2):
 
 def _trim_template(template_card, take_points=10):
     """Trim the template card so that the number of x-values to compute is much smaller"""
-    card_info = OperatorCard.from_dict(
-        safe_load(template_card.read_text(encoding="utf-8"))
-    )
+    raw_card = safe_load(template_card.read_text(encoding="utf-8"))
+    raw_card["init"] = (raw_card["mu0"], 4)
+    del raw_card["mu0"]
+    card_info = OperatorCard.from_dict(raw_card)
     original_x = card_info.xgrid
     size = len(original_x.raw)
     skip = int(size / take_points)
diff --git a/docs/source/overview/prerequisites.rst b/docs/source/overview/prerequisites.rst
index 05802e1c..f1ab340b 100644
--- a/docs/source/overview/prerequisites.rst
+++ b/docs/source/overview/prerequisites.rst
@@ -94,7 +94,9 @@ An example is the following::
     polarized: False
     time_like: False
     n_integration_cores: 1
-  mu0: 1.65
+  init: 
+    - 1.65
+    - 4
   mugrid:
     - - 50.0
       - 5
diff --git a/pyproject.toml b/pyproject.toml
index 90c9f08d..91b958eb 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -27,8 +27,8 @@ packages = [{ include = "pineko", from = "src" }]
 
 [tool.poetry.dependencies]
 python = ">=3.9,<3.13"
-eko = "^0.14.2"
-pineappl = "^0.8.2"
+eko = "^0.15.0"
+pineappl = "^0.8.6"
 PyYAML = "^6.0"
 numpy = "^1.21.0"
 pandas = "^2.1"
diff --git a/src/pineko/check.py b/src/pineko/check.py
index 3aebfc39..a6f6b96a 100644
--- a/src/pineko/check.py
+++ b/src/pineko/check.py
@@ -76,15 +76,19 @@ def in1d(a, b, rtol=1e-05, atol=1e-08):
     )
 
 
-def check_grid_and_eko_compatible(pineappl_grid, operators, xif, max_as, max_al):
-    """Check whether the EKO operators and the PineAPPL grid are compatible.
+def check_grid_and_eko_compatible(
+    pineappl_grid, eko_xgrid, eko_mu2, xif, max_as, max_al
+):
+    """Check whether the EKO operator and the PineAPPL grid are compatible.
 
     Parameters
     ----------
     pineappl_grid : pineappl.grid.Grid
         grid
-    operators : eko.EKO
-        operators
+    eko_xgrid : list
+        eko interpolation xgrid
+    eko_mu2: float
+        eko mu2 scale
     xif : float
         factorization scale variation
     max_as: int
@@ -95,7 +99,7 @@ def check_grid_and_eko_compatible(pineappl_grid, operators, xif, max_as, max_al)
     Raises
     ------
     ValueError
-        If the operators and the grid are not compatible.
+        If the operator and the grid are not compatible.
     """
     order_mask = pineappl.grid.Order.create_mask(
         pineappl_grid.orders(), max_as, max_al, True
@@ -104,16 +108,12 @@ def check_grid_and_eko_compatible(pineappl_grid, operators, xif, max_as, max_al)
     x_grid = evol_info.x1
     muf2_grid = evol_info.fac1
     # Q2 grid
-    if not np.all(
-        in1d(np.unique(list(operators.mu2grid)), xif * xif * np.array(muf2_grid))
-    ):
+    if not np.all(in1d(np.array([eko_mu2]), xif * xif * np.array(muf2_grid))):
         raise ValueError(
             "Q2 grid in pineappl grid and eko operator are NOT compatible!"
         )
     # x-grid
-    if not np.all(
-        in1d(np.unique(operators.bases.targetgrid.tolist()), np.array(x_grid))
-    ):
+    if not np.all(in1d(np.unique(eko_xgrid), np.array(x_grid))):
         raise ValueError("x grid in pineappl grid and eko operator are NOT compatible!")
 
 
diff --git a/src/pineko/cli/check.py b/src/pineko/cli/check.py
index a515c017..32a903cd 100644
--- a/src/pineko/cli/check.py
+++ b/src/pineko/cli/check.py
@@ -37,14 +37,16 @@ def sub_compatibility(grid_path, operator_path, xif, max_as, max_al):
 
     """
     pineappl_grid = pineappl.grid.Grid.read(grid_path)
+    pineappl_grid.optimize()
     with eko.EKO.read(operator_path) as operators:
-        try:
-            check.check_grid_and_eko_compatible(
-                pineappl_grid, operators, xif, max_as, max_al
-            )
-            rich.print("[green]Success:[/] grids are compatible")
-        except ValueError as e:
-            rich.print("[red]Error:[/]", e)
+        for (q2, _), _ in operators.items():
+            try:
+                check.check_grid_and_eko_compatible(
+                    pineappl_grid, operators.xgrid.tolist(), q2, xif, max_as, max_al
+                )
+            except ValueError as e:
+                rich.print("[red]Error:[/]", e)
+    rich.print("[green]Success:[/] grids and eko are compatible.")
 
 
 @dataclass
diff --git a/src/pineko/cli/convolve.py b/src/pineko/cli/convolve.py
index 6244ae35..e31921d4 100644
--- a/src/pineko/cli/convolve.py
+++ b/src/pineko/cli/convolve.py
@@ -67,6 +67,7 @@ def subcommand(
     PDF is an optional PDF set compatible with the EKO to compare grid and FK table.
     """
     grid = pineappl.grid.Grid.read(grid_path)
+    grid.optimize()
     n_ekos = len(op_paths)
     with eko.EKO.edit(op_paths[0]) as operators1:
         rich.print(
diff --git a/src/pineko/evolve.py b/src/pineko/evolve.py
index ca6c1784..20c1b406 100644
--- a/src/pineko/evolve.py
+++ b/src/pineko/evolve.py
@@ -16,6 +16,8 @@
 import rich.panel
 import yaml
 from eko import basis_rotation
+from eko.interpolation import XGrid
+from eko.io import manipulate
 from eko.io.types import ScaleVariationsMethod
 from eko.matchings import Atlas, nf_default
 from eko.quantities import heavy_quarks
@@ -124,9 +126,11 @@ def write_operator_card_from_file(
     if not pathlib.Path(pineappl_path).exists():
         raise FileNotFoundError(pineappl_path)
     pineappl_grid = pineappl.grid.Grid.read(pineappl_path)
+    pineappl_grid.optimize()
     default_card = yaml.safe_load(
         pathlib.Path(default_card_path).read_text(encoding="utf-8")
     )
+
     return write_operator_card(pineappl_grid, default_card, card_path, tcard)
 
 
@@ -201,10 +205,12 @@ def write_operator_card(pineappl_grid, default_card, card_path, tcard):
     # update scale variation method
     operators_card["configs"]["scvar_method"] = sv_method
 
-    # Make sure that we are using the theory Q0 and fail if the template has a different one
-    operators_card["mu0"] = tcard["Q0"]
-    if default_card.get("mu0") is not None and default_card["mu0"] != tcard["Q0"]:
-        raise ValueError("Template declares a value of Q0 different from theory")
+    operators_card["init"] = (tcard["Q0"], tcard["nf0"])
+    if default_card.get("init") is not None and default_card["init"] != [
+        tcard["Q0"],
+        tcard["nf0"],
+    ]:
+        raise ValueError("Template declares a value of Q0, nf0 different from theory")
 
     q2_grid = (xif * xif * muf2_grid).tolist()
     masses = np.array([tcard["mc"], tcard["mb"], tcard["mt"]]) ** 2
@@ -353,40 +359,19 @@ def evolve_grid(
 
     evol_info = grid.evolve_info(order_mask)
     x_grid = evol_info.x1
+    if "integrability_version" in grid.key_values():
+        x_grid = np.append(x_grid, 1.0)
+
     mur2_grid = evol_info.ren1
     xif = 1.0 if operators1.operator_card.configs.scvar_method is not None else xif
     tcard = operators1.theory_card
     opcard = operators1.operator_card
-    # rotate the targetgrid
-    if "integrability_version" in grid.key_values():
-        x_grid = np.append(x_grid, 1.0)
-
-    def xgrid_reshape(full_operator):
-        """Reinterpolate operators on output and/or input grids."""
-        eko.io.manipulate.xgrid_reshape(
-            full_operator, targetgrid=eko.interpolation.XGrid(x_grid)
-        )
-        check.check_grid_and_eko_compatible(grid, full_operator, xif, max_as, max_al)
-        # rotate to evolution (if doable and necessary)
-        if np.allclose(full_operator.bases.inputpids, br.flavor_basis_pids):
-            eko.io.manipulate.to_evol(full_operator)
-        # Here we are checking if the EKO contains the rotation matrix (flavor to evol)
-        elif not np.allclose(
-            full_operator.bases.inputpids, br.rotate_flavor_to_evolution
-        ):
-            raise ValueError("The EKO is neither in flavor nor in evolution basis.")
-
-    xgrid_reshape(operators1)
-    if operators2 is not None:
-        xgrid_reshape(operators2)
 
     # PineAPPL wants alpha_s = 4*pi*a_s
     # remember that we already accounted for xif in the opcard generation
     evmod = eko.couplings.couplings_mod_ev(opcard.configs.evolution_method)
     # Couplings ask for the square of the masses
     thresholds_ratios = np.power(tcard.heavy.matching_ratios, 2.0)
-    for q in range(tcard.couplings.max_num_flavs + 1, 6 + 1):
-        thresholds_ratios[q - 4] = np.inf
     sc = eko.couplings.Couplings(
         tcard.couplings,
         tcard.order,
@@ -398,25 +383,31 @@ def xgrid_reshape(full_operator):
     # To compute the alphas values we are first reverting the factorization scale shift
     # and then obtaining the renormalization scale using xir.
     ren_grid2 = xir * xir * mur2_grid
+    nfgrid = [x[1] for x in operators1.operator_card.mugrid]
     alphas_values = [
-        4.0
-        * np.pi
-        * sc.a_s(
-            mur2,
-        )
-        for mur2 in ren_grid2
+        4.0 * np.pi * sc.a_s(mur2, nf_to=nf) for mur2, nf in zip(ren_grid2, nfgrid)
     ]
 
     def prepare(operator):
         """Match the raw operator with its relevant metadata."""
         for (q2, _), op in operator.items():
+            # reshape the x-grid output
+            op = manipulate.xgrid_reshape(
+                op,
+                operator.xgrid,
+                opcard.configs.interpolation_polynomial_degree,
+                targetgrid=XGrid(x_grid),
+            )
+            # rotate the input to evolution basis
+            op = manipulate.to_evol(op, source=True)
+            check.check_grid_and_eko_compatible(grid, x_grid, q2, xif, max_as, max_al)
             info = PyOperatorSliceInfo(
                 fac0=operator.mu20,
-                x0=operator.bases.inputgrid.raw,
+                x0=operator.xgrid.tolist(),
                 pids0=basis_rotation.evol_basis_pids,
                 fac1=q2,
-                x1=operator.bases.targetgrid.raw,
-                pids1=operator.bases.targetpids,
+                x1=x_grid.tolist(),
+                pids1=basis_rotation.flavor_basis_pids,
                 pid_basis=PyPidBasis.Evol,
             )
             yield (info, op.operator)
diff --git a/src/pineko/theory.py b/src/pineko/theory.py
index fa58d169..a9580728 100644
--- a/src/pineko/theory.py
+++ b/src/pineko/theory.py
@@ -11,7 +11,6 @@
 import time
 
 import eko
-import eko.io.legacy
 import numpy as np
 import pineappl
 import rich
diff --git a/tests/test_evolve.py b/tests/test_evolve.py
index f529b80b..8c2fc653 100644
--- a/tests/test_evolve.py
+++ b/tests/test_evolve.py
@@ -52,20 +52,20 @@ def test_write_operator_card_q0(tmp_path):
     g = FakePine()
     t = copy.deepcopy(default_card)
     o = copy.deepcopy(example.raw_operator())
-    # 1. Same Q0 and mu0, all ok
+    # 1. Same Q0 and init, all ok
     t["Q0"] = 5.0
-    o["mu0"] = 5.0
+    o["init"] = [5.0, 4]
     _xs, _mu2s = pineko.evolve.write_operator_card(g, o, p, t)
     with open(p, encoding="utf8") as f:
         oo = yaml.safe_load(f)
-    np.testing.assert_allclose(oo["mu0"], t["Q0"])
+    np.testing.assert_allclose(oo["init"][0], t["Q0"])
     # 2. Q0 only in theory, all ok
-    o.pop("mu0")
+    o.pop("init")
     _xs, _mu2s = pineko.evolve.write_operator_card(g, o, p, t)
     with open(p, encoding="utf8") as f:
         oo = yaml.safe_load(f)
-    np.testing.assert_allclose(oo["mu0"], t["Q0"])
+    np.testing.assert_allclose(oo["init"][0], t["Q0"])
     # 3. op is different, raises error
-    o["mu0"] = 11.0
+    o["init"] = [11.0, 3]
     with pytest.raises(ValueError):
         _xs, _mu2s = pineko.evolve.write_operator_card(g, o, p, t)