Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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 53 : class BF_Legendre : public BasisFunctions {
97 : bool scaled_;
98 : virtual void setupUniformIntegrals();
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;
103 : };
104 :
105 :
106 6505 : PLUMED_REGISTER_ACTION(BF_Legendre,"BF_LEGENDRE")
107 :
108 :
109 54 : void BF_Legendre::registerKeywords(Keywords& keys) {
110 54 : BasisFunctions::registerKeywords(keys);
111 162 : keys.addFlag("SCALED",false,"Scale the polynomials such that they are orthonormal to 1.");
112 54 : }
113 :
114 53 : BF_Legendre::BF_Legendre(const ActionOptions&ao):
115 : PLUMED_VES_BASISFUNCTIONS_INIT(ao),
116 53 : scaled_(false)
117 : {
118 159 : parseFlag("SCALED",scaled_); addKeywordToList("SCALED",scaled_);
119 53 : setNumberOfBasisFunctions(getOrder()+1);
120 159 : setIntrinsicInterval("-1.0","+1.0");
121 : setNonPeriodic();
122 : setIntervalBounded();
123 106 : setType("Legendre");
124 106 : setDescription("Legendre polynomials");
125 106 : setLabelPrefix("L");
126 53 : setupBF();
127 53 : checkRead();
128 53 : }
129 :
130 :
131 1464140 : 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 1464140 : inside_range=true;
135 1464140 : argT=translateArgument(arg, inside_range);
136 1464140 : std::vector<double> derivsT(derivs.size());
137 : //
138 1464140 : values[0]=1.0;
139 1464140 : derivsT[0]=0.0;
140 1464140 : derivs[0]=0.0;
141 1464140 : values[1]=argT;
142 1464140 : derivsT[1]=1.0;
143 1464140 : derivs[1]=intervalDerivf();
144 26811744 : for(unsigned int i=1; i < getOrder(); i++) {
145 12673802 : double io = static_cast<double>(i);
146 50695208 : values[i+1] = ((2.0*io+1.0)/(io+1.0))*argT*values[i] - (io/(io+1.0))*values[i-1];
147 50695208 : derivsT[i+1] = ((2.0*io+1.0)/(io+1.0))*(values[i]+argT*derivsT[i])-(io/(io+1.0))*derivsT[i-1];
148 25347604 : derivs[i+1] = intervalDerivf()*derivsT[i+1];
149 : }
150 1464140 : if(scaled_) {
151 : // L0 is also scaled!
152 14848896 : for(unsigned int i=0; i<values.size(); i++) {
153 4644704 : double io = static_cast<double>(i);
154 4644704 : double sf = sqrt(io+0.5);
155 4644704 : values[i] *= sf;
156 4644704 : derivs[i] *= sf;
157 : }
158 : }
159 1714981 : if(!inside_range) {for(unsigned int i=0; i<derivs.size(); i++) {derivs[i]=0.0;}}
160 1464140 : }
161 :
162 :
163 51 : void BF_Legendre::setupUniformIntegrals() {
164 : setAllUniformIntegralsToZero();
165 : double L0_int = 1.0;
166 51 : if(scaled_) {L0_int = sqrt(0.5);}
167 : setUniformIntegral(0,L0_int);
168 51 : }
169 :
170 :
171 : }
172 4839 : }
|
