Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "ManyRestraintsBase.h"
      23             : #include "core/ActionRegister.h"
      24             : 
      25             : namespace PLMD {
      26             : namespace manyrestraints {
      27             : 
      28             : //+PLUMEDOC MCOLVARB LWALLS
      29             : /*
      30             : Add \ref LOWER_WALLS restraints on all the multicolvar values
      31             : 
      32             : This action takes the set of values calculated by the colvar specified by label in the DATA
      33             : keyword and places a restraint on each quantity, \f$x\f$, with the following functional form:
      34             : 
      35             : \f$
      36             :   k((x-a+o)/s)^e
      37             : \f$
      38             : 
      39             : \f$k\f$ (KAPPA) is an energy constant in internal unit of the code, \f$s\f$ (EPS) a rescaling factor and
      40             : \f$e\f$ (EXP) the exponent determining the power law. By default: EXP = 2, EPS = 1.0, OFF = 0.
      41             : 
      42             : \par Examples
      43             : 
      44             : The following set of commands can be used to stop any of the 800 atoms in group A from moving more than 2.46425 nm
      45             : in the z direction from atom 34137.  This is done by adding a lower wall on the z-distance between all the atoms in
      46             : group A and the position of 34137.
      47             : 
      48             : \plumedfile
      49             : l: ZDISTANCES GROUPA=1-800 GROUPB=34137 NOPBC
      50             : LWALLS DATA=l AT=2.46465 KAPPA=150.0 EXP=2 EPS=1 OFFSET=0 LABEL=lwall
      51             : \endplumedfile
      52             : 
      53             : 
      54             : */
      55             : //+ENDPLUMEDOC
      56             : 
      57             : 
      58             : class LWalls : public ManyRestraintsBase {
      59             : private:
      60             :   double at;
      61             :   double kappa;
      62             :   double exp;
      63             :   double eps;
      64             :   double offset;
      65             : public:
      66             :   static void registerKeywords( Keywords& keys );
      67             :   explicit LWalls( const ActionOptions& );
      68             :   double calcPotential( const double& val, double& df ) const override;
      69             : };
      70             : 
      71       10419 : PLUMED_REGISTER_ACTION(LWalls,"LWALLS")
      72             : 
      73           1 : void LWalls::registerKeywords( Keywords& keys ) {
      74           1 :   ManyRestraintsBase::registerKeywords( keys );
      75           2 :   keys.add("compulsory","AT","the radius of the sphere");
      76           2 :   keys.add("compulsory","KAPPA","the force constant for the wall.  The k_i in the expression for a wall.");
      77           2 :   keys.add("compulsory","OFFSET","0.0","the offset for the start of the wall.  The o_i in the expression for a wall.");
      78           2 :   keys.add("compulsory","EXP","2.0","the powers for the walls.  The e_i in the expression for a wall.");
      79           2 :   keys.add("compulsory","EPS","1.0","the values for s_i in the expression for a wall");
      80           1 : }
      81             : 
      82           0 : LWalls::LWalls(const ActionOptions& ao):
      83             :   Action(ao),
      84           0 :   ManyRestraintsBase(ao)
      85             : {
      86           0 :   parse("AT",at);
      87           0 :   parse("OFFSET",offset);
      88           0 :   parse("EPS",eps);
      89           0 :   parse("EXP",exp);
      90           0 :   parse("KAPPA",kappa);
      91           0 :   checkRead();
      92           0 : }
      93             : 
      94           0 : double LWalls::calcPotential( const double& val, double& df ) const {
      95           0 :   double uscale = (val - at + offset)/eps;
      96           0 :   if( uscale < 0. ) {
      97           0 :     double power = pow( uscale, exp );
      98           0 :     df = ( kappa / eps ) * exp * power / uscale;
      99             : 
     100           0 :     return kappa*power;
     101             :   }
     102             : 
     103             :   return 0.0;
     104             : }
     105             : 
     106             : }
     107             : }
     108             :