diff --git a/include/tadah/models/cutoffs.h b/include/tadah/models/cutoffs.h
index 882f0c99e37ae34e3350c52e199c0ab975280d34..9d43af6bd44b6249fc6fa82acb5bdf2cf2c8d3bf 100644
--- a/include/tadah/models/cutoffs.h
+++ b/include/tadah/models/cutoffs.h
@@ -18,16 +18,33 @@ class Cut_Base {
         void test_rcut(const double r);
 };
 
-/** @brief Dummy cutoff function
+/**
+ * @class Cut_Dummy
+ * @brief Represents a basic cutoff function with a sharp transition at the cutoff radius.
+ *
+ * The `Cut_Dummy` class implements a simple cutoff function defined by:
  * \f[
  * f_c(r) =
- * \begin{equation}
  * \begin{cases}
- *  1 & \text{if  } r \leq r_c\\
- *  0 & \text{otherwise}\\
+ * 1, & \text{if } r \leq r_{\text{cut}} \\
+ * 0, & \text{if } r > r_{\text{cut}}
  * \end{cases}
- * \end{equation}
  * \f]
+ *
+ * and its derivative:
+ * \f[
+ * f_c'(r) = 0
+ * \f]
+ *
+ * where:
+ * - \f$ r \f$ is the radial distance.
+ * - \f$ r_{\text{cut}} \f$ is the cutoff radius.
+ *
+ * **Characteristics:**
+ * - The function value is constant (1) within the cutoff radius and zero beyond it.
+ * - The derivative of the function is zero everywhere except at \f$ r = r_{\text{cut}} \f$ where it is undefined due to the discontinuity.
+ *
+ * **Note:** Since the function is discontinuous at \f$ r = r_{\text{cut}} \f$ its derivative is not defined at that point. In practical computations, the derivative function `calc_prime` returns `0.0` for all \f$ r \f$
  */
 class Cut_Dummy : public Cut_Base {
     private:
@@ -44,22 +61,25 @@ class Cut_Dummy : public Cut_Base {
         double calc_prime(double r);
 };
 
-/** @brief Cos cutoff function
+/**
+ * @class Cut_Cos
+ * @brief Cosine cutoff function.
  *
+ * The `Cut_Cos` class implements a smooth cosine cutoff function defined by:
  * \f[
  * f_c(r) =
- * \begin{equation}
  * \begin{cases}
- *  \frac{1}{2}\big[ \cos\big( \frac{\pi r}{r_c} \big)+1 \big] & \text{if  } r \leq r_c\\
- *  0 & \text{otherwise}\\
+ *  \dfrac{1}{2}\left[ \cos\left( \dfrac{\pi r}{r_c} \right) + 1 \right], & \text{if } r \leq r_c \\
+ *  0, & \text{if } r > r_c
  * \end{cases}
- * \end{equation}
  * \f]
  *
- * <div class="csl-entry">Behler, J., Parrinello, M. (2007).
- * Generalized neural-network representation of high-dimensional
- * potential-energy surfaces. <i>Physical Review Letters</i>,
- * <i>98</i>(14), 146401. https://doi.org/10.1103/PhysRevLett.98.146401</div>
+ * This function smoothly transitions from 1 to 0 over the interval \f$ 0 \leq r \leq r_c \f$.
+ * It is commonly used in molecular simulations to smoothly truncate interactions without introducing discontinuities in the potential energy or its derivatives.
+ *
+ * **Reference:**
+ *
+ * Behler, J., & Parrinello, M. (2007). *Generalized neural-network representation of high-dimensional potential-energy surfaces*. Physical Review Letters, 98(14), 146401. [DOI:10.1103/PhysRevLett.98.146401](https://doi.org/10.1103/PhysRevLett.98.146401)
  */
 class Cut_Cos : public Cut_Base {
     private:
@@ -76,21 +96,24 @@ class Cut_Cos : public Cut_Base {
         double calc_prime(double r);
 };
 
-/** @brief Tanh cutoff function.
+/**
+ * @class Cut_Tanh
+ * @brief Hyperbolic tangent cutoff function.
+ *
+ * The `Cut_Tanh` class implements a smooth cutoff function using the hyperbolic tangent, defined by:
  * \f[
  * f_c(r) =
- * \begin{equation}
- *  \begin{cases}
- *   \tanh^3\big( 1 -\frac{r}{r_c} \big)  & \text{if  } r \leq r_c\\
- *   0 & \text{otherwise}\\
- *  \end{cases}
- * \end{equation}
+ * \begin{cases}
+ *  \tanh^3\left( 1 - \dfrac{r}{r_c} \right), & \text{if } r \leq r_c \\
+ *  0, & \text{if } r > r_c
+ * \end{cases}
  * \f]
  *
- * <div class="csl-entry">Behler, J. (2011). Atom-centered symmetry
- * functions for constructing high-dimensional neural network potentials.
- * <i>J. Chem. Phys.</i>, <i>134</i>(7), 074106.
- * https://doi.org/10.1063/1.3553717</div>
+ * This function smoothly transitions from 1 to 0 over the interval \f$ 0 \leq r \leq r_c \f$, with the transition shape controlled by the cubic power of the hyperbolic tangent.
+ *
+ * **Reference:**
+ *
+ * Behler, J. (2011). Atom-centered symmetry functions for constructing high-dimensional neural network potentials. *Journal of Chemical Physics*, 134(7), 074106. [DOI:10.1063/1.3553717](https://doi.org/10.1063/1.3553717)
  */
 class Cut_Tanh : public Cut_Base {
     private:
@@ -107,31 +130,40 @@ class Cut_Tanh : public Cut_Base {
         double calc_prime(double r);
 };
 
-/** @brief Polynomial-2 cutoff function.
+/**
+ * @class Cut_PolyS
+ * @brief Polynomial cutoff function ensuring smoothness up to the second derivative.
  *
+ * The `Cut_PolyS` class implements a polynomial cutoff function defined as:
+
  * \f[
  * f_c(r) =
- * \begin{equation}
- *  \begin{cases}
- *   1 & \text{if  } r \leq (r_c-1)\\
- *   r^3(r(15-6r)-10)+1  & \text{if  } (r_c-1) < r \leq r_c\\
- *   0 & \text{otherwise}\\
- *  \end{cases}
- * \end{equation}
+ * \begin{cases}
+ * 1, & \text{if } r \leq r_{\text{inner}} \\
+ * \left( r - r_c + 1 \right )^3 \left[ r \left( 15 - 6 r \right ) - 10 \right ] + 1, & \text{if } r_{\text{inner}} < r \leq r_c \\
+ * 0, & \text{if } r > r_c
+ * \end{cases}
  * \f]
 
- *<div class="csl-entry">Singraber, A., Rg Behler, J., Dellago, C. (2019).
- Library-Based LAMMPS Implementation of High-Dimensional
- Neural Network Potentials. https://doi.org/10.1021/acs.jctc.8b00770</div>
+ * where:
+ * - \f$ r \f$ is the radial distance.
+ * - \f$ r_c \f$ is the cutoff radius.
+ * - \f$ r_{\text{inner}} = r_c - 1 \f$ is the inner cutoff radius.
+
+ * This function provides a smooth transition from 1 to 0 between \f$ r_{\text{inner}} \f$ and \f$ r_c \f$, ensuring continuity and smoothness in the potential and its derivatives up to the second order.
+
+ * **Reference:**
+
+ * Singraber, A., Behler, J., & Dellago, C. (2019). *Library-Based LAMMPS Implementation of High-Dimensional Neural Network Potentials*. Journal of Chemical Theory and Computation, 15(1), 182–190. [DOI:10.1021/acs.jctc.8b00770](https://doi.org/10.1021/acs.jctc.8b00770)
  */
-class Cut_Poly2 : public Cut_Base {
+class Cut_PolyS : public Cut_Base {
     private:
-        std::string lab = "Cut_Poly2";
+        std::string lab = "Cut_PolyS";
         double rcut, rcut_sq, rcut_inv;
         double rcut_inner;
     public:
-        Cut_Poly2();
-        Cut_Poly2(double rcut);
+        Cut_PolyS();
+        Cut_PolyS(double rcut);
         std::string label() ;
         void set_rcut(const double r);
         double get_rcut();
@@ -139,14 +171,94 @@ class Cut_Poly2 : public Cut_Base {
         double calc(double r);
         double calc_prime(double r);
 };
-class Cut_Cos_S : public Cut_Base {
+/**
+ * @class Cut_Poly
+ * @brief Represents a 5th order polynomial cutoff function with zero first and second derivatives at the cutoff radius.
+ *
+ * This class implements a smooth cutoff function defined by:
+ * \f[
+ * f(r) =
+ * \begin{cases}
+ * 1 - 10 \left( \dfrac{r}{r_{\text{cut}}} \right)^3 + 15 \left( \dfrac{r}{r_{\text{cut}}} \right)^4 - 6 \left( \dfrac{r}{r_{\text{cut}}} \right)^5, & \text{for } r \leq r_{\text{cut}} \\
+ * 0, & \text{for } r > r_{\text{cut}}
+ * \end{cases}
+ * \f]
+ * and its derivative:
+ * \f[
+ * f'(r) =
+ * \begin{cases}
+ * \left( -30 \left( \dfrac{r}{r_{\text{cut}}} \right)^2 + 60 \left( \dfrac{r}{r_{\text{cut}}} \right)^3 - 30 \left( \dfrac{r}{r_{\text{cut}}} \right)^4 \right) \dfrac{1}{r_{\text{cut}}}, & \text{for } r \leq r_{\text{cut}} \\
+ * 0, & \text{for } r > r_{\text{cut}}
+ * \end{cases}
+ * \f]
+ * where:
+ * - \f$ r \f$ is the radial distance.
+ * - \f$ r_{\text{cut}} \f$ is the cutoff radius.
+ *
+ * The function smoothly transitions from 1 to 0 within the cutoff radius \f$ r_{\text{cut}} \f$, ensuring that the function and its first and second derivatives are zero at \f$ r \geq r_{\text{cut}} \f$.
+ *
+ */
+class Cut_Poly : public Cut_Base {
     private:
-        std::string lab = "Cut_Cos_S";
+        std::string lab = "Cut_Poly";
+        double rcut, rcut_sq, rcut_inv;
+    public:
+        Cut_Poly();
+        Cut_Poly(double rcut);
+        std::string label() ;
+        void set_rcut(const double r);
+        double get_rcut();
+        double get_rcut_sq();
+        double calc(double r);
+        double calc_prime(double r);
+};
+/**
+ * @class Cut_CosS
+ * @brief Smooth cosine cutoff function with a smoothing interval.
+ *
+ * The `Cut_CosS` class implements a smooth cosine cutoff function that transitions smoothly from 1 to 0 over a specified smoothing interval between the inner cutoff radius \f$ r_{\text{inner}} \f$ and the outer cutoff radius \f$ r_{\text{cut}} \f$ The function is defined as:
+ *
+ * \f[
+ * f(r) =
+ * \begin{cases}
+ * 1, & \text{if } r \leq r_{\text{inner}} \\
+ * \dfrac{1}{2}\left[1 + \cos\left( \dfrac{\pi (r - r_{\text{inner}})}{r_{\text{cut}} - r_{\text{inner}}} \right) \right], & \text{if } r_{\text{inner}} < r < r_{\text{cut}} \\
+ * 0, & \text{if } r \geq r_{\text{cut}}
+ * \end{cases}
+ * \f]
+ *
+ * and its derivative:
+ *
+ * \f[
+ * f'(r) =
+ * \begin{cases}
+ * 0, & \text{if } r \leq r_{\text{inner}} \text{ or } r \geq r_{\text{cut}} \\
+ * -\dfrac{\pi}{2 (r_{\text{cut}} - r_{\text{inner}})} \sin\left( \dfrac{\pi (r - r_{\text{inner}})}{r_{\text{cut}} - r_{\text{inner}}} \right), & \text{if } r_{\text{inner}} < r < r_{\text{cut}}
+ * \end{cases}
+ * \f]
+ *
+ * where:
+ *
+ * - \f$ r \f$ is the radial distance.
+ * - \f$ r_{\text{cut}} \f$ is the outer cutoff radius (`rcut`).
+ * - \f$ r_{\text{inner}} \f$ is the inner cutoff radius (`rcut_inner`), defined as \f$ r_{\text{inner}} = r_{\text{cut}} - 1.0 \f$.
+ *
+ * **Characteristics:**
+ *
+ * - For \f$ r \leq r_{\text{inner}} \f$ the function \f$ f(r) = 1 \f$ and \f$ f'(r) = 0 \f$
+ * - For \f$ r_{\text{inner}} < r < r_{\text{cut}} \f$, the function transitions smoothly from 1 to 0.
+ * - For \f$ r \geq r_{\text{cut}} \f$, the function \f$ f(r) = 0 \f$ and \f$ f'(r) = 0 \f$.
+ *
+ * **Note:** The function ensures continuity and smoothness in simulations or calculations where a smooth cutoff is required.
+ */
+class Cut_CosS : public Cut_Base {
+    private:
+        std::string lab = "Cut_CosS";
         double rcut, rcut_sq, rcut_inv;
         double rcut_inner;
     public:
-        Cut_Cos_S();
-        Cut_Cos_S(double rcut);
+        Cut_CosS();
+        Cut_CosS(double rcut);
         std::string label() ;
         void set_rcut(const double r);
         double get_rcut();
@@ -154,5 +266,42 @@ class Cut_Cos_S : public Cut_Base {
         double calc(double r);
         double calc_prime(double r);
 };
-//template<> inline Registry<Cut_Base,double>::Map Registry<Cut_Base,double>::registry{};
+/**
+ * @class Cut_PT
+ * @brief Represents a smooth cutoff function using hyperbolic tangent smoothing.
+ *
+ * This class implements a smooth cutoff function defined by:
+ * \f[
+ * f(r) = \dfrac{1}{2} + \dfrac{1}{2} \tanh\left( \dfrac{1}{2} \left( \dfrac{r_c}{r + \dfrac{r_c}{2}} + \dfrac{r_c}{r - r_c} \right) \right)
+ * \f]
+ * and its derivative:
+ * \f[
+ * f'(r) = -\dfrac{r_c}{4} \left( \dfrac{1}{\cosh\left( \dfrac{r_c}{2} \left( \dfrac{1}{r + \dfrac{r_c}{2}} + \dfrac{1}{r - r_c} \right) \right)} \right)^2 \left( \dfrac{1}{\left( r + \dfrac{r_c}{2} \right)^2} + \dfrac{1}{\left( r - r_c \right)^2} \right)
+ * \f]
+ * where:
+ * - \f$ r \f$ is the radial distance.
+ * - \f$ r_c \f$ is the cutoff radius.
+ *
+ * The cutoff function smoothly transitions from 1 to 0 between \f$ r = -\dfrac{r_c}{2} \f$ and \f$ r = r_c \f$.
+ *
+ */
+class Cut_PT : public Cut_Base {
+    private:
+        std::string lab="Cut_PT";
+        double rcut;
+        double rcut_sq;
+        double rcut_inv;
+        double rcut_inner;
+
+    public:
+        Cut_PT();
+        Cut_PT(double rcut_value);
+        std::string label();
+        void set_rcut(const double r);
+        double get_rcut();
+        double get_rcut_sq();
+        double calc(double r);
+        double calc_prime(double r);
+};
+
 #endif
diff --git a/include/tadah/models/descriptors/d_basis_functions.h b/include/tadah/models/descriptors/d_basis_functions.h
index 09d468380c9013e6d68c0ddf762605de64e712ff..dcf14f4e205fcb8516584d3f861f227119a31ba8 100644
--- a/include/tadah/models/descriptors/d_basis_functions.h
+++ b/include/tadah/models/descriptors/d_basis_functions.h
@@ -5,101 +5,104 @@
 
 // GAUSSIANS
 inline double G(double r, double eta, double miu) {
-    return exp(-eta*(r-miu)*(r-miu));
+  return exp(-eta*(r-miu)*(r-miu));
 }
 inline double G(double r, double eta, double miu, double f) {
-    return exp(-eta*(r-miu)*(r-miu))*f;
+  return exp(-eta*(r-miu)*(r-miu))*f;
 }
 inline double dG(double r, double eta, double miu) {
-    return -2.0*eta*(r-miu)*exp(-eta*(r-miu)*(r-miu));
+  return -2.0*eta*(r-miu)*exp(-eta*(r-miu)*(r-miu));
 }
 inline double dG(double r, double eta, double miu, double f, double fp) {
-    return exp(-eta*(r-miu)*(r-miu) )*(-2.0*f*eta*(r-miu)+fp);
+  return exp(-eta*(r-miu)*(r-miu) )*(-2.0*f*eta*(r-miu)+fp);
 }
 
 // BLIPS
 inline double B(double x_c)
-    // x_c = eta*(rij-r_s)
-    // return B(x_c)
+  // r_s = miu
+  // x_c = eta*(rij-r_s)
 {
-    if (-2.0 < x_c && x_c <= -1.0) {
-        double b = 2.0 + x_c;
-        return 0.25 * b*b*b;
-    }
-    else if (-1.0 < x_c && x_c <= 1.0) {
-        double x2 = x_c*x_c;
-        double x3 = x2*x_c;
-        return 1.0 - 1.5*x2 + 0.75*std::fabs(x3);
-    }
-    else if (1.0 < x_c && x_c < 2.0) {
-        double b = 2.0 - x_c;
-        return 0.25 * b*b*b;
-    }
-    else {
-        return 0.0;
-    }
+  if (-2.0 < x_c && x_c <= -1.0) {
+    double b = 2.0 + x_c;
+    return 0.25 * b*b*b;
+  }
+  else if (-1.0 < x_c && x_c <= 1.0) {
+    double x2 = x_c*x_c;
+    double x3 = x2*x_c;
+    return 1.0 - 1.5*x2 + 0.75*std::fabs(x3);
+  }
+  else if (1.0 < x_c && x_c < 2.0) {
+    double b = 2.0 - x_c;
+    return 0.25 * b*b*b;
+  }
+  else {
+    return 0.0;
+  }
 }
 inline double B(double r, double eta, double miu) {
   return B(eta*(r-miu));
 }
 inline double B(double x_c, double f)
-    // x_c = eta*(rij-r_s)
-    // return B(x_c)*f
+  // x_c = eta*(rij-r_s)
+  // return B(x_c)*f
 {
-    if (-2.0 < x_c && x_c <= -1.0) {
-        double b = 2.0 + x_c;
-        return 0.25 * b*b*b *f;
-    }
-    else if (-1.0 < x_c && x_c <= 1.0) {
-        double x2 = x_c*x_c;
-        double x3 = x2*x_c;
-        return f*(1.0 - 1.5*x2 + 0.75*std::fabs(x3));
-    }
-    else if (1.0 < x_c && x_c < 2.0) {
-        double b = 2.0 - x_c;
-        return f*(0.25 * b*b*b);
-    }
-    else {
-        return 0.0;
-    }
+  if (-2.0 < x_c && x_c <= -1.0) {
+    double b = 2.0 + x_c;
+    return 0.25 * b*b*b *f;
+  }
+  else if (-1.0 < x_c && x_c <= 1.0) {
+    double x2 = x_c*x_c;
+    double x3 = x2*x_c;
+    return f*(1.0 - 1.5*x2 + 0.75*std::fabs(x3));
+  }
+  else if (1.0 < x_c && x_c < 2.0) {
+    double b = 2.0 - x_c;
+    return f*(0.25 * b*b*b);
+  }
+  else {
+    return 0.0;
+  }
 }
 inline double dB(double x_c, double eta)
-    // def: x_c = eta*(rij-r_s)
-    // return d/dr_ij B(x_c) = eta*d/dx_c B(x_c)
+  // def: x_c = eta*(rij-r_s)
+  // return d/dr_ij B(x_c) = eta*d/dx_c B(x_c)
 {
-    if (-2.0 < x_c && x_c <= -1.0) {
-        double d = 2.0+x_c;
-        return 0.75*eta*d*d;
-    }
-    else if (-1.0 < x_c && x_c <= 1.0) {
-        return -3.0*eta*x_c + 2.25*eta*x_c*std::fabs(x_c);
-    }
-    else if (1.0 < x_c && x_c < 2.0) {
-        double d = 2.0-x_c;
-        return -0.75*eta*d*d;
-    }
-    else {
-        return 0.0;
-    }
+  if (-2.0 < x_c && x_c <= -1.0) {
+    double d = 2.0+x_c;
+    return 0.75*eta*d*d;
+  }
+  else if (-1.0 < x_c && x_c <= 1.0) {
+    return -3.0*eta*x_c + 2.25*eta*x_c*std::fabs(x_c);
+  }
+  else if (1.0 < x_c && x_c < 2.0) {
+    double d = 2.0-x_c;
+    return -0.75*eta*d*d;
+  }
+  else {
+    return 0.0;
+  }
 }
 inline double dB(double x_c, double eta, double f, double fp)
-    // def: x_c = eta*(rij-r_s)
-    // return d/dr_ij f(r_ij)*B(x_c)
+  // def: x_c = eta*(rij-r_s)
+  // return d/dr_ij f(r_ij)*B(x_c)
 {
-    if (-2.0 < x_c && x_c <= -1.0) {
-        double d = 2.0+x_c;
-        return f*(0.75*eta*d*d)+fp*(0.25 * d*d*d);
-    }
-    else if (-1.0 < x_c && x_c <= 1.0) {
-        return f*(-3.0*eta*x_c + 2.25*eta*x_c*std::fabs(x_c))
-            +fp*(1.0 - 1.5*x_c*x_c + 0.75*std::fabs(x_c*x_c*x_c));
-    }
-    else if (1.0 < x_c && x_c < 2.0) {
-        double d = 2.0-x_c;
-        return f*(-0.75*eta*d*d)+fp*(0.25 * d*d*d);
-    }
-    else {
-        return 0.0;
-    }
+  if (-2.0 < x_c && x_c <= -1.0) {
+    double d = 2.0+x_c;
+    return f*(0.75*eta*d*d)+fp*(0.25 * d*d*d);
+  }
+  else if (-1.0 < x_c && x_c <= 1.0) {
+    return f*(-3.0*eta*x_c + 2.25*eta*x_c*std::fabs(x_c))
+      +fp*(1.0 - 1.5*x_c*x_c + 0.75*std::fabs(x_c*x_c*x_c));
+  }
+  else if (1.0 < x_c && x_c < 2.0) {
+    double d = 2.0-x_c;
+    return f*(-0.75*eta*d*d)+fp*(0.25 * d*d*d);
+  }
+  else {
+    return 0.0;
+  }
+}
+inline double dB(double r, double eta, double miu) {
+  return dB(eta*(r-miu), eta);
 }
 #endif
diff --git a/include/tadah/models/linear_regressor.h b/include/tadah/models/linear_regressor.h
index bcb9eddd24393d03283b3f5b35ce662d3bf7fdf6..ce111bb3d052cec94acb17a5f7ff5e5dde8e60d9 100644
--- a/include/tadah/models/linear_regressor.h
+++ b/include/tadah/models/linear_regressor.h
@@ -35,9 +35,10 @@ class LinearRegressor {
 
       int verbose = config.get<int>("VERBOSE");
       double lambda = config.get<double>("LAMBDA");
+      double rcond = config.size("LAMBDA")==2 ? config.get<double>("LAMBDA",1) : 1e-8;
 
       if (lambda == 0) {
-        OLS::solve(Phi, T, weights);
+        OLS::solve(Phi, T, weights, rcond);
       } else {
         double alpha = config.get<double>("ALPHA");
         double beta = config.get<double>("BETA");
diff --git a/include/tadah/models/m_krr_train.h b/include/tadah/models/m_krr_train.h
index d4e79033377c3051f1cc9c3df054ff2ea439378a..f5b1c8ed3b7494d4d388976d2960c1b4cf27386b 100644
--- a/include/tadah/models/m_krr_train.h
+++ b/include/tadah/models/m_krr_train.h
@@ -155,6 +155,7 @@ public:
         if (M_KRR_Core<Kern>::is_verbose()) std::cout << "Matrix condition number: " << condition_number(K) << std::endl;
 
         if (M_KRR_Core<Kern>::is_verbose()) std::cout << "Solving..." << std::flush;
+        //double rcond = config.template size("LAMBDA")==2 ? config.template get<double>("LAMBDA",1) : 1e-8;
         weights = solve_posv(K, T);
         if (M_KRR_Core<Kern>::is_verbose()) std::cout << "Done" << std::endl;
 
diff --git a/include/tadah/models/ols.h b/include/tadah/models/ols.h
index c507b3466bf2a95c5c5aab317d612793ce7de22f..b95c774060570667a6f072bc4ff26589292ab43f 100644
--- a/include/tadah/models/ols.h
+++ b/include/tadah/models/ols.h
@@ -26,7 +26,7 @@ class OLS {
      * @param weights Output vector containing the computed weights.
      */
     template <typename M, typename V>
-    static void solve(M &A, V &B, V &weights) {
+    static void solve(M &A, V &B, V &weights, double rcond) {
         // Resize B if necessary to match A's column count.
         if (B.size() < A.cols())
             B.resize(A.cols());
@@ -40,7 +40,6 @@ class OLS {
         double *b = B.ptr();
         int ldb = std::max(m, n);
         double *s = new double[std::min(m, n)]; // Singular values
-        double rcond = 1e-8; // Condition for singularity
         int rank;
         double *work;
         int lwork = -1; // Workspace query
diff --git a/include/tadah/models/ridge_regression.h b/include/tadah/models/ridge_regression.h
index 4eca71146254cd3603064ee002cb7af7edd547e4..f2d6691ddea72382c9b9392c9fb04a48f7d8d31d 100644
--- a/include/tadah/models/ridge_regression.h
+++ b/include/tadah/models/ridge_regression.h
@@ -46,7 +46,7 @@ class RidgeRegression {
      * @param lambda Regularization parameter.
      */
     template <typename V, typename W>
-    static void solve(const SVD &svd, V b, W &weights, const double lambda) {
+    static void solve(const SVD &svd, V b, W &weights, double lambda) {
       double *U = svd.getU(); // Matrix U from SVD (m x m)
       double *s = svd.getS(); // Singular values (as a vector)
       double *VT = svd.getVT(); // Matrix V^T from SVD (n x n)
diff --git a/src/cutoffs.cpp b/src/cutoffs.cpp
index b818cdab76a7d70389bed0b258fca95b2de6751c..7b17adbe57251e4ff0881e4e4765f99d9945626a 100644
--- a/src/cutoffs.cpp
+++ b/src/cutoffs.cpp
@@ -5,10 +5,12 @@
 
 template<> CONFIG::Registry<Cut_Base,double>::Map CONFIG::Registry<Cut_Base,double>::registry{};
 CONFIG::Registry<Cut_Base,double>::Register<Cut_Cos> Cut_Cos_1( "Cut_Cos" );
-CONFIG::Registry<Cut_Base,double>::Register<Cut_Cos_S> Cut_Cos_S_1( "Cut_Cos_S" );
-CONFIG::Registry<Cut_Base,double>::Register<Cut_Tanh> Cut_Tanh_1( "Cut_Tanh" );
-CONFIG::Registry<Cut_Base,double>::Register<Cut_Poly2> Cut_Poly_1( "Cut_Poly2" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_CosS> Cut_Cos_2( "Cut_CosS" );
 CONFIG::Registry<Cut_Base,double>::Register<Cut_Dummy> Cut_Dummy_1( "Cut_Dummy" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_Poly> Cut_Poly_1( "Cut_Poly" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_PolyS> Cut_Poly_2( "Cut_PolyS" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_PT> Cut_PT_1( "Cut_PT" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_Tanh> Cut_Tanh_1( "Cut_Tanh" );
 
 
 Cut_Base::~Cut_Base() {}
@@ -91,7 +93,6 @@ Cut_Tanh::Cut_Tanh() {}
 Cut_Tanh::Cut_Tanh(double rcut)
 {
     set_rcut(rcut);
-    test_rcut(rcut);
 }
 
 std::string Cut_Tanh::label() {
@@ -125,81 +126,81 @@ double Cut_Tanh::calc_prime(double r) {
     return -3.0*t*t*rcut_inv*s*s;
 }
 
-Cut_Poly2::Cut_Poly2() {}
-Cut_Poly2::Cut_Poly2(double rcut)
+Cut_PolyS::Cut_PolyS() {}
+Cut_PolyS::Cut_PolyS(double rcut)
 {
     if (rcut<=1.0)
         throw std::runtime_error(
-                "Cut_Poly2 cutoff distance cannot be smaller than 1.0.");
+                "Cut_PolyS cutoff distance cannot be smaller than 1.0.");
     set_rcut(rcut);
     test_rcut(rcut);
 }
 
-std::string Cut_Poly2::label() {
+std::string Cut_PolyS::label() {
     return lab;
 }
 
-void Cut_Poly2::set_rcut(const double r) {
+void Cut_PolyS::set_rcut(const double r) {
     test_rcut(r);
     rcut=r;
     rcut_sq=r*r;
     rcut_inv= r<=0 ? 0.0 : 1.0/r;
     rcut_inner=rcut-1.0;
 }
-double Cut_Poly2::get_rcut() {
+double Cut_PolyS::get_rcut() {
     return rcut;
 }
 
-double Cut_Poly2::get_rcut_sq() {
+double Cut_PolyS::get_rcut_sq() {
     return rcut_sq;
 }
 
-double Cut_Poly2::calc(double r) {
+double Cut_PolyS::calc(double r) {
     if (r>=rcut) return 0.0;
     else if (r<= rcut_inner) return 1.0;
     double rs=r-rcut_inner;
     return rs*rs*rs*(rs*(15.0-6.0*rs)-10.0)+1;
 }
-double Cut_Poly2::calc_prime(double r) {
+double Cut_PolyS::calc_prime(double r) {
     if (r>=rcut) return 0.0;
     else if (r<= rcut_inner) return 0.0;
     double rs=r-rcut_inner;
     return -30.0*(rs-1.0)*(rs-1.0)*rs*rs;
 }
 
-Cut_Cos_S::Cut_Cos_S() {}
-Cut_Cos_S::Cut_Cos_S(double rcut)
+Cut_CosS::Cut_CosS() {}
+Cut_CosS::Cut_CosS(double rcut)
 {
     set_rcut(rcut);
     test_rcut(rcut);
 }
 
-std::string Cut_Cos_S::label() {
+std::string Cut_CosS::label() {
     return lab;
 }
 
-void Cut_Cos_S::set_rcut(const double r) {
+void Cut_CosS::set_rcut(const double r) {
     test_rcut(r);
     rcut=r;
     rcut_sq=r*r;
     rcut_inv= r<=0 ? 0.0 : 1.0/r;
     rcut_inner=rcut-1.0;
 }
-double Cut_Cos_S::get_rcut() {
+double Cut_CosS::get_rcut() {
     return rcut;
 }
 
-double Cut_Cos_S::get_rcut_sq() {
+double Cut_CosS::get_rcut_sq() {
     return rcut_sq;
 }
 
-double Cut_Cos_S::calc(double r) {
+double Cut_CosS::calc(double r) {
   if (r>=rcut) return 0.0;
   else if (r<= rcut_inner) return 1.0;
   double rs = (r-rcut_inner)/(rcut-rcut_inner);
   return 0.5*(1+std::cos(M_PI*rs));
 }
-double Cut_Cos_S::calc_prime(double r) {
+double Cut_CosS::calc_prime(double r) {
   if (r>=rcut || r<= rcut_inner) return 0.0;
   else {
     double rs = (r - rcut_inner) / (rcut - rcut_inner);
@@ -207,3 +208,85 @@ double Cut_Cos_S::calc_prime(double r) {
     return -0.5 * M_PI * std::sin(M_PI * rs) * drs_dr;
   }
 }
+
+Cut_PT::Cut_PT() {}
+Cut_PT::Cut_PT(double rcut) {
+  set_rcut(rcut);
+  test_rcut(rcut);
+}
+std::string Cut_PT::label() {
+  return lab;
+}
+
+// Set the outer cutoff radius
+void Cut_PT::set_rcut(const double r) {
+    test_rcut(r);
+    rcut=r;
+    rcut_sq=r*r;
+    rcut_inv= r<=0 ? 0.0 : 1.0/r;
+    rcut_inner=-rcut/2;
+}
+double Cut_PT::get_rcut() {
+  return rcut;
+}
+double Cut_PT::get_rcut_sq() {
+  return rcut_sq;
+}
+double Cut_PT::calc(double r) {
+  if (r>=rcut) return 0.0;
+  double u = 0.5 * ((rcut / (r - rcut_inner)) + (rcut / (r - rcut)));
+  return 0.5 + 0.5 * tanh(u);
+}
+double Cut_PT::calc_prime(double r) {
+  if (r>=rcut) return 0.0;
+  double u = 0.5 * ((rcut / (r - rcut_inner)) + (rcut / (r - rcut)));
+  double cosh_u = cosh(u);
+  double sech2_u = 1.0 / (cosh_u * cosh_u);
+
+  double term1 = 1.0 / ((r - rcut_inner) * (r - rcut_inner));
+  double term2 = 1.0 / ((r - rcut) * (r - rcut));
+  double sum_terms = term1 + term2;
+
+  return - (rcut / 4.0) * sech2_u * sum_terms;
+}
+Cut_Poly::Cut_Poly() {}
+Cut_Poly::Cut_Poly(double rcut_value) {
+  set_rcut(rcut_value);
+}
+std::string Cut_Poly::label() {
+    return lab;
+}
+void Cut_Poly::set_rcut(const double r) {
+    test_rcut(r);
+    rcut = r;
+    rcut_sq = r * r;
+    rcut_inv= r<=0 ? 0.0 : 1.0/r;
+}
+double Cut_Poly::get_rcut() {
+    return rcut;
+}
+double Cut_Poly::get_rcut_sq() {
+    return rcut_sq;
+}
+double Cut_Poly::calc(double r) {
+    if (r >= rcut) {
+        return 0.0;
+    } else {
+        double x = r * rcut_inv; // x = r / rcut
+        double x3 = x * x * x;
+        double x4 = x3 * x;
+        double x5 = x4 * x;
+        return 1.0 - 10.0 * x3 + 15.0 * x4 - 6.0 * x5;
+    }
+}
+double Cut_Poly::calc_prime(double r) {
+    if (r >= rcut) {
+        return 0.0;
+    } else {
+        double x = r * rcut_inv; // x = r / rcut
+        double x2 = x * x;
+        double x3 = x2 * x;
+        double x4 = x3 * x;
+        return (-30.0 * x2 + 60.0 * x3 - 30.0 * x4) * rcut_inv;
+    }
+}
diff --git a/src/d3_base.cpp b/src/d3_base.cpp
index fe5fa6932ef38b78dac3fd4d1bdf9f5f990629df..5c8eec24e88307c9d7da32411ffa7e68200c893b 100644
--- a/src/d3_base.cpp
+++ b/src/d3_base.cpp
@@ -1,6 +1,6 @@
 #include <tadah/models/descriptors/d3/d3_base.h>
 //template struct Registry<D3_Base, Config &>;
-std::vector<std::string> D3_Base::get_init_atoms(Config &c) {
+std::vector<std::string> D3_Base::get_init_atoms(Config &) {
   return {};
 }
 
diff --git a/src/kern_base.cpp b/src/kern_base.cpp
index 1a431efa800ebe9bbf85bfd11f4f5b1c712a4cfe..25a15f442fbecc8e698b40710743f3dcb87c709e 100644
--- a/src/kern_base.cpp
+++ b/src/kern_base.cpp
@@ -36,18 +36,18 @@ force_type Kern_Base::fpredict(const t_type &kweights, const fd_type &fdij,
   }
   return v;
 }
-void Kern_Base::read_basis_from_config(const Config &config, Matrix &b) {
+void Kern_Base::read_basis_from_config(const Config &config, Matrix &basis) {
   // storage is in column-major order
   if (!(config.exist("SBASIS") && config.exist("BASIS"))) return;
   size_t sbasis = config.get<size_t>("SBASIS");
   size_t vlen = config("BASIS").size()/sbasis;
   std::vector<double> temp(sbasis*vlen);
   config.get<std::vector<double>>("BASIS",temp);
-  b.resize(vlen,sbasis);  // aeds are in column-major order
+  basis.resize(vlen,sbasis);  // aeds are in column-major order
   size_t ii=0;
   for (size_t i=0;i<sbasis; ++i) {
     for (size_t j=0;j<vlen; ++j) {
-      b(j,i) = temp[ii++];
+      basis(j,i) = temp[ii++];
     }
   }
 }
diff --git a/tests/test_cutoffs.cpp b/tests/test_cutoffs.cpp
index aaac9d0561436505138407296c7b83e4475939c8..1127f2871276b77e4248a3a5693be8c8a371d0db 100644
--- a/tests/test_cutoffs.cpp
+++ b/tests/test_cutoffs.cpp
@@ -119,13 +119,13 @@ TEST_CASE( "Testing Cutoffs: Cut_Tanh", "[Cut_Tanh]" ) {
 
     delete c2b;
 }
-TEST_CASE( "Testing Cutoffs: Cut_Poly2", "[Cut_Poly2]" ) {
+TEST_CASE( "Testing Cutoffs: Cut_PolyS", "[Cut_PolyS]" ) {
 
-    REQUIRE_NOTHROW(Cut_Poly2());
+    REQUIRE_NOTHROW(Cut_PolyS());
     double rcut2b = 6.2;
     double rcut2bsq = rcut2b*rcut2b;
-    using Cut = Cut_Poly2;
-    std::string cuttype = "Cut_Poly2";
+    using Cut = Cut_PolyS;
+    std::string cuttype = "Cut_PolyS";
     Cut_Base *c2b = new Cut( rcut2b );
 
     REQUIRE( c2b->label() == cuttype );
@@ -162,13 +162,13 @@ TEST_CASE( "Testing Cutoffs: Cut_Poly2", "[Cut_Poly2]" ) {
                Catch::Matchers::WithinAbs(-1.728, 1e-12));
     delete c2b;
 }
-TEST_CASE( "Testing Cutoffs: Cut_Cos_s", "[Cut_Cos_S]" ) {
+TEST_CASE( "Testing Cutoffs: Cut_CosS", "[Cut_CosS]" ) {
 
-    REQUIRE_NOTHROW(Cut_Cos_S());
+    REQUIRE_NOTHROW(Cut_CosS());
     double rcut2b = 6.2;
     double rcut2bsq = rcut2b*rcut2b;
-    using Cut = Cut_Cos_S;
-    std::string cuttype = "Cut_Cos_S";
+    using Cut = Cut_CosS;
+    std::string cuttype = "Cut_CosS";
     Cut_Base *c2b = new Cut( rcut2b );
 
     REQUIRE( c2b->label() == cuttype );
@@ -205,3 +205,46 @@ TEST_CASE( "Testing Cutoffs: Cut_Cos_s", "[Cut_Cos_S]" ) {
                Catch::Matchers::WithinAbs(-1.4939160824, 1e-10));
     delete c2b;
 }
+TEST_CASE( "Testing Cutoffs: Cut_PT", "[Cut_PT]" ) {
+
+    REQUIRE_NOTHROW(Cut_PT());
+    double rcut2b = 6.2;
+    double rcut2bsq = rcut2b*rcut2b;
+    using Cut = Cut_PT;
+    std::string cuttype = "Cut_PT";
+    Cut_Base *c2b = new Cut( rcut2b );
+
+    REQUIRE( c2b->label() == cuttype );
+
+    REQUIRE( c2b->calc(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( c2b->calc_prime(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut()-rcut2b)<std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut_sq()-rcut2bsq)<std::numeric_limits<double>::min() );
+
+    // cutoff cannot be negative
+    double temp = -0.1;
+    REQUIRE_THROWS(Cut( temp ));
+    REQUIRE_THROWS_AS(c2b->set_rcut(temp), std::runtime_error);
+    REQUIRE_THROWS_AS(c2b->test_rcut(temp), std::runtime_error);
+
+    // recheck after resetting cutoff
+    rcut2b=3.4;
+    rcut2bsq=rcut2b*rcut2b;
+    c2b->set_rcut(100000);
+    c2b->set_rcut(rcut2b);
+    REQUIRE( c2b->calc(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( c2b->calc_prime(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut()-rcut2b)<std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut_sq()-rcut2bsq)<std::numeric_limits<double>::min() );
+
+  REQUIRE_THAT(c2b->calc(2.0),
+               Catch::Matchers::WithinAbs(0.180990296913, 1e-12));
+  REQUIRE_THAT(c2b->calc_prime(2.0),
+               Catch::Matchers::WithinAbs(-0.293953123268, 1e-12));
+
+  REQUIRE_THAT(c2b->calc(3.0),
+               Catch::Matchers::WithinAbs(0.000419261765536, 1e-10));
+  REQUIRE_THAT(c2b->calc_prime(3.0),
+               Catch::Matchers::WithinAbs(-0.00897008113779, 1e-10));
+    delete c2b;
+}
diff --git a/tests/test_d2.cpp b/tests/test_d2.cpp
index 84faae1792def1f4ed44f0dec4f6747da653d05a..d703644e9c4f5d7c37e2e45d55a321d3977deabd 100644
--- a/tests/test_d2.cpp
+++ b/tests/test_d2.cpp
@@ -75,7 +75,7 @@ void test_D2(Config &c) {
       test_d2<D,Cut_Dummy>(c,r);
       test_d2<D,Cut_Cos>(c,r);
       test_d2<D,Cut_Tanh>(c,r);
-      test_d2<D,Cut_Poly2>(c,r);
+      test_d2<D,Cut_PolyS>(c,r);
     }
   }
   Config c2; // INIT2B=false
diff --git a/tests/test_dm.cpp b/tests/test_dm.cpp
index 1800b4c1e71cb3c58f56ddea67cc212bb0111df0..1e93f0418e0c095ed4259cf2d2b9c1e20889a5a8 100644
--- a/tests/test_dm.cpp
+++ b/tests/test_dm.cpp
@@ -315,7 +315,7 @@ TEST_CASE_METHOD(TestFixture, "Testing DM_EAD 1", "[DM]") {
     test_case<DM_EAD, Cut_Dummy>(config_mb, "DM_EAD", 3 * 2);
     test_case<DM_EAD, Cut_Cos>(config_mb, "DM_EAD", 3 * 2);
     test_case<DM_EAD, Cut_Tanh>(config_mb, "DM_EAD", 3 * 2);
-    test_case<DM_EAD, Cut_Poly2>(config_mb, "DM_EAD", 3 * 2);
+    test_case<DM_EAD, Cut_PolyS>(config_mb, "DM_EAD", 3 * 2);
 }
 
 TEST_CASE_METHOD(TestFixture, "Testing DM_EAD 2", "[DM]") {
@@ -332,7 +332,7 @@ TEST_CASE_METHOD(TestFixture, "Testing DM_EAD 2", "[DM]") {
     test_case<DM_EAD, Cut_Dummy>(config_mb, "DM_EAD", 4 * 3);
     test_case<DM_EAD, Cut_Cos>(config_mb, "DM_EAD", 4 * 3);
     test_case<DM_EAD, Cut_Tanh>(config_mb, "DM_EAD", 4 * 3);
-    test_case<DM_EAD, Cut_Poly2>(config_mb, "DM_EAD", 4 * 3);
+    test_case<DM_EAD, Cut_PolyS>(config_mb, "DM_EAD", 4 * 3);
 }
 
 TEST_CASE_METHOD(TestFixture, "Testing DM_Blip 1", "[DM]") {
@@ -349,7 +349,7 @@ TEST_CASE_METHOD(TestFixture, "Testing DM_Blip 1", "[DM]") {
     test_case<DM_Blip, Cut_Dummy>(config_mb, "DM_Blip", 3 * 2);
     test_case<DM_Blip, Cut_Cos>(config_mb, "DM_Blip", 3 * 2);
     test_case<DM_Blip, Cut_Tanh>(config_mb, "DM_Blip", 3 * 2);
-    test_case<DM_Blip, Cut_Poly2>(config_mb, "DM_Blip", 3 * 2);
+    test_case<DM_Blip, Cut_PolyS>(config_mb, "DM_Blip", 3 * 2);
 }
 
 TEST_CASE_METHOD(TestFixture, "Testing DM_Blip 2", "[DM]") {
@@ -366,7 +366,7 @@ TEST_CASE_METHOD(TestFixture, "Testing DM_Blip 2", "[DM]") {
     test_case<DM_Blip, Cut_Dummy>(config_mb, "DM_Blip", 4 * 3);
     test_case<DM_Blip, Cut_Cos>(config_mb, "DM_Blip", 4 * 3);
     test_case<DM_Blip, Cut_Tanh>(config_mb, "DM_Blip", 4 * 3);
-    test_case<DM_Blip, Cut_Poly2>(config_mb, "DM_Blip", 4 * 3);
+    test_case<DM_Blip, Cut_PolyS>(config_mb, "DM_Blip", 4 * 3);
 }
 
 TEST_CASE_METHOD(TestFixture, "Testing DM_EAM 1", "[DM]") {
diff --git a/tests/test_factory_cutoffs.cpp b/tests/test_factory_cutoffs.cpp
index 2b19927b0722a5797634de870585a235bfedd7cd..5ac5b8db8fff30d55efb7eeca9bf511c481aaec4 100644
--- a/tests/test_factory_cutoffs.cpp
+++ b/tests/test_factory_cutoffs.cpp
@@ -2,40 +2,60 @@
 #include <limits>
 #include <stdexcept>
 #include <tadah/core/registry.h>
+#include <tadah/core/utils.h>
 #include <tadah/models/cutoffs.h>
 #include <string>
 
 TEST_CASE( "Testing Factory: Cutoffs", "[factory_cutoffs]" ) {
 
-    double rcut2b = 6.2;
-    double rcut2bsq = rcut2b*rcut2b;
-
-    for (auto c:CONFIG::Registry<Cut_Base,double>::registry) {
-        std::string cuttype = c.first;
-        Cut_Base *c2b = CONFIG::factory<Cut_Base,double>( cuttype, rcut2b );
-
-        REQUIRE( c2b->label() == cuttype );
-
-        REQUIRE( c2b->calc(rcut2b) < std::numeric_limits<double>::min() );
-        REQUIRE( c2b->calc_prime(rcut2b) < std::numeric_limits<double>::min() );
-        REQUIRE( std::abs(c2b->get_rcut()-rcut2b)<std::numeric_limits<double>::min() );
-        REQUIRE( std::abs(c2b->get_rcut_sq()-rcut2bsq)<std::numeric_limits<double>::min() );
-
-        // cutoff cannot be negative
-        double temp = -0.1;
-        REQUIRE_THROWS(CONFIG::factory<Cut_Base,double>( cuttype, temp ));
-        REQUIRE_THROWS_AS(c2b->set_rcut(temp), std::runtime_error);
-        REQUIRE_THROWS_AS(c2b->test_rcut(temp), std::runtime_error);
-
-        // recheck after resetting cutoff
-        rcut2b=3.4;
-        rcut2bsq=rcut2b*rcut2b;
-        c2b->set_rcut(100000);
-        c2b->set_rcut(rcut2b);
-        REQUIRE( c2b->calc(rcut2b) < std::numeric_limits<double>::min() );
-        REQUIRE( c2b->calc_prime(rcut2b) < std::numeric_limits<double>::min() );
-        REQUIRE( std::abs(c2b->get_rcut()-rcut2b)<std::numeric_limits<double>::min() );
-        REQUIRE( std::abs(c2b->get_rcut_sq()-rcut2bsq)<std::numeric_limits<double>::min() );
-        if (c2b) delete c2b;
-    };
+  double rcut2b = 6.2;
+  double rcut2bsq = rcut2b*rcut2b;
+
+  for (auto c:CONFIG::Registry<Cut_Base,double>::registry) {
+    std::string cuttype = c.first;
+    Cut_Base *c2b = CONFIG::factory<Cut_Base,double>( cuttype, rcut2b );
+
+    REQUIRE( c2b->label() == cuttype );
+
+    REQUIRE( c2b->calc(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( c2b->calc_prime(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut()-rcut2b)<std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut_sq()-rcut2bsq)<std::numeric_limits<double>::min() );
+
+    // cutoff cannot be negative
+    double temp = -0.1;
+    REQUIRE_THROWS(CONFIG::factory<Cut_Base,double>( cuttype, temp ));
+    REQUIRE_THROWS_AS(c2b->set_rcut(temp), std::runtime_error);
+    REQUIRE_THROWS_AS(c2b->test_rcut(temp), std::runtime_error);
+
+    // recheck after resetting cutoff
+    rcut2b=3.4;
+    rcut2bsq=rcut2b*rcut2b;
+    c2b->set_rcut(100000);
+    c2b->set_rcut(rcut2b);
+    REQUIRE( c2b->calc(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( c2b->calc_prime(rcut2b) < std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut()-rcut2b)<std::numeric_limits<double>::min() );
+    REQUIRE( std::abs(c2b->get_rcut_sq()-rcut2bsq)<std::numeric_limits<double>::min() );
+    if (c2b) delete c2b;
+  };
+}
+TEST_CASE( "Testing Factory: Cutoff's Derivatives", "[factory_cutoffs]" ) {
+
+  double rcut2b = 6.2;
+  std::vector<double> rs = {0.0, 1.0, 3.3, 4.234, rcut2b};
+
+  for (auto c:CONFIG::Registry<Cut_Base,double>::registry) {
+    std::string cuttype = c.first;
+    if (cuttype == "Cut_Dummy") continue; // no central difference for this one
+    Cut_Base *c2b = CONFIG::factory<Cut_Base,double>( cuttype, rcut2b );
+
+    REQUIRE( c2b->label() == cuttype );
+    for (const double r : rs) {
+      double prime_approx = central_difference([&](double x){ return c2b->calc(x); }, r, 1e-8);
+      REQUIRE( std::abs(c2b->calc_prime(r)-prime_approx)<1e-6);
+    }
+
+    if (c2b) delete c2b;
+  };
 }
diff --git a/tests/test_ols.cpp b/tests/test_ols.cpp
index 8e6e63ef051166372c659f795660ea477a1b22d4..f84c6b08d78400a6db196a83180f40cf8315b2fd 100644
--- a/tests/test_ols.cpp
+++ b/tests/test_ols.cpp
@@ -14,7 +14,7 @@ TEST_CASE("Testing OLS") {
   aed_type w(3);
   aed_type b2=b1;
 
-  OLS::solve(Phi1,b1,w);
+  OLS::solve(Phi1,b1,w, 1e-8);
   aed_type p= Phi2*w;
   REQUIRE_THAT(p[0], Catch::Matchers::WithinRel(b2[0]));
   REQUIRE_THAT(p[1], Catch::Matchers::WithinRel(b2[1]));