Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2017 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed 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 : plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "FunctionTemplateBase.h"
23 : #include "FunctionShortcut.h"
24 : #include "FunctionOfScalar.h"
25 : #include "FunctionOfVector.h"
26 : #include "core/ActionRegister.h"
27 : #include <array>
28 : #include <cmath>
29 :
30 : namespace PLMD {
31 : namespace function {
32 :
33 : //+PLUMEDOC FUNCTION BESSEL
34 : /*
35 : Calculate the value of a Bessel function.
36 :
37 : This action can be used to evaluate [Bessel functions](https://en.wikipedia.org/wiki/Bessel_function).
38 : This was action is needed because the modified Bessel function of the first kind of order 0 is needed in
39 : order to calculate the [von-Misses distribution](https://en.wikipedia.org/wiki/Von_Mises_distribution)
40 : that is used in the implementation of [KERNEL](KERNEL.md). You can thus only use this function to evaluate the
41 : modified Bessel function of the first kind of order 0. You can evaluate the value of this Bessel function
42 : at $x=1$ by using the following input:
43 :
44 : ```plumed
45 : c: CONSTANT VALUE=1
46 : b: BESSEL ARG=c ORDER=0
47 : ```
48 :
49 : Notice that you can also use a vector in the input for this action as illustrated below:
50 :
51 : ```plumed
52 : c: CONSTANT VALUES=1,1,5,2
53 : b: BESSEL ARG=c ORDER=0
54 : ```
55 :
56 : The value output by BESSEL in this case is also a vector with three components. These components give the values
57 : of the Bessel function at 1.0, 1.5 and 2.0 respectively.
58 :
59 : */
60 : //+ENDPLUMEDOC
61 :
62 : //+PLUMEDOC FUNCTION BESSEL_SCALAR
63 : /*
64 : Calculate the value of a Bessel function.
65 :
66 : \par Examples
67 :
68 : */
69 : //+ENDPLUMEDOC
70 :
71 : //+PLUMEDOC FUNCTION BESSEL_VECTOR
72 : /*
73 : Calculate the bessel function for all the elements in a vector
74 :
75 : \par Examples
76 :
77 : */
78 : //+ENDPLUMEDOC
79 :
80 :
81 48 : class Bessel : public FunctionTemplateBase {
82 : private:
83 : unsigned order;
84 : // Cheb coefficient for range [0,8]
85 : static std::vector<double> A;
86 : // Cheb coefficient for range (8,inf)
87 : static std::vector<double> B;
88 : double chbevl(double x,std::vector<double>& array) const; // sub copied from scipy in C
89 : void fill_coefficients();
90 : public:
91 0 : bool derivativesImplemented() override {
92 0 : return false;
93 : }
94 : void registerKeywords( Keywords& keys ) override;
95 : void read( ActionWithArguments* action ) override;
96 : void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
97 : };
98 :
99 : typedef FunctionShortcut<Bessel> BesselShortcut;
100 : PLUMED_REGISTER_ACTION(BesselShortcut,"BESSEL")
101 : typedef FunctionOfScalar<Bessel> ScalarBessel;
102 : PLUMED_REGISTER_ACTION(ScalarBessel,"BESSEL_SCALAR")
103 : typedef FunctionOfVector<Bessel> VectorBessel;
104 : PLUMED_REGISTER_ACTION(VectorBessel,"BESSEL_VECTOR")
105 :
106 34 : void Bessel::registerKeywords(Keywords& keys) {
107 34 : keys.add("compulsory","ORDER","0","the order of Bessel function to use. Can only be zero at the moment.");
108 68 : keys.setValueDescription("scalar/vector","the value of the bessel function");
109 34 : }
110 :
111 7 : void Bessel::read( ActionWithArguments* action ) {
112 7 : if( action->getNumberOfArguments()!=1 ) {
113 0 : action->error("should only be one argument to bessel actions");
114 : }
115 7 : if( action->getPntrToArgument(0)->isPeriodic() ) {
116 0 : action->error("cannot use this function on periodic functions");
117 : }
118 7 : action->parse("ORDER",order);
119 7 : action->log.printf(" computing %dth order bessel function \n", order );
120 7 : if( order!=0 ) {
121 0 : action->error("only zero order bessel function is implemented");
122 : }
123 7 : fill_coefficients();
124 7 : }
125 :
126 14 : double Bessel::chbevl(double x,std::vector<double>& array) const {
127 : double b0, b1, b2;
128 14 : int n = array.size();
129 :
130 14 : b0 = array[0];
131 : b1 = 0.0;
132 : b2 = b1 ;
133 :
134 375 : for(int index = 1 ; index < n ; index++) {
135 : b2 = b1;
136 : b1 = b0;
137 361 : b0 = x * b1 - b2 + array[index];
138 : }
139 14 : return (0.5 * (b0 - b2));
140 : }
141 :
142 :
143 14 : void Bessel::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
144 : plumed_dbg_assert( args.size()==1 );
145 14 : if( order==0 ) {
146 14 : double x = fabs(args[0]);
147 14 : if (x <= 8.0) {
148 5 : double y = (x / 2.0) - 2.0;
149 5 : vals[0] = chbevl(y, A) ;
150 5 : return;
151 : }
152 9 : vals[0] = chbevl(32.0 / x - 2.0, B) / sqrt(x) ;
153 : } else {
154 0 : plumed_error();
155 : }
156 : }
157 :
158 : std::vector<double> Bessel::A;
159 : std::vector<double> Bessel::B;
160 :
161 7 : void Bessel::fill_coefficients() {
162 7 : A.resize(30);
163 14 : A = {-4.41534164647933937950E-18,
164 : 3.33079451882223809783E-17,
165 : -2.43127984654795469359E-16,
166 : 1.71539128555513303061E-15,
167 : -1.16853328779934516808E-14,
168 : 7.67618549860493561688E-14,
169 : -4.85644678311192946090E-13,
170 : 2.95505266312963983461E-12,
171 : -1.72682629144155570723E-11,
172 : 9.67580903537323691224E-11,
173 : -5.18979560163526290666E-10,
174 : 2.65982372468238665035E-9,
175 : -1.30002500998624804212E-8,
176 : 6.04699502254191894932E-8,
177 : -2.67079385394061173391E-7,
178 : 1.11738753912010371815E-6,
179 : -4.41673835845875056359E-6,
180 : 1.64484480707288970893E-5,
181 : -5.75419501008210370398E-5,
182 : 1.88502885095841655729E-4,
183 : -5.76375574538582365885E-4,
184 : 1.63947561694133579842E-3,
185 : -4.32430999505057594430E-3,
186 : 1.05464603945949983183E-2,
187 : -2.37374148058994688156E-2,
188 : 4.93052842396707084878E-2,
189 : -9.49010970480476444210E-2,
190 : 1.71620901522208775349E-1,
191 : -3.04682672343198398683E-1,
192 : 6.76795274409476084995E-1
193 7 : };
194 7 : B.resize(25);
195 14 : B = {-7.23318048787475395456E-18,
196 : -4.83050448594418207126E-18,
197 : 4.46562142029675999901E-17,
198 : 3.46122286769746109310E-17,
199 : -2.82762398051658348494E-16,
200 : -3.42548561967721913462E-16,
201 : 1.77256013305652638360E-15,
202 : 3.81168066935262242075E-15,
203 : -9.55484669882830764870E-15,
204 : -4.15056934728722208663E-14,
205 : 1.54008621752140982691E-14,
206 : 3.85277838274214270114E-13,
207 : 7.18012445138366623367E-13,
208 : -1.79417853150680611778E-12,
209 : -1.32158118404477131188E-11,
210 : -3.14991652796324136454E-11,
211 : 1.18891471078464383424E-11,
212 : 4.94060238822496958910E-10,
213 : 3.39623202570838634515E-9,
214 : 2.26666899049817806459E-8,
215 : 2.04891858946906374183E-7,
216 : 2.89137052083475648297E-6,
217 : 6.88975834691682398426E-5,
218 : 3.36911647825569408990E-3,
219 : 8.04490411014108831608E-1
220 7 : };
221 7 : }
222 :
223 : }
224 : }
225 :
226 :
|
