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 : 27 : 28 : namespace PLMD { 29 : namespace ves { 30 : 31 : //+PLUMEDOC VES_BASISF BF_LEGENDRE 32 : /* 33 : Legendre polynomials basis functions. 34 : 35 : Use as basis functions [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) 36 : \f$P_{n}(x)\f$ defined on a bounded interval. 37 : You need to provide the interval \f$[a,b]\f$ 38 : on which the basis functions are to be used, and the order of the 39 : expansion \f$N\f$ (i.e. the highest order polynomial used). 40 : The total number of basis functions is \f$N+1\f$ as the constant \f$P_{0}(x)=1\f$ 41 : is also included. 42 : These basis functions should not be used for periodic CVs. 43 : 44 : Intrinsically the Legendre polynomials are defined on the interval \f$[-1,1]\f$. 45 : A variable \f$t\f$ in the interval \f$[a,b]\f$ is transformed to a variable \f$x\f$ 46 : in the intrinsic interval \f$[-1,1]\f$ by using the transform function 47 : \f[ 48 : x(t) = \frac{t-(a+b)/2} 49 : {(b-a)/2} 50 : \f] 51 : 52 : The Legendre polynomials are given by the recurrence relation 53 : \f{align}{ 54 : P_{0}(x) &= 1 \\ 55 : P_{1}(x) &= x \\ 56 : P_{n+1}(x) &= \frac{2n+1}{n+1} \, x \, P_{n}(x) - \frac{n}{n+1} \, P_{n-1}(x) 57 : \f} 58 : 59 : The first 6 polynomials are shown below 60 : \image html ves_basisf-legendre.png 61 : 62 : The Legendre polynomial are orthogonal over the interval \f$[-1,1]\f$ 63 : \f[ 64 : \int_{-1}^{1} dx \, P_{n}(x)\, P_{m}(x) = \frac{2}{2n+1} \delta_{n,m} 65 : \f] 66 : By using the SCALED keyword the polynomials are scaled by a factor of 67 : \f$ \sqrt{\frac{2n+1}{2}}\f$ such that they are orthonormal to 1. 68 : 69 : 70 : From the above equation it follows that integral of the basis functions 71 : over the uniform target distribution \f$p_{\mathrm{u}}(x)\f$ are given by 72 : \f[ 73 : \int_{-1}^{1} dx \, P_{n}(x) p_{\mathrm{u}}(x) = \delta_{n,0}, 74 : \f] 75 : and thus always zero except for the constant \f$P_{0}(x)=1\f$. 76 : 77 : 78 : For further mathematical properties of the Legendre polynomials see for example 79 : the [Wikipedia page](https://en.wikipedia.org/wiki/Legendre_polynomials). 80 : 81 : \par Examples 82 : 83 : Here we employ a Legendre expansion of order 20 over the interval -4.0 to 8.0. 84 : This results in a total number of 21 basis functions. 85 : The label used to identify the basis function action can then be 86 : referenced later on in the input file. 87 : \plumedfile 88 : bf_leg: BF_LEGENDRE MINIMUM=-4.0 MAXIMUM=8.0 ORDER=20 89 : \endplumedfile 90 : 91 : \par Examples 92 : 93 : */ 94 : //+ENDPLUMEDOC 95 : 96 : class BF_Legendre : public BasisFunctions { 97 : bool scaled_; 98 : void setupUniformIntegrals() override; 99 : public: 100 : static void registerKeywords(Keywords&); 101 : explicit BF_Legendre(const ActionOptions&); 102 : void getAllValues(const double, double&, bool&, std::vector<double>&, std::vector<double>&) const override; 103 : }; 104 : 105 : 106 : PLUMED_REGISTER_ACTION(BF_Legendre,"BF_LEGENDRE") 107 : 108 : 109 63 : void BF_Legendre::registerKeywords(Keywords& keys) { 110 63 : BasisFunctions::registerKeywords(keys); 111 126 : keys.addFlag("SCALED",false,"Scale the polynomials such that they are orthonormal to 1."); 112 63 : } 113 : 114 61 : BF_Legendre::BF_Legendre(const ActionOptions&ao): 115 : PLUMED_VES_BASISFUNCTIONS_INIT(ao), 116 61 : scaled_(false) 117 : { 118 183 : parseFlag("SCALED",scaled_); addKeywordToList("SCALED",scaled_); 119 61 : setNumberOfBasisFunctions(getOrder()+1); 120 122 : setIntrinsicInterval("-1.0","+1.0"); 121 : setNonPeriodic(); 122 : setIntervalBounded(); 123 61 : setType("Legendre"); 124 61 : setDescription("Legendre polynomials"); 125 61 : setLabelPrefix("L"); 126 61 : setupBF(); 127 61 : checkRead(); 128 61 : } 129 : 130 : 131 1625057 : void BF_Legendre::getAllValues(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const { 132 : // plumed_assert(values.size()==numberOfBasisFunctions()); 133 : // plumed_assert(derivs.size()==numberOfBasisFunctions()); 134 1625057 : inside_range=true; 135 1625057 : argT=translateArgument(arg, inside_range); 136 1625057 : std::vector<double> derivsT(derivs.size()); 137 : // 138 1625057 : values[0]=1.0; 139 1625057 : derivsT[0]=0.0; 140 1625057 : derivs[0]=0.0; 141 1625057 : values[1]=argT; 142 1625057 : derivsT[1]=1.0; 143 1625057 : derivs[1]=intervalDerivf(); 144 15470320 : for(unsigned int i=1; i < getOrder(); i++) { 145 13845263 : double io = static_cast<double>(i); 146 13845263 : values[i+1] = ((2.0*io+1.0)/(io+1.0))*argT*values[i] - (io/(io+1.0))*values[i-1]; 147 13845263 : derivsT[i+1] = ((2.0*io+1.0)/(io+1.0))*(values[i]+argT*derivsT[i])-(io/(io+1.0))*derivsT[i-1]; 148 13845263 : derivs[i+1] = intervalDerivf()*derivsT[i+1]; 149 : } 150 1625057 : if(scaled_) { 151 : // L0 is also scaled! 152 5088852 : for(unsigned int i=0; i<values.size(); i++) { 153 4632062 : double io = static_cast<double>(i); 154 4632062 : double sf = sqrt(io+0.5); 155 4632062 : values[i] *= sf; 156 4632062 : derivs[i] *= sf; 157 : } 158 : } 159 1756541 : if(!inside_range) {for(unsigned int i=0; i<derivs.size(); i++) {derivs[i]=0.0;}} 160 1625057 : } 161 : 162 : 163 59 : void BF_Legendre::setupUniformIntegrals() { 164 59 : setAllUniformIntegralsToZero(); 165 : double L0_int = 1.0; 166 59 : if(scaled_) {L0_int = sqrt(0.5);} 167 : setUniformIntegral(0,L0_int); 168 59 : } 169 : 170 : 171 : } 172 : }