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