Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 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 "SwitchingFunction.h"
23 : #include "Tools.h"
24 : #include "Keywords.h"
25 : #include "OpenMP.h"
26 : #include <cmath>
27 : #include <vector>
28 : #include <limits>
29 : #include <algorithm>
30 : #include <optional>
31 :
32 : namespace PLMD {
33 :
34 : //+PLUMEDOC INTERNAL switchingfunction
35 : /*
36 : Functions that measure whether values are less than a certain quantity.
37 :
38 : Switching functions \f$s(r)\f$ take a minimum of one input parameter \f$r_0\f$.
39 : For \f$r \le d_0 \quad s(r)=1.0\f$ while for \f$r > d_0\f$ the function decays smoothly to 0.
40 : The various switching functions available in PLUMED differ in terms of how this decay is performed.
41 :
42 : Where there is an accepted convention in the literature (e.g. \ref COORDINATION) on the form of the
43 : switching function we use the convention as the default. However, the flexibility to use different
44 : switching functions is always present generally through a single keyword. This keyword generally
45 : takes an input with the following form:
46 :
47 : \verbatim
48 : KEYWORD={TYPE <list of parameters>}
49 : \endverbatim
50 :
51 : The following table contains a list of the various switching functions that are available in PLUMED 2
52 : together with an example input.
53 :
54 : <table align=center frame=void width=95%% cellpadding=5%%>
55 : <tr>
56 : <td> TYPE </td> <td> FUNCTION </td> <td> EXAMPLE INPUT </td> <td> DEFAULT PARAMETERS </td>
57 : </tr> <tr> <td>RATIONAL </td> <td>
58 : \f$
59 : s(r)=\frac{ 1 - \left(\frac{ r - d_0 }{ r_0 }\right)^{n} }{ 1 - \left(\frac{ r - d_0 }{ r_0 }\right)^{m} }
60 : \f$
61 : </td> <td>
62 : {RATIONAL R_0=\f$r_0\f$ D_0=\f$d_0\f$ NN=\f$n\f$ MM=\f$m\f$}
63 : </td> <td> \f$d_0=0.0\f$, \f$n=6\f$, \f$m=2n\f$ </td>
64 : </tr> <tr>
65 : <td> EXP </td> <td>
66 : \f$
67 : s(r)=\exp\left(-\frac{ r - d_0 }{ r_0 }\right)
68 : \f$
69 : </td> <td>
70 : {EXP R_0=\f$r_0\f$ D_0=\f$d_0\f$}
71 : </td> <td> \f$d_0=0.0\f$ </td>
72 : </tr> <tr>
73 : <td> GAUSSIAN </td> <td>
74 : \f$
75 : s(r)=\exp\left(-\frac{ (r - d_0)^2 }{ 2r_0^2 }\right)
76 : \f$
77 : </td> <td>
78 : {GAUSSIAN R_0=\f$r_0\f$ D_0=\f$d_0\f$}
79 : </td> <td> \f$d_0=0.0\f$ </td>
80 : </tr> <tr>
81 : <td> SMAP </td> <td>
82 : \f$
83 : s(r) = \left[ 1 + ( 2^{a/b} -1 )\left( \frac{r-d_0}{r_0} \right)^a \right]^{-b/a}
84 : \f$
85 : </td> <td>
86 : {SMAP R_0=\f$r_0\f$ D_0=\f$d_0\f$ A=\f$a\f$ B=\f$b\f$}
87 : </td> <td> \f$d_0=0.0\f$ </td>
88 : </tr> <tr>
89 : <td> Q </td> <td>
90 : \f$
91 : s(r) = \frac{1}{1 + \exp(\beta(r_{ij} - \lambda r_{ij}^0))}
92 : \f$
93 : </td> <td>
94 : {Q REF=\f$r_{ij}^0\f$ BETA=\f$\beta\f$ LAMBDA=\f$\lambda\f$ }
95 : </td> <td> \f$\lambda=1.8\f$, \f$\beta=50 nm^-1\f$ (all-atom)<br/>\f$\lambda=1.5\f$, \f$\beta=50 nm^-1\f$ (coarse-grained) </td>
96 : </tr> <tr>
97 : <td> CUBIC </td> <td>
98 : \f$
99 : s(r) = (y-1)^2(1+2y) \qquad \textrm{where} \quad y = \frac{r - r_1}{r_0-r_1}
100 : \f$
101 : </td> <td>
102 : {CUBIC D_0=\f$r_1\f$ D_MAX=\f$r_0\f$}
103 : </td> <td> </td>
104 : </tr> <tr>
105 : <td> TANH </td> <td>
106 : \f$
107 : s(r) = 1 - \tanh\left( \frac{ r - d_0 }{ r_0 } \right)
108 : \f$
109 : </td> <td>
110 : {TANH R_0=\f$r_0\f$ D_0=\f$d_0\f$}
111 : </td> <td> </td>
112 : </tr> <tr>
113 : <td> COSINUS </td> <td>
114 : \f$s(r) =\left\{\begin{array}{ll}
115 : 1 & \mathrm{if } r \leq d_0 \\
116 : 0.5 \left( \cos ( \frac{ r - d_0 }{ r_0 } \pi ) + 1 \right) & \mathrm{if } d_0 < r\leq d_0 + r_0 \\
117 : 0 & \mathrm{if } r < d_0 + r_0
118 : \end{array}\right.
119 : \f$
120 : </td> <td>
121 : {COSINUS R_0=\f$r_0\f$ D_0=\f$d_0\f$}
122 : </td> <td> </td>
123 : </tr> <tr>
124 : <td> CUSTOM </td> <td>
125 : \f$
126 : s(r) = FUNC
127 : \f$
128 : </td> <td>
129 : {CUSTOM FUNC=1/(1+x^6) R_0=\f$r_0\f$ D_0=\f$d_0\f$}
130 : </td> <td> </td>
131 : </tr>
132 : </table>
133 :
134 : Notice that most commonly used rational functions are better optimized and might run faster.
135 :
136 : Notice that for backward compatibility we allow using `MATHEVAL` instead of `CUSTOM`.
137 : Also notice that if the a `CUSTOM` switching function only depends on even powers of `x` it can be
138 : made faster by using `x2` as a variable. For instance
139 : \verbatim
140 : {CUSTOM FUNC=1/(1+x2^3) R_0=0.3}
141 : \endverbatim
142 : is equivalent to
143 : \verbatim
144 : {CUSTOM FUNC=1/(1+x^6) R_0=0.3}
145 : \endverbatim
146 : but runs faster. The reason is that there is an expensive square root calculation that can be optimized out.
147 :
148 :
149 : \attention
150 : With the default implementation CUSTOM is slower than other functions
151 : (e.g., it is slower than an equivalent RATIONAL function by approximately a factor 2).
152 : Checkout page \ref Lepton to see how to improve its performance.
153 :
154 : For all the switching functions in the above table one can also specify a further (optional) parameter using the parameter
155 : keyword D_MAX to assert that for \f$r>d_{\textrm{max}}\f$ the switching function can be assumed equal to zero.
156 : In this case the function is brought smoothly to zero by stretching and shifting it.
157 : \verbatim
158 : KEYWORD={RATIONAL R_0=1 D_MAX=3}
159 : \endverbatim
160 : the resulting switching function will be
161 : \f$
162 : s(r) = \frac{s'(r)-s'(d_{max})}{s'(0)-s'(d_{max})}
163 : \f$
164 : where
165 : \f$
166 : s'(r)=\frac{1-r^6}{1-r^{12}}
167 : \f$
168 : Since PLUMED 2.2 this is the default. The old behavior (no stretching) can be obtained with the
169 : NOSTRETCH flag. The NOSTRETCH keyword is only provided for backward compatibility and might be
170 : removed in the future. Similarly, the STRETCH keyword is still allowed but has no effect.
171 :
172 : Notice that switching functions defined with the simplified syntax are never stretched
173 : for backward compatibility. This might change in the future.
174 :
175 : */
176 : //+ENDPLUMEDOC
177 :
178 : namespace switchContainers {
179 :
180 1702 : baseSwitch::baseSwitch(double D0,double DMAX, double R0, std::string_view name)
181 1702 : : d0(D0),
182 1702 : dmax(DMAX),
183 1702 : dmax_2([](const double d) {
184 1702 : if(d<std::sqrt(std::numeric_limits<double>::max())) {
185 290 : return d*d;
186 : } else {
187 : return std::numeric_limits<double>::max();
188 : }
189 : }(dmax)),
190 1702 : invr0(1.0/R0),
191 1702 : invr0_2(invr0*invr0),
192 1992 : mytype(name) {}
193 :
194 1702 : baseSwitch::~baseSwitch()=default;
195 :
196 162833480 : double baseSwitch::calculate(const double distance, double& dfunc) const {
197 : double res = 0.0;//RVO!
198 162833480 : dfunc = 0.0;
199 162833480 : if(distance <= dmax) {
200 : res = 1.0;
201 156015847 : const double rdist = (distance-d0)*invr0;
202 156015847 : if(rdist > 0.0) {
203 59653094 : res = function(rdist,dfunc);
204 : //the following comments came from the original
205 : // this is for the chain rule (derivative of rdist):
206 59653094 : dfunc *= invr0;
207 : // for any future switching functions, be aware that multiplying invr0 is only
208 : // correct for functions of rdist = (r-d0)/r0.
209 :
210 : // this is because calculate() sets dfunc to the derivative divided times the
211 : // distance.
212 : // (I think this is misleading and I would like to modify it - GB)
213 59653094 : dfunc /= distance;
214 : }
215 156015847 : res=res*stretch+shift;
216 156015847 : dfunc*=stretch;
217 : }
218 162833480 : return res;
219 : }
220 :
221 31818704 : double baseSwitch::calculateSqr(double distance2,double&dfunc) const {
222 31818704 : double res= calculate(std::sqrt(distance2),dfunc);//RVO!
223 31818704 : return res;
224 : }
225 8 : double baseSwitch::get_d0() const {
226 8 : return d0;
227 : }
228 1534 : double baseSwitch::get_r0() const {
229 1534 : return 1.0/invr0;
230 : }
231 536580542 : double baseSwitch::get_dmax() const {
232 536580542 : return dmax;
233 : }
234 49030642 : double baseSwitch::get_dmax2() const {
235 49030642 : return dmax_2;
236 : }
237 1502 : std::string baseSwitch::description() const {
238 1502 : std::ostringstream ostr;
239 1502 : ostr<<get_r0()
240 : <<". Using "
241 : << mytype
242 3004 : <<" switching function with parameters d0="<< d0
243 3004 : << specificDescription();
244 1502 : return ostr.str();
245 1502 : }
246 150 : std::string baseSwitch::specificDescription() const {
247 150 : return "";
248 : }
249 239 : void baseSwitch::setupStretch() {
250 239 : if(dmax!=std::numeric_limits<double>::max()) {
251 239 : stretch=1.0;
252 239 : shift=0.0;
253 : double dummy;
254 239 : double s0=calculate(0.0,dummy);
255 239 : double sd=calculate(dmax,dummy);
256 239 : stretch=1.0/(s0-sd);
257 239 : shift=-sd*stretch;
258 : }
259 239 : }
260 0 : void baseSwitch::removeStretch() {
261 0 : stretch=1.0;
262 0 : shift=0.0;
263 0 : }
264 : template<int N, std::enable_if_t< (N >0), bool> = true, std::enable_if_t< (N %2 == 0), bool> = true>
265 : class fixedRational :public baseSwitch {
266 263 : std::string specificDescription() const override {
267 263 : std::ostringstream ostr;
268 263 : ostr << " nn=" << N << " mm=" <<N*2;
269 263 : return ostr.str();
270 263 : }
271 : public:
272 284 : fixedRational(double D0,double DMAX, double R0)
273 284 : :baseSwitch(D0,DMAX,R0,"rational") {}
274 :
275 : template <int POW>
276 1382 : static inline double doRational(const double rdist, double&dfunc, double result=0.0) {
277 : const double rNdist=Tools::fastpow<POW-1>(rdist);
278 27485057 : result=1.0/(1.0+rNdist*rdist);
279 27485057 : dfunc = -POW*rNdist*result*result;
280 1382 : return result;
281 : }
282 :
283 16154945 : inline double function(double rdist,double&dfunc) const override {
284 : //preRes and preDfunc are passed already set
285 1382 : dfunc=0.0;
286 1382 : double result = doRational<N>(rdist,dfunc);
287 16154945 : return result;
288 : }
289 :
290 11475870 : double calculateSqr(double distance2,double&dfunc) const override {
291 : double result=0.0;
292 11475870 : dfunc=0.0;
293 11475870 : if(distance2 <= dmax_2) {
294 11330112 : const double rdist = distance2*invr0_2;
295 : result = doRational<N/2>(rdist,dfunc);
296 11330112 : dfunc*=2*invr0_2;
297 : // stretch:
298 11330112 : result=result*stretch+shift;
299 11330112 : dfunc*=stretch;
300 : }
301 11475870 : return result;
302 :
303 : }
304 : };
305 :
306 : //these enums are useful for clarifying the settings in the factory
307 : //and the code is autodocumented ;)
308 : enum class rationalPow:bool {standard, fast};
309 : enum class rationalForm:bool {standard, simplified};
310 :
311 : template<rationalPow isFast, rationalForm nis2m>
312 : class rational : public baseSwitch {
313 : protected:
314 : const int nn=6;
315 : const int mm=12;
316 : const double preRes;
317 : const double preDfunc;
318 : const double preSecDev;
319 : const int nnf;
320 : const int mmf;
321 : const double preDfuncF;
322 : const double preSecDevF;
323 : //I am using PLMD::epsilon to be certain to call the one defined in Tools.h
324 : static constexpr double moreThanOne=1.0+5.0e10*PLMD::epsilon;
325 : static constexpr double lessThanOne=1.0-5.0e10*PLMD::epsilon;
326 :
327 177 : std::string specificDescription() const override {
328 177 : std::ostringstream ostr;
329 177 : ostr << " nn=" << nn << " mm=" <<mm;
330 177 : return ostr.str();
331 177 : }
332 : public:
333 202 : rational(double D0,double DMAX, double R0, int N, int M)
334 : :baseSwitch(D0,DMAX,R0,"rational"),
335 202 : nn(N),
336 202 : mm([](int m,int n) {
337 202 : if (m==0) {
338 91 : return n*2;
339 : } else {
340 : return m;
341 : }
342 : }(M,N)),
343 202 : preRes(static_cast<double>(nn)/mm),
344 202 : preDfunc(0.5*nn*(nn-mm)/static_cast<double>(mm)),
345 : //wolfram <3:lim_(x->1) d^2/(dx^2) (1 - x^N)/(1 - x^M) = (N (M^2 - 3 M (-1 + N) + N (-3 + 2 N)))/(6 M)
346 202 : preSecDev ((nn * (mm * mm - 3.0* mm * (-1 + nn ) + nn *(-3 + 2* nn )))/(6.0* mm )),
347 202 : nnf(nn/2),
348 202 : mmf(mm/2),
349 202 : preDfuncF(0.5*nnf*(nnf-mmf)/static_cast<double>(mmf)),
350 202 : preSecDevF((nnf* (mmf*mmf - 3.0* mmf* (-1 + nnf) + nnf*(-3 + 2* nnf)))/(6.0* mmf)) {}
351 :
352 18240750 : static inline double doRational(const double rdist, double&dfunc,double secDev, const int N,
353 : const int M,double result=0.0) {
354 : //the result and dfunc are assigned in the drivers for doRational
355 : //if(rdist>(1.0-100.0*epsilon) && rdist<(1.0+100.0*epsilon)) {
356 : //result=preRes;
357 : //dfunc=preDfunc;
358 : //} else {
359 : if constexpr (nis2m==rationalForm::simplified) {
360 2114004 : const double rNdist=Tools::fastpow(rdist,N-1);
361 2114004 : result=1.0/(1.0+rNdist*rdist);
362 2114004 : dfunc = -N*rNdist*result*result;
363 : } else {
364 16126746 : if(!((rdist > lessThanOne) && (rdist < moreThanOne))) {
365 16126734 : const double rNdist=Tools::fastpow(rdist,N-1);
366 16126734 : const double rMdist=Tools::fastpow(rdist,M-1);
367 16126734 : const double num = 1.0-rNdist*rdist;
368 16126734 : const double iden = 1.0/(1.0-rMdist*rdist);
369 16126734 : result = num*iden;
370 16126734 : dfunc = ((M*result*rMdist)-(N*rNdist))*iden;
371 16126734 : } else {
372 : //here I imply that the correct initialized are being passed to doRational
373 12 : const double x =(rdist-1.0);
374 12 : result = result+ x * ( dfunc + 0.5 * x * secDev);
375 12 : dfunc = dfunc + x * secDev;
376 : }
377 : }
378 18240750 : return result;
379 : }
380 18240684 : inline double function(double rdist,double&dfunc) const override {
381 : //preRes and preDfunc are passed already set
382 18240684 : dfunc=preDfunc;
383 18240684 : double result = doRational(rdist,dfunc,preSecDev,nn,mm,preRes);
384 18240684 : return result;
385 : }
386 :
387 3408419 : double calculateSqr(double distance2,double&dfunc) const override {
388 : if constexpr (isFast==rationalPow::fast) {
389 : double result=0.0;
390 80 : dfunc=0.0;
391 80 : if(distance2 <= dmax_2) {
392 66 : const double rdist = distance2*invr0_2;
393 66 : dfunc=preDfuncF;
394 66 : result = doRational(rdist,dfunc,preSecDevF,nnf,mmf,preRes);
395 66 : dfunc*=2*invr0_2;
396 : // stretch:
397 66 : result=result*stretch+shift;
398 66 : dfunc*=stretch;
399 : }
400 80 : return result;
401 : } else {
402 3408339 : double res= calculate(std::sqrt(distance2),dfunc);//RVO!
403 3408339 : return res;
404 : }
405 : }
406 : };
407 :
408 :
409 : template<int EXP,std::enable_if_t< (EXP %2 == 0), bool> = true>
410 1087 : std::optional<std::unique_ptr<baseSwitch>> fixedRationalFactory(double D0,double DMAX, double R0, int N) {
411 : if constexpr (EXP == 0) {
412 0 : return std::nullopt;
413 : } else {
414 1087 : if (N==EXP) {
415 284 : return PLMD::Tools::make_unique<switchContainers::fixedRational<EXP>>(D0,DMAX,R0);
416 : } else {
417 803 : return fixedRationalFactory<EXP-2>(D0,DMAX,R0,N);
418 : }
419 : }
420 : }
421 :
422 : std::unique_ptr<baseSwitch>
423 486 : rationalFactory(double D0,double DMAX, double R0, int N, int M) {
424 486 : bool fast = N%2==0 && M%2==0 && D0==0.0;
425 : //if (M==0) M will automatically became 2*NN
426 : constexpr int highestPrecompiledPower=12;
427 : //precompiled rational
428 486 : if(((2*N)==M || M == 0) && fast && N<=highestPrecompiledPower) {
429 284 : auto tmp = fixedRationalFactory<highestPrecompiledPower>(D0,DMAX,R0,N);
430 284 : if(tmp) {
431 : return std::move(*tmp);
432 : }
433 : //else continue with the at runtime implementation
434 : }
435 : //template<bool isFast, bool n2m>
436 : //class rational : public baseSwitch
437 202 : if(2*N==M || M == 0) {
438 134 : if(fast) {
439 : //fast rational
440 : return PLMD::Tools::make_unique<switchContainers::rational<
441 0 : rationalPow::fast,rationalForm::simplified>>(D0,DMAX,R0,N,M);
442 : }
443 : return PLMD::Tools::make_unique<switchContainers::rational<
444 134 : rationalPow::standard,rationalForm::simplified>>(D0,DMAX,R0,N,M);
445 : }
446 68 : if(fast) {
447 : //fast rational
448 : return PLMD::Tools::make_unique<switchContainers::rational<
449 63 : rationalPow::fast,rationalForm::standard>>(D0,DMAX,R0,N,M);
450 : }
451 : return PLMD::Tools::make_unique<switchContainers::rational<
452 5 : rationalPow::standard,rationalForm::standard>>(D0,DMAX,R0,N,M);
453 : }
454 : //function =
455 :
456 : class exponentialSwitch: public baseSwitch {
457 : public:
458 77 : exponentialSwitch(double D0, double DMAX, double R0)
459 77 : :baseSwitch(D0,DMAX,R0,"exponential") {}
460 : protected:
461 2404268 : inline double function(const double rdist,double&dfunc) const override {
462 2404268 : double result = std::exp(-rdist);
463 2404268 : dfunc=-result;
464 2404268 : return result;
465 : }
466 : };
467 :
468 : class gaussianSwitch: public baseSwitch {
469 : public:
470 68 : gaussianSwitch(double D0, double DMAX, double R0)
471 68 : :baseSwitch(D0,DMAX,R0,"gaussian") {}
472 : protected:
473 279665 : inline double function(const double rdist,double&dfunc) const override {
474 279665 : double result = std::exp(-0.5*rdist*rdist);
475 279665 : dfunc=-rdist*result;
476 279665 : return result;
477 : }
478 : };
479 :
480 : class fastGaussianSwitch: public baseSwitch {
481 : public:
482 116 : fastGaussianSwitch(double /*D0*/, double DMAX, double /*R0*/)
483 116 : :baseSwitch(0.0,DMAX,1.0,"fastgaussian") {}
484 : protected:
485 14 : inline double function(const double rdist,double&dfunc) const override {
486 14 : double result = std::exp(-0.5*rdist*rdist);
487 14 : dfunc=-rdist*result;
488 14 : return result;
489 : }
490 38317832 : inline double calculateSqr(double distance2,double&dfunc) const override {
491 : double result = 0.0;
492 38317832 : dfunc = 0.0;
493 38317832 : if (distance2 <dmax_2) {
494 : result=1.0;
495 :
496 38317818 : if(distance2>0.0) {
497 38317754 : result = exp(-0.5*distance2);
498 38317754 : dfunc = -result;
499 : // stretch:
500 38317754 : result=result*stretch+shift;
501 38317754 : dfunc*=stretch;
502 : }
503 : }
504 38317832 : return result;
505 : }
506 : };
507 :
508 : class smapSwitch: public baseSwitch {
509 : const int a=0;
510 : const int b=0;
511 : const double c=0.0;
512 : const double d=0.0;
513 : protected:
514 15 : std::string specificDescription() const override {
515 15 : std::ostringstream ostr;
516 15 : ostr<<" a="<<a<<" b="<<b;
517 15 : return ostr.str();
518 15 : }
519 : public:
520 17 : smapSwitch(double D0, double DMAX, double R0, int A, int B)
521 17 : :baseSwitch(D0,DMAX,R0,"smap"),
522 17 : a(A),
523 17 : b(B),
524 17 : c(std::pow(2., static_cast<double>(a)/static_cast<double>(b) ) - 1.0),
525 17 : d(-static_cast<double>(b) / static_cast<double>(a)) {}
526 : protected:
527 21911351 : inline double function(const double rdist,double&dfunc) const override {
528 :
529 21911351 : const double sx=c*Tools::fastpow( rdist, a );
530 21911351 : double result=std::pow( 1.0 + sx, d );
531 21911351 : dfunc=-b*sx/rdist*result/(1.0+sx);
532 21911351 : return result;
533 : }
534 : };
535 :
536 : class cubicSwitch: public baseSwitch {
537 : protected:
538 15 : std::string specificDescription() const override {
539 15 : std::ostringstream ostr;
540 15 : ostr<<" dmax="<<dmax;
541 15 : return ostr.str();
542 15 : }
543 : public:
544 17 : cubicSwitch(double D0, double DMAX)
545 17 : :baseSwitch(D0,DMAX,DMAX-D0,"cubic") {
546 : //this operation should be already done!!
547 : // R0 = dmax - d0;
548 : // invr0 = 1/R0;
549 : // invr0_2 = invr0*invr0;
550 17 : }
551 17 : ~cubicSwitch()=default;
552 : protected:
553 127277 : inline double function(const double rdist,double&dfunc) const override {
554 127277 : const double tmp1 = rdist - 1.0;
555 127277 : const double tmp2 = 1.0+2.0*rdist;
556 : //double result = tmp1*tmp1*tmp2;
557 127277 : dfunc = 2*tmp1*tmp2 + 2*tmp1*tmp1;
558 127277 : return tmp1*tmp1*tmp2;
559 : }
560 : };
561 :
562 : class tanhSwitch: public baseSwitch {
563 : public:
564 6 : tanhSwitch(double D0, double DMAX, double R0)
565 6 : :baseSwitch(D0,DMAX,R0,"tanh") {}
566 : protected:
567 12743 : inline double function(const double rdist,double&dfunc) const override {
568 12743 : const double tmp1 = std::tanh(rdist);
569 : //was dfunc=-(1-tmp1*tmp1);
570 12743 : dfunc = tmp1 * tmp1 - 1.0;
571 : //return result;
572 12743 : return 1.0 - tmp1;
573 : }
574 : };
575 :
576 : class cosinusSwitch: public baseSwitch {
577 : public:
578 5 : cosinusSwitch(double D0, double DMAX, double R0)
579 5 : :baseSwitch(D0,DMAX,R0,"cosinus") {}
580 : protected:
581 522147 : inline double function(const double rdist,double&dfunc) const override {
582 : double result = 0.0;
583 522147 : dfunc=0.0;
584 522147 : if(rdist<=1.0) {
585 : // rdist = (r-r1)/(r2-r1) ; 0.0<=rdist<=1.0 if r1 <= r <=r2; (r2-r1)/(r2-r1)=1
586 227036 : double rdistPI = rdist * PLMD::pi;
587 227036 : result = 0.5 * (std::cos ( rdistPI ) + 1.0);
588 227036 : dfunc = -0.5 * PLMD::pi * std::sin ( rdistPI );
589 : }
590 522147 : return result;
591 : }
592 : };
593 :
594 : class nativeqSwitch: public baseSwitch {
595 : double beta = 50.0; // nm-1
596 : double lambda = 1.8; // unitless
597 : double ref=0.0;
598 : protected:
599 864 : std::string specificDescription() const override {
600 864 : std::ostringstream ostr;
601 864 : ostr<<" beta="<<beta<<" lambda="<<lambda<<" ref="<<ref;
602 864 : return ostr.str();
603 864 : }
604 0 : inline double function(const double rdist,double&dfunc) const override {
605 0 : return 0.0;
606 : }
607 : public:
608 : nativeqSwitch(double D0, double DMAX, double R0, double BETA, double LAMBDA,double REF)
609 866 : : baseSwitch(D0,DMAX,R0,"nativeq"),beta(BETA),lambda(LAMBDA),ref(REF) {}
610 292966 : double calculate(const double distance, double& dfunc) const override {
611 : double res = 0.0;//RVO!
612 292966 : dfunc = 0.0;
613 292966 : if(distance<=dmax) {
614 : res = 1.0;
615 292946 : if(distance > d0) {
616 292934 : const double rdist = beta*(distance - lambda * ref);
617 292934 : double exprdist=std::exp(rdist);
618 292934 : res=1.0/(1.0+exprdist);
619 : /*2.9
620 : //need to see if this (5op+assign)
621 : //double exprmdist=1.0 + exprdist;
622 : //dfunc = - (beta *exprdist)/(exprmdist*exprmdist);
623 : //or this (5op but 2 divisions) is faster
624 : dfunc = - beta /(exprdist+ 2.0 +1.0/exprdist);
625 : //this cames from - beta * exprdist/(exprdist*exprdist+ 2.0 *exprdist +1.0)
626 : //dfunc *= invr0;
627 : dfunc /= distance;
628 : */
629 : //2.10
630 292934 : dfunc = - beta /(exprdist+ 2.0 +1.0/exprdist) /distance;
631 :
632 292934 : dfunc*=stretch;
633 : }
634 292946 : res=res*stretch+shift;
635 : }
636 292966 : return res;
637 : }
638 : };
639 :
640 : class leptonSwitch: public baseSwitch {
641 : /// Lepton expression.
642 134 : class funcAndDeriv {
643 : lepton::CompiledExpression expression;
644 : lepton::CompiledExpression deriv;
645 : double* varRef=nullptr;
646 : double* varDevRef=nullptr;
647 : public:
648 44 : funcAndDeriv(const std::string &func) {
649 44 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
650 44 : expression=pe.createCompiledExpression();
651 46 : std::string arg="x";
652 :
653 : {
654 44 : auto vars=expression.getVariables();
655 44 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
656 44 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
657 :
658 44 : if(found_x2) {
659 : arg="x2";
660 : }
661 44 : if (vars.size()==0) {
662 : // this is necessary since in some cases lepton thinks a variable is not present even though it is present
663 : // e.g. func=0*x
664 0 : varRef=nullptr;
665 44 : } else if(vars.size()==1 && (found_x || found_x2)) {
666 42 : varRef=&expression.getVariableReference(arg);
667 : } else {
668 4 : plumed_error()
669 : <<"Please declare a function with only ONE argument that can only be x or x2. Your function is: "
670 4 : << func;
671 : }
672 : }
673 :
674 86 : lepton::ParsedExpression ped=lepton::Parser::parse(func).differentiate(arg).optimize(lepton::Constants());
675 42 : deriv=ped.createCompiledExpression();
676 : {
677 42 : auto vars=expression.getVariables();
678 42 : if (vars.size()==0) {
679 0 : varDevRef=nullptr;
680 : } else {
681 42 : varDevRef=&deriv.getVariableReference(arg);
682 : }
683 : }
684 :
685 46 : }
686 92 : funcAndDeriv (const funcAndDeriv& other):
687 92 : expression(other.expression),
688 92 : deriv(other.deriv) {
689 92 : std::string arg="x";
690 :
691 : {
692 92 : auto vars=expression.getVariables();
693 92 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
694 92 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
695 :
696 92 : if(found_x2) {
697 : arg="x2";
698 : }
699 92 : if (vars.size()==0) {
700 0 : varRef=nullptr;
701 92 : } else if(vars.size()==1 && (found_x || found_x2)) {
702 92 : varRef=&expression.getVariableReference(arg);
703 : }// UB: I assume that the function is already correct
704 : }
705 :
706 : {
707 92 : auto vars=expression.getVariables();
708 92 : if (vars.size()==0) {
709 0 : varDevRef=nullptr;
710 : } else {
711 92 : varDevRef=&deriv.getVariableReference(arg);
712 : }
713 : }
714 92 : }
715 :
716 : funcAndDeriv& operator= (const funcAndDeriv& other) {
717 : if(this != &other) {
718 : expression = other.expression;
719 : deriv = other.deriv;
720 : std::string arg="x";
721 :
722 : {
723 : auto vars=expression.getVariables();
724 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
725 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
726 :
727 : if(found_x2) {
728 : arg="x2";
729 : }
730 : if (vars.size()==0) {
731 : varRef=nullptr;
732 : } else if(vars.size()==1 && (found_x || found_x2)) {
733 : varRef=&expression.getVariableReference(arg);
734 : }// UB: I assume that the function is already correct
735 : }
736 :
737 : {
738 : auto vars=expression.getVariables();
739 : if (vars.size()==0) {
740 : varDevRef=nullptr;
741 : } else {
742 : varDevRef=&deriv.getVariableReference(arg);
743 : }
744 : }
745 : }
746 : return *this;
747 : }
748 :
749 6515577 : std::pair<double,double> operator()(double const x) const {
750 : //FAQ: why this works? this thing is const and you are modifying things!
751 : //Actually I am modifying something that is pointed at, not my pointers,
752 : //so I am not mutating the state of this!
753 6515577 : if(varRef) {
754 6515577 : *varRef=x;
755 : }
756 6515577 : if(varDevRef) {
757 6515577 : *varDevRef=x;
758 : }
759 : return std::make_pair(
760 6515577 : expression.evaluate(),
761 6515577 : deriv.evaluate());
762 : }
763 :
764 : auto& getVariables() const {
765 42 : return expression.getVariables();
766 : }
767 : auto& getVariables_derivative() const {
768 : return deriv.getVariables();
769 : }
770 : };
771 : /// Function for lepton
772 : std::string lepton_func;
773 : /// \warning Since lepton::CompiledExpression is mutable, a vector is necessary for multithreading!
774 : std::vector <funcAndDeriv> expressions{};
775 : /// Set to true if lepton only uses x2
776 : bool leptonx2=false;
777 : protected:
778 18 : std::string specificDescription() const override {
779 18 : std::ostringstream ostr;
780 18 : ostr<<" func=" << lepton_func;
781 18 : return ostr.str();
782 18 : }
783 0 : inline double function(const double,double&) const override {
784 0 : return 0.0;
785 : }
786 : public:
787 44 : leptonSwitch(double D0, double DMAX, double R0, const std::string & func)
788 44 : :baseSwitch(D0,DMAX,R0,"lepton"),
789 44 : lepton_func(func),
790 86 : expressions (OpenMP::getNumThreads(), lepton_func) {
791 : //this is a bit odd, but it works
792 : auto vars=expressions[0].getVariables();
793 42 : leptonx2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
794 44 : }
795 :
796 5878300 : double calculate(const double distance,double&dfunc) const override {
797 5878300 : double res = 0.0;//RVO!
798 5878300 : dfunc = 0.0;
799 5878300 : if(leptonx2) {
800 2 : res= calculateSqr(distance*distance,dfunc);
801 : } else {
802 5878298 : if(distance<=dmax) {
803 5573465 : res = 1.0;
804 5573465 : const double rdist = (distance-d0)*invr0;
805 5573465 : if(rdist > 0.0) {
806 5267475 : const unsigned t=OpenMP::getThreadNum();
807 5267475 : plumed_assert(t<expressions.size());
808 5267475 : std::tie(res,dfunc) = expressions[t](rdist);
809 5267475 : dfunc *= invr0;
810 5267475 : dfunc /= distance;
811 : }
812 5573465 : res=res*stretch+shift;
813 5573465 : dfunc*=stretch;
814 : }
815 : }
816 5878300 : return res;
817 : }
818 :
819 7126130 : double calculateSqr(const double distance2,double&dfunc) const override {
820 7126130 : double result =0.0;
821 7126130 : dfunc=0.0;
822 7126130 : if(leptonx2) {
823 1248110 : if(distance2<=dmax_2) {
824 1248102 : const unsigned t=OpenMP::getThreadNum();
825 1248102 : const double rdist_2 = distance2*invr0_2;
826 1248102 : plumed_assert(t<expressions.size());
827 1248102 : std::tie(result,dfunc) = expressions[t](rdist_2);
828 : // chain rule:
829 1248102 : dfunc*=2*invr0_2;
830 : // stretch:
831 1248102 : result=result*stretch+shift;
832 1248102 : dfunc*=stretch;
833 : }
834 : } else {
835 5878020 : result = calculate(std::sqrt(distance2),dfunc);
836 : }
837 7126130 : return result;
838 : }
839 : };
840 : } // namespace switchContainers
841 :
842 0 : void SwitchingFunction::registerKeywords( Keywords& keys ) {
843 0 : keys.add("compulsory","R_0","the value of R_0 in the switching function");
844 0 : keys.add("compulsory","D_0","0.0","the value of D_0 in the switching function");
845 0 : keys.add("optional","D_MAX","the value at which the switching function can be assumed equal to zero");
846 0 : keys.add("compulsory","NN","6","the value of n in the switching function (only needed for TYPE=RATIONAL)");
847 0 : keys.add("compulsory","MM","0","the value of m in the switching function (only needed for TYPE=RATIONAL); 0 implies 2*NN");
848 0 : keys.add("compulsory","A","the value of a in the switching function (only needed for TYPE=SMAP)");
849 0 : keys.add("compulsory","B","the value of b in the switching function (only needed for TYPE=SMAP)");
850 0 : }
851 :
852 1629 : void SwitchingFunction::set(const std::string & definition,std::string& errormsg) {
853 1629 : std::vector<std::string> data=Tools::getWords(definition);
854 : #define CHECKandPARSE(datastring,keyword,variable,errormsg) \
855 : if(Tools::findKeyword(datastring,keyword) && !Tools::parse(datastring,keyword,variable))\
856 : errormsg="could not parse " keyword; //adiacent strings are automagically concatenated
857 : #define REQUIREDPARSE(datastring,keyword,variable,errormsg) \
858 : if(!Tools::parse(datastring,keyword,variable))\
859 : errormsg=keyword " is required for " + name ; //adiacent strings are automagically concatenated
860 :
861 1629 : if( data.size()<1 ) {
862 : errormsg="missing all input for switching function";
863 : return;
864 : }
865 1629 : std::string name=data[0];
866 : data.erase(data.begin());
867 1629 : double r0=0.0;
868 1629 : double d0=0.0;
869 1629 : double dmax=std::numeric_limits<double>::max();
870 1629 : init=true;
871 2097 : CHECKandPARSE(data,"D_0",d0,errormsg);
872 2061 : CHECKandPARSE(data,"D_MAX",dmax,errormsg);
873 :
874 1629 : bool dostretch=false;
875 1629 : Tools::parseFlag(data,"STRETCH",dostretch); // this is ignored now
876 1629 : dostretch=true;
877 1629 : bool dontstretch=false;
878 1629 : Tools::parseFlag(data,"NOSTRETCH",dontstretch); // this is ignored now
879 1629 : if(dontstretch) {
880 199 : dostretch=false;
881 : }
882 1629 : if(name=="CUBIC") {
883 : //cubic is the only switch type that only uses d0 and dmax
884 17 : function = PLMD::Tools::make_unique<switchContainers::cubicSwitch>(d0,dmax);
885 : } else {
886 3224 : REQUIREDPARSE(data,"R_0",r0,errormsg);
887 1612 : if(name=="RATIONAL") {
888 412 : int nn=6;
889 412 : int mm=0;
890 680 : CHECKandPARSE(data,"NN",nn,errormsg);
891 670 : CHECKandPARSE(data,"MM",mm,errormsg);
892 824 : function = switchContainers::rationalFactory(d0,dmax,r0,nn,mm);
893 1200 : } else if(name=="SMAP") {
894 17 : int a=0;
895 17 : int b=0;
896 : //in the original a and b are "default=0",
897 : //but you divide by a and b during the initialization!
898 : //better an error message than an UB, so no default
899 34 : REQUIREDPARSE(data,"A",a,errormsg);
900 34 : REQUIREDPARSE(data,"B",b,errormsg);
901 17 : function = PLMD::Tools::make_unique<switchContainers::smapSwitch>(d0,dmax,r0,a,b);
902 1183 : } else if(name=="Q") {
903 866 : double beta = 50.0; // nm-1
904 866 : double lambda = 1.8; // unitless
905 : double ref;
906 2598 : CHECKandPARSE(data,"BETA",beta,errormsg);
907 2598 : CHECKandPARSE(data,"LAMBDA",lambda,errormsg);
908 1732 : REQUIREDPARSE(data,"REF",ref,errormsg);
909 : //the original error message was not standard
910 : // if(!Tools::parse(data,"REF",ref))
911 : // errormsg="REF (reference distaance) is required for native Q";
912 866 : function = PLMD::Tools::make_unique<switchContainers::nativeqSwitch>(d0,dmax,r0,beta,lambda,ref);
913 317 : } else if(name=="EXP") {
914 77 : function = PLMD::Tools::make_unique<switchContainers::exponentialSwitch>(d0,dmax,r0);
915 240 : } else if(name=="GAUSSIAN") {
916 184 : if ( r0==1.0 && d0==0.0 ) {
917 116 : function = PLMD::Tools::make_unique<switchContainers::fastGaussianSwitch>(d0,dmax,r0);
918 : } else {
919 68 : function = PLMD::Tools::make_unique<switchContainers::gaussianSwitch>(d0,dmax,r0);
920 : }
921 56 : } else if(name=="TANH") {
922 6 : function = PLMD::Tools::make_unique<switchContainers::tanhSwitch>(d0,dmax,r0);
923 50 : } else if(name=="COSINUS") {
924 5 : function = PLMD::Tools::make_unique<switchContainers::cosinusSwitch>(d0,dmax,r0);
925 87 : } else if((name=="MATHEVAL" || name=="CUSTOM")) {
926 : std::string func;
927 88 : Tools::parse(data,"FUNC",func);
928 42 : function = PLMD::Tools::make_unique<switchContainers::leptonSwitch>(d0,dmax,r0,func);
929 : } else {
930 2 : errormsg="cannot understand switching function type '"+name+"'";
931 : }
932 : }
933 : #undef CHECKandPARSE
934 : #undef REQUIREDPARSE
935 :
936 1627 : if( !data.empty() ) {
937 : errormsg="found the following rogue keywords in switching function input : ";
938 0 : for(unsigned i=0; i<data.size(); ++i) {
939 2 : errormsg = errormsg + data[i] + " ";
940 : }
941 : }
942 :
943 1627 : if(dostretch && dmax!=std::numeric_limits<double>::max()) {
944 165 : function->setupStretch();
945 : }
946 1629 : }
947 :
948 1502 : std::string SwitchingFunction::description() const {
949 : // if this error is necessary, something went wrong in the constructor
950 : // plumed_merror("Unknown switching function type");
951 1502 : return function->description();
952 : }
953 :
954 92146953 : double SwitchingFunction::calculateSqr(double distance2,double&dfunc)const {
955 92146953 : return function -> calculateSqr(distance2, dfunc);
956 : }
957 :
958 127899205 : double SwitchingFunction::calculate(double distance,double&dfunc)const {
959 127899205 : plumed_massert(init,"you are trying to use an unset SwitchingFunction");
960 127899205 : double result=function->calculate(distance,dfunc);
961 127899205 : return result;
962 : }
963 :
964 74 : void SwitchingFunction::set(const int nn,int mm, const double r0, const double d0) {
965 74 : init=true;
966 74 : if(mm == 0) {
967 70 : mm = 2*nn;
968 : }
969 74 : double dmax=d0+r0*std::pow(0.00001,1./(nn-mm));
970 148 : function = switchContainers::rationalFactory(d0,dmax,r0,nn,mm);
971 74 : function->setupStretch();
972 74 : }
973 :
974 32 : double SwitchingFunction::get_r0() const {
975 32 : return function->get_r0();
976 : }
977 :
978 8 : double SwitchingFunction::get_d0() const {
979 8 : return function->get_d0();
980 : }
981 :
982 536580542 : double SwitchingFunction::get_dmax() const {
983 536580542 : return function->get_dmax();
984 : }
985 :
986 49030642 : double SwitchingFunction::get_dmax2() const {
987 49030642 : return function->get_dmax2();
988 : }
989 :
990 : }// Namespace PLMD
|
