From 0dc4bd448175da4db2232f404ebb1e4194421784 Mon Sep 17 00:00:00 2001
From: mkirsz <s1351949@sms.ed.ac.uk>
Date: Fri, 7 Feb 2025 22:52:25 +0000
Subject: [PATCH] New Cos_PT cutoff
---
include/tadah/models/cutoffs.h | 40 +++++++++++++++++++++++++++++++
src/cutoffs.cpp | 42 +++++++++++++++++++++++++++++++++
tests/test_cutoffs.cpp | 43 ++++++++++++++++++++++++++++++++++
3 files changed, 125 insertions(+)
diff --git a/include/tadah/models/cutoffs.h b/include/tadah/models/cutoffs.h
index 882f0c9..a961c6f 100644
--- a/include/tadah/models/cutoffs.h
+++ b/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
diff --git a/src/cutoffs.cpp b/src/cutoffs.cpp
index b818cda..72ebb6f 100644
--- a/src/cutoffs.cpp
+++ b/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;
+}
diff --git a/tests/test_cutoffs.cpp b/tests/test_cutoffs.cpp
index aaac9d0..b90be3b 100644
--- a/tests/test_cutoffs.cpp
+++ b/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;
+}
--
GitLab