diff --git a/src/core/fft_engine.cpp b/src/core/fft_engine.cpp
index f0bafa4..7073382 100644
--- a/src/core/fft_engine.cpp
+++ b/src/core/fft_engine.cpp
@@ -1,96 +1,89 @@
 /**
 
  *  @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::~FFTEngine() {
-  for (auto&& tup : plans) {
-    fftw::destroy(tup.second.first);
-    fftw::destroy(tup.second.second);
-  }
-}
-
 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};
 
   // Flags does not contain FFTW_ESTIMATE | FFTW_UNALIGNED => allocate in & out
-  if (not (flags & FFTW_ESTIMATE and flags & FFTW_UNALIGNED)) {
+  if (not(flags & FFTW_ESTIMATE and flags & FFTW_UNALIGNED)) {
     const auto size =
         std::accumulate(n.begin(), n.end(), 1, std::multiplies<void>()) *
         components;
     auto hermit_n = GridHermitian<Real, 1>::hermitianDimensions(n);
     const auto hsize = std::accumulate(hermit_n.begin(), hermit_n.end(), 1,
                                        std::multiplies<void>()) *
                        components;
     in.ptr = fftw::helper<Real>::alloc_real(size);
     out.ptr = fftw::helper<Real>::alloc_complex(hsize);
   }
 
-  auto forward =
+  fftw::plan<Real> forward{
       fftw::plan_many_forward(rank, n.data(), howmany, in, inembed, istride,
-                              idist, out, onembed, ostride, odist, flags);
-  auto backward =
+                              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(forward, backward);
+                               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 184b298..2ab66d8 100644
--- a/src/core/fft_engine.hh
+++ b/src/core/fft_engine.hh
@@ -1,149 +1,148 @@
 /**
  *  @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::helper<Real>::plan, fftw::helper<Real>::plan>;
+  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() {}
-  ~FFTEngine();
 
   /// 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);
 
 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
   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
diff --git a/src/core/fftw_interface.hh b/src/core/fftw_interface.hh
index 69b03ad..cc76c0e 100644
--- a/src/core/fftw_interface.hh
+++ b/src/core/fftw_interface.hh
@@ -1,176 +1,194 @@
 /**
  *  @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_INTERFACE
 #define FFTW_INTERFACE
 /* -------------------------------------------------------------------------- */
 #include <cstddef>
+#include <utility>
+
 #include <fftw3.h>
 
 namespace fftw {
 template <typename T>
 struct helper;
 
 template <>
 struct helper<double> {
   using complex = fftw_complex;
   using plan = fftw_plan;
 
   static auto alloc_real(std::size_t size) { return fftw_alloc_real(size); }
   static auto alloc_complex(std::size_t size) {
     return fftw_alloc_complex(size);
   }
 };
 
 #ifdef LONG_PRECISION
 template <>
 struct helper<long double> {
   using complex = fftwl_complex;
   using plan = fftwl_plan;
 
   static auto alloc_real(std::size_t size) { return fftwl_alloc_real(size); }
   static auto alloc_complex(std::size_t size) {
     return fftwl_alloc_complex(size);
   }
 };
 #endif
 
 template <typename T>
 inline auto free(T* ptr) {
   fftw_free(ptr);
 }
 
 inline auto destroy(fftw_plan plan) { fftw_destroy_plan(plan); }
 #ifdef LONG_PRECISION
 inline auto destroy(fftwl_plan plan) { fftwl_destroy_plan(plan); }
 #endif
 
+/// Holder type for fftw plans
 template <typename T>
 struct plan {
   typename helper<T>::plan _plan;
 
-  ~plan() { fftw::destroy(_plan); }
+  /// Create from plan
+  explicit plan(typename helper<T>::plan _plan = nullptr) : _plan(_plan) {}
+  /// Move constructor to avoid accidental plan destruction
+  plan(plan&& o) : _plan(std::exchange(o._plan, nullptr)) {}
+  /// Move operator
+  plan& operator=(plan&& o) {
+    _plan = std::exchange(o._plan, nullptr);
+    return *this;
+  }
+  /// Destroy plan
+  ~plan() {
+    if (_plan)
+      fftw::destroy(_plan);
+  }
+
+  /// For seamless use with fftw api
   operator typename helper<T>::plan() const { return _plan; }
 };
 
 /* -------------------------------------------------------------------------- */
 
 inline auto plan_many_forward(int rank, const int* n, int howmany, double* in,
                               const int* inembed, int istride, int idist,
                               fftw_complex* out, const int* onembed,
                               int ostride, int odist, unsigned flags) {
   return fftw_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist,
                                 out, onembed, ostride, odist, flags);
 }
 
 inline auto plan_many_backward(int rank, const int* n, int howmany,
                                fftw_complex* in, const int* inembed,
                                int istride, int idist, double* out,
                                const int* onembed, int ostride, int odist,
                                unsigned flags) {
   return fftw_plan_many_dft_c2r(rank, n, howmany, in, inembed, istride, idist,
                                 out, onembed, ostride, odist, flags);
 }
 
 inline auto plan_1d_forward(int n, double* in, fftw_complex* out,
                             unsigned flags) {
   return fftw_plan_dft_r2c_1d(n, in, out, flags);
 }
 
 inline auto plan_1d_backward(int n, fftw_complex* in, double* out,
                              unsigned flags) {
   return fftw_plan_dft_c2r_1d(n, in, out, flags);
 }
 
 inline auto plan_2d_forward(int n0, int n1, double* in, fftw_complex* out,
                             unsigned flags) {
   return fftw_plan_dft_r2c_2d(n0, n1, in, out, flags);
 }
 
 inline auto plan_2d_backward(int n0, int n1, fftw_complex* out, double* in,
                              unsigned flags) {
   return fftw_plan_dft_c2r_2d(n0, n1, out, in, flags);
 }
 
 inline auto execute(fftw_plan plan) { fftw_execute(plan); }
 inline auto execute(fftw_plan plan, double* in, fftw_complex* out) {
   fftw_execute_dft_r2c(plan, in, out);
 }
 
 inline auto execute(fftw_plan plan, fftw_complex* in, double* out) {
   fftw_execute_dft_c2r(plan, in, out);
 }
 
 /* -------------------------------------------------------------------------- */
 
 #ifdef LONG_PRECISION
 inline auto plan_many_forward(int rank, const int* n, int howmany,
                               long double* in, const int* inembed, int istride,
                               int idist, fftwl_complex* out, const int* onembed,
                               int ostride, int odist, unsigned flags) {
   return fftwl_plan_many_dft_r2c(rank, n, howmany, in, inembed, istride, idist,
                                  out, onembed, ostride, odist, flags);
 }
 
 inline auto plan_many_backward(int rank, const int* n, int howmany,
                                fftwl_complex* in, const int* inembed,
                                int istride, int idist, long double* out,
                                const int* onembed, int ostride, int odist,
                                unsigned flags) {
   return fftwl_plan_many_dft_c2r(rank, n, howmany, in, inembed, istride, idist,
                                  out, onembed, ostride, odist, flags);
 }
 
 inline auto plan_1d_forward(int n, long double* in, fftwl_complex* out,
                             unsigned flags) {
   return fftwl_plan_dft_r2c_1d(n, in, out, flags);
 }
 
 inline auto plan_1d_backward(int n, fftwl_complex* in, long double* out,
                              unsigned flags) {
   return fftwl_plan_dft_c2r_1d(n, in, out, flags);
 }
 
 inline auto plan_2d_forward(int n0, int n1, long double* in, fftwl_complex* out,
                             unsigned flags) {
   return fftwl_plan_dft_r2c_2d(n0, n1, in, out, flags);
 }
 
 inline auto plan_2d_backward(int n0, int n1, fftwl_complex* out,
                              long double* in, unsigned flags) {
   return fftwl_plan_dft_c2r_2d(n0, n1, out, in, flags);
 }
 
 inline auto execute(fftwl_plan plan) { fftwl_execute(plan); }
 inline auto execute(fftwl_plan plan, long double* in, fftwl_complex* out) {
   fftwl_execute_dft_r2c(plan, in, out);
 }
 
 inline auto execute(fftwl_plan plan, fftwl_complex* in, long double* out) {
   fftwl_execute_dft_c2r(plan, in, out);
 }
 
 #endif  // LONG_PRECISION
 
 }  // namespace fftw
 
 #endif  // FFTW_INTERFACE
diff --git a/src/core/grid_hermitian.hh b/src/core/grid_hermitian.hh
index 252e534..69854ec 100644
--- a/src/core/grid_hermitian.hh
+++ b/src/core/grid_hermitian.hh
@@ -1,168 +1,162 @@
 /**
  *  @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 __GRID_HERMITIAN_HH__
 #define __GRID_HERMITIAN_HH__
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include "tamaas.hh"
 #include <complex>
+#include <type_traits>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 static complex<Real> dummy;
 
 /**
  * @brief Multi-dimensional, multi-component herimitian array
  *
  * This class represents an array of hermitian data, meaning it has dimensions
  * of: n1 * n2 * n3 * ... * (nx / 2 + 1)
  *
  * However, it acts as a fully dimensioned array, returning a dummy reference
  * for data outside its real dimension, allowing one to write full (and
  * inefficient) loops without really worrying about the reduced dimension.
  *
  * It would however be better to just use the true dimensions of the surface
  * for all intents and purposes, as it saves computation time.
  */
 template <typename T, UInt dim>
 class GridHermitian : public Grid<complex<T>, dim> {
 public:
   using value_type = complex<T>;
   /* ------------------------------------------------------------------------ */
   /* Constructors */
   /* ------------------------------------------------------------------------ */
 public:
   GridHermitian() = default;
   GridHermitian(const GridHermitian& o) = default;
   GridHermitian(GridHermitian&& o) noexcept = default;
   using Grid<value_type, dim>::Grid;
   using Grid<complex<T>, dim>::operator=;
   /* ------------------------------------------------------------------------ */
   /* Access operators */
   /* ------------------------------------------------------------------------ */
   template <typename... T1>
   inline complex<T>& operator()(T1... args);
   template <typename... T1>
   inline const complex<T>& operator()(T1... args) const;
 
-  static std::array<UInt, dim>
-  hermitianDimensions(std::array<UInt, dim> n) {
+  static std::array<UInt, dim> hermitianDimensions(std::array<UInt, dim> n) {
     n[dim - 1] /= 2;
     n[dim - 1] += 1;
     return n;
   }
 
-  static std::vector<UInt>
-  hermitianDimensions(std::vector<UInt> n) {
-    n.back() /= 2;
-    n.back() += 1;
-    return n;
-  }
-
-  static std::vector<Int>
-  hermitianDimensions(std::vector<Int> n) {
+  template <typename Int,
+            typename = std::enable_if_t<std::is_arithmetic<Int>::value>>
+  static std::vector<Int> hermitianDimensions(std::vector<Int> n) {
     n.back() /= 2;
     n.back() += 1;
     return n;
   }
 
 private:
   template <typename... T1>
   inline void packTuple(UInt* t, UInt i, T1... args) const;
   template <typename... T1>
   inline void packTuple(UInt* t, UInt i) const;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline function definitions                                                */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline void GridHermitian<T, dim>::packTuple(UInt* t, UInt i,
                                              T1... args) const {
   *t = i;
   packTuple(t + 1, args...);
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline void GridHermitian<T, dim>::packTuple(UInt* t, UInt i) const {
   *t = i;
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline complex<T>& GridHermitian<T, dim>::operator()(T1... args) {
   constexpr UInt nargs = sizeof...(T1);
   static_assert(nargs <= dim + 1, "Too many arguments");
   UInt last_index = std::get<nargs - 1>(std::forward_as_tuple(args...));
   if (nargs != 1 && last_index >= this->n[dim - 1]) {
     std::array<UInt, dim + 1> tuple = {{0}};
     packTuple(tuple.data(), args...);
 
     for (UInt i = 0; i < dim; i++) {
       if (tuple[i] && i != dim - 1)
         tuple[i] = this->n[i] - tuple[i];
       else if (tuple[i] && i == dim - 1)
         tuple[i] = 2 * (this->n[i] - 1) - tuple[i];
     }
 
     dummy = conj(this->Grid<complex<T>, dim>::operator()(tuple));
     return dummy;
   }
 
   else
     return Grid<complex<T>, dim>::operator()(args...);
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const complex<T>& GridHermitian<T, dim>::operator()(T1... args) const {
   constexpr UInt nargs = sizeof...(T1);
   static_assert(nargs <= dim + 1, "Too many arguments");
   UInt last_index = std::get<nargs - 1>(std::forward_as_tuple(args...));
   if (nargs != 1 && last_index >= this->n[dim - 1]) {
     std::array<UInt, dim + 1> tuple = {{0}};
     packTuple(tuple.data(), args...);
 
     for (UInt i = 0; i < dim; i++) {
       if (tuple[i] && i != dim - 1)
         tuple[i] = this->n[i] - tuple[i];
       else if (tuple[i] && i == dim - 1)
         tuple[i] = 2 * (this->n[i] - 1) - tuple[i];
     }
 
     dummy = conj(this->Grid<complex<T>, dim>::operator()(tuple));
     return dummy;
   }
 
   else
     return Grid<complex<T>, dim>::operator()(args...);
 }
 
 }  // namespace tamaas
 
 #endif  // __GRID_HERMITIAN_HH__