diff --git a/examples/appl2pine/appl2pine.cpp b/examples/appl2pine/appl2pine.cpp
index 6266b7e8..702e3016 100644
--- a/examples/appl2pine/appl2pine.cpp
+++ b/examples/appl2pine/appl2pine.cpp
@@ -444,9 +444,9 @@ pineappl_grid* convert_grid(appl::grid& grid, bool reweight)
     auto const& results = grid.vconvolute(evolvepdf, evolvepdf, alphaspdf, 1);
     std::vector<double> other_results(results.size());
 
-    pineappl_grid_convolute(
+    pineappl_grid_convolute_with_one(
         grids.at(0),
-        xfx,
+        2212,
         xfx,
         alphas,
         pdf.get(),
diff --git a/examples/cpp/dyaa.cpp b/examples/cpp/dyaa.cpp
index 7597f4e3..2089d942 100644
--- a/examples/cpp/dyaa.cpp
+++ b/examples/cpp/dyaa.cpp
@@ -137,7 +137,7 @@ int main() {
     };
 
     std::vector<double> dxsec(24);
-    pineappl_grid_convolute(grid, xfx, xfx, alphas, pdf, nullptr,
+    pineappl_grid_convolute_with_one(grid, 2212, xfx, alphas, pdf, nullptr,
         nullptr, 1.0, 1.0, dxsec.data());
 
     // print the results
diff --git a/examples/fnlo2pine/fnlo2pine.cpp b/examples/fnlo2pine/fnlo2pine.cpp
index 70dd1802..12d3601a 100644
--- a/examples/fnlo2pine/fnlo2pine.cpp
+++ b/examples/fnlo2pine/fnlo2pine.cpp
@@ -550,8 +550,6 @@ int main(int argc, char* argv[])
 
     std::vector<pineappl_grid*> grids;
 
-    std::size_t pdfs = 0;
-
     for (auto const id : ids)
     {
         auto coeff_table = file.GetCoeffTable(id);
@@ -561,7 +559,6 @@ int main(int argc, char* argv[])
 
         if (converted != nullptr)
         {
-            pdfs = converted->GetNPDF();
             grids.push_back(convert_coeff_add_fix(converted, file, bins, alpha));
         }
         else
@@ -570,7 +567,6 @@ int main(int argc, char* argv[])
 
             if (converted != nullptr)
             {
-                pdfs = converted->GetNPDF();
                 auto const mur_ff = file.GetMuRFunctionalForm();
                 auto const muf_ff = file.GetMuFFunctionalForm();
 
@@ -612,10 +608,10 @@ int main(int argc, char* argv[])
     auto const& results = file.GetCrossSection();
     std::vector<double> other_results(results.size());
 
-    pineappl_grid_convolute(
+    pineappl_grid_convolute_with_one(
         grids.at(0),
+        2212,
         xfx,
-        pdfs == 2 ? xfx : xfx2,
         alphas,
         pdf.get(),
         nullptr,
diff --git a/examples/python/dyaa.py b/examples/python/dyaa.py
index 3eb82a7b..6811add4 100755
--- a/examples/python/dyaa.py
+++ b/examples/python/dyaa.py
@@ -157,8 +157,9 @@ def main(calls, pdfname, filename):
         import lhapdf  # pylint: disable=import-outside-toplevel
 
         pdf = lhapdf.mkPDF(pdfname, 0)
+        pdg_id = int(pdf.set().get_entry('Particle'))
         # perform convolution
-        dxsec = grid.convolute(pdf.xfxQ2, pdf.xfxQ2, pdf.alphasQ2)
+        dxsec = grid.convolute_with_one(pdg_id, pdf.xfxQ2, pdf.alphasQ2)
         for i in range(len(dxsec)):
             print(f"{bins[i]:.1f} {bins[i + 1]:.1f} {dxsec[i]:.3e}")
 
diff --git a/pineappl/Cargo.toml b/pineappl/Cargo.toml
index 93130fef..831a3b6c 100644
--- a/pineappl/Cargo.toml
+++ b/pineappl/Cargo.toml
@@ -11,14 +11,13 @@ categories = ["science"]
 description = "PineAPPL is not an extension of APPLgrid"
 
 [dependencies]
-anyhow = "1.0"
 arrayvec = "0.5"
 bincode = "1.0"
 enum_dispatch = "0.3.5"
 float-cmp = "0.8"
 git-version = "0.3"
 itertools = "0.10"
-lz-fear = "0.1"
+lz4_flex = "0.8.2"
 ndarray = { features = ["serde"], version = "0.15.3" }
 rustc-hash = "1.1.0"
 serde = { features = ["derive"], version = "1.0" }
@@ -26,6 +25,7 @@ thiserror = "1.0"
 indicatif = "0.15.0"
 
 [dev-dependencies]
+anyhow = "1.0"
 lhapdf = "0.1.8"
 rand = { default-features = false, version = "0.7" }
 rand_pcg = { default-features = false, version = "0.2" }
diff --git a/pineappl/src/fk_table.rs b/pineappl/src/fk_table.rs
index bb721a08..6291d980 100644
--- a/pineappl/src/fk_table.rs
+++ b/pineappl/src/fk_table.rs
@@ -1,6 +1,7 @@
 //! Provides the [`FkTable`] type.
 
-use super::grid::{Grid, Order};
+use super::grid::{Grid, GridError, Order};
+use super::lumi::LumiCache;
 use super::subgrid::Subgrid;
 use ndarray::Array4;
 use std::convert::TryFrom;
@@ -139,22 +140,28 @@ impl FkTable {
     /// # Errors
     ///
     /// TODO
-    pub fn write(&self, writer: impl Write) -> anyhow::Result<()> {
+    pub fn write(&self, writer: impl Write) -> Result<(), GridError> {
         self.grid.write(writer)
     }
 
+    /// Propagate `write_lz4` to `Grid`.
+    ///
+    /// # Errors
+    ///
+    /// See [`Grid::write_lz4`].
+    pub fn write_lz4(&self, writer: impl Write) -> Result<(), GridError> {
+        self.grid.write_lz4(writer)
+    }
+
     /// Propagate convolute to grid
     pub fn convolute(
         &self,
-        xfx1: &dyn Fn(i32, f64, f64) -> f64,
-        xfx2: &dyn Fn(i32, f64, f64) -> f64,
-        order_mask: &[bool],
+        lumi_cache: &mut LumiCache,
         bin_indices: &[usize],
         lumi_mask: &[bool],
-        xi: &[(f64, f64)],
     ) -> Vec<f64> {
         self.grid
-            .convolute(xfx1, xfx2, &|_| 1.0, order_mask, bin_indices, lumi_mask, xi)
+            .convolute(lumi_cache, &[], bin_indices, lumi_mask, &[(1.0, 1.0)])
     }
 }
 
diff --git a/pineappl/src/grid.rs b/pineappl/src/grid.rs
index 665e5107..11e83d54 100644
--- a/pineappl/src/grid.rs
+++ b/pineappl/src/grid.rs
@@ -15,7 +15,7 @@ use float_cmp::approx_eq;
 use git_version::git_version;
 use indicatif::{ProgressBar, ProgressStyle};
 use itertools::Itertools;
-use lz_fear::{framed::DecompressionError::WrongMagic, LZ4FrameReader};
+use lz4_flex::frame::{self, FrameDecoder, FrameEncoder};
 use ndarray::{s, Array3, Array5, Dimension};
 use rustc_hash::FxHashMap;
 use serde::{Deserialize, Serialize};
@@ -23,7 +23,7 @@ use std::borrow::Cow;
 use std::cmp::Ordering;
 use std::collections::HashMap;
 use std::convert::{TryFrom, TryInto};
-use std::io::{Read, Seek, SeekFrom, Write};
+use std::io::{self, BufReader, BufWriter, Read, Seek, Write};
 use std::mem;
 use std::ops::Range;
 use std::ptr;
@@ -119,6 +119,9 @@ pub enum GridError {
     /// Returned when failed to write a Grid.
     #[error(transparent)]
     WriteFailure(bincode::Error),
+    /// Returned while performing IO operations.
+    #[error(transparent)]
+    IoFailure(io::Error),
 }
 
 #[derive(Clone, Deserialize, Serialize)]
@@ -336,177 +339,17 @@ impl Grid {
         Cow::Borrowed(self.lumi())
     }
 
-    /// Performs a convolution of the contained subgrids with the given PDFs, `xfx1` for the first
-    /// parton and `xfx2` for the second parton, `alphas` for the evaluation of the strong
-    /// coupling. The parameters `order_mask` and `lumi_mask` can be used to selectively enable
-    /// perturbative orders and luminosities; they must either be empty (everything enabled) or as
-    /// large as the orders and luminosity function, respectively. If the corresponding entry is
-    /// `true` the order/luminosity is enable, `false` disables the entry. The tuple `xi` can be
-    /// used to independently vary the renormalization (first element) and factorization scale
-    /// (second element) from their central value `(1.0, 1.0)`.
+    /// Perform a convolution using the PDFs and strong coupling in `lumi_cache`, and only
+    /// selecting only the orders, bins and luminosities corresponding to `order_mask`,
+    /// `bin_indices` and `lumi_mask`. A variation of the scales
+    /// is performed using the factors in `xi`; the first factor varies the renormalization scale,
+    /// the second the factorization scale. Note that for the variation to be trusted all non-zero
+    /// log-grids must be contained.
     ///
     /// # Panics
     ///
     /// TODO
     pub fn convolute(
-        &self,
-        xfx1: &dyn Fn(i32, f64, f64) -> f64,
-        xfx2: &dyn Fn(i32, f64, f64) -> f64,
-        alphas: &dyn Fn(f64) -> f64,
-        order_mask: &[bool],
-        bin_indices: &[usize],
-        lumi_mask: &[bool],
-        xi: &[(f64, f64)],
-    ) -> Vec<f64> {
-        let bin_indices = if bin_indices.is_empty() {
-            (0..self.bin_limits.bins()).collect()
-        } else {
-            bin_indices.to_vec()
-        };
-        let mut bins: Vec<f64> = vec![0.0; bin_indices.len() * xi.len()];
-        let bin_sizes = self.bin_info().normalizations();
-
-        // prepare pdf and alpha caches
-        let mut pdf_cache1 = FxHashMap::default();
-        let mut pdf_cache2 = FxHashMap::default();
-        let mut alphas_cache = FxHashMap::default();
-        let mut last_xif = 0.0;
-
-        let (mut mu2_grid, mut x1_grid, mut x2_grid) = self
-            .subgrids
-            .iter()
-            .find(|subgrid| !subgrid.is_empty())
-            .map_or_else(
-                || (Cow::default(), Cow::default(), Cow::default()),
-                |grid| (grid.mu2_grid(), grid.x1_grid(), grid.x2_grid()),
-            );
-        let use_cache = !mu2_grid.is_empty() && !x1_grid.is_empty() && !x2_grid.is_empty();
-        let two_caches = !ptr::eq(&xfx1, &xfx2);
-
-        let mut xir_values: Vec<_> = xi.iter().map(|xi| xi.0).collect();
-        xir_values.sort_by(|lhs, rhs| lhs.partial_cmp(rhs).unwrap());
-        xir_values.dedup();
-        let xir_indices: Vec<_> = xi
-            .iter()
-            .map(|xi| xir_values.iter().position(|xir| xi.0 == *xir).unwrap())
-            .collect();
-        let self_lumi = self.pdg_lumi();
-
-        // iterate over the elements of `xi` and a corresponding index, but sorted using the
-        // factorisation value of `xi`
-        for (l, (&(xir, xif), &xir_index)) in xi
-            .iter()
-            .zip(xir_indices.iter())
-            .enumerate()
-            .sorted_by(|lhs, rhs| (lhs.1).1.partial_cmp((rhs.1).1).unwrap())
-        {
-            // whenever the value `xif` changes we can clear the PDF cache
-            if xif != last_xif {
-                pdf_cache1.clear();
-                pdf_cache2.clear();
-                last_xif = xif;
-            }
-            // iterate subgrids
-            for ((i, j, k), subgrid) in self.subgrids.indexed_iter() {
-                let order = &self.orders[i];
-                // can we skip?
-                if ((order.logxir > 0) && (xir == 1.0)) || ((order.logxif > 0) && (xif == 1.0)) {
-                    continue;
-                }
-
-                if (!order_mask.is_empty() && !order_mask[i])
-                    || (!lumi_mask.is_empty() && !lumi_mask[k])
-                {
-                    continue;
-                }
-
-                let bin_index = match bin_indices.iter().position(|&index| index == j) {
-                    Some(i) => i,
-                    None => continue,
-                };
-
-                let lumi_entry = &self_lumi[k];
-
-                let mut value = if subgrid.is_empty() {
-                    0.0
-                } else if use_cache {
-                    // Now we actually have to do something
-                    let new_mu2_grid = subgrid.mu2_grid();
-                    let new_x1_grid = subgrid.x1_grid();
-                    let new_x2_grid = subgrid.x2_grid();
-                    let mu2_grid_changed = new_mu2_grid != mu2_grid;
-
-                    if mu2_grid_changed {
-                        mu2_grid = new_mu2_grid;
-                        alphas_cache.clear();
-                    }
-
-                    if mu2_grid_changed || (new_x1_grid != x1_grid) || (new_x2_grid != x2_grid) {
-                        x1_grid = new_x1_grid;
-                        x2_grid = new_x2_grid;
-                        pdf_cache1.clear();
-                        pdf_cache2.clear();
-                    }
-
-                    // pass convolute down, using caching
-                    subgrid.convolute(&x1_grid, &x2_grid, &mu2_grid, &mut |ix1, ix2, imu2| {
-                        let x1 = x1_grid[ix1];
-                        let x2 = x2_grid[ix2];
-                        let mu2 = &mu2_grid[imu2];
-                        let muf2 = xif * xif * mu2.fac;
-
-                        let mut lumi = 0.0;
-
-                        for entry in lumi_entry.entry() {
-                            let xfx1 = *pdf_cache1
-                                .entry((entry.0, ix1, imu2))
-                                .or_insert_with(|| xfx1(entry.0, x1, muf2));
-                            let xfx2 = if two_caches {
-                                *pdf_cache2
-                                    .entry((entry.1, ix2, imu2))
-                                    .or_insert_with(|| xfx2(entry.1, x2, muf2))
-                            } else {
-                                *pdf_cache1
-                                    .entry((entry.1, ix2, imu2))
-                                    .or_insert_with(|| xfx2(entry.1, x2, muf2))
-                            };
-                            lumi += xfx1 * xfx2 * entry.2 / (x1 * x2);
-                        }
-
-                        let alphas = alphas_cache
-                            .entry(xir_values.len() * imu2 + xir_index)
-                            .or_insert_with(|| alphas(xir * xir * mu2.ren));
-
-                        lumi *= alphas.powi(order.alphas.try_into().unwrap());
-                        lumi
-                    })
-                } else {
-                    todo!();
-                };
-
-                // enhance with logarithms
-                if order.logxir > 0 {
-                    value *= (xir * xir).ln().powi(order.logxir.try_into().unwrap());
-                }
-
-                if order.logxif > 0 {
-                    value *= (xif * xif).ln().powi(order.logxif.try_into().unwrap());
-                }
-
-                // fill in ...
-                bins[l + xi.len() * bin_index] += value / bin_sizes[j];
-            }
-        }
-
-        bins
-    }
-
-    /// TODO
-    ///
-    /// # Panics
-    ///
-    /// TODO
-    pub fn convolute2(
         &self,
         lumi_cache: &mut LumiCache,
         order_mask: &[bool],
@@ -588,6 +431,8 @@ impl Grid {
         bins
     }
 
+    // TODO: migrate to use a LumiCache
+
     /// Convolutes a single subgrid `(order, bin, lumi)` with the PDFs strong coupling given by
     /// `xfx1`, `xfx2` and `alphas`. The convolution result is fully differentially, such that the
     /// axes of the result correspond to the values given by the subgrid `q2`, `x1` and `x2` grid
@@ -701,29 +546,58 @@ impl Grid {
         }
     }
 
-    /// Constructs a `Grid` by deserializing it from `reader`. Reading is not buffered.
+    /// Construct a `Grid` by deserializing it from `reader`. Reading is buffered.
     ///
     /// # Errors
     ///
     /// If reading from the compressed or uncompressed stream fails an error is returned.
-    pub fn read(mut reader: impl Read + Seek) -> anyhow::Result<Self> {
-        match LZ4FrameReader::new(&mut reader) {
-            Ok(reader) => Ok(bincode::deserialize_from(reader.into_read())?),
-            Err(WrongMagic(_)) => {
-                reader.seek(SeekFrom::Start(0))?;
-                Ok(bincode::deserialize_from(reader)?)
+    pub fn read(mut reader: impl Read + Seek) -> Result<Self, GridError> {
+        let mut result = bincode::deserialize_from(FrameDecoder::new(&mut reader));
+
+        // check if there was a `WrongMagicNumber` error, in which case `reader` isn't LZ4
+        // compressed
+        let uncompressed = result.as_ref().err().map_or(false, |e| {
+            if let bincode::ErrorKind::Io(io_error) = e.as_ref() {
+                let wrong_magic_number = io_error
+                    .get_ref()
+                    .and_then(|error| error.downcast_ref::<frame::Error>())
+                    .map_or(false, |frame_error| {
+                        matches!(frame_error, frame::Error::WrongMagicNumber)
+                    });
+
+                if wrong_magic_number {
+                    return true;
+                }
             }
-            Err(e) => Err(anyhow::Error::new(e)),
+
+            false
+        });
+
+        if uncompressed {
+            // go back to the start and try without compression
+            reader.rewind().map_err(GridError::IoFailure)?;
+            result = bincode::deserialize_from(BufReader::new(reader));
         }
+
+        result.map_err(GridError::ReadFailure)
     }
 
-    /// Serializes `self` into `writer`. Writing is not buffered.
+    /// Serializes `self` into `writer`. Writing is buffered.
     ///
     /// # Errors
     ///
     /// If writing fails an error is returned.
-    pub fn write(&self, writer: impl Write) -> anyhow::Result<()> {
-        Ok(bincode::serialize_into(writer, self)?)
+    pub fn write(&self, writer: impl Write) -> Result<(), GridError> {
+        bincode::serialize_into(BufWriter::new(writer), self).map_err(GridError::WriteFailure)
+    }
+
+    /// Serializes `self` into `writer`, using LZ4 compression. Writing is buffered.
+    ///
+    /// # Errors
+    ///
+    /// If writing or compression fails an error is returned.
+    pub fn write_lz4(&self, writer: impl Write) -> Result<(), GridError> {
+        bincode::serialize_into(FrameEncoder::new(writer), self).map_err(GridError::WriteFailure)
     }
 
     /// Fills the grid with events for the parton momentum fractions `x1` and `x2`, the scale `q2`,
@@ -1273,7 +1147,12 @@ impl Grid {
     ///
     /// Panics if the parameters do not match with the given grid.
     #[must_use]
-    pub fn convolute_eko(&self, operator: Array5<f64>, eko_info: EkoInfo) -> Option<FkTable> {
+    pub fn convolute_eko(
+        &self,
+        operator: Array5<f64>,
+        eko_info: EkoInfo,
+        order_mask: &[bool],
+    ) -> Option<FkTable> {
         // Check operator layout
         let dim = operator.shape();
 
@@ -1308,23 +1187,23 @@ impl Grid {
         };
 
         let pids1 = if has_pdf1 {
-            eko_info.grid_axes.pids.to_vec()
+            eko_info.grid_axes.pids.clone()
         } else {
             vec![initial_state_1]
         };
         let pids2 = if has_pdf2 {
-            eko_info.grid_axes.pids.to_vec()
+            eko_info.grid_axes.pids.clone()
         } else {
             vec![initial_state_2]
         };
         // create target luminosities
         let tgt_pids1 = if has_pdf1 {
-            eko_info.target_pids.to_vec()
+            eko_info.target_pids.clone()
         } else {
             vec![initial_state_1]
         };
         let tgt_pids2 = if has_pdf2 {
-            eko_info.target_pids.to_vec()
+            eko_info.target_pids.clone()
         } else {
             vec![initial_state_2]
         };
@@ -1404,6 +1283,10 @@ impl Grid {
                 // iterate over the source grid orders and add all of them together into
                 // `src_array`, using the right powers of alphas
                 for (order, powers) in self.orders.iter().enumerate() {
+                    if !order_mask.is_empty() && !order_mask[order] {
+                        continue;
+                    }
+
                     let logs = if (eko_info.xir, eko_info.xif) == (1.0, 1.0) {
                         if (powers.logxir > 0) || (powers.logxif > 0) {
                             continue;
@@ -1446,10 +1329,12 @@ impl Grid {
                             .muf2_grid
                             .iter()
                             .position(|&q2| q2 == eko_info.xir * eko_info.xir * scale)
-                            .expect(&format!(
-                                "Couldn't find q2: {:?} with xir: {:?} and muf2_grid: {:?}",
-                                scale, eko_info.xir, eko_info.grid_axes.muf2_grid
-                            ));
+                            .unwrap_or_else(|| {
+                                panic!(
+                                    "Couldn't find q2: {:?} with xir: {:?} and muf2_grid: {:?}",
+                                    scale, eko_info.xir, eko_info.grid_axes.muf2_grid
+                                )
+                            });
 
                         let ix1 = if invert_x && has_pdf1 {
                             eko_info.grid_axes.x_grid.len() - ix1 - 1
@@ -1979,7 +1864,7 @@ mod tests {
                 .collect(),
         )
         .unwrap();
-        let fk = grid.convolute_eko(operator, eko_info).unwrap();
+        let fk = grid.convolute_eko(operator, eko_info, &[]).unwrap();
 
         assert_eq!(fk.bins(), 1);
     }
diff --git a/pineappl/src/lumi.rs b/pineappl/src/lumi.rs
index 5507cb61..57a8a503 100644
--- a/pineappl/src/lumi.rs
+++ b/pineappl/src/lumi.rs
@@ -105,7 +105,7 @@ impl LumiEntry {
     }
 }
 
-/// Helper macro to quickly generate a LumiEntry at compile time.
+/// Helper macro to quickly generate a `LumiEntry` at compile time.
 ///
 /// # Examples
 ///
diff --git a/pineappl/tests/drell_yan_lo.rs b/pineappl/tests/drell_yan_lo.rs
index 2277c759..3af72044 100644
--- a/pineappl/tests/drell_yan_lo.rs
+++ b/pineappl/tests/drell_yan_lo.rs
@@ -2,6 +2,7 @@ use float_cmp::approx_eq;
 use lhapdf::Pdf;
 use pineappl::bin::BinRemapper;
 use pineappl::grid::{Grid, Ntuple, Order};
+use pineappl::lumi::LumiCache;
 use pineappl::lumi_entry;
 use pineappl::subgrid::{ExtraSubgridParams, SubgridParams};
 use rand::Rng;
@@ -176,8 +177,8 @@ fn dy_aa_lagrange_subgrid_static() -> anyhow::Result<()> {
     assert!(lhapdf::available_pdf_sets().iter().any(|x| x == &pdf_set));
 
     let pdf = Pdf::with_setname_and_member(&pdf_set, 0);
-    let xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
-    let alphas = |_| 0.0;
+    let mut xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
+    let mut alphas = |_| 0.0;
 
     // check `read` and `write`
     let mut file = Cursor::new(Vec::new());
@@ -190,7 +191,8 @@ fn dy_aa_lagrange_subgrid_static() -> anyhow::Result<()> {
     grid.scale_by_order(10.0, 0.5, 10.0, 10.0, 1.0);
     grid.scale_by_order(10.0, 1.0, 10.0, 10.0, 4.0);
 
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let mut lumi_cache = LumiCache::with_one(2212, &mut xfx, &mut alphas);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
     let reference = vec![
         5.29438499470369e-1,
         5.407794857747981e-1,
@@ -236,7 +238,7 @@ fn dy_aa_lagrange_subgrid_static() -> anyhow::Result<()> {
             .collect::<Vec<(f64, f64)>>(),
     )?)?;
 
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
 
     // results are slightly different because of the static scale detection - the interpolation
     // error in the Q^2 dimension is removed
@@ -295,8 +297,8 @@ fn dy_aa_lagrange_subgrid_dynamic() -> anyhow::Result<()> {
     assert!(lhapdf::available_pdf_sets().iter().any(|x| x == &pdf_set));
 
     let pdf = Pdf::with_setname_and_member(&pdf_set, 0);
-    let xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
-    let alphas = |_| 0.0;
+    let mut xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
+    let mut alphas = |_| 0.0;
 
     // check `read` and `write`
     let mut file = Cursor::new(Vec::new());
@@ -309,7 +311,8 @@ fn dy_aa_lagrange_subgrid_dynamic() -> anyhow::Result<()> {
     grid.scale_by_order(10.0, 0.5, 10.0, 10.0, 1.0);
     grid.scale_by_order(10.0, 1.0, 10.0, 10.0, 4.0);
 
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let mut lumi_cache = LumiCache::with_one(2212, &mut xfx, &mut alphas);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
     let reference = vec![
         5.093090431949207e-1,
         5.191668797562395e-1,
@@ -345,7 +348,7 @@ fn dy_aa_lagrange_subgrid_dynamic() -> anyhow::Result<()> {
     grid.optimize();
 
     // check that the results are still the same
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
 
     for (result, reference) in bins.iter().zip(reference.iter()) {
         assert!(approx_eq!(f64, *result, *reference, ulps = 16));
@@ -375,8 +378,8 @@ fn dy_aa_lagrange_subgrid_v2_dynamic() -> anyhow::Result<()> {
     assert!(lhapdf::available_pdf_sets().iter().any(|x| x == &pdf_set));
 
     let pdf = Pdf::with_setname_and_member(&pdf_set, 0);
-    let xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
-    let alphas = |_| 0.0;
+    let mut xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
+    let mut alphas = |_| 0.0;
 
     // check `read` and `write`
     let mut file = Cursor::new(Vec::new());
@@ -389,7 +392,8 @@ fn dy_aa_lagrange_subgrid_v2_dynamic() -> anyhow::Result<()> {
     grid.scale_by_order(10.0, 0.5, 10.0, 10.0, 1.0);
     grid.scale_by_order(10.0, 1.0, 10.0, 10.0, 4.0);
 
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let mut lumi_cache = LumiCache::with_one(2212, &mut xfx, &mut alphas);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
     let reference = vec![
         5.093090431949207e-1,
         5.191668797562395e-1,
@@ -425,7 +429,7 @@ fn dy_aa_lagrange_subgrid_v2_dynamic() -> anyhow::Result<()> {
     //grid.optimize();
 
     // check that the results are still the same
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
 
     for (result, reference) in bins.iter().zip(reference.iter()) {
         assert!(approx_eq!(f64, *result, *reference, ulps = 16));
@@ -455,8 +459,8 @@ fn dy_aa_lagrange_sparse_subgrid_dynamic() -> anyhow::Result<()> {
     assert!(lhapdf::available_pdf_sets().iter().any(|x| x == &pdf_set));
 
     let pdf = Pdf::with_setname_and_member(&pdf_set, 0);
-    let xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
-    let alphas = |_| 0.0;
+    let mut xfx = |id, x, q2| pdf.xfx_q2(id, x, q2);
+    let mut alphas = |_| 0.0;
 
     // check `read` and `write`
     let mut file = Cursor::new(Vec::new());
@@ -469,7 +473,8 @@ fn dy_aa_lagrange_sparse_subgrid_dynamic() -> anyhow::Result<()> {
     grid.scale_by_order(10.0, 0.5, 10.0, 10.0, 1.0);
     grid.scale_by_order(10.0, 1.0, 10.0, 10.0, 4.0);
 
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let mut lumi_cache = LumiCache::with_one(2212, &mut xfx, &mut alphas);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
     let reference = vec![
         5.093090431949207e-1,
         5.191668797562395e-1,
@@ -505,7 +510,7 @@ fn dy_aa_lagrange_sparse_subgrid_dynamic() -> anyhow::Result<()> {
     grid.optimize();
 
     // check that the results are still the same
-    let bins = grid.convolute(&xfx, &xfx, &alphas, &[], &[], &[], &[(1.0, 1.0)]);
+    let bins = grid.convolute(&mut lumi_cache, &[], &[], &[], &[(1.0, 1.0)]);
 
     for (result, reference) in bins.iter().zip(reference.iter()) {
         assert!(approx_eq!(f64, *result, *reference, ulps = 16));
diff --git a/pineappl_capi/src/lib.rs b/pineappl_capi/src/lib.rs
index 8c09c697..149d115b 100644
--- a/pineappl_capi/src/lib.rs
+++ b/pineappl_capi/src/lib.rs
@@ -62,7 +62,7 @@ use pineappl::bin::BinRemapper;
 use pineappl::empty_subgrid::EmptySubgridV1;
 use pineappl::grid::{Grid, Ntuple, Order};
 use pineappl::import_only_subgrid::ImportOnlySubgridV2;
-use pineappl::lumi::LumiEntry;
+use pineappl::lumi::{LumiCache, LumiEntry};
 use pineappl::sparse_array3::SparseArray3;
 use pineappl::subgrid::{ExtraSubgridParams, Mu2, Subgrid, SubgridParams};
 use std::collections::HashMap;
@@ -248,10 +248,69 @@ pub unsafe extern "C" fn pineappl_grid_bin_limits_right(
     slice::from_raw_parts_mut(right, limits.len()).copy_from_slice(&limits);
 }
 
-/// Convolutes the specified grid with the PDFs `xfx1` and `xfx2` and strong coupling `alphas`.
-/// These functions must evaluate the PDFs for the given `x` and `q2` for the parton with the given
-/// PDG id, `pdg_id`, and return the result. Note that the value must be the PDF multiplied with
-/// its argument `x`. The value of the pointer `state` provided to these functions is the same one
+/// Convolutes the specified grid with the PDF `xfx`, which is the PDF for a hadron with the PDG id
+/// `pdg_id`, and strong coupling `alphas`. These functions must evaluate the PDFs for the given
+/// `x` and `q2` for the parton with the given PDG id, `pdg_id`, and return the result. Note that
+/// the value must be the PDF multiplied with its argument `x`. The value of the pointer `state`
+/// provided to these functions is the same one given to this function. The parameter `order_mask`
+/// must be as long as there are perturbative orders contained in `grid` and is used to selectively
+/// disable (`false`) or enable (`true`) individual orders. If `order_mask` is set to `NULL`, all
+/// orders are active. The parameter `lumi_mask` can be used similarly, but must be as long as the
+/// luminosity function `grid` was created with has entries, or `NULL` to enable all luminosities.
+/// The values `xi_ren` and `xi_fac` can be used to vary the renormalization and factorization from
+/// its central value, which corresponds to `1.0`. After convolution of the grid with the PDFs the
+/// differential cross section for each bin is written into `results`.
+///
+/// # Safety
+///
+/// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer,
+/// this function is not safe to call. The function pointers `xfx1`, `xfx2`, and `alphas` must not
+/// be null pointers and point to valid functions. The parameters `order_mask` and `lumi_mask` must
+/// either be null pointers or point to arrays that are as long as `grid` has orders and lumi
+/// entries, respectively. Finally, `results` must be as long as `grid` has bins.
+#[no_mangle]
+pub unsafe extern "C" fn pineappl_grid_convolute_with_one(
+    grid: *const Grid,
+    pdg_id: i32,
+    xfx: extern "C" fn(pdg_id: i32, x: f64, q2: f64, state: *mut c_void) -> f64,
+    alphas: extern "C" fn(q2: f64, state: *mut c_void) -> f64,
+    state: *mut c_void,
+    order_mask: *const bool,
+    lumi_mask: *const bool,
+    xi_ren: f64,
+    xi_fac: f64,
+    results: *mut f64,
+) {
+    let grid = &*grid;
+    let mut pdf = |id, x, q2| xfx(id, x, q2, state);
+    let mut als = |q2| alphas(q2, state);
+    let order_mask = if order_mask.is_null() {
+        vec![]
+    } else {
+        slice::from_raw_parts(order_mask, grid.orders().len()).to_vec()
+    };
+    let lumi_mask = if lumi_mask.is_null() {
+        vec![]
+    } else {
+        slice::from_raw_parts(lumi_mask, grid.lumi().len()).to_vec()
+    };
+    let results = slice::from_raw_parts_mut(results, grid.bin_info().bins());
+    let mut lumi_cache = LumiCache::with_one(pdg_id, &mut pdf, &mut als);
+
+    results.copy_from_slice(&grid.convolute(
+        &mut lumi_cache,
+        &order_mask,
+        &[],
+        &lumi_mask,
+        &[(xi_ren, xi_fac)],
+    ));
+}
+
+/// Convolutes the specified grid with the PDFs `xfx1` and `xfx2`, which are the PDFs of hadrons
+/// with PDG ids `pdg_id1` and `pdg_id2`, respectively, and strong coupling `alphas`. These
+/// functions must evaluate the PDFs for the given `x` and `q2` for the parton with the given PDG
+/// id, `pdg_id`, and return the result. Note that the value must be the PDF multiplied with its
+/// argument `x`. The value of the pointer `state` provided to these functions is the same one
 /// given to this function. The parameter `order_mask` must be as long as there are perturbative
 /// orders contained in `grid` and is used to selectively disable (`false`) or enable (`true`)
 /// individual orders. If `order_mask` is set to `NULL`, all orders are active. The parameter
@@ -264,14 +323,16 @@ pub unsafe extern "C" fn pineappl_grid_bin_limits_right(
 /// # Safety
 ///
 /// If `grid` does not point to a valid `Grid` object, for example when `grid` is the null pointer,
-/// this function is not safe to call. The function pointers `xfx1`, `xfx2`, and `alphas` must be
-/// null pointers and point to valid functions. The parameters `order_mask` and `lumi_mask` must
+/// this function is not safe to call. The function pointers `xfx1`, `xfx2`, and `alphas` must not
+/// be null pointers and point to valid functions. The parameters `order_mask` and `lumi_mask` must
 /// either be null pointers or point to arrays that are as long as `grid` has orders and lumi
 /// entries, respectively. Finally, `results` must be as long as `grid` has bins.
 #[no_mangle]
-pub unsafe extern "C" fn pineappl_grid_convolute(
+pub unsafe extern "C" fn pineappl_grid_convolute_with_two(
     grid: *const Grid,
+    pdg_id1: i32,
     xfx1: extern "C" fn(pdg_id: i32, x: f64, q2: f64, state: *mut c_void) -> f64,
+    pdg_id2: i32,
     xfx2: extern "C" fn(pdg_id: i32, x: f64, q2: f64, state: *mut c_void) -> f64,
     alphas: extern "C" fn(q2: f64, state: *mut c_void) -> f64,
     state: *mut c_void,
@@ -282,9 +343,9 @@ pub unsafe extern "C" fn pineappl_grid_convolute(
     results: *mut f64,
 ) {
     let grid = &*grid;
-    let pdf1 = |id, x, q2| xfx1(id, x, q2, state);
-    let pdf2 = |id, x, q2| xfx2(id, x, q2, state);
-    let als = |q2| alphas(q2, state);
+    let mut pdf1 = |id, x, q2| xfx1(id, x, q2, state);
+    let mut pdf2 = |id, x, q2| xfx2(id, x, q2, state);
+    let mut als = |q2| alphas(q2, state);
     let order_mask = if order_mask.is_null() {
         vec![]
     } else {
@@ -296,11 +357,10 @@ pub unsafe extern "C" fn pineappl_grid_convolute(
         slice::from_raw_parts(lumi_mask, grid.lumi().len()).to_vec()
     };
     let results = slice::from_raw_parts_mut(results, grid.bin_info().bins());
+    let mut lumi_cache = LumiCache::with_two(pdg_id1, &mut pdf1, pdg_id2, &mut pdf2, &mut als);
 
     results.copy_from_slice(&grid.convolute(
-        &pdf1,
-        &pdf2,
-        &als,
+        &mut lumi_cache,
         &order_mask,
         &[],
         &lumi_mask,
diff --git a/pineappl_cli/src/helpers.rs b/pineappl_cli/src/helpers.rs
index 0ca2ae2d..4b62bd35 100644
--- a/pineappl_cli/src/helpers.rs
+++ b/pineappl_cli/src/helpers.rs
@@ -6,25 +6,28 @@ use pineappl::lumi::LumiCache;
 use prettytable::format::{FormatBuilder, LinePosition, LineSeparator};
 use prettytable::Table;
 use std::fs::{File, OpenOptions};
-use std::io::{BufReader, BufWriter};
+use std::path::Path;
 
 pub const ONE_SIGMA: f64 = 68.268_949_213_708_58;
 
 pub fn read_grid(input: &str) -> Result<Grid> {
-    Grid::read(BufReader::new(
-        File::open(input).context(format!("unable to open '{}'", input))?,
-    ))
-    .context(format!("unable to read '{}'", input))
+    Grid::read(File::open(input).context(format!("unable to open '{}'", input))?)
+        .context(format!("unable to read '{}'", input))
 }
 
 pub fn write_grid(output: &str, grid: &Grid) -> Result<()> {
-    grid.write(BufWriter::new(
-        OpenOptions::new()
-            .write(true)
-            .create_new(true)
-            .open(output)
-            .context(format!("unable to write '{}'", output))?,
-    ))
+    let path = Path::new(output);
+    let file = OpenOptions::new()
+        .write(true)
+        .create_new(true)
+        .open(output)
+        .context(format!("unable to write '{}'", output))?;
+
+    if path.extension().map_or(false, |ext| ext == "lz4") {
+        Ok(grid.write_lz4(file)?)
+    } else {
+        Ok(grid.write(file)?)
+    }
 }
 
 pub fn create_table() -> Table {
@@ -101,7 +104,7 @@ pub fn convolute(
     let mut alphas = |q2| lhapdf.alphas_q2(q2);
     let mut cache = LumiCache::with_one(pdf_pdg_id, &mut pdf, &mut alphas);
 
-    grid.convolute2(&mut cache, &orders, bins, lumis, &SCALES_VECTOR[0..scales])
+    grid.convolute(&mut cache, &orders, bins, lumis, &SCALES_VECTOR[0..scales])
 }
 
 pub fn convolute_subgrid(
diff --git a/pineappl_py/pineappl/grid.py b/pineappl_py/pineappl/grid.py
index 0a6158c0..7b0c499a 100644
--- a/pineappl_py/pineappl/grid.py
+++ b/pineappl_py/pineappl/grid.py
@@ -61,6 +61,49 @@ def create(cls, lumi, orders, bin_limits, subgrid_params):
         orders = [o.raw for o in orders]
         return cls(PyGrid(lumi, orders, bin_limits, subgrid_params.raw))
 
+    def subgrid(self, order, bin_, lumi):
+        """
+        Retrieve the subgrid at the given position.
+
+        Convenience wrapper for :meth:`pineappl.pineappl.PyGrid.set_subgrid()`.
+
+        Parameters
+        ----------
+            order : int
+                index of order
+            bin_ : int
+                index of bin
+            lumi : int
+                index of luminosity
+
+        Returns
+        -------
+            subgrid : Subgrid
+                subgrid content
+        """
+        return self.raw.subgrid(order, bin_, lumi)
+
+    def __getitem__(self, key):
+        """
+        Retrieve the subgrid at the given position.
+
+        Syntactic sugar for :meth:`subgrid`
+
+        Parameters
+        ----------
+            key : (int, int, int)
+                a 3-element integers tuple, consisting in `(order, bin, lumi)`
+
+        Returns
+        -------
+            subgrid : Subgrid
+                subgrid content
+        """
+        if len(key) != 3:
+            raise ValueError("A tuple with `(order, bin, lumi)` is required as key.")
+
+        return self.subgrid(*key)
+
     def set_subgrid(self, order, bin_, lumi, subgrid):
         """
         Set the subgrid at the given position.
@@ -80,6 +123,24 @@ def set_subgrid(self, order, bin_, lumi, subgrid):
         """
         self.raw.set_subgrid(order, bin_, lumi, subgrid.into())
 
+    def __setitem__(self, key, subgrid):
+        """
+        Set the subgrid at the given position.
+
+        Syntactic sugar for :meth:`set_subgrid`
+
+        Parameters
+        ----------
+            key : (int, int, int)
+                a 3-element integers tuple, consisting in `(order, bin, lumi)`
+            subgrid : ImportOnlySubgridV1
+                subgrid content
+        """
+        if len(key) != 3:
+            raise ValueError("A tuple with `(order, bin, lumi)` is required as key.")
+
+        self.set_subgrid(*key, subgrid)
+
     def set_remapper(self, remapper):
         """
         Set the normalizations.
@@ -106,10 +167,10 @@ def orders(self):
         """
         return [Order(*pyorder.as_tuple()) for pyorder in self.raw.orders()]
 
-    def convolute(
+    def convolute_with_one(
         self,
-        xfx1,
-        xfx2,
+        pdg_id,
+        xfx,
         alphas,
         order_mask=(),
         bin_indices=(),
@@ -121,10 +182,10 @@ def convolute(
 
         Parameters
         ----------
-            xfx1 : callable
-                lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x_1`-grid
-            xfx2 : callable
-                lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x_2`-grid
+            pdg_id : int
+                PDG Monte Carlo ID of the hadronic particle `xfx` is the PDF for
+            xfx : callable
+                lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x`-grid
             alphas : callable
                 lhapdf like callable with arguments `Q2` returning :math:`\alpha_s`
             order_mask : list(bool)
@@ -151,11 +212,11 @@ def convolute(
                 the scale variation)
 
         """
-        return self.raw.convolute(
-            xfx1, xfx2, alphas, order_mask, bin_indices, lumi_mask, xi
+        return self.raw.convolute_with_one(
+            pdg_id, xfx, alphas, order_mask, bin_indices, lumi_mask, xi
         )
 
-    def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids"):
+    def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids", order_mask=()):
         """
         Create an FKTable with the EKO.
 
@@ -168,6 +229,9 @@ def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids"):
             lumi_id_types : str
                 kind of lumi types (e.g. "pdg_mc_ids" for flavor basis, "evol"
                 for evolution basis)
+            order_mask : list(bool)
+                Mask for selecting specific orders. The value `True` means the corresponding order
+                is included. An empty list corresponds to all orders being enabled.
 
         Returns
         ------
@@ -191,6 +255,7 @@ def convolute_eko(self, operators, lumi_id_types="pdg_mc_ids"):
                 operator_grid.flatten(),
                 operator_grid.shape,
                 lumi_id_types,
+                order_mask,
             )
         )
 
diff --git a/pineappl_py/src/fk_table.rs b/pineappl_py/src/fk_table.rs
index 9fe28a65..062f9510 100644
--- a/pineappl_py/src/fk_table.rs
+++ b/pineappl_py/src/fk_table.rs
@@ -1,6 +1,7 @@
 use numpy::{IntoPyArray, PyArray4};
 use pineappl::fk_table::FkTable;
 use pineappl::grid::Grid;
+use pineappl::lumi::LumiCache;
 use pyo3::prelude::*;
 use std::convert::TryFrom;
 use std::fs::File;
@@ -95,23 +96,19 @@ impl PyFkTable {
     ///
     /// Parameters
     /// ----------
-    ///     xfx1 : callable
-    ///         lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x_1`-grid
-    ///     xfx2 : callable
-    ///         lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x_2`-grid
+    ///     pdg_id : integer
+    ///         PDG Monte Carlo ID of the hadronic particle `xfx` is the PDF for
+    ///     xfx : callable
+    ///         lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x`-grid
     ///
     /// Returns
     /// -------
     ///     list(float) :
     ///         cross sections for all bins
-    pub fn convolute(&self, xfx1: &PyAny, xfx2: &PyAny) -> Vec<f64> {
-        self.fk_table.convolute(
-            &|id, x, q2| f64::extract(xfx1.call1((id, x, q2)).unwrap()).unwrap(),
-            &|id, x, q2| f64::extract(xfx2.call1((id, x, q2)).unwrap()).unwrap(),
-            &[],
-            &[],
-            &[],
-            &[(1.0, 1.0)],
-        )
+    pub fn convolute_with_one(&self, pdg_id: i32, xfx: &PyAny) -> Vec<f64> {
+        let mut xfx = |id, x, q2| f64::extract(xfx.call1((id, x, q2)).unwrap()).unwrap();
+        let mut alphas = |_| 1.0;
+        let mut lumi_cache = LumiCache::with_one(pdg_id, &mut xfx, &mut alphas);
+        self.fk_table.convolute(&mut lumi_cache, &[], &[])
     }
 }
diff --git a/pineappl_py/src/grid.rs b/pineappl_py/src/grid.rs
index 417d63c7..5527b747 100644
--- a/pineappl_py/src/grid.rs
+++ b/pineappl_py/src/grid.rs
@@ -1,4 +1,5 @@
 use pineappl::grid::{EkoInfo, Grid, GridAxes, Ntuple, Order};
+use pineappl::lumi::LumiCache;
 
 use super::bin::PyBinRemapper;
 use super::fk_table::PyFkTable;
@@ -7,9 +8,11 @@ use super::subgrid::{PySubgridEnum, PySubgridParams};
 
 use ndarray::{Array, Ix5};
 
+use std::collections::HashMap;
 use std::fs::File;
 use std::io::BufReader;
 
+use pyo3::exceptions::PyValueError;
 use pyo3::prelude::*;
 
 /// PyO3 wrapper to :rustdoc:`pineappl::grid::Order <grid/struct.Order.html>`
@@ -124,6 +127,19 @@ impl PyGrid {
         );
     }
 
+    /// Get metadata values stored in the grid.
+    ///
+    ///
+    /// Parameters
+    /// ----------
+    ///     key : str
+    ///         key
+    ///     value : str
+    ///         value
+    pub fn key_values(&self) -> HashMap<String, String> {
+        self.grid.key_values().unwrap().clone()
+    }
+
     /// Set a metadata key-value pair in the grid.
     ///
     /// **Usage:** `yadism`
@@ -138,6 +154,15 @@ impl PyGrid {
         self.grid.set_key_value(key, value);
     }
 
+    /// Retrieve a subgrid.
+    ///
+    /// **Usage:** `yadism`
+    pub fn subgrid(&self, order: usize, bin: usize, lumi: usize) -> PySubgridEnum {
+        PySubgridEnum {
+            subgrid_enum: self.grid.subgrid(order, bin, lumi).clone(),
+        }
+    }
+
     /// Set a subgrid.
     ///
     /// **Usage:** `yadism`
@@ -185,10 +210,10 @@ impl PyGrid {
     ///
     /// Parameters
     /// ----------
-    ///     xfx1 : callable
-    ///         lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x_1`-grid
-    ///     xfx2 : callable
-    ///         lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x_2`-grid
+    ///     pdg_id : int
+    ///         PDG Monte Carlo ID of the hadronic particle `xfx` is the PDF for
+    ///     xfx : callable
+    ///         lhapdf like callable with arguments `pid, x, Q2` returning x*pdf for :math:`x`-grid
     ///     alphas : callable
     ///         lhapdf like callable with arguments `Q2` returning :math:`\alpha_s`
     ///     order_mask : list(bool)
@@ -213,25 +238,21 @@ impl PyGrid {
     ///     list(float) :
     ///         cross sections for all bins, for each scale-variation tuple (first all bins, then
     ///         the scale variation)
-    pub fn convolute(
+    pub fn convolute_with_one(
         &self,
-        xfx1: &PyAny,
-        xfx2: &PyAny,
+        pdg_id: i32,
+        xfx: &PyAny,
         alphas: &PyAny,
         order_mask: Vec<bool>,
         bin_indices: Vec<usize>,
         lumi_mask: Vec<bool>,
         xi: Vec<(f64, f64)>,
     ) -> Vec<f64> {
-        self.grid.convolute(
-            &|id, x, q2| f64::extract(xfx1.call1((id, x, q2)).unwrap()).unwrap(),
-            &|id, x, q2| f64::extract(xfx2.call1((id, x, q2)).unwrap()).unwrap(),
-            &|q2| f64::extract(alphas.call1((q2,)).unwrap()).unwrap(),
-            &order_mask,
-            &bin_indices,
-            &lumi_mask,
-            &xi,
-        )
+        let mut xfx = |id, x, q2| f64::extract(xfx.call1((id, x, q2)).unwrap()).unwrap();
+        let mut alphas = |q2| f64::extract(alphas.call1((q2,)).unwrap()).unwrap();
+        let mut lumi_cache = LumiCache::with_one(pdg_id, &mut xfx, &mut alphas);
+        self.grid
+            .convolute(&mut lumi_cache, &order_mask, &bin_indices, &lumi_mask, &xi)
     }
 
     /// Convolute with with an evolution operator.
@@ -277,6 +298,7 @@ impl PyGrid {
         operator_flattened: Vec<f64>,
         operator_shape: Vec<usize>,
         lumi_id_types: String,
+        order_mask: Vec<bool>,
     ) -> PyFkTable {
         let operator = Array::from_shape_vec(operator_shape, operator_flattened).unwrap();
         let eko_info = EkoInfo {
@@ -296,7 +318,11 @@ impl PyGrid {
 
         let evolved_grid = self
             .grid
-            .convolute_eko(operator.into_dimensionality::<Ix5>().unwrap(), eko_info)
+            .convolute_eko(
+                operator.into_dimensionality::<Ix5>().unwrap(),
+                eko_info,
+                &order_mask,
+            )
             .unwrap();
         PyFkTable {
             fk_table: evolved_grid,
@@ -333,6 +359,16 @@ impl PyGrid {
         self.grid.write(File::create(path).unwrap()).unwrap();
     }
 
+    /// Write grid to compressed file.
+    ///
+    /// Parameters
+    /// ----------
+    ///     path : str
+    ///         file path
+    pub fn write_lz4(&self, path: &str) {
+        self.grid.write_lz4(File::create(path).unwrap()).unwrap();
+    }
+
     /// Optimize grid content.
     ///
     /// **Usage:** `yadism`
@@ -341,9 +377,12 @@ impl PyGrid {
     }
 
     /// Optimize grid content.
-    // pub fn merge(&mut self, mut other: Self) {
-    // self.grid.merge(other.grid);
-    // }
+    pub fn merge_from_file(&mut self, path: &str) -> PyResult<()> {
+        match self.grid.merge(Self::read(path).grid) {
+            Ok(()) => Ok(()),
+            Err(x) => Err(PyValueError::new_err(format!("{:?}", x))),
+        }
+    }
 
     /// Extract the number of dimensions for bins.
     ///
diff --git a/pineappl_py/tests/test_grid.py b/pineappl_py/tests/test_grid.py
index 4c7cccb5..1c7c23b2 100644
--- a/pineappl_py/tests/test_grid.py
+++ b/pineappl_py/tests/test_grid.py
@@ -80,7 +80,7 @@ def test_bins(self):
         np.testing.assert_allclose(g.bin_left(1), [2,3])
         np.testing.assert_allclose(g.bin_right(1), [3,5])
 
-    def test_convolute(self):
+    def test_convolute_with_one(self):
         g = self.fake_grid()
 
         # DIS grid
@@ -93,9 +93,9 @@ def test_convolute(self):
             [1.0],
         )
         g.set_subgrid(0, 0, 0, subgrid)
-        np.testing.assert_allclose(g.convolute(lambda pid,x,q2: 0., lambda pid,x,q2: 0., lambda q2: 0.), [0.]*2)
-        np.testing.assert_allclose(g.convolute(lambda pid,x,q2: 1, lambda pid,x,q2: 1., lambda q2: 1.), [5e6/9999,0.])
-        np.testing.assert_allclose(g.convolute(lambda pid,x,q2: 1, lambda pid,x,q2: 1., lambda q2: 2.), [2**3 * 5e6/9999,0.])
+        np.testing.assert_allclose(g.convolute_with_one(2212, lambda pid,x,q2: 0., lambda q2: 0.), [0.]*2)
+        np.testing.assert_allclose(g.convolute_with_one(2212, lambda pid,x,q2: 1, lambda q2: 1.), [5e6/9999,0.])
+        np.testing.assert_allclose(g.convolute_with_one(2212, lambda pid,x,q2: 1, lambda q2: 2.), [2**3 * 5e6/9999,0.])
 
     def test_axes(self):
         g = self.fake_grid()