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 "BasisFunctions.h"
      24             : 
      25             : #include "core/ActionRegister.h"
      26             : #include "core/PlumedMain.h"
      27             : 
      28             : namespace PLMD {
      29             : namespace ves {
      30             : 
      31             : //+PLUMEDOC VES_BASISF BF_GAUSSIANS
      32             : /*
      33             : Gaussian basis functions.
      34             : 
      35             : \attention
      36             : __These basis functions do not form orthogonal bases. We recommend using wavelets (\ref BF_WAVELETS) instead that do form orthogonal bases__.
      37             : 
      38             : Basis functions given by Gaussian distributions with shifted centers defined on a
      39             : bounded interval. See \cite ValssonPampel_Wavelets_2022 for full details.
      40             : 
      41             : You need to provide the interval \f$[a,b]\f$ on which the bias is to be
      42             : expanded.
      43             : The ORDER keyword of the basis set \f$N\f$ determines the number of equally sized
      44             : sub-intervalls to be used.
      45             : On the borders of each of these sub-intervalls the mean \f$\mu\f$ of a Gaussian
      46             : basis function is placed:
      47             : \f{align}{
      48             :   \mu_i = a + (i-1) \frac{b-a}{N}
      49             : \f}
      50             : 
      51             : The total number of basis functions is \f$N+4\f$ as the constant
      52             : \f$f_{0}(x)=1\f$, as well as two additional Gaussians at the Boundaries are also included.
      53             : 
      54             : The basis functions are given by
      55             : \f{align}{
      56             :   f_0(x) &= 1 \\
      57             :   f_i(x) &= \exp\left(-\frac{{\left(x-\mu_i\right)}^2}{2\sigma^2}\right)
      58             : \f}
      59             : 
      60             : When the Gaussians are used for a periodic CV (with the PERIODIC keyword),
      61             : the sub-intervals are chosen in the same way, but only \f$N+1\f$ functions
      62             : are required to fill it (the ones at the boundary coincide and the ones outside
      63             : can be omitted).
      64             : 
      65             : It is possible to specify the width \f$\sigma\f$ (i.e. the standard deviation)
      66             : of the Gaussians using the WIDTH keyword.
      67             : By default it is set to the sub-intervall length.
      68             : It was found that performance can be typically improved with a smaller value (around 75 % of the sub-interval length), although a too small overlap will prevent the basis set from working correctly at all.
      69             : 
      70             : The optimization procedure then adjusts the heigths of the individual Gaussians.
      71             : To avoid 'blind' optimization of the basis functions outside the currently sampled area, it is often beneficial to use the OPTIMIZATION_THRESHOLD keyword of the VES_LINEAR_EXPANSION (set it to a small value, e.g. 1e-6)
      72             : 
      73             : As an example two adjacent basis functions (with the mentioned width choice of 75% of the sub-interval length) can be seen below.
      74             : The full basis consists of shifted Gaussians in the full specified interval.
      75             : 
      76             : \image html ves_basisf-gaussians.png
      77             : 
      78             : 
      79             : \par Examples
      80             : 
      81             : The bias is expanded with Gaussian functions in the intervall from 0.0 to
      82             : 10.0 using order 20.
      83             : This results in 24 basis functions.
      84             : 
      85             : \plumedfile
      86             : bfG: BF_GAUSSIANS MINIMUM=0.0 MAXIMUM=10.0 ORDER=20
      87             : \endplumedfile
      88             : 
      89             : Because it was not specified, the width of the Gaussians is by default
      90             : set to the sub-intervall length, i.e.\ \f$\sigma=0.5\f$.
      91             : To e.g. enhance the overlap between neighbouring basis functions, it can be
      92             : specified explicitely:
      93             : 
      94             : \plumedfile
      95             : bfG: BF_GAUSSIANS MINIMUM=0.0 MAXIMUM=10.0 ORDER=20 WIDTH=0.7
      96             : \endplumedfile
      97             : 
      98             : */
      99             : //+ENDPLUMEDOC
     100             : 
     101             : class BF_Gaussians : public BasisFunctions {
     102             :   /// one over width of the Gaussians
     103             :   double inv_sigma_;
     104             :   /// positions of the centers
     105             :   std::vector<double> centers_;
     106             :   void setupLabels() override;
     107             : public:
     108             :   static void registerKeywords( Keywords&);
     109             :   explicit BF_Gaussians(const ActionOptions&);
     110             :   void getAllValues(const double, double&, bool&, std::vector<double>&, std::vector<double>&) const override;
     111             : };
     112             : 
     113             : 
     114             : PLUMED_REGISTER_ACTION(BF_Gaussians,"BF_GAUSSIANS")
     115             : 
     116             : 
     117           5 : void BF_Gaussians::registerKeywords(Keywords& keys) {
     118           5 :   BasisFunctions::registerKeywords(keys);
     119          10 :   keys.add("optional","WIDTH","The width (i.e. standart deviation) of the Gaussian functions. By default it is equal to the sub-intervall size.");
     120          10 :   keys.addFlag("PERIODIC", false, "Use periodic version of basis set.");
     121           5 :   keys.remove("NUMERICAL_INTEGRALS");
     122           5 : }
     123             : 
     124           3 : BF_Gaussians::BF_Gaussians(const ActionOptions&ao):
     125           3 :   PLUMED_VES_BASISFUNCTIONS_INIT(ao)
     126             : {
     127           3 :   log.printf("  Gaussian basis functions, see and cite ");
     128           6 :   log << plumed.cite("Pampel and Valsson, J. Chem. Theory Comput. 18, 4127-4141 (2022) - DOI:10.1021/acs.jctc.2c00197");
     129             : 
     130           3 :   setIntrinsicInterval(intervalMin(),intervalMax());
     131             : 
     132           3 :   double width = (intervalMax()-intervalMin()) / getOrder();
     133           3 :   parse("WIDTH",width);
     134           3 :   if(width <= 0.0) {plumed_merror("WIDTH should be larger than 0");}
     135           4 :   if(width != (intervalMax()-intervalMin())/getOrder()) {addKeywordToList("WIDTH",width);}
     136           3 :   inv_sigma_ = 1/(width);
     137             : 
     138           3 :   bool periodic = false;
     139           3 :   parseFlag("PERIODIC",periodic);
     140           4 :   if (periodic) {addKeywordToList("PERIODIC",periodic);}
     141             : 
     142             :   // 1 constant, getOrder() on interval, 1 (left) + 2 (right) at boundaries if not periodic
     143           3 :   unsigned int num_BFs = periodic ? getOrder()+1U : getOrder()+4U;
     144           3 :   setNumberOfBasisFunctions(num_BFs);
     145             : 
     146           3 :   centers_.push_back(0.0); // constant one
     147          69 :   for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
     148          66 :     centers_.push_back(intervalMin()+(static_cast<int>(i) - 1 - static_cast<int>(!periodic))*(intervalMax()-intervalMin())/getOrder());
     149             :   }
     150           3 :   periodic ? setPeriodic() : setNonPeriodic();
     151             :   setIntervalBounded();
     152           3 :   setType("gaussian_functions");
     153           3 :   setDescription("Gaussian functions with shifted centers that are being optimized in their height");
     154           3 :   setupBF();
     155           3 :   log.printf("   width: %f\n",width);
     156           3 :   checkRead();
     157           3 : }
     158             : 
     159             : 
     160        9236 : void BF_Gaussians::getAllValues(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
     161        9236 :   inside_range=true;
     162        9236 :   argT=checkIfArgumentInsideInterval(arg,inside_range);
     163        9236 :   values[0]=1.0;
     164        9236 :   derivs[0]=0.0;
     165      212043 :   for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
     166      202807 :     double dist = argT - centers_[i];
     167      202807 :     if(arePeriodic()) { // wrap around similar to MetaD
     168       64140 :       dist /= intervalRange();
     169       64140 :       dist = Tools::pbc(dist);
     170       64140 :       dist *= intervalRange();
     171             :     }
     172      202807 :     values[i] = exp(-0.5*pow(dist*inv_sigma_,2.0));
     173      202807 :     derivs[i] = -values[i] * (dist)*pow(inv_sigma_,2.0);
     174             :   }
     175       11156 :   if(!inside_range) {for (auto& d: derivs) {d=0.0;}}
     176        9236 : }
     177             : 
     178             : 
     179             : // label according to position of mean
     180           3 : void BF_Gaussians::setupLabels() {
     181           3 :   setLabel(0,"const");
     182          69 :   for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
     183          66 :     std::string is; Tools::convert(centers_[i],is);
     184         132 :     setLabel(i,"m="+is);
     185             :   }
     186           3 : }
     187             : 
     188             : }
     189             : }