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             : 
      29             : namespace PLMD {
      30             : namespace ves {
      31             : 
      32             : //+PLUMEDOC VES_BASISF BF_CUBIC_B_SPLINES
      33             : /*
      34             : Cubic B spline basis functions.
      35             : 
      36             : !!! attention ""
      37             : 
      38             :     __These basis functions do not form orthogonal bases. We recommend using wavelets ([BF_WAVELETS](BF_WAVELETS.md)) instead that do for orthogonal bases__.
      39             : 
      40             : A basis using cubic B spline functions according to the first paper cited below. See the second paper cited below for full details.
      41             : 
      42             : The mathematical expression of the individual splines is given by
      43             : $$
      44             : \begin{aligned}
      45             :   h\left(x\right) =
      46             :   \begin{cases}
      47             :     \left(2 - \lvert x \rvert\right)^3, & 1 \leq \lvert x \rvert \leq 2\\
      48             :     4 - 6\lvert x \rvert^2 + 3 \lvert x \rvert^3,\qquad & \lvert x \rvert \leq 1\\
      49             :     0, & \text{elsewhere}.
      50             :   \end{cases}
      51             : \end{aligned}
      52             : $$
      53             : 
      54             : The full basis consists of equidistant splines at positions $\mu_i$ which are optimized in their height:
      55             : $$
      56             :   f_i\left(x\right) = h\left(\frac{x-\mu_i}{\sigma}\right)
      57             : $$
      58             : 
      59             : Note that the distance between individual splines cannot be chosen freely but is equal to the width: $\mu_{i+1} = \mu_{i} + \sigma$.
      60             : 
      61             : 
      62             : The ORDER keyword of the basis set determines the number of equally sized sub-intervalls to be used.
      63             : On the borders of each of these sub-intervalls the mean $\mu_i$ of a spline function is placed.
      64             : 
      65             : The total number of basis functions is $\text{ORDER}+4$ as the constant $f_{0}(x)=1$, as well as the two splines with means just outside the interval are also included.
      66             : 
      67             : As an example two adjacent basis functions can be seen below.
      68             : The full basis consists of shifted splines in the full specified interval.
      69             : 
      70             : ![A graph illustrating the cubic B spline basisfunctions](figures/ves_basisf-splines.png)
      71             : 
      72             : When the splines are used for a periodic CV (with the PERIODIC keyword),
      73             : the sub-intervals are chosen in the same way, but only $\text{ORDER}+1$ functions
      74             : are required to fill it (the ones at the boundary coincide and the ones outside
      75             : can be omitted).
      76             : 
      77             : 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)
      78             : 
      79             : ## Examples
      80             : 
      81             : The bias is expanded with cubic B splines in the intervall from 0.0 to 10.0 specifying an order of 20.
      82             : This results in 24 basis functions.
      83             : 
      84             : ```plumed
      85             : bf: BF_CUBIC_B_SPLINES MINIMUM=0.0 MAXIMUM=10.0 ORDER=20
      86             : ```
      87             : 
      88             : */
      89             : //+ENDPLUMEDOC
      90             : 
      91             : class BF_CubicBsplines : public BasisFunctions {
      92             :   double spacing_;
      93             :   double inv_spacing_;
      94             :   double inv_normfactor_;
      95             :   double spline(const double, double&) const;
      96             : public:
      97             :   static void registerKeywords( Keywords&);
      98             :   explicit BF_CubicBsplines(const ActionOptions&);
      99             :   void getAllValues(const double, double&, bool&, std::vector<double>&, std::vector<double>&) const;
     100             : };
     101             : 
     102             : 
     103             : PLUMED_REGISTER_ACTION(BF_CubicBsplines,"BF_CUBIC_B_SPLINES")
     104             : 
     105             : // See DOI 10.1007/s10614-007-9092-4 for more information;
     106             : 
     107             : 
     108           5 : void BF_CubicBsplines::registerKeywords(Keywords& keys) {
     109           5 :   BasisFunctions::registerKeywords(keys);
     110           5 :   keys.add("optional","NORMALIZATION","the normalization factor that is used to normalize the basis functions by dividing the values. By default it is 2.");
     111           5 :   keys.addFlag("PERIODIC", false, "Use periodic version of basis set.");
     112           5 :   keys.remove("NUMERICAL_INTEGRALS");
     113           5 :   keys.addDOI("10.1007/s10614-007-9092-4");
     114           5 :   keys.addDOI("10.1021/acs.jctc.2c00197");
     115           5 : }
     116             : 
     117           3 : BF_CubicBsplines::BF_CubicBsplines(const ActionOptions&ao):
     118           3 :   PLUMED_VES_BASISFUNCTIONS_INIT(ao) {
     119           3 :   log.printf("  Cubic B spline basis functions, see and cite ");
     120           6 :   log << plumed.cite("Pampel and Valsson, J. Chem. Theory Comput. 18, 4127-4141 (2022) - DOI:10.1021/acs.jctc.2c00197");
     121             : 
     122           3 :   setIntrinsicInterval(intervalMin(),intervalMax());
     123           3 :   spacing_=(intervalMax()-intervalMin())/static_cast<double>(getOrder());
     124           3 :   inv_spacing_ = 1.0/spacing_;
     125             : 
     126           3 :   bool periodic = false;
     127           3 :   parseFlag("PERIODIC",periodic);
     128           3 :   if (periodic) {
     129           2 :     addKeywordToList("PERIODIC",periodic);
     130             :   }
     131             : 
     132             :   // 1 constant, getOrder() on interval, 1 (left) + 2 (right) at boundaries if not periodic
     133           3 :   unsigned int num_BFs = periodic ? getOrder()+1U : getOrder()+4U;
     134           3 :   setNumberOfBasisFunctions(num_BFs);
     135             : 
     136           3 :   double normfactor_=2.0;
     137           3 :   parse("NORMALIZATION",normfactor_);
     138           3 :   if(normfactor_!=2.0) {
     139           0 :     addKeywordToList("NORMALIZATION",normfactor_);
     140             :   }
     141           3 :   inv_normfactor_=1.0/normfactor_;
     142             : 
     143           3 :   periodic ? setPeriodic() : setNonPeriodic();
     144             :   setIntervalBounded();
     145           3 :   setType("splines_2nd-order");
     146           3 :   setDescription("Cubic B-splines (2nd order splines)");
     147           3 :   setLabelPrefix("S");
     148           3 :   setupBF();
     149           3 :   log.printf("   normalization factor: %f\n",normfactor_);
     150           3 :   checkRead();
     151           3 : }
     152             : 
     153             : 
     154        9236 : void BF_CubicBsplines::getAllValues(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
     155             :   // plumed_assert(values.size()==numberOfBasisFunctions());
     156             :   // plumed_assert(derivs.size()==numberOfBasisFunctions());
     157        9236 :   inside_range=true;
     158        9236 :   argT=checkIfArgumentInsideInterval(arg,inside_range);
     159             :   //
     160        9236 :   values[0]=1.0;
     161        9236 :   derivs[0]=0.0;
     162             :   //
     163      212043 :   for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
     164      202807 :     double argx = ((argT-intervalMin())*inv_spacing_) - (static_cast<double>(i) - 2.0);
     165      202807 :     if(arePeriodic()) { // periodic range of argx is [-intervalRange/spacing,+intervalRange/spacing]
     166       64140 :       argx /= intervalRange()*inv_spacing_;
     167       64140 :       argx = Tools::pbc(argx);
     168       64140 :       argx *= (intervalRange()*inv_spacing_);
     169             :     }
     170      202807 :     values[i] = spline(argx, derivs[i]);
     171      202807 :     derivs[i] *= inv_spacing_;
     172             :   }
     173        9236 :   if(!inside_range) {
     174        2000 :     for(unsigned int i=0; i<derivs.size(); i++) {
     175        1920 :       derivs[i]=0.0;
     176             :     }
     177             :   }
     178        9236 : }
     179             : 
     180             : 
     181      202807 : double BF_CubicBsplines::spline(const double arg, double& deriv) const {
     182             :   double value=0.0;
     183             :   double x=arg;
     184             :   // derivative of abs(x);
     185             :   double dx = 1.0;
     186      202807 :   if(x < 0) {
     187      101208 :     x=-x;
     188             :     dx = -1.0;
     189             :   }
     190             :   //
     191      202807 :   if(x > 2) {
     192             :     value=0.0;
     193      165556 :     deriv=0.0;
     194       37251 :   } else if(x >= 1) {
     195       18896 :     value = ((2.0-x)*(2.0-x)*(2.0-x));
     196       18896 :     deriv = dx*(-3.0*(2.0-x)*(2.0-x));
     197             :     // value=((2.0-x)*(2.0-x)*(2.0-x))/6.0;
     198             :     // deriv=-x*x*(2.0-x)*(2.0-x);
     199             :   } else {
     200       18355 :     value = 4.0-6.0*x*x+3.0*x*x*x;
     201       18355 :     deriv = dx*(-12.0*x+9.0*x*x);
     202             :     // value=x*x*x*0.5-x*x+2.0/3.0;
     203             :     // deriv=(3.0/2.0)*x*x-2.0*x;
     204             :   }
     205      202807 :   value *= inv_normfactor_;
     206      202807 :   deriv *= inv_normfactor_;
     207      202807 :   return value;
     208             : }
     209             : 
     210             : 
     211             : }
     212             : }