diff --git a/CHANGELOG.md b/CHANGELOG.md
index efbeae3..23dabbc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,8 +2,6 @@
 
 ## [Unreleased]
 
-## [Version 2.0.0-alpha1] - 2019-04-XY
-
 We have introduced a number of breaking changes, motivated by the need to
 modernize the library.
 
@@ -14,6 +12,10 @@ modernize the library.
 - **BREAKING** We require a fully compliant C++14 compiler.
 - **BREAKING** We require CMake 3.9
 - Project dependencies are now managed as a superbuild.
+- The implementation of the `RadialSolution` for the second order ODE
+  associated with the spherical diffuse Green's function is less heavily
+  `template`-d.
+- The dependency on Boost.Odeint has been dropped.
 
 ### Removed
 
diff --git a/src/green/InterfacesImpl.hpp b/src/green/InterfacesImpl.hpp
new file mode 100644
index 0000000..94b614d
--- /dev/null
+++ b/src/green/InterfacesImpl.hpp
@@ -0,0 +1,277 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2019 Roberto Di Remigio, Luca Frediani and contributors.
+ *
+ * This file is part of PCMSolver.
+ *
+ * PCMSolver is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PCMSolver 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to the
+ * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
+ */
+
+#pragma once
+
+#include <array>
+#include <cmath>
+#include <fstream>
+#include <functional>
+#include <iomanip>
+#include <tuple>
+#include <vector>
+
+#include <Eigen/Core>
+
+#include "utils/MathUtils.hpp"
+#include "utils/RungeKutta.hpp"
+#include "utils/SplineFunction.hpp"
+
+/*! \file InterfacesImpl.hpp */
+
+namespace pcm {
+namespace green {
+namespace detail {
+/*! \brief Abstract class for a system of ordinary differential equations
+ *  \tparam Order The order of the ordinary differential equation
+ *  \author Roberto Di Remigio
+ *  \date 2018
+ */
+template <size_t Order = 1> class ODESystem {
+public:
+  typedef std::array<double, Order> StateType;
+  size_t ODEorder() const { return Order; }
+  void operator()(const StateType & f, StateType & dfdx, const double t) const {
+    RHS(f, dfdx, t);
+  }
+  virtual ~ODESystem() {}
+
+private:
+  virtual void RHS(const StateType & f, StateType & dfdx, const double t) const = 0;
+};
+
+/*! \typedef ProfileEvaluator
+ *  \brief sort of a function pointer to the dielectric profile evaluation function
+ */
+typedef std::function<std::tuple<double, double>(const double)> ProfileEvaluator;
+
+/*! \class LnTransformedRadial
+ *  \brief system of ln-transformed first-order radial differential equations
+ *  \author Roberto Di Remigio
+ *  \date 2015
+ *
+ *  Provides a handle to the system of differential equations for the integrator.
+ *  The dielectric profile comes in as a boost::function object.
+ */
+class LnTransformedRadial final : public pcm::green::detail::ODESystem<2> {
+public:
+  /*! Type of the state vector of the ODE */
+  typedef pcm::green::detail::ODESystem<2>::StateType StateType;
+  /*! Constructor from profile evaluator and angular momentum */
+  LnTransformedRadial(const ProfileEvaluator & e, int lval) : eval_(e), l_(lval) {}
+
+private:
+  /*! Dielectric profile function and derivative evaluation */
+  const ProfileEvaluator eval_;
+  /*! Angular momentum */
+  const int l_;
+  /*! \brief Provides a functor for the evaluation of the system of first-order ODEs.
+   *  \param[in] rho state vector holding the function and its first derivative
+   *  \param[out] drhodr state vector holding the first and second derivative
+   *  \param[in] y logarithmic position on the integration grid
+   *
+   *  The second-order ODE and the system of first-order ODEs
+   *  are reported in the manuscript.
+   */
+  virtual void RHS(const StateType & rho, StateType & drhodr, const double y) const {
+    // Evaluate the dielectric profile
+    double eps = 0.0, epsPrime = 0.0;
+    std::tie(eps, epsPrime) = eval_(std::exp(y));
+    if (utils::numericalZero(eps))
+      PCMSOLVER_ERROR("Division by zero!");
+    double gamma_epsilon = std::exp(y) * epsPrime / eps;
+    // System of equations is defined here
+    drhodr[0] = rho[1];
+    drhodr[1] = -rho[1] * (rho[1] + 1.0 + gamma_epsilon) + l_ * (l_ + 1);
+  }
+};
+} // namespace detail
+
+using detail::ProfileEvaluator;
+
+/*! \class RadialFunction
+ *  \brief represents solutions to the radial 2nd order ODE
+ *  \author Roberto Di Remigio
+ *  \date 2018
+ *  \tparam ODE system of 1st order ODEs replacing the 2nd order ODE
+ *
+ * The sign of the integration interval, i.e. the sign of the difference
+ * between the minimum and the maximum of the interval, is used to discriminate
+ * between the two independent solutions (zeta-type and omega-type) of the
+ * differential equation:
+ *   - When the sign is positive (y_min_ < y_max_), we are integrating **forwards**
+ * for the zeta-type solution.
+ *   - When the sign is negative (y_max_ < y_min_), we are integrating **backwards**
+ * for the omega-type solution.
+ *
+ * This is based on an idea Luca had in Wuhan/Chongqing for the planar
+ * diffuse interface.
+ * With respect to the previous, much more heavily template-d implementation,
+ * it introduces a bit more if-else if-else branching in the function_impl and
+ * derivative_impl functions.
+ * This is largely offset by the better readability and reduced code duplication.
+ */
+template <typename ODE> class RadialFunction final {
+public:
+  RadialFunction() : L_(0), y_min_(0.0), y_max_(0.0), y_sign_(1) {}
+  RadialFunction(int l,
+                 double ymin,
+                 double ymax,
+                 double ystep,
+                 const ProfileEvaluator & eval)
+      : L_(l),
+        y_min_(ymin),
+        y_max_(ymax),
+        y_sign_(pcm::utils::sign(y_max_ - y_min_)) {
+    compute(ystep, eval);
+  }
+  std::tuple<double, double> operator()(double point) const {
+    return std::make_tuple(function_impl(point), derivative_impl(point));
+  }
+  friend std::ostream & operator<<(std::ostream & os, RadialFunction & obj) {
+    for (size_t i = 0; i < obj.function_[0].size(); ++i) {
+      // clang-format off
+      os << std::fixed << std::left << std::setprecision(14)
+         << obj.function_[0][i] << "    "
+         << obj.function_[1][i] << "    "
+         << obj.function_[2][i] << std::endl;
+      // clang-format on
+    }
+    return os;
+  }
+
+private:
+  typedef typename ODE::StateType StateType;
+  /*! Angular momentum of the solution */
+  int L_;
+  /*! Lower bound of the integration interval */
+  double y_min_;
+  /*! Upper bound of the integration interval */
+  double y_max_;
+  /*! The sign of the integration interval */
+  double y_sign_;
+  /*! The actual data: grid, function value and first derivative values */
+  std::array<std::vector<double>, 3> function_;
+  /*! Reports progress of differential equation integrator */
+  void push_back(const StateType & x, double y) {
+    function_[0].push_back(y);
+    function_[1].push_back(x[0]);
+    function_[2].push_back(x[1]);
+  }
+  /*! \brief Calculates radial solution
+   *  \param[in] step ODE integrator step
+   *  \param[in] eval dielectric profile evaluator function object
+   *  \return the number of integration steps
+   *
+   *  This function discriminates between the first (zeta-type), i.e. the one
+   *  with r^l behavior, and the second (omega-type) radial solution, i.e. the
+   *  one with r^(-l-1) behavior, based on the sign of the integration interval
+   *  y_sign_.
+   */
+  size_t compute(const double step, const ProfileEvaluator & eval) {
+    ODE system(eval, L_);
+    // Set initial conditions
+    StateType init;
+    if (y_sign_ > 0.0) { // zeta-type solution
+      init[0] = y_sign_ * L_ * y_min_;
+      init[1] = y_sign_ * L_;
+    } else { // omega-type solution
+      init[0] = y_sign_ * (L_ + 1) * y_min_;
+      init[1] = y_sign_ * (L_ + 1);
+    }
+    pcm::utils::RungeKutta4<StateType> stepper;
+    size_t nSteps =
+        pcm::utils::integrate_const(stepper,
+                                    system,
+                                    init,
+                                    y_min_,
+                                    y_max_,
+                                    y_sign_ * step,
+                                    std::bind(&RadialFunction<ODE>::push_back,
+                                              this,
+                                              std::placeholders::_1,
+                                              std::placeholders::_2));
+
+    // clang-format off
+    // Reverse order of function_ if omega-type solution was computed
+    // this ensures that they are in ascending order, as later expected by
+    // function_impl and derivative_impl
+    if (y_sign_ < 0.0) {
+      for (auto & comp : function_) {
+        std::reverse(comp.begin(), comp.end());
+      }
+      // clang-format on
+    }
+    return nSteps;
+  }
+  /*! \brief Returns value of function at given point
+   *  \param[in] point evaluation point
+   *
+   *  We first check if point is below y_min_, if yes we use
+   *  the asymptotic form.
+   */
+  double function_impl(double point) const {
+    double val = 0.0;
+    if (point <= y_min_ && y_sign_ > 0.0) { // Asymptotic zeta-type solution
+      val = y_sign_ * L_ * point;
+    } else if (point >= y_min_ && y_sign_ < 0.0) { // Asymptotic omega-type solution
+      val = y_sign_ * (L_ + 1) * point;
+    } else {
+      val = utils::splineInterpolation(point, function_[0], function_[1]);
+    }
+    return val;
+  }
+  /*! \brief Returns value of 1st derivative of function at given point
+   *  \param[in] point evaluation point
+   *
+   *  Below y_min_, the asymptotic form is used. Otherwise we interpolate.
+   */
+  double derivative_impl(double point) const {
+    double val = 0.0;
+    if (point <= y_min_ && y_sign_ > 0.0) { // Asymptotic zeta-type solution
+      val = y_sign_ * L_;
+    } else if (point >= y_min_ && y_sign_ < 0.0) { // Asymptotic omega-type solution
+      val = y_sign_ * (L_ + 1);
+    } else {
+      val = utils::splineInterpolation(point, function_[0], function_[2]);
+    }
+    return val;
+  }
+};
+
+/*! \brief Write contents of a RadialFunction to file
+ *  \param[in] f solution to the radial 2nd order ODE to print
+ *  \param[in] fname name of the file
+ *  \author Roberto Di Remigio
+ *  \date 2018
+ *  \tparam ODE system of 1st order ODEs replacing the 2nd order ODE
+ */
+template <typename ODE>
+void writeToFile(RadialFunction<ODE> & f, const std::string & fname) {
+  std::ofstream fout;
+  fout.open(fname.c_str());
+  fout << f << std::endl;
+  fout.close();
+}
+} // namespace green
+} // namespace pcm
diff --git a/src/green/SphericalDiffuse.cpp b/src/green/SphericalDiffuse.cpp
new file mode 100644
index 0000000..805fa99
--- /dev/null
+++ b/src/green/SphericalDiffuse.cpp
@@ -0,0 +1,404 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2019 Roberto Di Remigio, Luca Frediani and contributors.
+ *
+ * This file is part of PCMSolver.
+ *
+ * PCMSolver is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PCMSolver 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to the
+ * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
+ */
+
+#include "SphericalDiffuse.hpp"
+
+#include <iosfwd>
+
+#include <Eigen/Core>
+
+// Has to be included here
+#include "InterfacesImpl.hpp"
+
+#include "GreenData.hpp"
+#include "GreensFunction.hpp"
+#include "cavity/Element.hpp"
+#include "utils/Legendre.hpp"
+#include "utils/MathUtils.hpp"
+
+namespace pcm {
+namespace green {
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::coefficientCoulomb(
+    const Eigen::Vector3d & source,
+    const Eigen::Vector3d & probe) const {
+  // Obtain coefficient for the separation of the Coulomb singularity
+  return this->coefficient_impl(source, probe);
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::Coulomb(
+    const Eigen::Vector3d & source,
+    const Eigen::Vector3d & probe) const {
+  double r12 = (source - probe).norm();
+
+  // Obtain coefficient for the separation of the Coulomb singularity
+  return (1.0 / (this->coefficient_impl(source, probe) * r12));
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::imagePotential(
+    const Eigen::Vector3d & source,
+    const Eigen::Vector3d & probe) const {
+  // Obtain coefficient for the separation of the Coulomb singularity
+  double Cr12 = this->coefficient_impl(source, probe);
+
+  double gr12 = 0.0;
+  for (int L = 1; L <= maxLGreen_; ++L) {
+    gr12 += this->imagePotentialComponent_impl(L, source, probe, Cr12);
+  }
+
+  return gr12;
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::coefficientCoulombDerivative(
+    const Eigen::Vector3d & direction,
+    const Eigen::Vector3d & p1,
+    const Eigen::Vector3d & p2) const {
+  return threePointStencil(
+      std::bind(&SphericalDiffuse<ProfilePolicy>::coefficientCoulomb,
+                this,
+                std::placeholders::_1,
+                std::placeholders::_2),
+      p2,
+      p1,
+      direction,
+      this->delta_);
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::CoulombDerivative(
+    const Eigen::Vector3d & direction,
+    const Eigen::Vector3d & p1,
+    const Eigen::Vector3d & p2) const {
+  return threePointStencil(std::bind(&SphericalDiffuse<ProfilePolicy>::Coulomb,
+                                     this,
+                                     std::placeholders::_1,
+                                     std::placeholders::_2),
+                           p2,
+                           p1,
+                           direction,
+                           this->delta_);
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::imagePotentialDerivative(
+    const Eigen::Vector3d & direction,
+    const Eigen::Vector3d & p1,
+    const Eigen::Vector3d & p2) const {
+  return threePointStencil(
+      std::bind(&SphericalDiffuse<ProfilePolicy>::imagePotential,
+                this,
+                std::placeholders::_1,
+                std::placeholders::_2),
+      p2,
+      p1,
+      direction,
+      this->delta_);
+}
+
+template <typename ProfilePolicy>
+std::tuple<double, double> SphericalDiffuse<ProfilePolicy>::epsilon(
+    const Eigen::Vector3d & point) const {
+  return this->profile_((point + this->origin_).norm());
+}
+
+template <typename ProfilePolicy>
+void SphericalDiffuse<ProfilePolicy>::toFile(const std::string & prefix) {
+  std::string tmp;
+  prefix.empty() ? tmp = prefix : tmp = prefix + "-";
+  writeToFile(zetaC_, tmp + "zetaC.dat");
+  writeToFile(omegaC_, tmp + "omegaC.dat");
+  for (int L = 1; L <= maxLGreen_; ++L) {
+    std::ostringstream Llabel;
+    Llabel << std::setw(4) << std::setfill('0') << L;
+    writeToFile(zeta_[L], tmp + "zeta_" + Llabel.str() + ".dat");
+    writeToFile(omega_[L], tmp + "omega_" + Llabel.str() + ".dat");
+  }
+}
+
+template <typename ProfilePolicy>
+Stencil SphericalDiffuse<ProfilePolicy>::operator()(Stencil * sp,
+                                                    Stencil * pp) const {
+  // Transfer raw arrays to Eigen vectors using the Map type
+  Eigen::Map<Eigen::Matrix<double, 3, 1>> source(sp), probe(pp);
+
+  // Obtain coefficient for the separation of the Coulomb singularity
+  double Cr12 = this->coefficient_impl(source, probe);
+
+  double gr12 = 0.0;
+  for (int L = 0; L <= maxLGreen_; ++L) {
+    gr12 += this->imagePotentialComponent_impl(L, source, probe, Cr12);
+  }
+  double r12 = (source - probe).norm();
+
+  double gf = (1.0 / (Cr12 * r12) + gr12);
+  return gf;
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::kernelD_impl(
+    const Eigen::Vector3d & direction,
+    const Eigen::Vector3d & p1,
+    const Eigen::Vector3d & p2) const {
+  double eps_r2 = 0.0;
+  // Shift p2 by origin_
+  std::tie(eps_r2, std::ignore) = this->epsilon(p2);
+
+  return (eps_r2 * this->derivativeProbe(direction, p1, p2));
+}
+
+template <typename ProfilePolicy>
+KernelS SphericalDiffuse<ProfilePolicy>::exportKernelS_impl() const {
+  return std::bind(&SphericalDiffuse<ProfilePolicy>::kernelS,
+                   *this,
+                   std::placeholders::_1,
+                   std::placeholders::_2);
+}
+
+template <typename ProfilePolicy>
+KernelD SphericalDiffuse<ProfilePolicy>::exportKernelD_impl() const {
+  return std::bind(&SphericalDiffuse<ProfilePolicy>::kernelD,
+                   *this,
+                   std::placeholders::_1,
+                   std::placeholders::_2,
+                   std::placeholders::_3);
+}
+
+template <typename DerivativeTraits>
+DerivativeProbe SphericalDiffuse<DerivativeTraits>::exportDerivativeProbe_impl()
+    const {
+  return std::bind(&SphericalDiffuse<DerivativeTraits>::derivativeProbe,
+                   *this,
+                   std::placeholders::_1,
+                   std::placeholders::_2,
+                   std::placeholders::_3);
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::singleLayer_impl(const Element & e,
+                                                         double factor) const {
+  // Diagonal of S inside the cavity
+  double Sii_I = detail::diagonalSi(e.area(), factor);
+  // "Diagonal" of Coulomb singularity separation coefficient
+  double coulomb_coeff = this->coefficientCoulomb(e.center(), e.center());
+  // "Diagonal" of the image Green's function
+  double image = this->imagePotential(e.center(), e.center());
+  return (Sii_I / coulomb_coeff + image);
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::doubleLayer_impl(const Element & e,
+                                                         double factor) const {
+  // Diagonal of S inside the cavity
+  double Sii_I = detail::diagonalSi(e.area(), factor);
+  // Diagonal of D inside the cavity
+  double Dii_I = detail::diagonalDi(e.area(), e.sphere().radius, factor);
+  // "Diagonal" of Coulomb singularity separation coefficient
+  double coulomb_coeff = this->coefficientCoulomb(e.center(), e.center());
+  // "Diagonal" of the directional derivative of the Coulomb singularity
+  // separation coefficient
+  double coeff_grad =
+      this->coefficientCoulombDerivative(e.normal(), e.center(), e.center()) /
+      std::pow(coulomb_coeff, 2);
+  // "Diagonal" of the directional derivative of the image Green's function
+  double image_grad =
+      this->imagePotentialDerivative(e.normal(), e.center(), e.center());
+
+  double eps_r2 = 0.0;
+  std::tie(eps_r2, std::ignore) = this->epsilon(e.center());
+
+  return (eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad));
+}
+
+template <typename ProfilePolicy>
+std::ostream & SphericalDiffuse<ProfilePolicy>::printObject(std::ostream & os) {
+  Eigen::IOFormat CleanFmt(Eigen::StreamPrecision, 0, ", ", "\n", "(", ")");
+  os << "Green's function type: spherical diffuse" << std::endl;
+  os << this->profile_ << std::endl;
+  os << "Sphere center        = " << this->origin_.transpose().format(CleanFmt)
+     << std::endl;
+  os << "Angular momentum (Green's function)    = " << this->maxLGreen_ << std::endl;
+  os << "Angular momentum (Coulomb coefficient) = " << this->maxLC_;
+  return os;
+}
+
+template <typename ProfilePolicy>
+void SphericalDiffuse<ProfilePolicy>::initSphericalDiffuse() {
+  using namespace detail;
+
+  // Parameters for the numerical solution of the radial differential equation
+  /*! Lower bound of the integration interval */
+  double r_0_ = 0.1;
+  double y_0_ = std::log(r_0_);
+  /*! Upper bound of the integration interval */
+  double r_infinity_ = this->profile_.upperLimit() + 30.0;
+  double y_infinity_ = std::log(r_infinity_);
+  double relative_width = this->profile_.relativeWidth();
+  /*! Time step between observer calls */
+  double step_ = 1.0e-2 * relative_width;
+
+  ProfileEvaluator eval_ =
+      std::bind(&ProfilePolicy::operator(), this->profile_, std::placeholders::_1);
+
+  zetaC_ =
+      RadialFunction<LnTransformedRadial>(maxLC_, y_0_, y_infinity_, step_, eval_);
+  omegaC_ =
+      RadialFunction<LnTransformedRadial>(maxLC_, y_infinity_, y_0_, step_, eval_);
+  zeta_.reserve(maxLGreen_ + 1);
+  omega_.reserve(maxLGreen_ + 1);
+  for (int L = 0; L <= maxLGreen_; ++L) {
+    zeta_.push_back(
+        RadialFunction<LnTransformedRadial>(L, y_0_, y_infinity_, step_, eval_));
+    omega_.push_back(
+        RadialFunction<LnTransformedRadial>(L, y_infinity_, y_0_, step_, eval_));
+  }
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::imagePotentialComponent_impl(
+    int L,
+    const Eigen::Vector3d & sp,
+    const Eigen::Vector3d & pp,
+    double Cr12) const {
+  Eigen::Vector3d sp_shift = sp + this->origin_;
+  Eigen::Vector3d pp_shift = pp + this->origin_;
+  double r1 = sp_shift.norm();
+  double r2 = pp_shift.norm();
+  double cos_gamma = sp_shift.dot(pp_shift) / (r1 * r2);
+  // Evaluate Legendre polynomial of order L
+  // First of all clean-up cos_gamma, Legendre polynomials
+  // are only defined for -1 <= x <= 1
+  if (utils::numericalZero(cos_gamma - 1))
+    cos_gamma = 1.0;
+  if (utils::numericalZero(cos_gamma + 1))
+    cos_gamma = -1.0;
+  double pl_x = Legendre::Pn<double>(L, cos_gamma);
+
+  /* Zeta and Omega are now on a logarithmic scale We need to pass the
+         arguments in the correct scale and then use the chain rule to get
+         the proper derivative */
+  double y1 = std::log(r1);
+  double y2 = std::log(r2);
+
+  /* Sample zeta_[L] */
+  double zeta1 = 0.0, zeta2 = 0.0, d_zeta2 = 0.0;
+  /* Value of zeta_[L] at point with index 1 */
+  std::tie(zeta1, std::ignore) = zeta_[L](y1);
+  /* Value of zeta_[L] and its first derivative at point with index 2 */
+  std::tie(zeta2, d_zeta2) = zeta_[L](y2);
+
+  /* Sample omega_[L] */
+  double omega1 = 0.0, omega2 = 0.0, d_omega2 = 0.0;
+  /* Value of omega_[L] at point with index 1 */
+  std::tie(omega1, std::ignore) = omega_[L](y1);
+  /* Value of omega_[L] and its first derivative at point with index 2 */
+  std::tie(omega2, d_omega2) = omega_[L](y2);
+
+  double eps_r2 = 0.0;
+  std::tie(eps_r2, std::ignore) = this->profile_(pp_shift.norm());
+
+  /* Evaluation of the Wronskian and the denominator */
+  double denominator = (d_zeta2 - d_omega2) * r2 * eps_r2;
+
+  double gr12 = 0.0;
+  if (r1 < r2) {
+    gr12 = std::exp(zeta1 - zeta2) * (2 * L + 1) / denominator;
+    double f_L = r1 / r2;
+    for (int i = 1; i < L; ++i) {
+      f_L *= r1 / r2;
+    }
+    gr12 = (gr12 - f_L / (r2 * Cr12)) * pl_x;
+  } else {
+    gr12 = std::exp(omega1 - omega2) * (2 * L + 1) / denominator;
+    double f_L = r2 / r1;
+    for (int i = 1; i < L; ++i) {
+      f_L *= r2 / r1;
+    }
+    gr12 = (gr12 - f_L / (r1 * Cr12)) * pl_x;
+  }
+  return gr12;
+}
+
+template <typename ProfilePolicy>
+double SphericalDiffuse<ProfilePolicy>::coefficient_impl(
+    const Eigen::Vector3d & sp,
+    const Eigen::Vector3d & pp) const {
+  double r1 = (sp + this->origin_).norm();
+  double r2 = (pp + this->origin_).norm();
+
+  double y1 = std::log(r1);
+  double y2 = std::log(r2);
+
+  /* Sample zetaC_ */
+  double zeta1 = 0.0, zeta2 = 0.0, d_zeta2 = 0.0;
+  /* Value of zetaC_ at point with index 1 */
+  std::tie(zeta1, std::ignore) = zetaC_(y1);
+  /* Value of zetaC_ and its first derivative at point with index 2 */
+  std::tie(zeta2, d_zeta2) = zetaC_(y2);
+
+  /* Sample omegaC_ */
+  double omega1 = 0.0, omega2 = 0.0, d_omega2 = 0.0;
+  /* Value of omegaC_ at point with index 1 */
+  std::tie(omega1, std::ignore) = omegaC_(y1);
+  /* Value of omegaC_ and its first derivative at point with index 2 */
+  std::tie(omega2, d_omega2) = omegaC_(y2);
+
+  double tmp = 0.0, coeff = 0.0;
+  double eps_r2 = 0.0;
+  std::tie(eps_r2, std::ignore) = this->profile_(r2);
+
+  /* Evaluation of the Wronskian and the denominator */
+  double denominator = (d_zeta2 - d_omega2) * r2 * eps_r2;
+
+  if (r1 < r2) {
+    double f_L = r1 / r2;
+    for (int i = 1; i < maxLC_; ++i) {
+      f_L *= r1 / r2;
+    }
+    tmp = std::exp(zeta1 - zeta2) * (2 * maxLC_ + 1) / denominator;
+    coeff = f_L / (tmp * r2);
+  } else {
+    double f_L = r2 / r1;
+    for (int i = 1; i < maxLC_; ++i) {
+      f_L *= r2 / r1;
+    }
+    tmp = std::exp(omega1 - omega2) * (2 * maxLC_ + 1) / denominator;
+    coeff = f_L / (tmp * r1);
+  }
+
+  return coeff;
+}
+
+using dielectric_profile::OneLayerTanh;
+template class SphericalDiffuse<OneLayerTanh>;
+
+using dielectric_profile::OneLayerErf;
+template class SphericalDiffuse<OneLayerErf>;
+
+using dielectric_profile::OneLayerLog;
+template class SphericalDiffuse<OneLayerLog>;
+
+} // namespace green
+} // namespace pcm
diff --git a/src/green/SphericalDiffuse.hpp b/src/green/SphericalDiffuse.hpp
new file mode 100644
index 0000000..0567328
--- /dev/null
+++ b/src/green/SphericalDiffuse.hpp
@@ -0,0 +1,272 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2019 Roberto Di Remigio, Luca Frediani and contributors.
+ *
+ * This file is part of PCMSolver.
+ *
+ * PCMSolver is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PCMSolver 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to the
+ * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
+ */
+
+#pragma once
+
+#include <iosfwd>
+
+#include <Eigen/Core>
+
+/*! \file SphericalDiffuse.hpp */
+
+namespace pcm {
+namespace cavity {
+class Element;
+} // namespace cavity
+} // namespace pcm
+
+#include "GreenData.hpp"
+#include "GreensFunction.hpp"
+#include "InterfacesImpl.hpp"
+#include "dielectric_profile/OneLayerErf.hpp"
+#include "dielectric_profile/OneLayerLog.hpp"
+#include "dielectric_profile/OneLayerTanh.hpp"
+
+namespace pcm {
+namespace green {
+/*! \class SphericalDiffuse
+ *  \brief Green's function for a diffuse interface with spherical symmetry
+ *  \author Hui Cao, Ville Weijo, Luca Frediani and Roberto Di Remigio
+ *  \date 2010-2015
+ *  \tparam ProfilePolicy functional form of the diffuse layer
+ *
+ *  The origin of the dielectric sphere can be changed by means of the constructor.
+ *  The solution of the differential equation defining the Green's function is
+ *  **always** performed assuming that the dielectric sphere is centered in the
+ * origin of the
+ *  coordinate system. Whenever the public methods are invoked to "sample" the
+ * Green's function
+ *  at a pair of points, a translation of the sampling points is performed first.
+ */
+template <typename ProfilePolicy = dielectric_profile::OneLayerLog>
+class SphericalDiffuse final : public GreensFunction<Stencil, ProfilePolicy> {
+public:
+  /*! Constructor for a one-layer interface
+   * \param[in] e1 left-side dielectric constant
+   * \param[in] e2 right-side dielectric constant
+   * \param[in] w width of the interface layer
+   * \param[in] c center of the diffuse layer
+   * \param[in] o center of the sphere
+   * \param[in] l maximum value of angular momentum
+   */
+  SphericalDiffuse(double e1,
+                   double e2,
+                   double w,
+                   double c,
+                   const Eigen::Vector3d & o,
+                   int l)
+      : GreensFunction<Stencil, ProfilePolicy>(ProfilePolicy(e1, e2, w, c)),
+        origin_(o),
+        maxLGreen_(l),
+        maxLC_(2 * l) {
+    initSphericalDiffuse();
+  }
+
+  virtual double permittivity() const override final {
+    PCMSOLVER_ERROR("permittivity() only implemented for uniform dielectrics");
+  }
+
+  /*! \brief Returns Coulomb singularity separation coefficient
+   *  \param[in] source location of the source charge
+   *  \param[in] probe location of the probe charge
+   */
+  double coefficientCoulomb(const Eigen::Vector3d & source,
+                            const Eigen::Vector3d & probe) const;
+
+  /*! \brief Returns singular part of the Green's function
+   *  \param[in] source location of the source charge
+   *  \param[in] probe location of the probe charge
+   */
+  double Coulomb(const Eigen::Vector3d & source,
+                 const Eigen::Vector3d & probe) const;
+
+  /*! \brief Returns non-singular part of the Green's function (image potential)
+   *  \param[in] source location of the source charge
+   *  \param[in] probe location of the probe charge
+   */
+  double imagePotential(const Eigen::Vector3d & source,
+                        const Eigen::Vector3d & probe) const;
+
+  /*! Returns value of the directional derivative of the
+   *  Coulomb singularity separation coefficient for the pair of points p1, p2:
+   *  \f$ \nabla_{\mathbf{p_2}}G(\mathbf{p}_1, \mathbf{p}_2)\cdot
+   *\mathbf{n}_{\mathbf{p}_2}\f$
+   *  Notice that this method returns the directional derivative with respect
+   *  to the probe point, thus assuming that the direction is relative to that
+   *point.
+   *
+   *  \param[in] direction the direction
+   *  \param[in]        p1 first point
+   *  \param[in]        p2 second point
+   */
+  double coefficientCoulombDerivative(const Eigen::Vector3d & direction,
+                                      const Eigen::Vector3d & p1,
+                                      const Eigen::Vector3d & p2) const;
+
+  /*! Returns value of the directional derivative of the
+   *  singular part of the Greens's function for the pair of points p1, p2:
+   *  \f$ \nabla_{\mathbf{p_2}}G(\mathbf{p}_1, \mathbf{p}_2)\cdot
+   *\mathbf{n}_{\mathbf{p}_2}\f$
+   *  Notice that this method returns the directional derivative with respect
+   *  to the probe point, thus assuming that the direction is relative to that
+   *point.
+   *
+   *  \param[in] direction the direction
+   *  \param[in]        p1 first point
+   *  \param[in]        p2 second point
+   */
+  double CoulombDerivative(const Eigen::Vector3d & direction,
+                           const Eigen::Vector3d & p1,
+                           const Eigen::Vector3d & p2) const;
+
+  /*! Returns value of the directional derivative of the
+   *  non-singular part (image potential) of the Greens's function for the pair of
+   *points p1, p2:
+   *  \f$ \nabla_{\mathbf{p_2}}G(\mathbf{p}_1, \mathbf{p}_2)\cdot
+   *\mathbf{n}_{\mathbf{p}_2}\f$
+   *  Notice that this method returns the directional derivative with respect
+   *  to the probe point, thus assuming that the direction is relative to that
+   *point.
+   *
+   *  \param[in] direction the direction
+   *  \param[in]        p1 first point
+   *  \param[in]        p2 second point
+   */
+  double imagePotentialDerivative(const Eigen::Vector3d & direction,
+                                  const Eigen::Vector3d & p1,
+                                  const Eigen::Vector3d & p2) const;
+
+  /*! Handle to the dielectric profile evaluation */
+  std::tuple<double, double> epsilon(const Eigen::Vector3d & point) const;
+  void toFile(const std::string & prefix = "");
+  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+private:
+  /*! Evaluates the Green's function given a pair of points
+   *  \param[in] sp the source point
+   *  \param[in] pp the probe point
+   *
+   *  \note This takes care of the origin shift
+   */
+  virtual Stencil operator()(Stencil * sp, Stencil * pp) const override;
+
+  /*! Returns value of the kernel of the \f$\mathcal{D}\f$ integral operator for the
+   * pair of points p1, p2:
+   *  \f$ [\boldsymbol{\varepsilon}\nabla_{\mathbf{p_2}}G(\mathbf{p}_1,
+   * \mathbf{p}_2)]\cdot \mathbf{n}_{\mathbf{p}_2}\f$
+   *  To obtain the kernel of the \f$\mathcal{D}^\dagger\f$ operator call this
+   * methods with \f$\mathbf{p}_1\f$
+   *  and \f$\mathbf{p}_2\f$ exchanged and with \f$\mathbf{n}_{\mathbf{p}_2} =
+   * \mathbf{n}_{\mathbf{p}_1}\f$
+   *  \param[in] direction the direction
+   *  \param[in]        p1 first point
+   *  \param[in]        p2 second point
+   */
+  virtual double kernelD_impl(const Eigen::Vector3d & direction,
+                              const Eigen::Vector3d & p1,
+                              const Eigen::Vector3d & p2) const override;
+
+  virtual KernelS exportKernelS_impl() const override;
+
+  virtual KernelD exportKernelD_impl() const override;
+
+  virtual DerivativeProbe exportDerivativeProbe_impl() const override;
+
+  virtual double singleLayer_impl(const Element & e, double factor) const override;
+
+  virtual double doubleLayer_impl(const Element & e, double factor) const override;
+
+  virtual std::ostream & printObject(std::ostream & os) override;
+  /*! This calculates all the components needed to evaluate the Green's function */
+  void initSphericalDiffuse();
+
+  /*! Center of the dielectric sphere */
+  Eigen::Vector3d origin_;
+
+  /*! @{ Parameters and functions for the calculation of the Green's function,
+   * including Coulomb singularity */
+  /*! Maximum angular momentum in the final summation over Legendre polynomials to
+   * obtain G */
+  int maxLGreen_;
+
+  /*! \brief First independent radial solution, used to build Green's function.
+   *  \note The vector has dimension maxLGreen_ and has r^l behavior
+   */
+  std::vector<RadialFunction<detail::LnTransformedRadial>> zeta_;
+
+  /*! \brief Second independent radial solution, used to build Green's function.
+   *  \note The vector has dimension maxLGreen_  and has r^(-l-1) behavior
+   */
+  std::vector<RadialFunction<detail::LnTransformedRadial>> omega_;
+
+  /*! \brief Returns L-th component of the radial part of the Green's function
+   *  \param[in] L  angular momentum
+   *  \param[in] sp source point
+   *  \param[in] pp probe point
+   *  \param[in] Cr12 Coulomb singularity separation coefficient
+   *  \note This function shifts the given source and probe points by the location
+   * of
+   * the
+   *  dielectric sphere.
+   */
+  double imagePotentialComponent_impl(int L,
+                                      const Eigen::Vector3d & sp,
+                                      const Eigen::Vector3d & pp,
+                                      double Cr12) const;
+  /*! @}*/
+  /*! @{ Parameters and functions for the calculation of the Coulomb singularity
+   * separation coefficient */
+  /*! Maximum angular momentum to obtain C(r, r'), needed to separate the Coulomb
+   * singularity */
+  int maxLC_; // = 2 * maxLGreen_;
+
+  /*! \brief First independent radial solution, used to build coefficient.
+   *  \note This is needed to separate the Coulomb singularity and has r^l behavior
+   */
+  RadialFunction<detail::LnTransformedRadial> zetaC_;
+
+  /*! \brief Second independent radial solution, used to build coefficient.
+   *  \note This is needed to separate the Coulomb singularity and has r^(-l-1)
+   * behavior
+   */
+  RadialFunction<detail::LnTransformedRadial> omegaC_;
+
+  /*! \brief Returns coefficient for the separation of the Coulomb singularity
+   *  \param[in] sp first point
+   *  \param[in] pp second point
+   *  \note This function shifts the given source and probe points by the location
+   * of
+   * the
+   *  dielectric sphere.
+   */
+  double coefficient_impl(const Eigen::Vector3d & sp,
+                          const Eigen::Vector3d & pp) const;
+  /*! @}*/
+};
+
+template <typename ProfilePolicy>
+IGreensFunction * createSphericalDiffuse(const GreenData & data) {
+  return new SphericalDiffuse<ProfilePolicy>(
+      data.epsilon1, data.epsilon2, data.width, data.center, data.origin, data.maxL);
+}
+} // namespace green
+} // namespace pcm
diff --git a/src/interface/Meddle.cpp b/src/interface/Meddle.cpp
index f8add79..ec00da5 100644
--- a/src/interface/Meddle.cpp
+++ b/src/interface/Meddle.cpp
@@ -442,7 +442,7 @@ void Meddle::initInput(int nr_nuclei,
   Eigen::Matrix3Xd centers = Eigen::Map<Eigen::Matrix3Xd>(coordinates, 3, nr_nuclei);
 
   if (input_.mode() != "EXPLICIT") {
-    Symmetry pg = buildGroup(
+    Symmetry pg(
         symmetry_info[0], symmetry_info[1], symmetry_info[2], symmetry_info[3]);
     input_.molecule(detail::initMolecule(input_, pg, nr_nuclei, chg, centers));
   }
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index bc75099..f58e200 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -22,6 +22,7 @@ list(APPEND headers_list
   Molecule.hpp
   PhysicalConstants.hpp
   QuadratureRules.hpp
+  RungeKutta.hpp
   Solvent.hpp
   Sphere.hpp
   SplineFunction.hpp
diff --git a/src/utils/MathUtils.hpp b/src/utils/MathUtils.hpp
index 24c8042..03b0538 100644
--- a/src/utils/MathUtils.hpp
+++ b/src/utils/MathUtils.hpp
@@ -44,48 +44,6 @@
 
 namespace pcm {
 namespace utils {
-/*! \fn inline int parity(std::bitset<nBits> bitrep)
- *  \param[in] bitrep a bitset
- *  \tparam nBits lenght of the input bitset
- *
- *  Calculate the parity of the bitset as defined by:
- *     bitrep[0] XOR bitrep[1] XOR ... XOR bitrep[nBits-1]
- */
-template <size_t nBits> inline int parity(std::bitset<nBits> bitrep) {
-  int parity = 0;
-  for (size_t i = 0; i < bitrep.size(); ++i) {
-    parity ^= bitrep[i];
-  }
-  return parity;
-}
-
-/*! \fn inline double parity(unsigned int i)
- *  \param[in] i an integer, usually an index for an irrep or a symmetry operation
- *
- * Returns parity of input integer.
- * The parity is defined as the result of using XOR on the bitrep
- * of the given integer. For example:
- *   2 -> 010 -> 0^1^0 = 1 -> -1.0
- *   6 -> 110 -> 1^1^0 = 0 ->  1.0
- *
- * It can also be interpreted as the action of a given operation on the
- * Cartesian axes:
- *      zyx         Parity
- *   0  000    E      1.0
- *   1  001   Oyz    -1.0
- *   2  010   Oxz    -1.0
- *   3  011   C2z     1.0
- *   4  100   Oxy    -1.0
- *   5  101   C2y     1.0
- *   6  110   C2x     1.0
- *   7  111    i     -1.0
- */
-inline double parity(unsigned int i) {
-  // Use a ternary if construct. If the bitset is odd return -1.0 Return +1.0
-  // otherwise.
-  return (parity(std::bitset<3>(i)) ? -1.0 : 1.0);
-}
-
 /*! \fn inline bool isZero(double value, double threshold)
  *  \param[in] value     the value to be checked
  *  \param[in] threshold the threshold
@@ -112,110 +70,7 @@ inline bool numericalZero(double value) { return (isZero(value, 1.0e-14)); }
  */
 template <typename T> inline int sign(T val) { return (T(0) < val) - (val < T(0)); }
 
-/*! \fn inline void symmetryBlocking(Eigen::MatrixXd & matrix, int cavitySize, int
- * ntsirr, int nr_irrep)
- *  \param[out] matrix the matrix to be block-diagonalized
- *  \param[in]  cavitySize the size of the cavity (size of the matrix)
- *  \param[in]  ntsirr     the size of the irreducible portion of the cavity (size of
- * the blocks)
- *  \param[in]  nr_irrep   the number of irreducible representations (number of
- * blocks)
- */
-inline void symmetryBlocking(Eigen::MatrixXd & matrix,
-                             int cavitySize,
-                             int ntsirr,
-                             int nr_irrep) {
-  // This function implements the simmetry-blocking of the PCM
-  // matrix due to point group symmetry as reported in:
-  // L. Frediani, R. Cammi, C. S. Pomelli, J. Tomasi and K. Ruud, J. Comput.Chem. 25,
-  // 375 (2003)
-  // u is the character table for the group (t in the paper)
-  Eigen::MatrixXd u = Eigen::MatrixXd::Zero(nr_irrep, nr_irrep);
-  for (int i = 0; i < nr_irrep; ++i) {
-    for (int j = 0; j < nr_irrep; ++j) {
-      u(i, j) = parity(i & j);
-    }
-  }
-  // Naming of indices:
-  //     a, b, c, d   run over the total size of the cavity (N)
-  //     i, j, k, l   run over the number of irreps (n)
-  //     p, q, r, s   run over the irreducible size of the cavity (N/n)
-  // Instead of forming U (T in the paper) and then perform the dense
-  // multiplication, we multiply block-by-block using just the u matrix.
-  //      matrix = U * matrix * Ut; U * Ut = Ut * U = id
-  // First half-transformation, i.e. first_half = matrix * Ut
-  Eigen::MatrixXd first_half = Eigen::MatrixXd::Zero(cavitySize, cavitySize);
-  for (int i = 0; i < nr_irrep; ++i) {
-    int ioff = i * ntsirr;
-    for (int k = 0; k < nr_irrep; ++k) {
-      int koff = k * ntsirr;
-      for (int j = 0; j < nr_irrep; ++j) {
-        int joff = j * ntsirr;
-        double ujk = u(j, k) / nr_irrep;
-        for (int p = 0; p < ntsirr; ++p) {
-          int a = ioff + p;
-          for (int q = 0; q < ntsirr; ++q) {
-            int b = joff + q;
-            int c = koff + q;
-            first_half(a, c) += matrix(a, b) * ujk;
-          }
-        }
-      }
-    }
-  }
-  // Second half-transformation, i.e. matrix = U * first_half
-  matrix.setZero(cavitySize, cavitySize);
-  for (int i = 0; i < nr_irrep; ++i) {
-    int ioff = i * ntsirr;
-    for (int k = 0; k < nr_irrep; ++k) {
-      int koff = k * ntsirr;
-      for (int j = 0; j < nr_irrep; ++j) {
-        int joff = j * ntsirr;
-        double uij = u(i, j);
-        for (int p = 0; p < ntsirr; ++p) {
-          int a = ioff + p;
-          int b = joff + p;
-          for (int q = 0; q < ntsirr; ++q) {
-            int c = koff + q;
-            matrix(a, c) += uij * first_half(b, c);
-          }
-        }
-      }
-    }
-  }
-  // Traverse the matrix and discard numerical zeros
-  for (int a = 0; a < cavitySize; ++a) {
-    for (int b = 0; b < cavitySize; ++b) {
-      if (numericalZero(matrix(a, b))) {
-        matrix(a, b) = 0.0;
-      }
-    }
-  }
-}
-
-/*! \fn inline void symmetryPacking(std::vector<Eigen::MatrixXd> & blockedMatrix,
- * const Eigen::MatrixXd & fullMatrix, int nrBlocks, int dimBlock)
- *  \param[out] blockedMatrix the result of packing fullMatrix
- *  \param[in]  fullMatrix the matrix to be packed
- *  \param[in]  dimBlock the dimension of the square blocks
- *  \param[in]  nrBlocks the number of square blocks
- */
-inline void symmetryPacking(std::vector<Eigen::MatrixXd> & blockedMatrix,
-                            const Eigen::MatrixXd & fullMatrix,
-                            int dimBlock,
-                            int nrBlocks) {
-  // This function packs the square block diagonal fullMatrix with nrBlocks of
-  // dimension dimBlock
-  // into a std::vector<Eigen::MatrixXd> containing nrBlocks square matrices of
-  // dimenion dimBlock.
-  int j = 0;
-  for (int i = 0; i < nrBlocks; ++i) {
-    blockedMatrix.push_back(fullMatrix.block(j, j, dimBlock, dimBlock));
-    j += dimBlock;
-  }
-}
-
-/*! \fn inline void hermitivitize(Eigen::MatrixBase<Derived> & obj_)
+/*! \brief Enforce Hermitian symmetry on a matrix
  *  \param[out] obj_ the Eigen object to be hermitivitized
  *  \tparam     Derived the numeric type of obj_ elements
  *
@@ -245,11 +100,11 @@ inline void hermitivitize(Eigen::MatrixBase<Derived> & obj_) {
  *rotation
  *  matrix will be:
  *  \f[
- *  	R = \begin{pmatrix}
- *  	      c_1c_3 - s_1c_2s_3 & -s_1c_3 - c_1c_2s_3 &  s_2s_3 \\
- *  	      c_1s_3 + s_1c_2c_3 & -s_1s_3 + c_1c_2c_3 & -s_2c_3 \\
- *  	      s_1s_2             & c_1s_2              &  c_2
- *  	    \end{pmatrix}
+ *    R = \begin{pmatrix}
+ *          c_1c_3 - s_1c_2s_3 & -s_1c_3 - c_1c_2s_3 &  s_2s_3 \\
+ *          c_1s_3 + s_1c_2c_3 & -s_1s_3 + c_1c_2c_3 & -s_2c_3 \\
+ *          s_1s_2             & c_1s_2              &  c_2
+ *        \end{pmatrix}
  *  \f]
  *  Eigen's geometry module is used to calculate the rotation matrix
  */
@@ -264,69 +119,6 @@ inline void eulerRotation(Eigen::Matrix3d & R_,
        Eigen::AngleAxis<double>(phi, Eigen::Vector3d::UnitZ());
 }
 
-/*! \brief Return value of function defined on grid at an arbitrary point
- *  \param[in] point     where the function has to be evaluated
- *  \param[in] grid      holds points on grid where function is known
- *  \param[in] function  holds known function values
- *
- *  This function finds the nearest values for the given point
- *  and performs a linear interpolation.
- *  \warning This function assumes that grid has already been sorted!
- */
-inline double linearInterpolation(const double point,
-                                  const std::vector<double> & grid,
-                                  const std::vector<double> & function) {
-  // Find nearest points on grid to the arbitrary point given
-  size_t index = std::distance(grid.begin(),
-                               std::lower_bound(grid.begin(), grid.end(), point)) -
-                 1;
-
-  // Parameters for the interpolating line
-  double y_1 = function[index], y_0 = function[index - 1];
-  double x_1 = grid[index], x_0 = grid[index - 1];
-  double m = (y_1 - y_0) / (x_1 - x_0);
-
-  return (m * (point - x_0) + y_0);
-}
-
-/*! \brief Return value of function defined on grid at an arbitrary point
- *  \param[in] point     where the function has to be evaluated
- *  \param[in] grid      holds points on grid where function is known
- *  \param[in] function  holds known function values
- *
- *  This function finds the nearest values for the given point
- *  and performs a cubic spline interpolation.
- *  \warning This function assumes that grid has already been sorted!
- */
-inline double splineInterpolation(const double point,
-                                  const std::vector<double> & grid,
-                                  const std::vector<double> & function) {
-  // Find nearest points on grid to the arbitrary point given
-  int index =
-      std::distance(grid.begin(), std::lower_bound(grid.begin(), grid.end(), point));
-
-  int imax = grid.size() - 1;
-  if (index <= 0)
-    index = 1;
-  if (index >= imax - 1)
-    index = imax - 2;
-
-  // Parameters for the interpolating spline
-  Eigen::Vector4d x = (Eigen::Vector4d() << grid[index - 1],
-                       grid[index],
-                       grid[index + 1],
-                       grid[index + 2])
-                          .finished();
-  Eigen::Vector4d y = (Eigen::Vector4d() << function[index - 1],
-                       function[index],
-                       function[index + 1],
-                       function[index + 2])
-                          .finished();
-  SplineFunction s(x, y);
-
-  return s(point);
-}
-
 /*! \brief Prints Eigen object (matrix or vector) to file
  *  \param[in] matrix Eigen object
  *  \param[in] fname  name of the file
diff --git a/src/utils/Molecule.cpp b/src/utils/Molecule.cpp
index bf21fa9..dab578f 100644
--- a/src/utils/Molecule.cpp
+++ b/src/utils/Molecule.cpp
@@ -56,7 +56,7 @@ Molecule::Molecule(int nat,
       atoms_(at),
       spheres_(sph) {
   rotor_ = findRotorType();
-  pointGroup_ = buildGroup(0, 0, 0, 0);
+  pointGroup_ = Symmetry(0, 0, 0, 0);
 }
 
 Molecule::Molecule(int nat,
@@ -74,7 +74,7 @@ Molecule::Molecule(int nat,
       atoms_(at),
       spheres_(sph) {
   rotor_ = findRotorType();
-  pointGroup_ = buildGroup(nr_gen, gen[0], gen[1], gen[2]);
+  pointGroup_ = Symmetry(nr_gen, gen);
 }
 
 Molecule::Molecule(int nat,
@@ -107,7 +107,7 @@ Molecule::Molecule(const std::vector<Sphere> & sph)
     atoms_.push_back(Atom("Dummy", "Du", charge, mass, mass, geometry_.col(i)));
   }
   rotor_ = findRotorType();
-  pointGroup_ = buildGroup(0, 0, 0, 0);
+  pointGroup_ = Symmetry(0, 0, 0, 0);
 }
 
 Molecule::Molecule(const Molecule & other) { *this = other; }
diff --git a/src/utils/Molecule.hpp b/src/utils/Molecule.hpp
index 8c9a9ec..6a3628e 100644
--- a/src/utils/Molecule.hpp
+++ b/src/utils/Molecule.hpp
@@ -46,6 +46,7 @@ class Element;
 namespace pcm {
 using utils::Atom;
 using utils::Sphere;
+using utils::Symmetry;
 
 enum rotorType { rtAsymmetric, rtSymmetric, rtSpherical, rtLinear, rtAtom };
 const std::string rotorTypeList[] = {"Asymmetric",
@@ -90,7 +91,7 @@ class Molecule {
    */
   Molecule() {
     rotor_ = rtAsymmetric;
-    pointGroup_ = buildGroup(0, 0, 0, 0);
+    pointGroup_ = Symmetry(0, 0, 0, 0);
   }
   /*! \brief Constructor from full molecular data
    *  \param[in] nat number of atoms
diff --git a/src/utils/RungeKutta.hpp b/src/utils/RungeKutta.hpp
new file mode 100644
index 0000000..4f0fb36
--- /dev/null
+++ b/src/utils/RungeKutta.hpp
@@ -0,0 +1,195 @@
+/*
+ * PCMSolver, an API for the Polarizable Continuum Model
+ * Copyright (C) 2018 Roberto Di Remigio, Luca Frediani and contributors.
+ *
+ * This file is part of PCMSolver.
+ *
+ * PCMSolver is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * PCMSolver 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 Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with PCMSolver.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ * For information on the complete list of contributors to the
+ * PCMSolver API, see: <http://pcmsolver.readthedocs.io/>
+ */
+
+#pragma once
+
+#include <limits>
+
+/*! \file RungeKutta.hpp */
+
+namespace pcm {
+namespace utils {
+namespace detail {
+/*! \return t1 < t2 if dt > 0 and t1 > t2 if dt < 0 with epsilon accuracy
+ *  \note This function is part of Boost.Odeint
+ */
+template <typename T> bool less_with_sign(T t1, T t2, T dt) {
+  if (dt > 0) { // return t1 < t2
+    return t2 - t1 > std::numeric_limits<T>::epsilon();
+  } else { // return t1 > t2
+    return t1 - t2 > std::numeric_limits<T>::epsilon();
+  }
+}
+
+/*! \return t1 <= t2 if dt > 0 and t1 => t2 if dt < 0 with epsilon accuracy
+ *  \note This function is part of Boost.Odeint
+ */
+template <typename T> bool less_eq_with_sign(T t1, T t2, T dt) {
+  if (dt > 0) {
+    return t1 - t2 <= std::numeric_limits<T>::epsilon();
+  } else {
+    return t2 - t1 <= std::numeric_limits<T>::epsilon();
+  }
+}
+} // namespace detail
+
+/*! \brief Integrates an ODE given a stepper.
+ *  \author Roberto Di Remigio
+ *  \date 2018
+ *  \tparam Stepper  integrate stepping function
+ *  \tparam System   the ODE system to integrate
+ *  \tparam Observer function reporting integration progress
+ *
+ *  The interface mimics the one of the Boost.Odeint package.
+ */
+template <typename Stepper, typename System, typename Observer>
+size_t integrate_const(const Stepper & stepper,
+                       const System & system,
+                       typename System::StateType & f,
+                       const double tStart,
+                       const double tEnd,
+                       const double dt,
+                       const Observer & observer) {
+  size_t nSteps = 0;
+  double t = tStart;
+  while (detail::less_eq_with_sign(t + dt, tEnd, dt)) {
+    observer(f, t);
+    stepper.doStep(system, f, t, dt);
+    ++nSteps;
+    t = tStart + nSteps * dt;
+  }
+  observer(f, t);
+  return nSteps;
+}
+
+/*! \brief 4th-order Runge-Kutta stepper.
+ *  \author Roberto Di Remigio
+ *  \date 2018
+ *  \tparam StateType the data structure used to hold the state of the ODE
+ *
+ *  The interface mimics the one of the Boost.Odeint package.
+ */
+template <typename StateType> struct RungeKutta4 final {
+  /*! \brief Implements the "classic" 4th-order Runge-Kutta stepper.
+   *  \author Roberto Di Remigio
+   *  \date 2018
+   *  \tparam System the ODE system to integrate
+   *
+   *  This function implements the "classic" Runge-Kutta stepper, with the
+   *  following Butcher Tableau:
+   *
+   *  \f[
+   *    \begin{array}{c|cccc}
+   *     0   & 0   & 0   & 0   & 0\\
+   *     1/2 & 1/2 & 0   & 0   & 0\\
+   *     1/2 & 0   & 1/2 & 0   & 0\\
+   *     1   & 0   & 0   & 1   & 0\\
+   *     \hline
+   *         & 1/6 & 1/3 & 1/3 & 1/6\\
+   *    \end{array}
+   *  \f]
+   */
+  template <typename System>
+  void doStep(const System & system,
+              StateType & f,
+              const double t,
+              const double dt) const {
+    StateType k1, k2, k3, k4;
+    StateType x1, x2, x3;
+
+    // Step 1
+    system(f, k1, t);
+    // Step2
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      x1[i] = f[i] + 0.5 * dt * k1[i];
+    }
+    system(x1, k2, t + 0.5 * dt);
+    // Step 3
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      x2[i] = f[i] + 0.5 * dt * k2[i];
+    }
+    system(x2, k3, t + 0.5 * dt);
+    // Step 4
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      x3[i] = f[i] + dt * k3[i];
+    }
+    system(x3, k4, t + dt);
+    // Update
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      f[i] += dt * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0;
+    }
+  }
+
+  /*! \brief Implements the 3/8 4th-order Runge-Kutta stepper.
+   *  \author Roberto Di Remigio
+   *  \date 2018
+   *  \tparam System the ODE system to integrate
+   *
+   *  This function implements the 3/8 Runge-Kutta stepper, with the
+   *  following Butcher Tableau:
+   *
+   *  \f[
+   *    \begin{array}{c|cccc}
+   *    0   & 0   & 0   & 0   & 0\\
+   *    1/3 & 1/3 & 0   & 0   & 0\\
+   *    2/3 & -1/3   & 1 & 0   & 0\\
+   *    1   & 1   & -1   & 1   & 0\\
+   *    \hline
+   *        & 1/8 & 3/8 & 3/8 & 1/8\\
+   *    \end{array}
+   *  \f]
+   */
+  template <typename System>
+  void doStep38(const System & system,
+                StateType & f,
+                const double t,
+                const double dt) const {
+    StateType k1{}, k2{}, k3{}, k4{};
+    StateType x1{}, x2{}, x3{};
+
+    // Step 1
+    system(f, k1, t);
+    // Step 2
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      x1[i] = f[i] + (dt / 3.0) * k1[i];
+    }
+    system(x1, k2, t + (dt / 3.0));
+    // Step 3
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      x2[i] = f[i] + dt * (-1.0 / 3.0 * k1[i] + k2[i]);
+    }
+    system(x2, k3, t + dt * (2.0 / 3.0));
+    // Step 4
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      x3[i] = f[i] + dt * (k1[i] - k2[i] + k3[i]);
+    }
+    system(x3, k4, t + dt);
+    // Update
+    for (size_t i = 0; i < system.ODEorder(); ++i) {
+      f[i] += dt * ((1.0 / 8.0) * k1[i] + (3.0 / 8.0) * k2[i] + (3.0 / 8.0) * k3[i] +
+                    (1.0 / 8.0) * k4[i]);
+    }
+  }
+};
+} // namespace utils
+} // namespace pcm
diff --git a/src/utils/RungeKutta4.cpp b/src/utils/RungeKutta4.cpp
deleted file mode 100644
index c35a405..0000000
--- a/src/utils/RungeKutta4.cpp
+++ /dev/null
@@ -1,134 +0,0 @@
-/*
- * The Runge-Kutta source code is Copyright(c) 2010 John Burkardt.
- *
- * This Runge-Kutta ODE solver is free software: you can redistribute it and/or
- * modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This Runge-Kutts ODE Solver 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with the code.  If not, see <http://www.gnu.org/licenses/>.
- *
- * See also:
- * https://people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html
- */
-
-#include "rk4.hpp"
-
-#include <cmath>
-#include <cstdlib>
-#include <iomanip>
-#include <iostream>
-
-namespace pcm {
-namespace utils {
-using namespace std;
-
-double rk4(double t0, double u0, double dt, double f(double t, double u)) {
-  double f0;
-  double f1;
-  double f2;
-  double f3;
-  double t1;
-  double t2;
-  double t3;
-  double u;
-  double u1;
-  double u2;
-  double u3;
-  //
-  //  Get four sample values of the derivative.
-  //
-  f0 = f(t0, u0);
-
-  t1 = t0 + dt / 2.0;
-  u1 = u0 + dt * f0 / 2.0;
-  f1 = f(t1, u1);
-
-  t2 = t0 + dt / 2.0;
-  u2 = u0 + dt * f1 / 2.0;
-  f2 = f(t2, u2);
-
-  t3 = t0 + dt;
-  u3 = u0 + dt * f2;
-  f3 = f(t3, u3);
-  //
-  //  Combine to estimate the solution at time T0 + DT.
-  //
-  u = u0 + dt * (f0 + 2.0 * f1 + 2.0 * f2 + f3) / 6.0;
-
-  return u;
-}
-
-double * rk4vec(double t0,
-                int m,
-                double u0[],
-                double dt,
-                double * f(double t, int m, double u[]))
-
-{
-  double * f0;
-  double * f1;
-  double * f2;
-  double * f3;
-  int i;
-  double t1;
-  double t2;
-  double t3;
-  double * u;
-  double * u1;
-  double * u2;
-  double * u3;
-  //
-  //  Get four sample values of the derivative.
-  //
-  f0 = f(t0, m, u0);
-
-  t1 = t0 + dt / 2.0;
-  u1 = new double[m];
-  for (i = 0; i < m; i++) {
-    u1[i] = u0[i] + dt * f0[i] / 2.0;
-  }
-  f1 = f(t1, m, u1);
-
-  t2 = t0 + dt / 2.0;
-  u2 = new double[m];
-  for (i = 0; i < m; i++) {
-    u2[i] = u0[i] + dt * f1[i] / 2.0;
-  }
-  f2 = f(t2, m, u2);
-
-  t3 = t0 + dt;
-  u3 = new double[m];
-  for (i = 0; i < m; i++) {
-    u3[i] = u0[i] + dt * f2[i];
-  }
-  f3 = f(t3, m, u3);
-  //
-  //  Combine them to estimate the solution.
-  //
-  u = new double[m];
-  for (i = 0; i < m; i++) {
-    u[i] = u0[i] + dt * (f0[i] + 2.0 * f1[i] + 2.0 * f2[i] + f3[i]) / 6.0;
-  }
-  //
-  //  Free memory.
-  //
-  delete[] f0;
-  delete[] f1;
-  delete[] f2;
-  delete[] f3;
-  delete[] u1;
-  delete[] u2;
-  delete[] u3;
-
-  return u;
-}
-} // namespace utils
-} // namespace pcm
diff --git a/src/utils/RungeKutta4.hpp b/src/utils/RungeKutta4.hpp
deleted file mode 100644
index 8759168..0000000
--- a/src/utils/RungeKutta4.hpp
+++ /dev/null
@@ -1,82 +0,0 @@
-/*
- * The Runge-Kutta source code is Copyright(c) 2013 John Burkardt.
- *
- * This Runge-Kutta ODE solver is free software: you can redistribute it and/or
- * modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This Runge-Kutts ODE Solver 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 Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with the code.  If not, see <http://www.gnu.org/licenses/>.
- *
- * See also:
- * https://people.sc.fsu.edu/~jburkardt/cpp_src/rk4/rk4.html
- */
-
-/*! \file RungeKutta4.hpp */
-
-namespace pcm {
-namespace utils {
-/*! Right-hand side for a scalar ordinary-differential equation */
-typedef std::function<double(double, double)> ScalarODERHS;
-
-/*! \brief  RK4 takes one Runge-Kutta step for a scalar ODE.
- *  \author John Burkardt, Roberto Di Remigio
- *  \param[in] t0 the current time.
- *  \param[in] u0 the solution estimate at the current time.
- *  \param[in] dt the time step.
- *  \param[in] F a function which evaluates the derivative, or right hand side of the
- * problem.
- *  \return the fourth-order Runge-Kutta solution estimate at time t0+dt.
- *
- *  It is assumed that an initial value problem, of the form
- *
- *    du/dt = f( t, u )
- *    u(t0) = u0
- *
- *  is being solved.
- *
- *  If the user can supply current values of t, u, a stepsize dt, and a
- *  function to evaluate the derivative, this function can compute the
- *  fourth-order Runge Kutta estimate to the solution at time t+dt.
- */
-double rk4(double t0, double u0, double dt, double f(double t, double u));
-
-/*! Right-hand side for a vector ordinary-differential equation */
-typedef std::function<double *(double, int, double)> VectorODERHS;
-
-/*!
- *  \brief Takes one Runge-Kutta step for a vector ODE.
- *  \author John Burkardt, Roberto Di Remigio
- *  \param[in] t0 the current time.
- *  \param[in] n the spatial dimension.
- *  \param[in] u0[M] the solution estimate at the current time.
- *  \param[in] dt the time step.
- *  \param[in] F a function which evaluates the derivative, or right hand side of the
- * problem.
- *  \return the fourth-order Runge-Kutta solution estimate at time t0+dt.
- *
- *  It is assumed that an initial value problem, of the form
- *
- *    du/dt = f ( t, u )
- *    u(t0) = u0
- *
- *  is being solved.
- *
- *  If the user can supply current values of t, u, a stepsize dt, and a
- *  function to evaluate the derivative, this function can compute the
- *  fourth-order Runge Kutta estimate to the solution at time t+dt.
- */
-double * rk4vec(double t0,
-                int n,
-                double u0[],
-                double dt,
-                double * f(double t, int n, double u[]));
-} // namespace utils
-} // namespace pcm
diff --git a/src/utils/SplineFunction.hpp b/src/utils/SplineFunction.hpp
index 79284c3..c9eba69 100644
--- a/src/utils/SplineFunction.hpp
+++ b/src/utils/SplineFunction.hpp
@@ -24,19 +24,22 @@
 #pragma once
 
 #include <algorithm>
+#include <vector>
 
 #include <Eigen/Core>
 #include <unsupported/Eigen/Splines>
 
-/*! \file SplineFunction.hpp
- *  \class SplineFunction
+/*! \file SplineFunction.hpp */
+
+namespace pcm {
+namespace utils {
+/*! \class SplineFunction
  *  \brief Spline interpolation of a function
  *  \author Roberto Di Remigio
  *  \date 2015
  *
  *  Taken from StackOverflow http://stackoverflow.com/a/29825204/2528668
  */
-
 class SplineFunction final {
 private:
   typedef Eigen::Spline<double, 1> CubicSpline;
@@ -84,3 +87,43 @@ class SplineFunction final {
         .transpose();
   }
 };
+
+/*! \brief Return value of function defined on grid at an arbitrary point
+ *  \param[in] point     where the function has to be evaluated
+ *  \param[in] grid      holds points on grid where function is known
+ *  \param[in] function  holds known function values
+ *
+ *  This function finds the nearest values for the given point
+ *  and performs a cubic spline interpolation.
+ *  \warning This function assumes that grid has already been sorted!
+ */
+inline double splineInterpolation(const double point,
+                                  const std::vector<double> & grid,
+                                  const std::vector<double> & function) {
+  // Find nearest points on grid to the arbitrary point given
+  int index =
+      std::distance(grid.begin(), std::lower_bound(grid.begin(), grid.end(), point));
+
+  int imax = grid.size() - 1;
+  if (index <= 0)
+    index = 1;
+  if (index >= imax - 1)
+    index = imax - 2;
+
+  // Parameters for the interpolating spline
+  Eigen::Vector4d x = (Eigen::Vector4d() << grid[index - 1],
+                       grid[index],
+                       grid[index + 1],
+                       grid[index + 2])
+                          .finished();
+  Eigen::Vector4d y = (Eigen::Vector4d() << function[index - 1],
+                       function[index],
+                       function[index + 1],
+                       function[index + 2])
+                          .finished();
+  SplineFunction s(x, y);
+
+  return s(point);
+}
+} // namespace utils
+} // namespace pcm
diff --git a/src/utils/Symmetry.cpp b/src/utils/Symmetry.cpp
index efb6098..090b538 100644
--- a/src/utils/Symmetry.cpp
+++ b/src/utils/Symmetry.cpp
@@ -23,25 +23,91 @@
 
 #include "Symmetry.hpp"
 
+#include <vector>
+
 #include <Eigen/Core>
 
-/* Indexing of symmetry operations and their mapping to a bitstring:
- *      zyx         Parity
- *   0  000    E      1.0
- *   1  001   Oyz    -1.0
- *   2  010   Oxz    -1.0
- *   3  011   C2z     1.0
- *   4  100   Oxy    -1.0
- *   5  101   C2y     1.0
- *   6  110   C2x     1.0
- *   7  111    i     -1.0
- */
+#include "MathUtils.hpp"
 
-Symmetry buildGroup(int _nr_gen, int _gen1, int _gen2, int _gen3) {
-  int gen[3];
-  gen[0] = _gen1;
-  gen[1] = _gen2;
-  gen[2] = _gen3;
+namespace pcm {
+namespace utils {
+void symmetryBlocking(Eigen::MatrixXd & matrix,
+                      int cavitySize,
+                      int ntsirr,
+                      int nr_irrep) {
+  // u is the character table for the group (t in the paper)
+  Eigen::MatrixXd u = Eigen::MatrixXd::Zero(nr_irrep, nr_irrep);
+  for (int i = 0; i < nr_irrep; ++i) {
+    for (int j = 0; j < nr_irrep; ++j) {
+      u(i, j) = parity(i & j);
+    }
+  }
+  // Naming of indices:
+  //     a, b, c, d   run over the total size of the cavity (N)
+  //     i, j, k, l   run over the number of irreps (n)
+  //     p, q, r, s   run over the irreducible size of the cavity (N/n)
+  // Instead of forming U (T in the paper) and then perform the dense
+  // multiplication, we multiply block-by-block using just the u matrix.
+  //      matrix = U * matrix * Ut; U * Ut = Ut * U = id
+  // First half-transformation, i.e. first_half = matrix * Ut
+  Eigen::MatrixXd first_half = Eigen::MatrixXd::Zero(cavitySize, cavitySize);
+  for (int i = 0; i < nr_irrep; ++i) {
+    int ioff = i * ntsirr;
+    for (int k = 0; k < nr_irrep; ++k) {
+      int koff = k * ntsirr;
+      for (int j = 0; j < nr_irrep; ++j) {
+        int joff = j * ntsirr;
+        double ujk = u(j, k) / nr_irrep;
+        for (int p = 0; p < ntsirr; ++p) {
+          int a = ioff + p;
+          for (int q = 0; q < ntsirr; ++q) {
+            int b = joff + q;
+            int c = koff + q;
+            first_half(a, c) += matrix(a, b) * ujk;
+          }
+        }
+      }
+    }
+  }
+  // Second half-transformation, i.e. matrix = U * first_half
+  matrix.setZero(cavitySize, cavitySize);
+  for (int i = 0; i < nr_irrep; ++i) {
+    int ioff = i * ntsirr;
+    for (int k = 0; k < nr_irrep; ++k) {
+      int koff = k * ntsirr;
+      for (int j = 0; j < nr_irrep; ++j) {
+        int joff = j * ntsirr;
+        double uij = u(i, j);
+        for (int p = 0; p < ntsirr; ++p) {
+          int a = ioff + p;
+          int b = joff + p;
+          for (int q = 0; q < ntsirr; ++q) {
+            int c = koff + q;
+            matrix(a, c) += uij * first_half(b, c);
+          }
+        }
+      }
+    }
+  }
+  // Traverse the matrix and discard numerical zeros
+  for (int a = 0; a < cavitySize; ++a) {
+    for (int b = 0; b < cavitySize; ++b) {
+      if (numericalZero(matrix(a, b))) {
+        matrix(a, b) = 0.0;
+      }
+    }
+  }
+}
 
-  return Symmetry(_nr_gen, gen);
+void symmetryPacking(std::vector<Eigen::MatrixXd> & blockedMatrix,
+                     const Eigen::MatrixXd & fullMatrix,
+                     int dimBlock,
+                     int nrBlocks) {
+  int j = 0;
+  for (int i = 0; i < nrBlocks; ++i) {
+    blockedMatrix.push_back(fullMatrix.block(j, j, dimBlock, dimBlock));
+    j += dimBlock;
+  }
 }
+} // namespace utils
+} // namespace pcm
diff --git a/src/utils/Symmetry.hpp b/src/utils/Symmetry.hpp
index 4d6aa15..e26bc26 100644
--- a/src/utils/Symmetry.hpp
+++ b/src/utils/Symmetry.hpp
@@ -24,57 +24,134 @@
 #pragma once
 
 #include <algorithm>
+#include <bitset>
 #include <cmath>
+#include <vector>
 
-/*! \file Symmetry.hpp
- *  \class Symmetry
+#include <Eigen/Core>
+
+/*! \file Symmetry.hpp */
+
+namespace pcm {
+namespace utils {
+/*! \class Symmetry
  *  \brief Contains very basic info about symmetry (only Abelian groups)
  *  \author Roberto Di Remigio
- *  \date 2014
+ *  \date 2014, 2018
+ *  \note C1 is built as Symmetry C1(0, 0, 0, 0);
  *
- *  Just a wrapper around a vector containing the generators of the group
+ *  Indexing of symmetry operations and their mapping to a bitstring:
+ *       zyx         Parity
+ *    0  000    E      1.0
+ *    1  001   Oyz    -1.0
+ *    2  010   Oxz    -1.0
+ *    3  011   C2z     1.0
+ *    4  100   Oxy    -1.0
+ *    5  101   C2y     1.0
+ *    6  110   C2x     1.0
+ *    7  111    i     -1.0
  */
+class Symmetry final {
+public:
+  Symmetry()
+      : nrGenerators_(0),
+        nrIrrep_(static_cast<int>(std::pow(2.0, nrGenerators_))),
+        generators_({{0, 0, 0}}) {}
+  Symmetry(int nr_gen, int gen1, int gen2, int gen3)
+      : nrGenerators_(nr_gen),
+        nrIrrep_(static_cast<int>(std::pow(2.0, nrGenerators_))),
+        generators_({{gen1, gen2, gen3}}) {}
+  Symmetry(int nr_gen, int gen[3])
+      : nrGenerators_(nr_gen),
+        nrIrrep_(static_cast<int>(std::pow(2.0, nrGenerators_))) {
+    generators_[0] = gen[0];
+    generators_[1] = gen[1];
+    generators_[2] = gen[2];
+  }
+
+  int nrGenerators() const { return nrGenerators_; }
+  int nrIrrep() const { return nrIrrep_; }
+  int generators(size_t i) const { return generators_[i]; }
 
-class Symmetry {
 private:
-  /*!
-   * Number of generators
-   */
+  /*! Number of generators */
   int nrGenerators_;
-  /*!
-   * Generators
-   */
-  int generators_[3];
-  /*!
-   * Number of irreps
-   */
+  /*! Number of irreps */
   int nrIrrep_;
+  /*! Generators */
+  std::array<int, 3> generators_;
+};
 
-public:
-  /*! \brief Default constructor sets up C1 point group
-   */
-  Symmetry() : nrGenerators_(0) {
-    std::fill(generators_, generators_ + 3, 0);
-    nrIrrep_ = int(std::pow(2.0, nrGenerators_));
-  }
-  Symmetry(int nr_gen, int gen[3]) : nrGenerators_(nr_gen) {
-    // Transfer the passed generators array into generators_
-    std::copy(gen, gen + nrGenerators_, generators_);
-    // We can now initialize the number of irreps
-    nrIrrep_ = int(std::pow(2.0, nrGenerators_));
-  }
-  Symmetry(const Symmetry & other)
-      : nrGenerators_(other.nrGenerators_), nrIrrep_(other.nrIrrep_) {
-    std::copy(other.generators_, other.generators_ + nrGenerators_, generators_);
+/*! \brief Transform matrix to block diagonal form by symmetry
+ *  \param[out] matrix the matrix to be block-diagonalized
+ *  \param[in]  cavitySize the size of the cavity (size of the matrix)
+ *  \param[in]  ntsirr     the size of the irreducible portion of the cavity (size of
+ * the blocks)
+ *  \param[in]  nr_irrep   the number of irreducible representations (number of
+ * blocks)
+ *
+ * This function implements the symmetry-blocking of the PCM matrix due to
+ * point group symmetry as reported in:
+ * L. Frediani, R. Cammi, C. S. Pomelli, J. Tomasi and K. Ruud, J. Comput.Chem. 25,
+ * 375 (2003)
+ */
+void symmetryBlocking(Eigen::MatrixXd & matrix,
+                      int cavitySize,
+                      int ntsirr,
+                      int nr_irrep);
+
+/*! \brief Packs symmetry blocked matrix, i.e. only stores non-zero blocks
+ *  \param[out] blockedMatrix the result of packing fullMatrix
+ *  \param[in]  fullMatrix the matrix to be packed
+ *  \param[in]  dimBlock the dimension of the square blocks
+ *  \param[in]  nrBlocks the number of square blocks
+ *
+ *  This function packs the square block diagonal fullMatrix with nrBlocks of
+ *  dimension dimBlock into a std::vector<Eigen::MatrixXd> containing nrBlocks
+ *  square matrices of dimenion dimBlock.
+ */
+void symmetryPacking(std::vector<Eigen::MatrixXd> & blockedMatrix,
+                     const Eigen::MatrixXd & fullMatrix,
+                     int dimBlock,
+                     int nrBlocks);
+
+/*! \brief Calculate the parity of the bitset
+ *  \param[in] bitrep a bitset
+ *  \tparam nBits length of the input bitset
+ *
+ *  Calculate the parity of the bitset as defined by:
+ *     bitrep[0] XOR bitrep[1] XOR ... XOR bitrep[nBits-1]
+ */
+template <size_t nBits> inline int parity(std::bitset<nBits> bitrep) {
+  int parity = 0;
+  for (size_t i = 0; i < bitrep.size(); ++i) {
+    parity ^= bitrep[i];
   }
-  ~Symmetry() {}
-  int nrGenerators() const { return nrGenerators_; }
-  int generators(int i) const { return generators_[i]; }
-  int nrIrrep() const { return nrIrrep_; }
-};
+  return parity;
+}
 
-/*! Builds Symmetry object.
+/*! \brief Returns parity of input integer.
+ *  \param[in] i an integer, usually an index for an irrep or a symmetry operation
+ *
+ * The parity is defined as the result of using XOR on the bitrep
+ * of the given integer. For example:
+ *   2 -> 010 -> 0^1^0 = 1 -> -1.0
+ *   6 -> 110 -> 1^1^0 = 0 ->  1.0
  *
- * \note C1 is built as Symmetry C1 = buildGroup(0, 0, 0, 0);
+ * It can also be interpreted as the action of a given operation on the
+ * Cartesian axes:
+ *      zyx         Parity
+ *   0  000    E      1.0
+ *   1  001   Oyz    -1.0
+ *   2  010   Oxz    -1.0
+ *   3  011   C2z     1.0
+ *   4  100   Oxy    -1.0
+ *   5  101   C2y     1.0
+ *   6  110   C2x     1.0
+ *   7  111    i     -1.0
  */
-Symmetry buildGroup(int _nr_gen, int _gen1, int _gen2, int _gen3);
+inline double parity(unsigned int i) {
+  return (parity(std::bitset<3>(i)) ? -1.0 : 1.0);
+}
+} // namespace utils
+} // namespace pcm
diff --git a/tests/TestingMolecules.hpp b/tests/TestingMolecules.hpp
index 38e4b3d..bb1e8e4 100644
--- a/tests/TestingMolecules.hpp
+++ b/tests/TestingMolecules.hpp
@@ -82,38 +82,38 @@ template <int group> Molecule dummy(double radius, const Eigen::Vector3d & cente
   Symmetry pGroup;
   switch (group) {
     case (pgC1):
-      pGroup = buildGroup(0, 0, 0, 0);
+      pGroup = Symmetry(0, 0, 0, 0);
       break;
     case (pgC2):
       // C2 as generated by C2z
-      pGroup = buildGroup(1, 3, 0, 0);
+      pGroup = Symmetry(1, 3, 0, 0);
       break;
     case (pgCs):
       // Cs as generated by Oyz
-      pGroup = buildGroup(1, 1, 0, 0);
+      pGroup = Symmetry(1, 1, 0, 0);
       break;
     case (pgCi):
       // Ci as generated by i
-      pGroup = buildGroup(1, 7, 0, 0);
+      pGroup = Symmetry(1, 7, 0, 0);
       break;
     case (pgD2):
       // D2 as generated by C2z and C2x
-      pGroup = buildGroup(2, 3, 6, 0);
+      pGroup = Symmetry(2, 3, 6, 0);
       break;
     case (pgC2v):
       // C2v as generated by Oyz and Oxz
-      pGroup = buildGroup(2, 1, 2, 0);
+      pGroup = Symmetry(2, 1, 2, 0);
       break;
     case (pgC2h):
       // C2h as generated by Oxy and i
-      pGroup = buildGroup(2, 4, 7, 0);
+      pGroup = Symmetry(2, 4, 7, 0);
       break;
     case (pgD2h):
       // D2h as generated by Oxy, Oxz and Oyz
-      pGroup = buildGroup(3, 4, 2, 1);
+      pGroup = Symmetry(3, 4, 2, 1);
       break;
     default:
-      pGroup = buildGroup(0, 0, 0, 0);
+      pGroup = Symmetry(0, 0, 0, 0);
       break;
   }
 
@@ -159,7 +159,7 @@ Molecule NH3() {
   spheres.push_back(sph4);
 
   // C1
-  Symmetry pGroup = buildGroup(0, 0, 0, 0);
+  Symmetry pGroup(0, 0, 0, 0);
 
   return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };
@@ -197,14 +197,14 @@ template <int group> Molecule H3() {
   Symmetry pGroup;
   switch (group) {
     case (pgC1):
-      pGroup = buildGroup(0, 0, 0, 0);
+      pGroup = Symmetry(0, 0, 0, 0);
       break;
     case (pgC2v):
       // C2v as generated by Oyz and Oxz
-      pGroup = buildGroup(2, 1, 2, 0);
+      pGroup = Symmetry(2, 1, 2, 0);
       break;
     default:
-      pGroup = buildGroup(0, 0, 0, 0);
+      pGroup = Symmetry(0, 0, 0, 0);
       break;
   }
 
@@ -235,7 +235,7 @@ Molecule H2() {
   spheres.push_back(sph2);
   spheres.push_back(sph3);
 
-  Symmetry pGroup = buildGroup(0, 0, 0, 0);
+  Symmetry pGroup(0, 0, 0, 0);
 
   return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };
@@ -274,38 +274,38 @@ template <int group> Molecule CO2() {
   Symmetry pGroup;
   switch (group) {
     case (pgC1):
-      pGroup = buildGroup(0, 0, 0, 0);
+      pGroup = Symmetry(0, 0, 0, 0);
       break;
     case (pgC2):
       // C2 as generated by C2z
-      pGroup = buildGroup(1, 3, 0, 0);
+      pGroup = Symmetry(1, 3, 0, 0);
       break;
     case (pgCs):
       // Cs as generated by Oyz
-      pGroup = buildGroup(1, 1, 0, 0);
+      pGroup = Symmetry(1, 1, 0, 0);
       break;
     case (pgCi):
       // Ci as generated by i
-      pGroup = buildGroup(1, 7, 0, 0);
+      pGroup = Symmetry(1, 7, 0, 0);
       break;
     case (pgD2):
       // D2 as generated by C2z and C2x
-      pGroup = buildGroup(2, 3, 6, 0);
+      pGroup = Symmetry(2, 3, 6, 0);
       break;
     case (pgC2v):
       // C2v as generated by Oyz and Oxz
-      pGroup = buildGroup(2, 1, 2, 0);
+      pGroup = Symmetry(2, 1, 2, 0);
       break;
     case (pgC2h):
       // C2h as generated by Oxy and i
-      pGroup = buildGroup(2, 4, 7, 0);
+      pGroup = Symmetry(2, 4, 7, 0);
       break;
     case (pgD2h):
       // D2h as generated by Oxy, Oxz and Oyz
-      pGroup = buildGroup(3, 4, 2, 1);
+      pGroup = Symmetry(3, 4, 2, 1);
       break;
     default:
-      pGroup = buildGroup(0, 0, 0, 0);
+      pGroup = Symmetry(0, 0, 0, 0);
       break;
   }
 
@@ -348,7 +348,7 @@ Molecule CH3() {
   spheres.push_back(sph4);
 
   // Cs as generated by Oxy
-  Symmetry pGroup = buildGroup(1, 4, 0, 0);
+  Symmetry pGroup(1, 4, 0, 0);
 
   return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };
@@ -399,7 +399,7 @@ Molecule C2H4() {
   spheres.push_back(sph6);
 
   // D2h as generated by Oxy, Oxz, Oyz
-  Symmetry pGroup = buildGroup(3, 4, 2, 1);
+  Symmetry pGroup(3, 4, 2, 1);
 
   return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };
@@ -496,7 +496,7 @@ Molecule C6H6() {
   spheres.push_back(sph11);
   spheres.push_back(sph12);
 
-  Symmetry pGroup = buildGroup(0, 0, 0, 0);
+  Symmetry pGroup(0, 0, 0, 0);
 
   return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };
@@ -533,7 +533,7 @@ Molecule H2O() {
   spheres.push_back(sph2);
   spheres.push_back(sph3);
 
-  Symmetry pGroup = buildGroup(0, 0, 0, 0);
+  Symmetry pGroup(0, 0, 0, 0);
 
   return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
 };