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             : \par Examples
      74             : 
      75             : Generally you only need to provide the the filename of the external grid
      76             : file.
      77             : \plumedfile
      78             : td: TD_GRID FILE=input-grid.data
      79             : \endplumedfile
      80             : 
      81             : The input grid is then specified using the usual format employed by PLUMED an example of which
      82             : is shown below:
      83             : 
      84             : \auxfile{input-grid.data}
      85             : #! FIELDS d1 external.bias der_d1
      86             : #! SET min_d1 1.14
      87             : #! SET max_d1 1.32
      88             : #! SET nbins_d1 6
      89             : #! SET periodic_d1 false
      90             :    1.1400   0.0031   0.1101
      91             :    1.1700   0.0086   0.2842
      92             :    1.2000   0.0222   0.6648
      93             :    1.2300   0.0521   1.4068
      94             :    1.2600   0.1120   2.6873
      95             :    1.2900   0.2199   4.6183
      96             :    1.3200   0.3948   7.1055
      97             : \endauxfile
      98             : 
      99             : */
     100             : //+ENDPLUMEDOC
     101             : 
     102             : 
     103             : class TD_Grid : public TargetDistribution {
     104             :   std::unique_ptr<GridBase> distGrid_;
     105             :   std::vector<double> minima_;
     106             :   std::vector<double> maxima_;
     107             :   std::vector<bool> periodic_;
     108             :   bool zero_outside_;
     109             :   double shift_;
     110             : public:
     111             :   static void registerKeywords( Keywords&);
     112             :   explicit TD_Grid(const ActionOptions& ao);
     113             :   double getValue(const std::vector<double>&) const override;
     114             : };
     115             : 
     116             : 
     117       10428 : PLUMED_REGISTER_ACTION(TD_Grid,"TD_GRID")
     118             : 
     119             : 
     120          10 : void TD_Grid::registerKeywords(Keywords& keys) {
     121          10 :   TargetDistribution::registerKeywords(keys);
     122          20 :   keys.add("compulsory","FILE","The name of the external grid file to be used as a target distribution.");
     123          20 :   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");
     124          20 :   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.");
     125          20 :   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.");
     126          10 :   keys.use("WELLTEMPERED_FACTOR");
     127          10 :   keys.use("SHIFT_TO_ZERO");
     128          10 : }
     129             : 
     130           9 : TD_Grid::TD_Grid(const ActionOptions& ao):
     131             :   PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
     132           9 :   minima_(0),
     133           9 :   maxima_(0),
     134           9 :   zero_outside_(false),
     135          18 :   shift_(0.0)
     136             : {
     137             :   std::string filename;
     138           9 :   parse("FILE",filename);
     139           9 :   parse("SHIFT",shift_);
     140           9 :   if(shift_!=0.0) {
     141           1 :     if(isTargetDistGridShiftedToZero()) {plumed_merror(getName() + ": using both SHIFT and SHIFT_TO_ZERO is not allowed.");}
     142             :   }
     143           9 :   parseFlag("ZERO_OUTSIDE",zero_outside_);
     144             : 
     145           9 :   bool do_not_normalize=false;
     146           9 :   parseFlag("DO_NOT_NORMALIZE",do_not_normalize);
     147           9 :   if(do_not_normalize && shift_!=0.0) {plumed_merror(getName() + ": using both SHIFT and DO_NOT_NORMALIZE is not allowed.");}
     148           9 :   if(!do_not_normalize) {setForcedNormalization();}
     149             : 
     150           9 :   checkRead();
     151             : 
     152             :   std::string gridlabel;
     153             :   std::vector<std::string> arglabels;
     154             :   std::vector<std::string> argmin;
     155             :   std::vector<std::string> argmax;
     156             :   std::vector<bool> argperiodic;
     157             :   std::vector<unsigned int> argnbins;
     158           9 :   bool has_deriv = false;
     159           9 :   unsigned int nargs = VesTools::getGridFileInfo(filename,gridlabel,arglabels,argmin,argmax,argperiodic,argnbins,has_deriv);
     160           9 :   if(nargs==0) {
     161           0 :     plumed_merror(getName() + ": problem in parsing information from grid file");
     162             :   }
     163             : 
     164           9 :   setDimension(arglabels.size());
     165           9 :   std::vector<std::unique_ptr<Value>> arguments(arglabels.size());
     166          22 :   for(unsigned int i=0; i < arglabels.size(); i++) {
     167          26 :     arguments[i]= Tools::make_unique<Value>(nullptr,arglabels[i],false);
     168          13 :     if(argperiodic[i]) {
     169           0 :       arguments[i]->setDomain(argmin[i],argmax[i]);
     170             :     }
     171             :     else {
     172          13 :       arguments[i]->setNotPeriodic();
     173             :     }
     174             :   }
     175             : 
     176           9 :   IFile gridfile; gridfile.open(filename);
     177           9 :   if(has_deriv) {
     178           0 :     distGrid_=GridBase::create(gridlabel,Tools::unique2raw(arguments),gridfile,false,true,true);
     179             :   }
     180             :   else {
     181          27 :     distGrid_=GridBase::create(gridlabel,Tools::unique2raw(arguments),gridfile,false,false,false);
     182             :   }
     183           9 :   gridfile.close();
     184             : 
     185             : 
     186           9 :   minima_.resize(getDimension());
     187           9 :   maxima_.resize(getDimension());
     188           9 :   periodic_.resize(getDimension());
     189          22 :   for (unsigned int i=0; i < getDimension(); i++) {
     190          13 :     Tools::convert(distGrid_->getMin()[i],minima_[i]);
     191          13 :     Tools::convert(distGrid_->getMax()[i],maxima_[i]);
     192          13 :     periodic_[i] = argperiodic[i];
     193          13 :     if(periodic_[i]) {maxima_[i]-=distGrid_->getDx()[i];}
     194             :   }
     195             : 
     196          27 : }
     197             : 
     198             : 
     199       49369 : double TD_Grid::getValue(const std::vector<double>& argument) const {
     200             :   double outside = 0.0;
     201       49369 :   std::vector<double> arg = argument;
     202      145566 :   for(unsigned int k=0; k<getDimension(); k++) {
     203       96533 :     if(zero_outside_ && (argument[k] < minima_[k] || argument[k] > maxima_[k])) {
     204             :       return outside;
     205             :     }
     206       96197 :     else if(argument[k] < minima_[k]) {
     207         168 :       arg[k] = minima_[k];
     208             :     }
     209       96029 :     else if(argument[k] > maxima_[k]) {
     210         168 :       arg[k] =maxima_[k];
     211             :     }
     212             :   }
     213       49033 :   return GridLinearInterpolation::getGridValueWithLinearInterpolation(distGrid_.get(),arg)+shift_;
     214             : }
     215             : 
     216             : 
     217             : }
     218             : }