Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2014-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 UWALLS
      29             : /*
      30             : Add \ref UPPER_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 a cluster composed of 20 atoms subliming.  The position of
      45             : the center of mass of the cluster is calculated by the \ref COM command labelled c1.  The \ref DISTANCES
      46             : command labelled d1 is then used to calculate the distance between each of the 20 atoms in the cluster
      47             : and the center of mass of the cluster.  These distances are then passed to the UWALLS command, which adds
      48             : a \ref UPPER_WALLS restraint on each of them and thereby prevents each of them from moving very far from the center
      49             : of mass of the cluster.
      50             : 
      51             : \plumedfile
      52             : COM ATOMS=1-20 LABEL=c1
      53             : DISTANCES GROUPA=c1 GROUPB=1-20 LABEL=d1
      54             : UWALLS DATA=d1 AT=2.5 KAPPA=0.2 LABEL=sr
      55             : \endplumedfile
      56             : 
      57             : 
      58             : */
      59             : //+ENDPLUMEDOC
      60             : 
      61             : 
      62             : class UWalls : public ManyRestraintsBase {
      63             : private:
      64             :   double at;
      65             :   double kappa;
      66             :   double exp;
      67             :   double eps;
      68             :   double offset;
      69             : public:
      70             :   static void registerKeywords( Keywords& keys );
      71             :   explicit UWalls( const ActionOptions& );
      72             :   double calcPotential( const double& val, double& df ) const override;
      73             : };
      74             : 
      75       10423 : PLUMED_REGISTER_ACTION(UWalls,"UWALLS")
      76             : 
      77           3 : void UWalls::registerKeywords( Keywords& keys ) {
      78           3 :   ManyRestraintsBase::registerKeywords( keys );
      79           6 :   keys.add("compulsory","AT","the radius of the sphere");
      80           6 :   keys.add("compulsory","KAPPA","the force constant for the wall.  The k_i in the expression for a wall.");
      81           6 :   keys.add("compulsory","OFFSET","0.0","the offset for the start of the wall.  The o_i in the expression for a wall.");
      82           6 :   keys.add("compulsory","EXP","2.0","the powers for the walls.  The e_i in the expression for a wall.");
      83           6 :   keys.add("compulsory","EPS","1.0","the values for s_i in the expression for a wall");
      84           3 : }
      85             : 
      86           2 : UWalls::UWalls(const ActionOptions& ao):
      87             :   Action(ao),
      88           2 :   ManyRestraintsBase(ao)
      89             : {
      90           2 :   parse("AT",at);
      91           2 :   parse("OFFSET",offset);
      92           2 :   parse("EPS",eps);
      93           2 :   parse("EXP",exp);
      94           2 :   parse("KAPPA",kappa);
      95           2 :   checkRead();
      96           2 : }
      97             : 
      98       14800 : double UWalls::calcPotential( const double& val, double& df ) const {
      99       14800 :   double uscale = (val - at + offset)/eps;
     100       14800 :   if( uscale > 0. ) {
     101        1036 :     double power = pow( uscale, exp );
     102        1036 :     df = ( kappa / eps ) * exp * power / uscale;
     103             : 
     104        1036 :     return kappa*power;
     105             :   }
     106             : 
     107             :   return 0.0;
     108             : }
     109             : 
     110             : }
     111             : }
     112             :