Implemented short-distance suppression of symmetry functions

.gitignore

@@ -1,3 +1,8 @@
 *.o
 compile_commands.json
 .clangd
+lib/
+bin/
+include/
+application/
+*.so

src/libnnp/CutoffFunction.cpp

@@ -84,6 +84,13 @@ void CutoffFunction::setCutoffType(CutoffType const cutoffType)
          dfPtr = &CutoffFunction:: dfCORE;
         fdfPtr = &CutoffFunction::fdfCORE;
     }
+    else if (cutoffType == CT_ICOS)
+    {
+          fPtr = &CutoffFunction::  fICOS;
+         dfPtr = &CutoffFunction:: dfICOS;
+        fdfPtr = &CutoffFunction::fdfICOS;
+    }
+
     else
     {
         throw invalid_argument("ERROR: Unknown cutoff type.\n");
@@ -111,6 +118,13 @@ void CutoffFunction::setCutoffParameter(double const alpha)
     return;
 }
 
+void CutoffFunction::setInnerCutoffParameter(double const icut_beta,
+					     double const icut_gamma)
+{
+  this->icut_beta = icut_beta;
+  this->icut_gamma = icut_gamma;
+  return;
+}
 double CutoffFunction::fCOS(double r) const
 {
     if (r < rci) return 1.0;
@@ -211,3 +225,64 @@ void CutoffFunction::fdfCORE(double r, double& fc, double& dfc) const
     dfc *= iw;
     return;
 }
+
+double CutoffFunction::fICOS(double r) const
+{
+  double const rii = this->icut_beta * rc; 
+  double const rio = this->icut_gamma * rc; 
+  if (r < rii)
+    return 0.0;
+  else if (r >= rii && r <= rio) {
+    double const x = (r - rii) / (rio - rii);
+    return 0.5 * (-cos(PI * x) + 1.0);
+  }
+  else if (r > rio && r < rci)
+    return 1.0;
+  else {
+    double const x = (r - rci) / (rc - rci);
+    return 0.5 * (cos(PI * x) + 1.0);
+  }
+}
+
+double CutoffFunction::dfICOS(double r) const
+{
+  double const rii = this->icut_beta * rc; 
+  double const rio = this->icut_gamma * rc; 
+  if (r < rii)
+    return 0.0;
+  else if (r >= rii && r <= rio) {
+    double const x = (r - rii) / (rio - rii);
+    return 0.5 * PI * sin(PI * x) / (rio - rii);
+  }
+  else if (r > rio && r < rci)
+    return 0.0;
+  else {
+    double const x = (r - rci) / (rc - rci);
+    return -0.5 * PI * sin(PI * x) / (rc - rci);
+  }
+}
+
+void CutoffFunction::fdfICOS(double r, double& fc, double& dfc) const
+{
+  double const rii = this->icut_beta * rc; 
+  double const rio = this->icut_gamma * rc;
+  if (r < rii) {
+    fc = 0.0;
+    dfc = 0.0;
+  }
+  else if (r >= rii && r <= rio) {
+    double const x = (r - rii) / (rio - rii);
+    fc = 0.5 * (-cos(PI * x) + 1.0);
+    dfc = 0.5 * PI * sin(PI * x) / (rio - rii);
+  }
+  else if (r > rio && r < rci){
+    fc = 1.0;
+    dfc = 0.0;
+  }
+  else if (r >= rci){
+    double const x = (r - rci) / (rc - rci);
+    fc = 0.5 * (cos(PI * x) + 1.0);
+    dfc = -0.5 * PI * sin(PI * x) / (rc - rci);
+  }
+  return;
+}

src/libnnp/CutoffFunction.h

@@ -68,7 +68,10 @@ class CutoffFunction
         CT_POLY3,
         /** @f$f(x) = (x(x((315 - 70x)x - 540) + 420) - 126)x^5 + 1@f$
          */
-        CT_POLY4
+        CT_POLY4,
+	/** @f$f(x) = 1/2(-\cos(2\pi x) + 1)@f$
+         */
+	CT_ICOS
     };
 
     /** Constructor, initializes to ´CT_HARD´.
@@ -99,6 +102,14 @@ class CutoffFunction
      * @param[in] alpha Width parameter @f$\alpha@f$.
      */
     void       setCutoffParameter(double const alpha);
+  /* Set parameter for inner cutoff
+   * 
+   * @param[in] icut_beta inner cutoff function parameter icut_beta
+   * @param[in] icut_gamma inner cutoff function parameter icut_gamma
+   */
+  void setInnerCutoffParameter(double const icut_beta,
+			       double const icut_gamma);
+
     /** Getter for #alpha.
      *
      * @return Cutoff parameter used.
@@ -145,6 +156,11 @@ class CutoffFunction
     double       iw;
     /// Core functions used by POLYN, if any.
     CoreFunction core;
+  /// inner cutoff function parameter beta
+  double icut_beta;
+  /// inner cutoff function parameter gamma
+  double icut_gamma;
+
     /// Function pointer to f.
     double       (CutoffFunction::*fPtr)(double r) const;
     /// Function pointer to df.
@@ -174,6 +190,10 @@ class CutoffFunction
     double   fCORE (double r) const;
     double  dfCORE (double r) const;
     void   fdfCORE (double r, double& fc, double& dfc) const;
+
+    double  fICOS  (double r) const;
+    double dfICOS  (double r) const;
+    void  fdfICOS  (double r, double& fc, double& dfc) const;
 };
 
 //////////////////////////////////

src/libnnp/Element.cpp

@@ -339,6 +339,23 @@ void Element::setCutoffFunction(CutoffFunction::CutoffType const cutoffType,
     return;
 }
 
+void Element::setCutoffFunction(CutoffFunction::CutoffType const cutoffType,
+                                double const                     cutoffAlpha,
+				double const icut_beta,
+				double const icut_gamma)
+{
+  setCutoffFunction(cutoffType, cutoffAlpha);
+    for (vector<SymFnc*>::const_iterator
+         it = symmetryFunctions.begin(); it != symmetryFunctions.end(); ++it)
+    {
+        SymFncBaseCutoff* sfcb = dynamic_cast<SymFncBaseCutoff*>(*it);
+        if (sfcb != nullptr)
+        {
+	  sfcb->setInnerCutoffFunction(icut_beta, icut_gamma);
+        }
+    }
+
+}
 void Element::setScalingNone() const
 {
     for (size_t i = 0; i < symmetryFunctions.size(); ++i)

src/libnnp/Element.h

@@ -124,6 +124,17 @@ class Element
     void                     setCutoffFunction(
                                  CutoffFunction::CutoffType const cutoffType,
                                  double const                     cutoffAlpha);
+    /** Set cutoff function for all symmetry functions.
+     *
+     * @param[in] cutoffType Type of cutoff function.
+     * @param[in] cutoffAlpha Cutoff parameter for all functions.
+     */
+    void                     setCutoffFunction(
+                                 CutoffFunction::CutoffType const cutoffType,
+                                 double const                     cutoffAlpha,
+				 double const icut_beta,
+				 double const icut_gamma);
+
     /** Set no scaling of symmetry function.
      *
      * Still scaling factors need to be initialized.

src/libnnp/Mode.cpp

@@ -421,6 +421,19 @@ void Mode::setupCutoff()
         log << "f(r) = 1\n";
         log << "WARNING: Hard cutoff used!\n";
     }
+    else if (cutoffType == CutoffFunction::CT_ICOS)
+    {
+        log << strpr("CutoffFunction::CT_ICOS (%d)\n", cutoffType);
+        log << "x := (r - rc * alpha) / (rc - rc * alpha)\n";
+        log << "f(x) = 1/2 * (-cos(2*pi*x) + 1)\n";
+
+	if (settings.keywordExists("inner_cutoff"))
+	  {
+	    vector<string> icut_args = split(settings["inner_cutoff"]);
+	    icut_beta  = atof(icut_args.at(0).c_str());
+	    icut_gamma = atof(icut_args.at(1).c_str());
+	  }
+    }
     else
     {
         throw invalid_argument("ERROR: Unknown cutoff type.\n");
@@ -474,8 +487,13 @@ void Mode::setupSymmetryFunctions()
     {
         if (normalize) it->changeLengthUnitSymmetryFunctions(convLength);
         it->sortSymmetryFunctions();
-        maxCutoffRadius = max(it->getMaxCutoffRadius(), maxCutoffRadius);
-        it->setCutoffFunction(cutoffType, cutoffAlpha);
+	maxCutoffRadius = max(it->getMaxCutoffRadius(), maxCutoffRadius);
+
+	if (settings.keywordExists("inner_cutoff")) {
+	  it->setCutoffFunction(cutoffType, cutoffAlpha, icut_beta, icut_gamma);
+	} else {
+	  it->setCutoffFunction(cutoffType, cutoffAlpha);
+	}
         log << strpr("Short range atomic symmetry functions element %2s :\n",
                      it->getSymbol().c_str());
         log << "--------------------------------------------------"

src/libnnp/Mode.h

@@ -527,6 +527,9 @@ class Mode
     CutoffFunction::CutoffType cutoffType;
     std::vector<Element>       elements;
 
+  double icut_beta;
+  double icut_gamma;
+
     /** Read in weights for a specific type of neural network.
      *
      * @param[in] type Actual network type to initialize ("short" or "charge").

src/libnnp/Settings.cpp

@@ -37,6 +37,7 @@ map<string, shared_ptr<Settings::Key>> const createKnownKeywordsMap()
     m["elements"                      ] = "";
     m["atom_energy"                   ] = "";
     m["cutoff_type"                   ] = "";
+    m["inner_cutoff"                  ] = ""; 
     m["symfunction_short"             ] = "";
     m["scale_symmetry_functions"      ] = "";
     m["scale_min_short"               ] = "";

src/libnnp/SymFncBaseCutoff.cpp

@@ -48,6 +48,15 @@ void SymFncBaseCutoff::setCutoffFunction(
     return;
 }
 
+void SymFncBaseCutoff::setInnerCutoffFunction(double const icut_beta,
+					      double const icut_gamma)
+{
+  this->icut_beta = icut_beta;
+  this->icut_gamma = icut_gamma;
+  fc.setInnerCutoffParameter(icut_beta, icut_gamma);
+}
+
+
 SymFncBaseCutoff::SymFncBaseCutoff(size_t type,
                                      ElementMap const& elementMap) :
     SymFnc(type, elementMap),

src/libnnp/SymFncBaseCutoff.h

@@ -45,6 +45,14 @@ class SymFncBaseCutoff : public SymFnc
     void                     setCutoffFunction(CutoffFunction::
                                           CutoffType cutoffType,
                                           double     cutoffAlpha);
+    /** Set cutoff function type and parameter.
+     *
+     * @param[in] icut_beta Inner cutoff function parameter beta
+     * @param[in] icut_gamma Inner cutoff function parameter gamma
+     */
+    void                     setInnerCutoffFunction(double icut_beta,
+						    double icut_gamma);
+
     /** Get private #cutoffAlpha member variable.
      */
     double                   getCutoffAlpha() const;
@@ -56,6 +64,9 @@ class SymFncBaseCutoff : public SymFnc
     CutoffFunction::
     CutoffType               getCutoffType() const;
 
+  double geticut_beta() const;
+  double geticut_gamma() const;
+
 protected:
     /// Cutoff parameter @f$\alpha@f$.
     double                     cutoffAlpha;
@@ -66,6 +77,8 @@ class SymFncBaseCutoff : public SymFnc
     /// Cutoff type used by this symmetry function.
     CutoffFunction::CutoffType cutoffType;
 
+  double icut_beta;
+  double icut_gamma;
     /** Constructor, initializes #type.
      */
     SymFncBaseCutoff(std::size_t type, ElementMap const&);
@@ -82,7 +95,11 @@ inline CutoffFunction::CutoffType SymFncBaseCutoff::getCutoffType() const
 {
     return cutoffType;
 }
+
+inline double SymFncBaseCutoff::geticut_beta() const { return icut_beta; }
+inline double SymFncBaseCutoff::geticut_gamma() const { return icut_gamma; }
 
 }
 
+
 #endif

src/libnnp/SymGrpBaseCutoff.h

@@ -37,6 +37,8 @@ class SymGrpBaseCutoff : public SymGrp
     double                     rc;
     /// Cutoff function parameter @f$\alpha@f$ (common feature).
     double                     cutoffAlpha;
+  double icut_beta;
+  double icut_gamma;
     /// Subtype string (specifies cutoff type) (common feature).
     std::string                subtype;
     /// Cutoff function used by this symmetry function group.

src/libnnp/SymGrpExpAngn.cpp

@@ -82,10 +82,13 @@ bool SymGrpExpAngn::addMember(SymFnc const* const symmetryFunction)
         e1          = sf->getE1();
         e2          = sf->getE2();
         convLength  = sf->getConvLength();
+	icut_beta   = sf->geticut_beta();
+	icut_gamma  = sf->geticut_gamma(); 
 
         fc.setCutoffType(cutoffType);
         fc.setCutoffRadius(rc);
         fc.setCutoffParameter(cutoffAlpha);
+	fc.setInnerCutoffParameter(icut_beta, icut_gamma);
     }
 
     if (sf->getCutoffType()  != cutoffType ) return false;

src/libnnp/SymGrpExpRad.cpp

@@ -85,10 +85,13 @@ bool SymGrpExpRad::addMember(SymFnc const* const symmetryFunction)
         rc          = sf->getRc();
         e1          = sf->getE1();
         convLength  = sf->getConvLength();
+	icut_beta   = sf->geticut_beta();
+	icut_gamma  = sf->geticut_gamma(); 
 
         fc.setCutoffType(cutoffType);
         fc.setCutoffRadius(rc);
         fc.setCutoffParameter(cutoffAlpha);
+	fc.setInnerCutoffParameter(icut_beta, icut_gamma);
     }
 
     if (sf->getCutoffType()  != cutoffType ) return false;

src/makefile

@@ -13,13 +13,13 @@ PROJECT_LIB=$(PROJECT_DIR)/lib
 # GENERAL SETTINGS
 ###############################################################################
 # Default compiler (gnu/intel/llvm).
-COMP=gnu
+COMP=intel
 
 # Default build mode.
 # Possible modes are "static", "shared" and "test".
 # Note: For target "pynnp" there is no "static" mode, "shared" is automatically
 # used instead.
-MODE=static
+MODE=shared
 
 # Installation directory for binaries.
 INSTALL_BIN=$(HOME)/local/bin

src/makefile.gnu

@@ -5,8 +5,8 @@
 ###############################################################################
 # Enter here paths to GSL or EIGEN if they are not in your standard include
 # path. DO NOT completely remove the entry, leave at least "./".
-PROJECT_GSL=./
-PROJECT_EIGEN=/usr/include/eigen3/
+PROJECT_GSL=.//jet/home/xyttyxyx/selfcompiled-programs/lib/
+PROJECT_EIGEN=/jet/home/xyttyxyx/selfcompiled-programs/compile/eigen-3.3.9/
 
 ###############################################################################
 # COMPILERS AND FLAGS

src/makefile.intel

@@ -12,7 +12,7 @@ PROJECT_EIGEN=/usr/include/eigen3/
 # COMPILERS AND FLAGS
 ###############################################################################
 PROJECT_CC=icpc
-PROJECT_MPICC=mpiicpc
+PROJECT_MPICC=mpicxx
 # OpenMP parallelization is disabled by default, add flag "-qopenmp" to enable.
 PROJECT_CFLAGS=-O3 -xHost -std=c++11 -ipo
 PROJECT_CFLAGS_MPI=-Wno-long-long
@@ -41,7 +41,7 @@ PROJECT_LDFLAGS_BLAS=-lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
 #PROJECT_OPTIONS+= -DN2P2_NO_TIME
 
 # Disable check for low number of neighbors.
-#PROJECT_OPTIONS+= -DN2P2_NO_NEIGH_CHECK
+PROJECT_OPTIONS+= -DN2P2_NO_NEIGH_CHECK
 
 # Use alternative (older) memory layout for symmetry function derivatives.
 #PROJECT_OPTIONS+= -DN2P2_FULL_SFD_MEMORY