Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2021 The VES code team
       3             :    (see the PEOPLE-VES file at the root of this folder for a list of names)
       4             : 
       5             :    See http://www.ves-code.org for more information.
       6             : 
       7             :    This file is part of VES code module.
       8             : 
       9             :    The VES code module 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             :    The VES code module 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 the VES code module.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : 
      23             : #include "TargetDistribution.h"
      24             : #include "GridIntegrationWeights.h"
      25             : #include "VesTools.h"
      26             : #include "GridLinearInterpolation.h"
      27             : 
      28             : #include "core/ActionRegister.h"
      29             : #include "tools/Grid.h"
      30             : #include "core/Value.h"
      31             : #include "tools/File.h"
      32             : 
      33             : 
      34             : namespace PLMD {
      35             : namespace ves {
      36             : 
      37             : 
      38             : //+PLUMEDOC VES_TARGETDIST TD_GRID
      39             : /*
      40             : Target distribution from an external grid file (static).
      41             : 
      42             : Using this keyword you can use a target distribution that is read from an
      43             : external grid file that is in the proper PLUMED file format. You do not to
      44             : give any information about the external grid file as all relevant information
      45             : should be automatically detected. It is assumed that the distribution read in
      46             : from the grid is a proper probability distribution, i.e. always non-negative
      47             : and can be normalized.
      48             : 
      49             : By default the target distribution from the external grid is always normalized
      50             : inside the code. You can disable this normalization by using DO_NOT_NORMALIZE
      51             : keyword. However, be warned that this will generally lead to the wrong
      52             : behavior if the distribution from the external grid is not properly
      53             : normalized to 1.
      54             : 
      55             : If the distribution from the external grid file has for some reason
      56             : negative values can you use the SHIFT keyword to shift the distribution
      57             : by a given value. Another option is to use
      58             : the SHIFT_TO_ZERO keyword to shift the minimum of the distribution to zero.
      59             : 
      60             : Note that the number of grid bins used in the external grid file do not have
      61             : to be the same as used in the bias or action where the target distribution is
      62             : employed as the code will employ a linear (or bilinear for two dimensions)
      63             : interpolation to calculate values. Currently only one or two dimensional grids
      64             : are supported.
      65             : 
      66             : It can happen that the intervals on which the target distribution is defined is
      67             : larger than the intervals covered by the external grid file. In this case the
      68             : default option is to consider the target distribution as continuous such that
      69             : values outside the boundary of the external grid file are the same as at
      70             : the boundary. This can be changed by using the ZERO_OUTSIDE keyword which
      71             : will make values outside to be taken as zero.
      72             : 
      73             : ## Examples
      74             : 
      75             : Generally you only need to provide the the filename of the external grid
      76             : file.
      77             : 
      78             : ```plumed
      79             : #SETTINGS INPUTFILES=regtest/ves/rt-td-grid/targetdist-input-grid-1d.data INPUTFILELINES=1-20
      80             : 
      81             : td: TD_GRID FILE=regtest/ves/rt-td-grid/targetdist-input-grid-1d.data
      82             : ```
      83             : 
      84             : The input grid is then specified using the usual format employed by PLUMED. You can see an example
      85             : by clicking on the name of the file in the input above.
      86             : 
      87             : */
      88             : //+ENDPLUMEDOC
      89             : 
      90             : 
      91             : class TD_Grid : public TargetDistribution {
      92             :   std::unique_ptr<GridBase> distGrid_;
      93             :   std::vector<double> minima_;
      94             :   std::vector<double> maxima_;
      95             :   std::vector<bool> periodic_;
      96             :   bool zero_outside_;
      97             :   double shift_;
      98             : public:
      99             :   static void registerKeywords( Keywords&);
     100             :   explicit TD_Grid(const ActionOptions& ao);
     101             :   double getValue(const std::vector<double>&) const override;
     102             : };
     103             : 
     104             : 
     105             : PLUMED_REGISTER_ACTION(TD_Grid,"TD_GRID")
     106             : 
     107             : 
     108          11 : void TD_Grid::registerKeywords(Keywords& keys) {
     109          11 :   TargetDistribution::registerKeywords(keys);
     110          11 :   keys.add("compulsory","FILE","The name of the external grid file to be used as a target distribution.");
     111          11 :   keys.add("optional","SHIFT","Shift the grid read in by some constant value. Due to normalization the final shift in the target distribution will generally not be the same as the value given here");
     112          11 :   keys.addFlag("ZERO_OUTSIDE",false,"By default the target distribution is continuous such that values outside the boundary of the external grid file are the same as at the boundary. This can be changed by using this flag which will make values outside to be taken as zero.");
     113          11 :   keys.addFlag("DO_NOT_NORMALIZE",false,"By default the target distribution from the external grid is always normalized inside the code. You can use this flag to disable this normalization. However, be warned that this will generally lead to the wrong behavior if the distribution from the external grid is not properly normalized to 1.");
     114          11 :   keys.use("WELLTEMPERED_FACTOR");
     115          11 :   keys.use("SHIFT_TO_ZERO");
     116          11 : }
     117             : 
     118           9 : TD_Grid::TD_Grid(const ActionOptions& ao):
     119             :   PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
     120           9 :   minima_(0),
     121           9 :   maxima_(0),
     122           9 :   zero_outside_(false),
     123          18 :   shift_(0.0) {
     124             :   std::string filename;
     125           9 :   parse("FILE",filename);
     126           9 :   parse("SHIFT",shift_);
     127           9 :   if(shift_!=0.0) {
     128           1 :     if(isTargetDistGridShiftedToZero()) {
     129           0 :       plumed_merror(getName() + ": using both SHIFT and SHIFT_TO_ZERO is not allowed.");
     130             :     }
     131             :   }
     132           9 :   parseFlag("ZERO_OUTSIDE",zero_outside_);
     133             : 
     134           9 :   bool do_not_normalize=false;
     135           9 :   parseFlag("DO_NOT_NORMALIZE",do_not_normalize);
     136           9 :   if(do_not_normalize && shift_!=0.0) {
     137           0 :     plumed_merror(getName() + ": using both SHIFT and DO_NOT_NORMALIZE is not allowed.");
     138             :   }
     139           9 :   if(!do_not_normalize) {
     140             :     setForcedNormalization();
     141             :   }
     142             : 
     143           9 :   checkRead();
     144             : 
     145             :   std::string gridlabel;
     146             :   std::vector<std::string> arglabels;
     147             :   std::vector<std::string> argmin;
     148             :   std::vector<std::string> argmax;
     149             :   std::vector<bool> argperiodic;
     150             :   std::vector<unsigned int> argnbins;
     151           9 :   bool has_deriv = false;
     152           9 :   unsigned int nargs = VesTools::getGridFileInfo(filename,gridlabel,arglabels,argmin,argmax,argperiodic,argnbins,has_deriv);
     153           9 :   if(nargs==0) {
     154           0 :     plumed_merror(getName() + ": problem in parsing information from grid file");
     155             :   }
     156             : 
     157           9 :   setDimension(arglabels.size());
     158           9 :   std::vector<std::unique_ptr<Value>> arguments(arglabels.size());
     159          22 :   for(unsigned int i=0; i < arglabels.size(); i++) {
     160          26 :     arguments[i]= Tools::make_unique<Value>(nullptr,arglabels[i],false);
     161          13 :     if(argperiodic[i]) {
     162           0 :       arguments[i]->setDomain(argmin[i],argmax[i]);
     163             :     } else {
     164          13 :       arguments[i]->setNotPeriodic();
     165             :     }
     166             :   }
     167             : 
     168           9 :   IFile gridfile;
     169           9 :   gridfile.open(filename);
     170           9 :   if(has_deriv) {
     171           0 :     distGrid_=GridBase::create(gridlabel,Tools::unique2raw(arguments),gridfile,false,true,true);
     172             :   } else {
     173          18 :     distGrid_=GridBase::create(gridlabel,Tools::unique2raw(arguments),gridfile,false,false,false);
     174             :   }
     175           9 :   gridfile.close();
     176             : 
     177             : 
     178           9 :   minima_.resize(getDimension());
     179           9 :   maxima_.resize(getDimension());
     180           9 :   periodic_.resize(getDimension());
     181          22 :   for (unsigned int i=0; i < getDimension(); i++) {
     182          13 :     Tools::convert(distGrid_->getMin()[i],minima_[i]);
     183          13 :     Tools::convert(distGrid_->getMax()[i],maxima_[i]);
     184          13 :     periodic_[i] = argperiodic[i];
     185          13 :     if(periodic_[i]) {
     186           0 :       maxima_[i]-=distGrid_->getDx()[i];
     187             :     }
     188             :   }
     189             : 
     190          27 : }
     191             : 
     192             : 
     193       49369 : double TD_Grid::getValue(const std::vector<double>& argument) const {
     194             :   double outside = 0.0;
     195       49369 :   std::vector<double> arg = argument;
     196      145566 :   for(unsigned int k=0; k<getDimension(); k++) {
     197       96533 :     if(zero_outside_ && (argument[k] < minima_[k] || argument[k] > maxima_[k])) {
     198             :       return outside;
     199       96197 :     } else if(argument[k] < minima_[k]) {
     200         168 :       arg[k] = minima_[k];
     201       96029 :     } else if(argument[k] > maxima_[k]) {
     202         168 :       arg[k] =maxima_[k];
     203             :     }
     204             :   }
     205       49033 :   return GridLinearInterpolation::getGridValueWithLinearInterpolation(distGrid_.get(),arg)+shift_;
     206             : }
     207             : 
     208             : 
     209             : }
     210             : }