From 97ffe36f9a11ae58f90b2485b3b350afa6470e8e Mon Sep 17 00:00:00 2001
From: mkirsz <s1351949@sms.ed.ac.uk>
Date: Sat, 8 Feb 2025 08:37:22 +0000
Subject: [PATCH] Poly1 cutoff

---
 include/tadah/models/cutoffs.h | 145 +++++++++++++++++++++++++--------
 src/cutoffs.cpp                |  45 +++++++++-
 2 files changed, 152 insertions(+), 38 deletions(-)

diff --git a/include/tadah/models/cutoffs.h b/include/tadah/models/cutoffs.h
index a961c6f..55ef5b0 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:
+ * - \( r \) is the radial distance.
+ * - \( r_{\text{cut}} \) 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 \( r = r_{\text{cut}} \), where it is undefined due to the discontinuity.
+ *
+ * **Note:** Since the function is discontinuous at \( r = r_{\text{cut}} \), its derivative is not defined at that point. In practical computations, the derivative function `calc_prime` returns `0.0` for all \( r \).
  */
 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 \( 0 \leq r \leq r_c \).
+ * 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 \( 0 \leq r \leq r_c \), with the transition shape controlled by the cubic power of the hyperbolic tangent. It is particularly useful in constructing atom-centered symmetry functions for high-dimensional neural network potentials.
+ *
+ * **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,22 +130,31 @@ class Cut_Tanh : public Cut_Base {
         double calc_prime(double r);
 };
 
-/** @brief Polynomial-2 cutoff function.
+/**
+ * @class Cut_Poly2
+ * @brief Polynomial cutoff function ensuring smoothness up to the second derivative.
  *
+ * The `Cut_Poly2` 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:
+ * - \( r \) is the radial distance.
+ * - \( r_c \) is the cutoff radius.
+ * - \( r_{\text{inner}} = r_c - 1 \) is the inner cutoff radius.
+
+ * This function provides a smooth transition from 1 to 0 between \( r_{\text{inner}} \) and \( r_c \), 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 {
     private:
@@ -139,6 +171,47 @@ class Cut_Poly2 : public Cut_Base {
         double calc(double r);
         double calc_prime(double r);
 };
+/**
+ * @class Cut_Poly1
+ * @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:
+ * - \( r \) is the radial distance.
+ * - \( r_{\text{cut}} \) is the cutoff radius.
+ *
+ * The function smoothly transitions from 1 to 0 within the cutoff radius \( r_{\text{cut}} \), ensuring that the function and its first and second derivatives are zero at \( r \geq r_{\text{cut}} \).
+ *
+ */
+class Cut_Poly1 : public Cut_Base {
+    private:
+        std::string lab = "Cut_Poly1";
+        double rcut, rcut_sq, rcut_inv;
+    public:
+        Cut_Poly1();
+        Cut_Poly1(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_Cos_S : public Cut_Base {
     private:
         std::string lab = "Cut_Cos_S";
diff --git a/src/cutoffs.cpp b/src/cutoffs.cpp
index 72ebb6f..82033cf 100644
--- a/src/cutoffs.cpp
+++ b/src/cutoffs.cpp
@@ -7,7 +7,8 @@ template<> CONFIG::Registry<Cut_Base,double>::Map CONFIG::Registry<Cut_Base,doub
 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_Poly1> Cut_Poly_1( "Cut_Poly1" );
+CONFIG::Registry<Cut_Base,double>::Register<Cut_Poly2> Cut_Poly_2( "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" );
 
@@ -92,7 +93,6 @@ Cut_Tanh::Cut_Tanh() {}
 Cut_Tanh::Cut_Tanh(double rcut)
 {
     set_rcut(rcut);
-    test_rcut(rcut);
 }
 
 std::string Cut_Tanh::label() {
@@ -249,3 +249,44 @@ double Cut_PT::calc_prime(double r) {
 
   return - (rcut / 4.0) * sech2_u * sum_terms;
 }
+Cut_Poly1::Cut_Poly1() {}
+Cut_Poly1::Cut_Poly1(double rcut_value) {
+  set_rcut(rcut_value);
+}
+std::string Cut_Poly1::label() {
+    return lab;
+}
+void Cut_Poly1::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_Poly1::get_rcut() {
+    return rcut;
+}
+double Cut_Poly1::get_rcut_sq() {
+    return rcut_sq;
+}
+double Cut_Poly1::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_Poly1::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;
+    }
+}
-- 
GitLab