Line data    Source code

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