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