Proposing: compacting lower and upper walls into a single file

src/bias/LWalls.cpp → src/bias/WallsScalar.cpp

@@ -25,6 +25,52 @@
 namespace PLMD {
 namespace bias {
 
+//+PLUMEDOC BIAS UPPER_WALLS
+/*
+Defines a wall for the value of one or more collective variables,
+ which limits the region of the phase space accessible during the simulation.
+
+The restraining potential starts acting on the system when the value of the CV is greater
+(in the case of UPPER_WALLS) or lower (in the case of LOWER_WALLS) than a certain limit \f$a_i\f$ (AT)
+minus an offset \f$o_i\f$ (OFFSET).
+The expression for the bias due to the wall is given by:
+
+\f$
+  \sum_i {k_i}((x_i-a_i+o_i)/s_i)^e_i
+\f$
+
+\f$k_i\f$ (KAPPA) is an energy constant in internal unit of the code, \f$s_i\f$ (EPS) a rescaling factor and
+\f$e_i\f$ (EXP) the exponent determining the power law. By default: EXP = 2, EPS = 1.0, OFFSET = 0.
+
+
+\par Examples
+
+The following input tells plumed to add both a lower and an upper walls on the distance
+between atoms 3 and 5 and the distance between atoms 2 and 4. The lower and upper limits
+are defined at different values. The strength of the walls is the same for the four cases.
+It also tells plumed to print the energy of the walls.
+\plumedfile
+DISTANCE ATOMS=3,5 LABEL=d1
+DISTANCE ATOMS=2,4 LABEL=d2
+UPPER_WALLS ARG=d1,d2 AT=1.0,1.5 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=uwall
+LOWER_WALLS ARG=d1,d2 AT=0.0,1.0 KAPPA=150.0,150.0 EXP=2,2 EPS=1,1 OFFSET=0,0 LABEL=lwall
+PRINT ARG=uwall.bias,lwall.bias
+\endplumedfile
+(See also \ref DISTANCE and \ref PRINT).
+
+*/
+//+ENDPLUMEDOC
+
+//+PLUMEDOC BIAS UPPER_WALLS_SCALAR
+/*
+Defines a wall for the value of one or more collective variables,
+ which limits the region of the phase space accessible during the simulation.
+
+\par Examples
+
+*/
+//+ENDPLUMEDOC
+
 //+PLUMEDOC BIAS LOWER_WALLS
 /*
 Defines a wall for the value of one or more collective variables,
@@ -71,22 +117,58 @@ Defines a wall for the value of one or more collective variables,
 */
 //+ENDPLUMEDOC
 
-class LWalls : public Bias {
+enum class WallType {UPPER,LOWER};
+
+template<WallType W>
+class WallsScalar : public Bias {
   std::vector<double> at;
   std::vector<double> kappa;
   std::vector<double> exp;
   std::vector<double> eps;
   std::vector<double> offset;
+  static constexpr auto force2="force2";
+  inline double walloff(double cv, double off) {
+    if constexpr(W==WallType::UPPER) {
+      return cv+off;
+    } else {
+      return cv-off;
+    }
+  }
+  inline double wallpow(double scale, double exponent) {
+    if constexpr(W==WallType::UPPER) {
+      return std::pow(scale,exponent);
+    } else {
+      return std::pow(-scale,exponent);
+    }
+  }
+  inline bool check(double scale) {
+    if constexpr(W==WallType::UPPER) {
+      return scale > 0.0;
+    } else {
+      return scale < 0.0;
+    }
+  }
 public:
-  explicit LWalls(const ActionOptions&);
+  explicit WallsScalar(const ActionOptions&);
   void calculate() override;
   static void registerKeywords(Keywords& keys);
 };
 
+using UWalls = WallsScalar<WallType::UPPER>;
+PLUMED_REGISTER_ACTION(UWalls,"UPPER_WALLS_SCALAR")
+
+using LWalls = WallsScalar<WallType::LOWER>;
 PLUMED_REGISTER_ACTION(LWalls,"LOWER_WALLS_SCALAR")
 
-void LWalls::registerKeywords(Keywords& keys) {
-  Bias::registerKeywords(keys); keys.setDisplayName("LOWER_WALLS");
+template<WallType W>
+void WallsScalar<W>::registerKeywords(Keywords& keys) {
+  Bias::registerKeywords(keys);
+  if constexpr(W==WallType::UPPER) {
+    keys.setDisplayName("UPPER_WALLS");
+  } else {
+    keys.setDisplayName("LOWER_WALLS");
+  }
+
   keys.use("ARG"); keys.add("hidden","NO_ACTION_LOG","suppresses printing from action on the log");
   keys.add("compulsory","AT","the positions of the wall. The a_i in the expression for a wall.");
   keys.add("compulsory","KAPPA","the force constant for the wall.  The k_i in the expression for a wall.");
@@ -96,7 +178,8 @@ void LWalls::registerKeywords(Keywords& keys) {
   keys.addOutputComponent("force2","default","the instantaneous value of the squared force due to this bias potential");
 }
 
-LWalls::LWalls(const ActionOptions&ao):
+template<WallType W>
+WallsScalar<W>::WallsScalar(const ActionOptions&ao):
   PLUMED_BIAS_INIT(ao),
   at(getNumberOfArguments(),0),
   kappa(getNumberOfArguments(),0.0),
@@ -128,30 +211,33 @@ LWalls::LWalls(const ActionOptions&ao):
   for(unsigned i=0; i<eps.size(); i++) log.printf(" %f",eps[i]);
   log.printf("\n");
 
-  addComponent("force2"); componentIsNotPeriodic("force2");
+  addComponent(force2); componentIsNotPeriodic(force2);
 }
 
-void LWalls::calculate() {
+
+
+template<WallType W>
+void WallsScalar<W>::calculate() {
   double ene = 0.0;
   double totf2 = 0.0;
   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     double f = 0.0;
     const double cv=difference(i,at[i],getArgument(i));
     const double off=offset[i];
     const double epsilon=eps[i];
-    const double lscale = (cv-off)/epsilon;
-    if( lscale < 0.) {
+    const double scale = walloff(cv,off)/epsilon;
+    if( check(scale) ) {
       const double k=kappa[i];
       const double exponent=exp[i];
-      double power = std::pow( lscale, exponent );
-      f = -( k / epsilon ) * exponent * power / lscale;
+      double power = wallpow( scale, exponent );
+      f = -( k / epsilon ) * exponent * power / scale;
       ene += k * power;
       totf2 += f * f;
     }
     setOutputForce(i,f);
   }
   setBias(ene);
-  getPntrToComponent("force2")->set(totf2);
+  getPntrToComponent(force2)->set(totf2);
 }
 
 }