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