I have put the switchingFunc in Tools and added a comment in Coordination

on how to call it. Would it be better to separate parsing tools and Math
tools?

src/ColvarCoordination.cpp

@@ -134,6 +134,10 @@ void ColvarCoordination::calculate()
   }
 
 // we should define switching functions and derivatives in Tools 
+// Carlo: I put an implementation of the switching function in Tools
+// It can be called as
+  double dfunc=0.;
+  ncoord += Tools::switchingFunc(distance.modulo(), nn, mm, r_0, d_0, &dfunc);
 
   const double rdist = (distance.modulo()-d_0)/r_0;
   double dfunc=0.;

src/Tools.cpp

@@ -158,6 +158,27 @@ void Tools::interpretRanges(std::vector<std::string>&s){
   s=news;
 }
 
-
-
-
+double Tools::switchingFunc(const double distance, const int nn, const int mm, const double r_0, const double d_0, double *dfunc){
+  double result=0.;
+  const double threshold=pow(0.00001,1./(nn-mm));
+  const double rdist = (distance-d_0)/r_0;
+
+  if(rdist<=0.){
+     result+=1.;
+     *dfunc=0.;
+  }else if(rdist>(1.-epsilon) && rdist<(1+epsilon)){
+     result+=nn/mm;
+     *dfunc=0.5*nn*(nn-mm)/mm;
+  }else if(rdist>threshold){
+     *dfunc=0.;
+  }else{
+     double rNdist = pow(rdist, nn-1);
+     double rMdist = pow(rdist, mm-1);
+     double num = 1.-rNdist*rdist;
+     double iden = 1./(1.-rMdist*rdist);
+     double func = num*iden;
+     result += func;
+     *dfunc = ((-nn*rNdist*iden)+(func*(iden*mm)*rMdist))/(distance*r_0);
+  }
+  return result;
+}

src/Tools.h

@@ -58,6 +58,8 @@ public:
   static bool parseFlag(std::vector<std::string>&line,const std::string&key,bool&val);
 /// Interpret atom ranges
   static void interpretRanges(std::vector<std::string>&);
+/// Switching function
+  static double switchingFunc(const double, const int, const int, const double, const double, double *);
 };
 
 template <class T>