diff --git a/src/core/fft_engine.hh b/src/core/fft_engine.hh
index 9c31695..7e0c266 100644
--- a/src/core/fft_engine.hh
+++ b/src/core/fft_engine.hh
@@ -1,136 +1,141 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef FFT_ENGINE_H
 #define FFT_ENGINE_H
 /* -------------------------------------------------------------------------- */
 #include "fftw_interface.hh"
 #include "grid.hh"
 #include "grid_base.hh"
 #include "grid_hermitian.hh"
+#include "partitioner.hh"
 #include <algorithm>
 #include <array>
 #include <map>
 #include <memory>
 #include <string>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 class FFTEngine {
 protected:
   using plan_t = std::pair<fftw::plan<Real>, fftw::plan<Real>>;
   using key_t = std::basic_string<UInt>;
   using complex_t = fftw::helper<Real>::complex;
 
 public:
   virtual ~FFTEngine() noexcept = default;
 
   /// Execute a forward plan on real and spectral 1D
   virtual void forward(Grid<Real, 1>& real,
                        GridHermitian<Real, 1>& spectral) = 0;
   /// Execute a forward plan on real and spectral 2D
   virtual void forward(Grid<Real, 2>& real,
                        GridHermitian<Real, 2>& spectral) = 0;
   /// Execute a backward plan on real and spectral 1D
   virtual void backward(Grid<Real, 1>& real,
                         GridHermitian<Real, 1>& spectral) = 0;
   /// Execute a backward plan on real and spectral 2D
   virtual void backward(Grid<Real, 2>& real,
                         GridHermitian<Real, 2>& spectral) = 0;
 
   /// Fill a grid with wavevector values in appropriate ordering
   template <typename T, UInt dim, bool hermitian>
   static Grid<T, dim> computeFrequencies(const std::array<UInt, dim>& sizes);
 
   /// Instanciate an appropriate FFTEngine subclass
   static std::unique_ptr<FFTEngine>
   makeEngine(unsigned int flags = FFTW_ESTIMATE);
 
 protected:
   /// Make a transform signature from a pair of grids
   template <UInt dim>
   static key_t make_key(Grid<Real, dim>& real,
                         GridHermitian<Real, dim>& spectral);
   /// Cast to FFTW complex type
   static auto cast(Complex* data) { return reinterpret_cast<complex_t*>(data); }
 };
 
 template <UInt dim>
 FFTEngine::key_t FFTEngine::make_key(Grid<Real, dim>& real,
                                      GridHermitian<Real, dim>& spectral) {
   if (real.getNbComponents() != spectral.getNbComponents())
     TAMAAS_EXCEPTION("Components do not match");
 
   auto hermitian_dims =
       GridHermitian<Real, dim>::hermitianDimensions(real.sizes());
   if (not std::equal(hermitian_dims.begin(), hermitian_dims.end(),
                      spectral.sizes().begin()))
     TAMAAS_EXCEPTION("Spectral grid does not have hermitian size");
 
   // Reserve space for dimensions + components + both last strides
   key_t key(real.getDimension() + 3, 0);
   std::copy_n(real.sizes().begin(), dim, key.begin());
   key[dim] = real.getNbComponents();
   key[dim + 1] = real.getStrides().back();
   key[dim + 2] = spectral.getStrides().back();
   return key;
 }
 
 template <typename T, UInt dim, bool hermitian>
-Grid<T, dim> FFTEngine::computeFrequencies(const std::array<UInt, dim>& sizes) {
+Grid<T, dim>
+FFTEngine::computeFrequencies(const std::array<UInt, dim>& global_sizes) {
   // If hermitian is true, we suppose the dimensions of freq are
   // reduced based on hermitian symetry and that it has dim components
 
-  auto& n = sizes;
-  Grid<T, dim> freq(n, dim);
+  Grid<T, dim> freq(Partitioner<dim>::local_size(global_sizes), dim);
   constexpr UInt dmax = dim - static_cast<UInt>(hermitian);
 
+  // For MPI
+  const auto& n = global_sizes;
+  const auto offset = Partitioner<dim>::local_offset(freq) / dim;
+
 #pragma omp parallel for  // to get rid of a compilation warning from nvcc
   for (UInt i = 0; i < freq.getNbPoints(); ++i) {
-    UInt index = i;
-    VectorProxy<T, dim> wavevector(freq(index * freq.getNbComponents()));
+    UInt index = i + offset;
+    VectorProxy<T, dim> wavevector(freq(i * freq.getNbComponents()));
     std::array<UInt, dim> tuple{{0}};
     /// Computing tuple from index
     for (Int d = dim - 1; d >= 0; d--) {
       tuple[d] = index % n[d];
       index -= tuple[d];
       index /= n[d];
     }
 
     if (hermitian)
       wavevector(dim - 1) = tuple[dim - 1];
 
     for (UInt d = 0; d < dmax; d++) {
       // Type conversion
       T td = tuple[d];
       T nd = n[d];
       T q = (tuple[d] < n[d] / 2) ? td : td - nd;
       wavevector(d) = q;
     }
   }
 
   return freq;
 }
 
 }  // namespace tamaas
 
 #endif
diff --git a/src/mpi/fftw_mpi_interface.hh b/src/mpi/fftw_mpi_interface.hh
index 716eead..3992d9d 100644
--- a/src/mpi/fftw_mpi_interface.hh
+++ b/src/mpi/fftw_mpi_interface.hh
@@ -1,65 +1,69 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef FFTW_MPI_INTERFACE
 #define FFTW_MPI_INTERFACE
 /* -------------------------------------------------------------------------- */
 #include "mpi_interface.hh"
 #include <cstdlib>
 #include <numeric>
 #include <tuple>
+#include <stdexcept>
 #ifdef USE_MPI
 #include <fftw3-mpi.h>
 #endif
 /* -------------------------------------------------------------------------- */
 namespace fftw {
 namespace mpi_dummy {
 inline void init() {}
 inline void cleanup() {}
 inline auto local_size_many(int rank, const std::ptrdiff_t* size,
                             std::ptrdiff_t howmany) {
   return std::make_tuple(howmany * std::accumulate(size, size + rank,
                                                    std::ptrdiff_t{1},
                                                    std::multiplies<void>()),
                          size[0], 0);
 }
 }  // namespace mpi_dummy
 
 #ifdef USE_MPI
 namespace mpi_impl {
 inline void init() { fftw_mpi_init(); }
 inline void cleanup() { fftw_mpi_cleanup(); }
 inline auto local_size_many(int rank, const std::ptrdiff_t* size,
                             std::ptrdiff_t howmany) {
+  if (rank < 2)
+    throw std::domain_error("FFTW-MPI cannot be used for 1D transforms");
+
   std::ptrdiff_t local_n0, local_n0_offset;
   auto res = fftw_mpi_local_size_many(
       rank, size, howmany, FFTW_MPI_DEFAULT_BLOCK, ::tamaas::mpi::comm::world,
       &local_n0, &local_n0_offset);
   return std::make_tuple(res, local_n0, local_n0_offset);
 }
 }  // namespace mpi_impl
 namespace mpi = mpi_impl;
 #else
 namespace mpi = mpi_dummy;
 #endif
 }  // namespace fftw
 #endif
diff --git a/tests/test_fftfreq.cpp b/tests/test_fftfreq.cpp
index 2749695..3479a8b 100644
--- a/tests/test_fftfreq.cpp
+++ b/tests/test_fftfreq.cpp
@@ -1,76 +1,98 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "fft_engine.hh"
 #include "grid.hh"
 #include "grid_hermitian.hh"
 #include "grid_view.hh"
+#include "mpi_interface.hh"
+#include "partitioner.hh"
 #include "test.hh"
 
 #include <pybind11/embed.h>
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 
 using namespace tamaas;
 namespace py = pybind11;
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTEngine, Frequencies1D) {
   std::array<UInt, 1> sizes = {{10}};
+
+  if (mpi::initialized()) {
+    auto throw_call = [](auto&& sizes) {
+      return FFTEngine::computeFrequencies<Real, 1, false>(sizes);
+    };
+    ASSERT_THROW(throw_call(sizes), std::domain_error);
+    return;
+  }
+
   auto freq = FFTEngine::computeFrequencies<Real, 1, false>(sizes);
 
   py::module fftfreq = py::module::import("fftfreq");
 
   std::vector<Real> reference(freq.dataSize());
   py::array py_arg(reference.size(), reference.data(), py::none());
 
   fftfreq.attr("frequencies1D")(py_arg);
 
-  ASSERT_TRUE(compare(reference, freq, AreFloatEqual()))
+  EXPECT_TRUE(compare(reference, freq, AreFloatEqual()))
       << "Non hermitian frequencies are wrong";
 
   auto hfreq = FFTEngine::computeFrequencies<Real, 1, true>(sizes);
   std::iota(reference.begin(), reference.end(), 0);
 
-  ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual()))
+  EXPECT_TRUE(compare(reference, hfreq, AreFloatEqual()))
       << "Hermitian frequencies are wrong";
 }
 
 /* -------------------------------------------------------------------------- */
 TEST(TestFFTEngine, Frequencies2D) {
   std::array<UInt, 2> sizes = {{10, 10}};
+
+  Grid<Real, 2> gathered(sizes, 2), hgathered(sizes, 2);
+
   auto freq = FFTEngine::computeFrequencies<Real, 2, false>(sizes);
+  auto hfreq = FFTEngine::computeFrequencies<Real, 2, true>(sizes);
+
+  mpi::gather(freq.getInternalData(), gathered.getInternalData(),
+              freq.dataSize());
+  mpi::gather(hfreq.getInternalData(), hgathered.getInternalData(),
+              hfreq.dataSize());
+
+  if (mpi::rank() != 0)
+    return;
 
   py::module fftfreq = py::module::import("fftfreq");
 
-  std::vector<Real> reference(freq.dataSize());
-  py::array py_arg({10, 10, 2}, reference.data(), py::none());
+  std::vector<Real> reference(gathered.dataSize());
+  py::array py_arg({sizes[0], sizes[1], 2u}, reference.data(), py::none());
 
   fftfreq.attr("frequencies2D")(py_arg);
-  ASSERT_TRUE(compare(reference, freq, AreFloatEqual()))
+  EXPECT_TRUE(compare(reference, gathered, AreFloatEqual()))
       << "Non hermitian frequencies are wrong";
 
-  auto hfreq = FFTEngine::computeFrequencies<Real, 2, true>(sizes);
   fftfreq.attr("hfrequencies2D")(py_arg);
-  ASSERT_TRUE(compare(reference, hfreq, AreFloatEqual()))
+  EXPECT_TRUE(compare(reference, hgathered, AreFloatEqual()))
       << "Hermitian frequencies are wrong";
 }
diff --git a/tests/test_integration.cpp b/tests/test_integration.cpp
index ff4726e..bf36c8a 100644
--- a/tests/test_integration.cpp
+++ b/tests/test_integration.cpp
@@ -1,164 +1,170 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-2020 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "fft_engine.hh"
 #include "integration/accumulator.hh"
+#include "mpi_interface.hh"
 #include "test.hh"
 #include <cmath>
 #include <expolit/expolit>
 /* -------------------------------------------------------------------------- */
 using namespace tamaas;
 namespace ex = expolit;
 
 static Real tol = 1e-13;
 
 class AccumulatorTest : public ::testing::Test {
 protected:
   void SetUp() override {
     nodal_values.resize(n);
 
     for (auto&& grid : nodal_values) {
       grid.setNbComponents(6);
       grid.resize({nh, nh});
       grid = 1;
     }
 
     accumulator.makeUniformMesh(nodal_values.size(), size);
     wavevectors = FFTEngine::computeFrequencies<Real, 2, true>({nh, nh});
   }
 
   UInt n = 10, nh = 10;
   Real size = 3;
   std::vector<GridHermitian<Real, 2>> nodal_values;
   Grid<Real, 2> wavevectors;
   Accumulator<model_type::volume_2d, SymMatrixProxy<Complex, 3>> accumulator;
 };
 
 TEST_F(AccumulatorTest, uniformMesh) {
   auto& pos = this->accumulator.nodePositions();
   std::vector<Real> sol(pos.size());
   std::iota(sol.begin(), sol.end(), 0);
   std::transform(sol.begin(), sol.end(), sol.begin(),
                  [&](auto& x) { return size * x / (sol.size() - 1); });
   ASSERT_TRUE(compare(sol, pos, AreFloatEqual())) << "Uniform mesh fail";
 }
 
 // Litteral type
 using Q = ex::Litteral<struct Q_>;
 
 TEST_F(AccumulatorTest, forward) {
+  if (mpi::size() != 1)
+    GTEST_SKIP();  // test needs work with MPI
+
   // Setup for layer checking
   std::vector<int> layers_seen, layers_to_see(n);
   std::iota(layers_to_see.begin(), layers_to_see.end(), 0);
 
   // Setup for integral checking
   constexpr Q q;
   constexpr auto z = ex::Polynomial<Real, 1>({0, 1});
   constexpr auto g0 = ex::exp(q * z);
   constexpr auto g1 = q * z * ex::exp(q * z);
 
   for (auto&& tuple :
        accumulator.forward(this->nodal_values, this->wavevectors)) {
     auto&& l = std::get<0>(tuple);
     auto&& xl_ = std::get<1>(tuple);
     auto&& acc_g0 = std::get<2>(tuple);
     auto&& acc_g1 = std::get<3>(tuple);
     layers_seen.push_back(l);
 
     // fundamental mode
     auto testing = acc_g0(0, 0, 0);
     ASSERT_NEAR(testing.real(), xl_, tol) << "Fundamental mode G0 fail";
     testing = acc_g1(0, 0, 0);
     ASSERT_NEAR(testing.real(), 0, tol) << "Fundamental mode G1 fail";
 
     // other modes
     testing = acc_g0(1, 0, 0);
     ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), tol) << "Symmetry fail";
     ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), tol) << "Symmetry fail";
 
     testing = acc_g0(1, 1, 0);
     Real sqrt_two = std::sqrt(Real(2));
     auto ref_value = ex::definite_integral(
         std::pair<Real, Real>(0, xl_),
         ex::substitute<Q::tag>(ex::Constant<Real>({sqrt_two}), g0));
     ASSERT_NEAR(testing.real(), ref_value, tol)
         << "Integration of exp(qy)*ϕ(y) fail";
 
     testing = acc_g1(1, 1, 0);
     ref_value = ex::definite_integral(
         std::pair<Real, Real>(0, xl_),
         ex::substitute<Q::tag>(ex::Constant<Real>({sqrt_two}), g1));
     ASSERT_NEAR(testing.real(), ref_value, tol)
         << "Integration of qy*exp(qy)*ϕ(y) fail";
   }
 
   ASSERT_TRUE(compare(layers_seen, layers_to_see))
       << "Did not see correct layers";
 }
 
 TEST_F(AccumulatorTest, backward) {
+  if (mpi::size() != 1)
+    GTEST_SKIP();
   // Setup for layer checking
   std::vector<int> layers_seen, layers_to_see(n);
   std::iota(layers_to_see.rbegin(), layers_to_see.rend(), 0);
 
   // Setup for integral checking
   constexpr Q q;
   constexpr auto z = ex::Polynomial<Real, 1>({0, 1});
   constexpr auto g0 = ex::exp(q * (z * (-1)));
   constexpr auto g1 = q * z * ex::exp(q * ((-1) * z));
 
   for (auto&& tuple : accumulator.backward(nodal_values, wavevectors)) {
     auto&& l = std::get<0>(tuple);
     auto&& xl_ = std::get<1>(tuple);
     auto&& acc_g0 = std::get<2>(tuple);
     auto&& acc_g1 = std::get<3>(tuple);
     layers_seen.push_back(l);
 
     // fundamental mode
     auto testing = acc_g0(0, 0, 0);
     ASSERT_NEAR(testing.real(), size - xl_, tol) << "Fundamental mode G0 fail";
     testing = acc_g1(0, 0, 0);
     ASSERT_NEAR(testing.real(), 0, tol) << "Fundamental mode G1 fail";
 
     // other modes
     testing = acc_g0(1, 0, 0);
     ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), tol) << "Symmetry fail";
     ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), tol) << "Symmetry fail";
 
     testing = acc_g0(1, 1, 0);
     auto ref_value = ex::definite_integral(
         std::make_pair(xl_, size),
         ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g0));
     ASSERT_NEAR(testing.real(), ref_value, tol)
         << "Integration of exp(-qy)*ϕ(y) fail";
 
     testing = acc_g1(1, 1, 0);
     ref_value = ex::definite_integral(
         std::make_pair(xl_, size),
         ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g1));
     ASSERT_NEAR(testing.real(), ref_value, tol)
         << "Integration of qy*exp(-qy)*ϕ(y) fail";
   }
 
   ASSERT_TRUE(compare(layers_seen, layers_to_see))
       << "Did not see correct layers";
 }