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             : 
      26             : #include "core/ActionRegister.h"
      27             : #include "tools/Grid.h"
      28             : 
      29             : #include "lepton/Lepton.h"
      30             : 
      31             : 
      32             : namespace PLMD {
      33             : namespace ves {
      34             : 
      35             : //+PLUMEDOC VES_TARGETDIST TD_CUSTOM
      36             : /*
      37             : Target distribution given by an arbitrary mathematical expression (static or dynamic).
      38             : 
      39             : Use as a target distribution the distribution defined by
      40             : \f[
      41             : p(\mathbf{s}) =
      42             : \frac{f(\mathbf{s})}{\int d\mathbf{s} \, f(\mathbf{s})}
      43             : \f]
      44             : where \f$f(\mathbf{s})\f$ is some arbitrary mathematical function that
      45             : is parsed by the lepton library.
      46             : 
      47             : The function \f$f(\mathbf{s})\f$ is given by the FUNCTION keywords by
      48             : using _s1_,_s2_,..., as variables for the arguments
      49             : \f$\mathbf{s}=(s_1,s_2,\ldots,s_d)\f$.
      50             : If one variable is not given the target distribution will be
      51             : taken as uniform in that argument.
      52             : 
      53             : It is also possible to include the free energy surface \f$F(\mathbf{s})\f$
      54             : in the target distribution by using the _FE_ variable. In this case the
      55             : target distribution is dynamic and needs to be updated with current
      56             : best estimate of \f$F(\mathbf{s})\f$, similarly as for the
      57             : \ref TD_WELLTEMPERED "well-tempered target distribution".
      58             : Furthermore, the inverse temperature \f$\beta = (k_{\mathrm{B}}T)^{-1}\f$ and
      59             : the thermal energy \f$k_{\mathrm{B}}T\f$ can be included
      60             : by using the _beta_ and \f$k_B T\f$ variables.
      61             : 
      62             : The target distribution will be automatically normalized over the region on
      63             : which it is defined on. Therefore, the function given in
      64             : FUNCTION needs to be non-negative and it must be possible to normalize the function. The
      65             : code will perform checks to make sure that this is indeed the case.
      66             : 
      67             : 
      68             : \par Examples
      69             : 
      70             : Here we use as shifted [Maxwell-Boltzmann distribution](https://en.wikipedia.org/wiki/Maxwell%E2%80%93Boltzmann_distribution)
      71             : as a target distribution in one-dimension.
      72             : Note that it is not need to include the normalization factor as the distribution will be
      73             : automatically normalized.
      74             : \plumedfile
      75             : TD_CUSTOM ...
      76             :  FUNCTION=(s1+20)^2*exp(-(s1+20)^2/(2*10.0^2))
      77             :  LABEL=td
      78             : ... TD_CUSTOM
      79             : \endplumedfile
      80             : 
      81             : Here we have a two dimensional target distribution where we
      82             : use a [generalized normal distribution](https://en.wikipedia.org/wiki/Generalized_normal_distribution)
      83             : for argument \f$s_2\f$ while the distribution for \f$s_1\f$ is taken as
      84             : uniform as the variable _s1_ is not included in the function.
      85             : \plumedfile
      86             : TD_CUSTOM ...
      87             :  FUNCTION=exp(-(abs(s2-20.0)/5.0)^4.0)
      88             :  LABEL=td
      89             : ... TD_CUSTOM
      90             : \endplumedfile
      91             : 
      92             : By using the _FE_ variable the target distribution can depend on
      93             : the free energy surface \f$F(\mathbf{s})\f$. For example,
      94             : the following input is identical to using \ref TD_WELLTEMPERED with
      95             : a bias factor of 10.
      96             : \plumedfile
      97             : TD_CUSTOM ...
      98             :  FUNCTION=exp(-(beta/10.0)*FE)
      99             :  LABEL=td
     100             : ... TD_CUSTOM
     101             : \endplumedfile
     102             : Here the inverse temperature is automatically obtained by using the _beta_
     103             : variable. It is also possible to use the \f$k_B T\f$ variable. The following
     104             : syntax will give the exact same results as the syntax above
     105             : \plumedfile
     106             : TD_CUSTOM ...
     107             :  FUNCTION=exp(-(1.0/(kBT*10.0))*FE)
     108             :  LABEL=td
     109             : ... TD_CUSTOM
     110             : \endplumedfile
     111             : 
     112             : 
     113             : */
     114             : //+ENDPLUMEDOC
     115             : 
     116             : class TD_Custom : public TargetDistribution {
     117             : private:
     118             :   void setupAdditionalGrids(const std::vector<Value*>&, const std::vector<std::string>&, const std::vector<std::string>&, const std::vector<unsigned int>&) override;
     119             :   //
     120             :   lepton::CompiledExpression expression;
     121             :   //
     122             :   std::vector<double*> cv_var_lepton_refs_;
     123             :   double* kbt_var_lepton_ref_;
     124             :   double* beta_var_lepton_ref_;
     125             :   double* fes_var_lepton_ref_;
     126             :   //
     127             :   std::vector<unsigned int> cv_var_idx_;
     128             :   std::vector<std::string> cv_var_str_;
     129             :   //
     130             :   std::string cv_var_prefix_str_;
     131             :   std::string fes_var_str_;
     132             :   std::string kbt_var_str_;
     133             :   std::string beta_var_str_;
     134             :   //
     135             :   bool use_fes_;
     136             :   bool use_kbt_;
     137             :   bool use_beta_;
     138             : public:
     139             :   static void registerKeywords( Keywords&);
     140             :   explicit TD_Custom(const ActionOptions& ao);
     141             :   void updateGrid() override;
     142             :   double getValue(const std::vector<double>&) const override;
     143          64 :   ~TD_Custom() {};
     144             : };
     145             : 
     146       10435 : PLUMED_REGISTER_ACTION(TD_Custom,"TD_CUSTOM")
     147             : 
     148             : 
     149          17 : void TD_Custom::registerKeywords(Keywords& keys) {
     150          17 :   TargetDistribution::registerKeywords(keys);
     151          34 :   keys.add("compulsory","FUNCTION","The function you wish to use for the target distribution where you should use the variables _s1_,_s2_,... for the arguments. You can also use the current estimate of the FES by using the variable _FE_ and the temperature by using the \\f$k_B T\\f$ and _beta_ variables.");
     152          17 :   keys.use("WELLTEMPERED_FACTOR");
     153          17 :   keys.use("SHIFT_TO_ZERO");
     154          17 : }
     155             : 
     156             : 
     157          16 : TD_Custom::TD_Custom(const ActionOptions& ao):
     158             :   PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
     159             : //
     160          16 :   cv_var_lepton_refs_(0,nullptr),
     161          16 :   kbt_var_lepton_ref_(nullptr),
     162          16 :   beta_var_lepton_ref_(nullptr),
     163          16 :   fes_var_lepton_ref_(nullptr),
     164             : //
     165          16 :   cv_var_idx_(0),
     166          16 :   cv_var_str_(0),
     167             : //
     168          16 :   cv_var_prefix_str_("s"),
     169          16 :   fes_var_str_("FE"),
     170          16 :   kbt_var_str_("kBT"),
     171          16 :   beta_var_str_("beta"),
     172             : //
     173          16 :   use_fes_(false),
     174          16 :   use_kbt_(false),
     175          32 :   use_beta_(false)
     176             : {
     177             :   std::string func_str;
     178          16 :   parse("FUNCTION",func_str);
     179          16 :   checkRead();
     180             :   //
     181             :   try {
     182          16 :     lepton::ParsedExpression pe=lepton::Parser::parse(func_str).optimize(lepton::Constants());
     183          16 :     log<<"  function as parsed by lepton: "<<pe<<"\n";
     184          16 :     expression=pe.createCompiledExpression();
     185             :   }
     186           0 :   catch(PLMD::lepton::Exception& exc) {
     187           0 :     plumed_merror("There was some problem in parsing the function "+func_str+" given in FUNCTION with lepton");
     188           0 :   }
     189             : 
     190          39 :   for(auto &p: expression.getVariables()) {
     191          23 :     std::string curr_var = p;
     192             :     unsigned int cv_idx;
     193          39 :     if(curr_var.substr(0,cv_var_prefix_str_.size())==cv_var_prefix_str_ && Tools::convertNoexcept(curr_var.substr(cv_var_prefix_str_.size()),cv_idx) && cv_idx>0) {
     194          16 :       cv_var_idx_.push_back(cv_idx-1);
     195             :     }
     196           7 :     else if(curr_var==fes_var_str_) {
     197           3 :       use_fes_=true;
     198             :       setDynamic();
     199             :       setFesGridNeeded();
     200             :     }
     201           4 :     else if(curr_var==kbt_var_str_) {
     202           2 :       use_kbt_=true;
     203             :     }
     204           2 :     else if(curr_var==beta_var_str_) {
     205           2 :       use_beta_=true;
     206             :     }
     207             :     else {
     208           0 :       plumed_merror(getName()+": problem with parsing formula with lepton, cannot recognise the variable "+curr_var);
     209             :     }
     210             :   }
     211             :   //
     212          16 :   std::sort(cv_var_idx_.begin(),cv_var_idx_.end());
     213          16 :   cv_var_str_.resize(cv_var_idx_.size());
     214          16 :   cv_var_lepton_refs_.resize(cv_var_str_.size());
     215          32 :   for(unsigned int j=0; j<cv_var_idx_.size(); j++) {
     216          16 :     std::string str1; Tools::convert(cv_var_idx_[j]+1,str1);
     217          32 :     cv_var_str_[j] = cv_var_prefix_str_+str1;
     218             :     try {
     219          16 :       cv_var_lepton_refs_[j] = &expression.getVariableReference(cv_var_str_[j]);
     220           0 :     } catch(PLMD::lepton::Exception& exc) {}
     221             :   }
     222             : 
     223          16 :   if(use_kbt_) {
     224             :     try {
     225           2 :       kbt_var_lepton_ref_ = &expression.getVariableReference(kbt_var_str_);
     226           0 :     } catch(PLMD::lepton::Exception& exc) {}
     227             :   }
     228          16 :   if(use_beta_) {
     229             :     try {
     230           2 :       beta_var_lepton_ref_ = &expression.getVariableReference(beta_var_str_);
     231           0 :     } catch(PLMD::lepton::Exception& exc) {}
     232             :   }
     233          16 :   if(use_fes_) {
     234             :     try {
     235           3 :       fes_var_lepton_ref_ = &expression.getVariableReference(fes_var_str_);
     236           0 :     } catch(PLMD::lepton::Exception& exc) {}
     237             :   }
     238             : 
     239          16 : }
     240             : 
     241             : 
     242          16 : void TD_Custom::setupAdditionalGrids(const std::vector<Value*>& arguments, const std::vector<std::string>& min, const std::vector<std::string>& max, const std::vector<unsigned int>& nbins) {
     243          16 :   if(cv_var_idx_.size()>0 && cv_var_idx_[cv_var_idx_.size()-1]>getDimension()) {
     244           0 :     plumed_merror(getName()+": mismatch between CVs given in FUNC and the dimension of the target distribution");
     245             :   }
     246          16 : }
     247             : 
     248             : 
     249           0 : double TD_Custom::getValue(const std::vector<double>& argument) const {
     250           0 :   plumed_merror("getValue not implemented for TD_Custom");
     251             :   return 0.0;
     252             : }
     253             : 
     254             : 
     255          46 : void TD_Custom::updateGrid() {
     256          46 :   if(use_fes_) {
     257          33 :     plumed_massert(getFesGridPntr()!=NULL,"the FES grid has to be linked to the free energy in the target distribution");
     258             :   }
     259          46 :   if(use_kbt_) {
     260          22 :     if(kbt_var_lepton_ref_) {*kbt_var_lepton_ref_= 1.0/getBeta();}
     261             :   }
     262          46 :   if(use_beta_) {
     263          22 :     if(beta_var_lepton_ref_) {*beta_var_lepton_ref_= getBeta();}
     264             :   }
     265             :   //
     266          92 :   std::vector<double> integration_weights = GridIntegrationWeights::getIntegrationWeights(getTargetDistGridPntr());
     267             :   double norm = 0.0;
     268             :   //
     269       40909 :   for(Grid::index_t l=0; l<targetDistGrid().getSize(); l++) {
     270       40863 :     std::vector<double> point = targetDistGrid().getPoint(l);
     271       99928 :     for(unsigned int k=0; k<cv_var_str_.size() ; k++) {
     272       59065 :       if(cv_var_lepton_refs_[k]) {*cv_var_lepton_refs_[k] = point[cv_var_idx_[k]];}
     273             :     }
     274       40863 :     if(use_fes_) {
     275        3300 :       if(fes_var_lepton_ref_) {*fes_var_lepton_ref_ = getFesGridPntr()->getValue(l);}
     276             :     }
     277       40863 :     double value = expression.evaluate();
     278             : 
     279       40863 :     if(value<0.0 && !isTargetDistGridShiftedToZero()) {plumed_merror(getName()+": The target distribution function gives negative values. You should change the definition of the function used for the target distribution to avoid this. You can also use the SHIFT_TO_ZERO keyword to avoid this problem.");}
     280       40863 :     targetDistGrid().setValue(l,value);
     281       40863 :     norm += integration_weights[l]*value;
     282       40863 :     logTargetDistGrid().setValue(l,-std::log(value));
     283             :   }
     284          46 :   if(norm>0.0) {
     285          45 :     targetDistGrid().scaleAllValuesAndDerivatives(1.0/norm);
     286             :   }
     287           1 :   else if(!isTargetDistGridShiftedToZero()) {
     288           0 :     plumed_merror(getName()+": The target distribution function cannot be normalized proberly. You should change the definition of the function used for the target distribution to avoid this. You can also use the SHIFT_TO_ZERO keyword to avoid this problem.");
     289             :   }
     290          46 :   logTargetDistGrid().setMinToZero();
     291          46 : }
     292             : 
     293             : 
     294             : }
     295             : }