New Cos_PT cutoff

include/tadah/models/cutoffs.h

@@ -154,5 +154,45 @@ class Cut_Cos_S : public Cut_Base {
         double calc(double r);
         double calc_prime(double r);
 };
+/**
+ * @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) = \frac{1}{2} + \frac{1}{2} \tanh\left( \frac{1}{2} \left( \frac{b}{r - a} + \frac{b}{r - R} \right) \right)
+ * \f]
+ * and its derivative:
+ * \f[
+ * f'(r) = -\frac{b}{4} \left( \frac{1}{\cosh\left( \dfrac{b}{2} \left( \dfrac{1}{r - a} + \dfrac{1}{r - R} \right) \right)} \right)^2 \left( \dfrac{1}{(r - a)^2} + \dfrac{1}{(r - R)^2} \right)
+ * \f]
+ * where:
+ * - \( r \) is the radial distance.
+ * - \( a \) is the inner cutoff radius (`rcut_inner`).
+ * - \( R \) is the outer cutoff radius (`rcut`).
+ * - \( b \) is the smoothing parameter.
+ *
+ * The function smoothly transitions from 0 to 1 between the inner cutoff radius \( a \) and the outer cutoff radius \( R \).
+ *
+ */
+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);
+};
+
 //template<> inline Registry<Cut_Base,double>::Map Registry<Cut_Base,double>::registry{};
 #endif

src/cutoffs.cpp

@@ -9,6 +9,7 @@ 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_Dummy> Cut_Dummy_1( "Cut_Dummy" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_PT> Cut_PT_1( "Cut_PT" );
 
 
 Cut_Base::~Cut_Base() {}
@@ -207,3 +208,44 @@ 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;
+}

tests/test_cutoffs.cpp

@@ -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;
+}