diff --git a/src/core/fft_engine.cpp b/src/core/fft_engine.cpp
index 345e459..7274bc0 100644
--- a/src/core/fft_engine.cpp
+++ b/src/core/fft_engine.cpp
@@ -1,89 +1,87 @@
 /**
 
  *  @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 <algorithm>
 #include <functional>
 #include <numeric>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 namespace detail {
 /// RAII helper for fftw_free
 template <typename T>
 struct fftw_ptr {
   T* ptr;
 
   ~fftw_ptr() noexcept {
     if (ptr)
       fftw::free(ptr);
   }
 
   operator T*() { return ptr; }
 };
 }  // namespace detail
 
 FFTEngine::plan_t& FFTEngine::getPlans(key_t key) {
   if (plans.find(key) != plans.end())
     return plans[key];
 
   const int rank = key.size() - 3;  // dimension of fft
   const auto components = key[rank];
 
   std::vector<int> n(rank);  // size of individual fft
   std::copy_n(key.begin(), rank, n.begin());
 
   const int howmany = components;  // how many fft to compute
   const int idist = 1,
             odist = 1;  // components are next to each other in memory
   const int istride = key[rank + 1] * components,
             ostride = key[rank + 2] * components;
   int *inembed = nullptr, *onembed = nullptr;  // row major
 
   detail::fftw_ptr<Real> in{nullptr};
-  detail::fftw_ptr<fftw::helper<Real>::complex> out{nullptr};
+  detail::fftw_ptr<complex_t> out{nullptr};
 
   // Flags does not contain FFTW_ESTIMATE | FFTW_UNALIGNED => allocate in & out
   if (not(_flags & FFTW_ESTIMATE and _flags & FFTW_UNALIGNED)) {
-    const auto size =
-        std::accumulate(n.begin(), n.end(), 1, std::multiplies<void>()) *
-        istride;
-    auto hermit_n = GridHermitian<Real, 1>::hermitianDimensions(n);
-    const auto hsize = std::accumulate(hermit_n.begin(), hermit_n.end(), 1,
-                                       std::multiplies<void>()) *
-                       ostride;
-    in.ptr = fftw::helper<Real>::alloc_real(size);
-    out.ptr = fftw::helper<Real>::alloc_complex(hsize);
+    const auto prod = [](auto&& n) {
+      return std::accumulate(n.cbegin(), n.cend(), 1, std::multiplies<void>());
+    };
+
+    const auto hermit_n = GridHermitian<Real, 1>::hermitianDimensions(n);
+    in.ptr = fftw::helper<Real>::alloc_real(prod(n) * istride);
+    out.ptr = fftw::helper<Real>::alloc_complex(prod(hermit_n) * ostride);
   }
 
   fftw::plan<Real> forward{
       fftw::plan_many_forward(rank, n.data(), howmany, in, inembed, istride,
                               idist, out, onembed, ostride, odist, _flags)};
   fftw::plan<Real> backward{
       fftw::plan_many_backward(rank, n.data(), howmany, out, onembed, ostride,
                                odist, in, inembed, istride, idist, _flags)};
   plans[key] = std::make_pair(std::move(forward), std::move(backward));
   return plans[key];
 }
 
 }  // namespace tamaas
diff --git a/src/core/fft_engine.hh b/src/core/fft_engine.hh
index ffaec6d..46158dc 100644
--- a/src/core/fft_engine.hh
+++ b/src/core/fft_engine.hh
@@ -1,150 +1,150 @@
 /**
  *  @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 <array>
 #include <map>
 #include <string>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 class FFTEngine {
   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:
   /// Initialize with flags
   explicit FFTEngine(unsigned int flags = FFTW_ESTIMATE)
       : _flags(flags), plans() {}
 
   /// Execute a forward plan on real and spectral
   template <UInt dim>
   void forward(Grid<Real, dim>& real, GridHermitian<Real, dim>& spectral);
   /// Execute a backward plan on real and spectral
   template <UInt dim>
   void backward(Grid<Real, dim>& real, GridHermitian<Real, dim>& spectral);
   /// 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);
   /// Get FFTW plan flags
   unsigned int flags() const { return _flags; }
 
 protected:
   /// Return the plans pair for a given transform signature
   plan_t& getPlans(key_t key);
 
   /// 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);
 
 private:
-  unsigned int _flags;             ///< FFTW flags
+  unsigned int _flags;            ///< FFTW flags
   std::map<key_t, plan_t> plans;  ///< plans corresponding to signatures
 };
 
 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);
   thrust::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 <UInt dim>
 void FFTEngine::forward(Grid<Real, dim>& real,
                         GridHermitian<Real, dim>& spectral) {
   auto& plans = getPlans(make_key(real, spectral));
   fftw::execute(plans.first, real.getInternalData(),
                 reinterpret_cast<complex_t*>(spectral.getInternalData()));
 }
 
 template <UInt dim>
 void FFTEngine::backward(Grid<Real, dim>& real,
                          GridHermitian<Real, dim>& spectral) {
   auto& plans = getPlans(make_key(real, spectral));
   fftw::execute(plans.second,
                 reinterpret_cast<complex_t*>(spectral.getInternalData()),
                 real.getInternalData());
   // Normalize
   real *= (1. / real.getNbPoints());
 }
 
 template <typename T, UInt dim, bool hermitian>
 Grid<T, dim> FFTEngine::computeFrequencies(const std::array<UInt, dim>& 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);
   constexpr UInt dmax = dim - static_cast<UInt>(hermitian);
 
 #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()));
     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