Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2012-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 "Bias.h"
      23             : #include "core/ActionRegister.h"
      24             : #include "tools/Grid.h"
      25             : #include "tools/Exception.h"
      26             : #include "tools/File.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace bias {
      30             : 
      31             : //+PLUMEDOC BIAS EXTERNAL
      32             : /*
      33             : Calculate a restraint that is defined on a grid that is read during start up
      34             : 
      35             : \par Examples
      36             : 
      37             : The following is an input for a calculation with an external potential that is
      38             : defined in the file bias.dat and that acts on the distance between atoms 3 and 5.
      39             : \plumedfile
      40             : DISTANCE ATOMS=3,5 LABEL=d1
      41             : EXTERNAL ARG=d1 FILE=bias.grid LABEL=external
      42             : \endplumedfile
      43             : 
      44             : The bias.grid will then look something like this:
      45             : \auxfile{bias.grid}
      46             : #! FIELDS d1 external.bias der_d1
      47             : #! SET min_d1 1.14
      48             : #! SET max_d1 1.32
      49             : #! SET nbins_d1 6
      50             : #! SET periodic_d1 false
      51             :    1.1400   0.0031   0.1101
      52             :    1.1700   0.0086   0.2842
      53             :    1.2000   0.0222   0.6648
      54             :    1.2300   0.0521   1.4068
      55             :    1.2600   0.1120   2.6873
      56             :    1.2900   0.2199   4.6183
      57             :    1.3200   0.3948   7.1055
      58             : \endauxfile
      59             : 
      60             : This should then be followed by the value of the potential and its derivative
      61             : at 100 equally spaced points along the distance between 0 and 1.
      62             : 
      63             : You can also include grids that are a function of more than one collective
      64             : variable.  For instance the following would be the input for an external
      65             : potential acting on two torsional angles:
      66             : \plumedfile
      67             : TORSION ATOMS=4,5,6,7 LABEL=t1
      68             : TORSION ATOMS=6,7,8,9 LABEL=t2
      69             : EXTERNAL ARG=t1,t2 FILE=bias2.grid LABEL=ext
      70             : \endplumedfile
      71             : 
      72             : The file bias2.grid for this calculation would need to look something like this:
      73             : \auxfile{bias2.grid}
      74             : #! FIELDS t1 t2 ext.bias der_t1 der_t2
      75             : #! SET min_t1 -pi
      76             : #! SET max_t1 pi
      77             : #! SET nbins_t1 3
      78             : #! SET periodic_t1 true
      79             : #! SET min_t2 -pi
      80             : #! SET max_t2 pi
      81             : #! SET nbins_t2 3
      82             : #! SET periodic_t2 true
      83             :  -3.141593 -3.141593 0.000000 -0.000000 -0.000000
      84             :  -1.047198 -3.141593 0.000000 0.000000 -0.000000
      85             :  1.047198 -3.141593 0.000000 -0.000000 -0.000000
      86             : 
      87             :  -3.141593 -1.047198 0.000000 -0.000000 0.000000
      88             :  -1.047198 -1.047198 0.007922 0.033185 0.033185
      89             :  1.047198 -1.047198 0.007922 -0.033185 0.033185
      90             : 
      91             :  -3.141593 1.047198 0.000000 -0.000000 -0.000000
      92             :  -1.047198 1.047198 0.007922 0.033185 -0.033185
      93             :  1.047198 1.047198 0.007922 -0.033185 -0.033185
      94             : \endauxfile
      95             : 
      96             : This would be then followed by 100 blocks of data.  In the first block of data the
      97             : value of t1 (the value in the first column) is kept fixed and the value of
      98             : the function is given at 100 equally spaced values for t2 between \f$-pi\f$ and \f$+pi\f$.  In the
      99             : second block of data t1 is fixed at \f$-pi + \frac{2pi}{100}\f$ and the value of the function is
     100             : given at 100 equally spaced values for t2 between \f$-pi\f$ and \f$+pi\f$. In the third block of
     101             : data the same is done but t1 is fixed at \f$-pi + \frac{4pi}{100}\f$ and so on until you get to
     102             : the one hundredth block of data where t1 is fixed at \f$+pi\f$.
     103             : 
     104             : Please note the order that the order of arguments in the plumed.dat file must be the same as
     105             : the order of arguments in the header of the grid file.
     106             : */
     107             : //+ENDPLUMEDOC
     108             : 
     109             : class External : public Bias {
     110             : 
     111             : private:
     112             :   std::unique_ptr<GridBase> BiasGrid_;
     113             :   double scale_;
     114             : 
     115             : public:
     116             :   explicit External(const ActionOptions&);
     117             :   void calculate() override;
     118             :   static void registerKeywords(Keywords& keys);
     119             : };
     120             : 
     121             : PLUMED_REGISTER_ACTION(External,"EXTERNAL")
     122             : 
     123           4 : void External::registerKeywords(Keywords& keys) {
     124           4 :   Bias::registerKeywords(keys);
     125           4 :   keys.use("ARG");
     126           8 :   keys.add("compulsory","FILE","the name of the file containing the external potential.");
     127           8 :   keys.addFlag("NOSPLINE",false,"specifies that no spline interpolation is to be used when calculating the energy and forces due to the external potential");
     128           8 :   keys.addFlag("SPARSE",false,"specifies that the external potential uses a sparse grid");
     129           8 :   keys.add("compulsory","SCALE","1.0","a factor that multiplies the external potential, useful to invert free energies");
     130           4 : }
     131             : 
     132           2 : External::External(const ActionOptions& ao):
     133           2 :   PLUMED_BIAS_INIT(ao)
     134             : {
     135             :   std::string filename;
     136           4 :   parse("FILE",filename);
     137           2 :   if( filename.length()==0 ) error("No external potential file was specified");
     138           2 :   bool sparsegrid=false;
     139           2 :   parseFlag("SPARSE",sparsegrid);
     140           2 :   bool nospline=false;
     141           2 :   parseFlag("NOSPLINE",nospline);
     142           2 :   bool spline=!nospline;
     143           3 :   parse("SCALE",scale_);
     144             : 
     145           2 :   checkRead();
     146             : 
     147           2 :   log.printf("  External potential from file %s\n",filename.c_str());
     148           2 :   log.printf("  Multiplied by %lf\n",scale_);
     149           2 :   if(spline) {log.printf("  External potential uses spline interpolation\n");}
     150           2 :   if(sparsegrid) {log.printf("  External potential uses sparse grid\n");}
     151             : 
     152             : // read grid
     153           2 :   IFile gridfile; gridfile.open(filename);
     154           1 :   std::string funcl=getLabel() + ".bias";
     155           2 :   BiasGrid_=GridBase::create(funcl,getArguments(),gridfile,sparsegrid,spline,true);
     156           1 :   if(BiasGrid_->getDimension()!=getNumberOfArguments()) error("mismatch between dimensionality of input grid and number of arguments");
     157           3 :   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     158           4 :     if( getPntrToArgument(i)->isPeriodic()!=BiasGrid_->getIsPeriodic()[i] ) error("periodicity mismatch between arguments and input bias");
     159             :   }
     160           6 : }
     161             : 
     162           5 : void External::calculate()
     163             : {
     164           5 :   unsigned ncv=getNumberOfArguments();
     165           5 :   std::vector<double> cv(ncv), der(ncv);
     166             : 
     167          15 :   for(unsigned i=0; i<ncv; ++i) {cv[i]=getArgument(i);}
     168             : 
     169           5 :   double ene=scale_*BiasGrid_->getValueAndDerivatives(cv,der);
     170             : 
     171           5 :   setBias(ene);
     172             : 
     173          15 :   for(unsigned i=0; i<ncv; ++i) {
     174          10 :     const double f=-scale_*der[i];
     175          10 :     setOutputForce(i,f);
     176             :   }
     177           5 : }
     178             : 
     179             : }
     180             : }