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 "FunctionSetup.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. Notice that the function that is computed is the
42 : exponentially scaled one, $e^{-|x|} I_0(x)$ (for example, it is equal to 0.4658 at $x=1$, where $I_0(1)=1.2661$),
43 : which is the quantity needed to normalize the von Mises kernels that are used in [KERNEL](KERNEL.md).
44 : Derivatives are not implemented, so this action cannot be used in a bias.
45 : You can evaluate the value of this Bessel function at $x=1$ by using the following input:
46 :
47 : ```plumed
48 : c: CONSTANT VALUE=1
49 : b: BESSEL ARG=c ORDER=0
50 : ```
51 :
52 : Notice that you can also use a vector in the input for this action as illustrated below:
53 :
54 : ```plumed
55 : c: CONSTANT VALUES=1,1.5,2
56 : b: BESSEL ARG=c ORDER=0
57 : ```
58 :
59 : The value output by BESSEL in this case is also a vector with three components. These components give the values
60 : of the Bessel function at 1.0, 1.5 and 2.0 respectively.
61 :
62 : ## The MASK keyword
63 :
64 : As this action only takes one argument in input the MASK keyword action works in the way that is described in the documentation for
65 : [LESS_THAN](LESS_THAN.md). This keyword is only used if the input argument is a vector. You pass to it a vector that contains ones
66 : and zeros and that has the same length as the vector that is passed in the input to ARG. If the $i$th element of the vector that was passed
67 : in the input to MASK is one then the corresponding element of the output vector is evaluated. If it is equal to zero then that element is not
68 : evaluated. Consequently in the example input shown below only the 1st and 4th element of the output vector are evaluated.
69 :
70 : ```plumed
71 : m: CONSTANT VALUES=1,0,0,1
72 : c: CONSTANT VALUES=1,1.5,2,5
73 : b: BESSEL ARG=c MASK=m ORDER=0
74 : ```
75 :
76 : */
77 : //+ENDPLUMEDOC
78 :
79 : class Bessel {
80 : static constexpr std::array<double,30> A{-4.41534164647933937950E-18,
81 : 3.33079451882223809783E-17,
82 : -2.43127984654795469359E-16,
83 : 1.71539128555513303061E-15,
84 : -1.16853328779934516808E-14,
85 : 7.67618549860493561688E-14,
86 : -4.85644678311192946090E-13,
87 : 2.95505266312963983461E-12,
88 : -1.72682629144155570723E-11,
89 : 9.67580903537323691224E-11,
90 : -5.18979560163526290666E-10,
91 : 2.65982372468238665035E-9,
92 : -1.30002500998624804212E-8,
93 : 6.04699502254191894932E-8,
94 : -2.67079385394061173391E-7,
95 : 1.11738753912010371815E-6,
96 : -4.41673835845875056359E-6,
97 : 1.64484480707288970893E-5,
98 : -5.75419501008210370398E-5,
99 : 1.88502885095841655729E-4,
100 : -5.76375574538582365885E-4,
101 : 1.63947561694133579842E-3,
102 : -4.32430999505057594430E-3,
103 : 1.05464603945949983183E-2,
104 : -2.37374148058994688156E-2,
105 : 4.93052842396707084878E-2,
106 : -9.49010970480476444210E-2,
107 : 1.71620901522208775349E-1,
108 : -3.04682672343198398683E-1,
109 : 6.76795274409476084995E-1};
110 : static constexpr std::array<double,25>B {-7.23318048787475395456E-18,
111 : -4.83050448594418207126E-18,
112 : 4.46562142029675999901E-17,
113 : 3.46122286769746109310E-17,
114 : -2.82762398051658348494E-16,
115 : -3.42548561967721913462E-16,
116 : 1.77256013305652638360E-15,
117 : 3.81168066935262242075E-15,
118 : -9.55484669882830764870E-15,
119 : -4.15056934728722208663E-14,
120 : 1.54008621752140982691E-14,
121 : 3.85277838274214270114E-13,
122 : 7.18012445138366623367E-13,
123 : -1.79417853150680611778E-12,
124 : -1.32158118404477131188E-11,
125 : -3.14991652796324136454E-11,
126 : 1.18891471078464383424E-11,
127 : 4.94060238822496958910E-10,
128 : 3.39623202570838634515E-9,
129 : 2.26666899049817806459E-8,
130 : 2.04891858946906374183E-7,
131 : 2.89137052083475648297E-6,
132 : 6.88975834691682398426E-5,
133 : 3.36911647825569408990E-3,
134 : 8.04490411014108831608E-1};
135 : public:
136 : unsigned order{0};
137 : static void registerKeywords( Keywords& keys );
138 : explicit Bessel();
139 : static void read( Bessel& func, ActionWithArguments* action, FunctionOptions& options );
140 : template <size_t N>
141 : static double chbevl(double x, const std::array<double,N>& array);
142 : static void calc( const Bessel& func, bool noderiv, View<const double> args, FunctionOutput& funcout );
143 : };
144 :
145 : typedef FunctionShortcut<Bessel> BesselShortcut;
146 : PLUMED_REGISTER_ACTION(BesselShortcut,"BESSEL")
147 : typedef FunctionOfScalar<Bessel> ScalarBessel;
148 : PLUMED_REGISTER_ACTION(ScalarBessel,"BESSEL_SCALAR")
149 : typedef FunctionOfVector<Bessel> VectorBessel;
150 : PLUMED_REGISTER_ACTION(VectorBessel,"BESSEL_VECTOR")
151 :
152 40 : void Bessel::registerKeywords(Keywords& keys) {
153 40 : keys.add("compulsory","ORDER","0","the order of Bessel function to use. Can only be zero at the moment.");
154 80 : keys.setValueDescription("scalar/vector","the value of the exponentially scaled modified bessel function of the first kind, exp(-|x|) I_0(x)");
155 40 : }
156 :
157 14 : Bessel::Bessel() {
158 14 : }
159 :
160 7 : void Bessel::read( Bessel& func, ActionWithArguments* action, FunctionOptions& options ) {
161 7 : if( action->getNumberOfArguments()!=1 ) {
162 0 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( action );
163 0 : if( !av || (av && action->getNumberOfArguments()-av->getNumberOfMasks()!=1) ) {
164 0 : action->error("should only be one argument to this action");
165 : }
166 : }
167 7 : if( action->getPntrToArgument(0)->isPeriodic() ) {
168 0 : action->error("cannot use this function on periodic functions");
169 : }
170 7 : action->parse("ORDER",func.order);
171 7 : action->log.printf(" computing %dth order bessel function \n", func.order );
172 7 : if( func.order!=0 ) {
173 0 : action->error("only zero order bessel function is implemented");
174 : }
175 7 : }
176 :
177 : template <size_t N>
178 21 : double Bessel::chbevl(const double x, const std::array<double,N>& array) {
179 : double b0, b1, b2;
180 :
181 21 : b0 = array[0];
182 : b1 = 0.0;
183 : b2 = b1 ;
184 :
185 585 : for(unsigned index = 1 ; index < N ; index++) {
186 : b2 = b1;
187 : b1 = b0;
188 564 : b0 = x * b1 - b2 + array[index];
189 : }
190 21 : return (0.5 * (b0 - b2));
191 : }
192 :
193 :
194 21 : void Bessel::calc( const Bessel& func,
195 : bool noderiv,
196 : const View<const double> args,
197 : FunctionOutput& funcout ) {
198 : // NDEBUG seems to be ignored by nvcc...
199 : plumed_dbg_assert( args.size()==1 );
200 21 : if( !noderiv ) {
201 0 : plumed_merror("derivatives have not been implemented for this function");
202 : }
203 : //I trust the parser into avoiding geting here with order!=0
204 21 : double x = fabs(args[0]);
205 21 : if (x <= 8.0) {
206 12 : double y = (x / 2.0) - 2.0;
207 12 : funcout.values[0] = chbevl(y, Bessel::A) ;
208 12 : return;
209 : }
210 9 : funcout.values[0] = chbevl(32.0 / x - 2.0, Bessel::B) / sqrt(x) ;
211 : ///TODO: higher order bessel functions
212 :
213 : }
214 :
215 : }
216 : }
217 :
218 :
|
