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 :
31 : #pragma GCC diagnostic error "-Wswitch"
32 : /*
33 : IMPORTANT NOTE FOR DEVELOPERS:
34 :
35 : If you add a new type of switching function in this file please add documentation for your new switching function type in function/LessThan.cpp
36 : */
37 :
38 : namespace PLMD {
39 : namespace switchContainers {
40 :
41 : std::string description(switchType type, const Data& data);
42 :
43 2443 : Switch::~Switch()=default;
44 :
45 : template <typename SF>
46 : class SwitchInterface :public Switch {
47 : switchType type;
48 : Data data;
49 : public:
50 2383 : explicit SwitchInterface(const std::pair <switchType,Data>& d):
51 2383 : type(d.first),
52 2383 : data(d.second) {}
53 129799979 : double calculate(double distance, double& dfunc) const override {
54 129799979 : auto [f,d] = SF::calculate(data,distance);
55 129799979 : dfunc=d;
56 129799979 : return f;
57 : }
58 85028956 : double calculateSqr(double distance2, double& dfunc) const override {
59 85028956 : auto [f,d] = SF::calculateSqr(data,distance2);
60 85028956 : dfunc=d;
61 85028956 : return f;
62 : }
63 289 : void setupStretch() override {
64 289 : if(data.dmax<std::numeric_limits<double>::max()) {
65 289 : data.stretch=1.0;
66 289 : data.shift=0.0;
67 289 : double s0=SF::calculate(data,0.0).first;
68 289 : double sd=SF::calculate(data,data.dmax).first;
69 289 : data.stretch=1.0/(s0-sd);
70 289 : data.shift=-sd*data.stretch;
71 : }
72 289 : }
73 589016459 : const Data& getData() const override {
74 589016459 : return data;
75 : }
76 322 : switchType getType() const override {
77 322 : return type;
78 : }
79 1937 : std::string description() const override {
80 1937 : return switchContainers::description(type,data);
81 : }
82 : };
83 :
84 2132 : Data Data::init(const double D0,const double DMAX, const double R0) {
85 2132 : Data data ;
86 2132 : data.d0=D0;
87 2132 : data.dmax=DMAX;
88 2132 : data.dmax_2=(data.dmax<std::sqrt(std::numeric_limits<double>::max()))
89 2132 : ? data.dmax*data.dmax
90 : : std::numeric_limits<double>::max();
91 2132 : data.invr0=1.0/R0;
92 2132 : data.invr0_2=data.invr0*data.invr0;
93 2132 : return data;
94 : }
95 :
96 0 : void Data::toACCDevice() const {
97 : #pragma acc enter data copyin(this[0:1], d0, dmax, dmax_2, invr0, invr0_2, \
98 : stretch, shift, nn, mm, preRes, preDfunc, preSecDev, \
99 : nnf, mmf, preDfuncF, preSecDevF, a, b, c, d, beta, \
100 : lambda, ref )
101 0 : }
102 0 : void Data::removeFromACCDevice() const {
103 : #pragma acc exit data delete(ref, lambda, beta, d, c, b, a, preSecDevF, preDfuncF, \
104 : mmf, nnf, preSecDev, preDfunc, preRes, mm, nn, shift, \
105 : stretch, invr0_2, invr0, dmax_2, dmax, d0, this[0:1])
106 0 : }
107 :
108 : std::string typeToString(switchType type);
109 : namespace switchContainersUtils {
110 : template<class, class = void>
111 : constexpr bool has_function_data = false;
112 :
113 : //this verifies that T has a method gatherForces_custom that can be called with this signature
114 : template<class T>
115 : constexpr bool has_function_data <T, std::void_t<
116 : decltype(T::function(
117 : std::declval<const Data& >(),
118 : std::declval<double >()
119 : ))> > = true;
120 : } //namespace switchContainersUtils
121 :
122 : using ValueDerivative=std::pair<double,double>;
123 :
124 : inline ValueDerivative applystretch(const Data&data,double distance,ValueDerivative in) {
125 62790085 : in.first=in.first*data.stretch+data.shift;
126 62790085 : in.second*=data.stretch;
127 62790085 : in.second *= data.invr0;
128 62790085 : in.second /= distance;
129 : return in;
130 : }
131 :
132 : /// container for the actual switching function used by PLMD::SwitchingFunction
133 : template<typename switching>
134 : struct baseSwitch {
135 163945014 : static ValueDerivative calculate(const Data&data, const double distance) {
136 163945014 : const double rdist = (distance-data.d0)*data.invr0;
137 : if constexpr (switchContainersUtils::has_function_data<switching>) {
138 155535467 : return (distance > data.dmax) ? ValueDerivative{0.0,0.0}
139 : :
140 155504479 : (rdist > 0.0) ?
141 58335639 : applystretch(data,distance,switching::function(data,rdist))
142 97168840 : : ValueDerivative{data.stretch+data.shift,0.0};
143 : } else {
144 8409547 : return (distance > data.dmax) ? ValueDerivative{0.0,0.0}
145 : :
146 7443219 : (rdist > 0.0) ?
147 4454446 : applystretch(data,distance,switching::function(rdist))
148 2988773 : : ValueDerivative{data.stretch+data.shift,0.0};
149 : }
150 : /* //not branchless (legacy) implementation
151 : double res = 0.0;
152 : double dfunc = 0.0;
153 :
154 : if(distance <= data.dmax) {
155 : res = 1.0;
156 : const double rdist = (distance-data.d0)*data.invr0;
157 : if(rdist > 0.0) {
158 : if constexpr (switchContainersUtils::has_function_data<switching>) {
159 : std::tie(res,dfunc) = switching::function(data,rdist);
160 : } else {
161 : std::tie(res,dfunc) = switching::function(rdist);
162 : }
163 : //the following comments came from the original
164 : // this is for the chain rule (derivative of rdist):
165 : dfunc *= data.invr0;
166 : // for any future switching functions, be aware that multiplying invr0 is only
167 : // correct for functions of rdist = (r-d0)/r0.
168 :
169 : // this is because calculate() sets dfunc to the derivative divided times the
170 : // distance.
171 : // (I think this is misleading and I would like to modify it - GB)
172 : dfunc /= distance;
173 : }
174 : res=res*data.stretch+data.shift;
175 : dfunc*=data.stretch;
176 : }
177 : return {res,dfunc};
178 : */
179 : }
180 :
181 30938593 : static ValueDerivative calculateSqr(const Data&data, double distance2) {
182 30938593 : return switching::calculate(data,std::sqrt(distance2));
183 : }
184 : };
185 : template<int N,
186 : std::enable_if_t< (N >0), bool> = true,
187 : std::enable_if_t< (N %2 == 0), bool> = true>
188 : struct fixedRational :public baseSwitch<fixedRational<N>> {
189 :
190 : template <int POW>
191 1574 : static inline ValueDerivative doRational(const double rdist, double dfunc=0.0, double result=0.0) {
192 : const double rNdist=Tools::fastpow<POW-1>(rdist);
193 28058762 : result=1.0/(1.0+rNdist*rdist);
194 1574 : dfunc = -POW*rNdist*result*result;
195 1574 : return {result,dfunc};
196 : }
197 :
198 16461419 : static inline ValueDerivative function(const Data&, double rdist) {
199 : //preRes and preDfunc are passed already set
200 1574 : return doRational<N>(rdist);
201 : }
202 :
203 11810403 : static ValueDerivative calculateSqr(const Data& data,double distance2) {
204 11810403 : double result=0.0;
205 11810403 : double dfunc=0.0;
206 11810403 : if(distance2 <= data.dmax_2) {
207 11595769 : const double rdist = distance2*data.invr0_2;
208 : std::tie(result,dfunc) = doRational<N/2>(rdist);
209 11595769 : dfunc*=2*data.invr0_2;
210 : // stretch:
211 11595769 : result=result*data.stretch+data.shift;
212 11595769 : dfunc*=data.stretch;
213 : }
214 11810403 : return {result,dfunc};
215 : }
216 : };
217 :
218 : //these enums are useful for clarifying the settings in the factory
219 : //and the code is autodocumented ;)
220 : enum class rationalPow:bool {standard, fast};
221 : enum class rationalForm:bool {standard, simplified};
222 : template<rationalPow isFast, rationalForm nis2m>
223 : struct rational : public baseSwitch<rational<isFast,nis2m>> {
224 : //I am using PLMD::epsilon to be certain to call the one defined in Tools.h
225 : static constexpr double moreThanOne=1.0+5.0e10*PLMD::epsilon;
226 : static constexpr double lessThanOne=1.0-5.0e10*PLMD::epsilon;
227 252 : static std::pair <switchType,Data> init(double D0,double DMAX, double R0, int N, int M) {
228 252 : auto data= Data::init(D0,DMAX,R0);
229 : data.nn=N;
230 252 : data.mm = (M==0) ? N*2 : M;
231 252 : data.preRes=static_cast<double>(data.nn)/data.mm;
232 252 : data.preDfunc=0.5*data.nn*(data.nn-data.mm)/static_cast<double>(data.mm);
233 : //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);
234 252 : data.preSecDev = (data.nn * (data.mm * data.mm - 3.0* data.mm * (-1 + data.nn ) + data.nn *(-3 + 2* data.nn )))/(6.0* data.mm );
235 252 : data.nnf = data.nn/2;
236 252 : data.mmf = data.mm/2;
237 252 : data.preDfuncF = 0.5*data.nnf*(data.nnf-data.mmf)/static_cast<double>(data.mmf);
238 252 : data.preSecDevF = (data.nnf* (data.mmf*data.mmf - 3.0* data.mmf* (-1 + data.nnf) + data.nnf*(-3 + 2* data.nnf)))/(6.0* data.mmf);
239 :
240 252 : const bool fast = N%2==0 && M%2==0 && D0==0.0;
241 252 : if(2*N==M || M == 0) {
242 155 : if(fast) {
243 : return {switchType::rationalSimpleFast,data};
244 : }
245 : return {switchType::rationalSimple,data};
246 : }
247 97 : if(fast) {
248 : return {switchType::rationalFast,data};
249 : }
250 : return {switchType::rational,data};
251 : }
252 :
253 19334191 : static inline ValueDerivative doRational(const double rdist,double secDev, const int N,
254 : const int M, double dfunc=0.0,double result=0.0) {
255 : //the result and dfunc are assigned in the drivers for doRational
256 : //if(rdist>(1.0-100.0*epsilon) && rdist<(1.0+100.0*epsilon)) {
257 : //result=preRes;
258 : //dfunc=preDfunc;
259 : //} else {
260 : if constexpr (nis2m==rationalForm::simplified) {
261 2986057 : const double rNdist=Tools::fastpow(rdist,N-1);
262 2986057 : result=1.0/(1.0+rNdist*rdist);
263 2986057 : dfunc = -N*rNdist*result*result;
264 : } else {
265 16348134 : if(!((rdist > lessThanOne) && (rdist < moreThanOne))) {
266 16348115 : const double rNdist=Tools::fastpow(rdist,N-1);
267 16348115 : const double rMdist=Tools::fastpow(rdist,M-1);
268 16348115 : const double num = 1.0-rNdist*rdist;
269 16348115 : const double iden = 1.0/(1.0-rMdist*rdist);
270 16348115 : result = num*iden;
271 16348115 : dfunc = ((M*result*rMdist)-(N*rNdist))*iden;
272 16348115 : } else {
273 : //here I imply that the correct initialized are being passed to doRational
274 19 : const double x =(rdist-1.0);
275 19 : result = result+ x * ( dfunc + 0.5 * x * secDev);
276 19 : dfunc = dfunc + x * secDev;
277 : }
278 : }
279 19334191 : return {result,dfunc};
280 : }
281 19334111 : static inline ValueDerivative function(const Data&data, double rdist) {
282 : //preRes and preDfunc are passed already set
283 19334111 : return doRational(rdist,data.preSecDev,data.nn,data.mm,data.preDfunc,data.preRes);
284 : }
285 :
286 3645002 : static ValueDerivative calculateSqr(const Data&data, double distance2) {
287 : if constexpr (isFast==rationalPow::fast) {
288 100 : double result=0.0;
289 100 : double dfunc=0.0;
290 100 : if(distance2 <= data.dmax_2) {
291 80 : const double rdist = distance2*data.invr0_2;
292 : std::tie(result,dfunc) =
293 80 : doRational(rdist,data.preSecDevF,data.nnf,data.mmf,data.preDfuncF,data.preRes);
294 80 : dfunc*=2*data.invr0_2;
295 : // stretch:
296 80 : result=result*data.stretch+data.shift;
297 80 : dfunc*=data.stretch;
298 : }
299 100 : return {result,dfunc};
300 : } else {
301 3644902 : return baseSwitch<rational<isFast,nis2m>>::calculate(data,std::sqrt(distance2));
302 : }
303 : }
304 : };
305 :
306 : std::unique_ptr<Switch>
307 616 : rationalFactory (double D0,double DMAX, double R0, int N, int M) {
308 616 : bool fast = N%2==0 && M%2==0 && D0==0.0;
309 : //if (M==0) M will automatically became 2*NN
310 : constexpr int highestPrecompiledPower=12;
311 : //precompiled rational
312 616 : if(((2*N)==M || M == 0) && fast && N<=highestPrecompiledPower) {
313 : #define FIXEDRATIONALENUM(x) case x: \
314 : return std::make_unique<SwitchInterface<fixedRational<x>>>( \
315 : std::pair <switchType,Data> {switchType::rationalfix##x,Data::init(D0,DMAX,R0)});
316 370 : switch (N) {
317 32 : FIXEDRATIONALENUM(12)
318 0 : FIXEDRATIONALENUM(10)
319 2 : FIXEDRATIONALENUM(8)
320 351 : FIXEDRATIONALENUM(6)
321 1 : FIXEDRATIONALENUM(4)
322 0 : FIXEDRATIONALENUM(2)
323 : default:
324 : break;
325 : }
326 : #undef FIXEDRATIONALENUM
327 : }
328 : #undef FIXEDRATIONALENUM
329 : //continue with the 'at runtime implementation'
330 246 : auto data = rational<rationalPow::standard,rationalForm::standard>::init(D0,DMAX,R0,N,M);
331 246 : if(2*N==M || M == 0) {
332 153 : if(fast) {
333 : //fast rational
334 : return PLMD::Tools::make_unique<SwitchInterface<rational<
335 0 : rationalPow::fast,rationalForm::simplified>>>(data);
336 : }
337 : return PLMD::Tools::make_unique<SwitchInterface<rational<
338 153 : rationalPow::standard,rationalForm::simplified>>>(data);
339 : }
340 93 : if(fast) {
341 : //fast rational
342 : return PLMD::Tools::make_unique<SwitchInterface<rational<
343 85 : rationalPow::fast,rationalForm::standard>>>(data);
344 : }
345 : return PLMD::Tools::make_unique<SwitchInterface<rational<
346 8 : rationalPow::standard,rationalForm::standard>>>(data);
347 : }
348 :
349 : std::pair <switchType,Data>
350 8 : initRational(double D0,double DMAX, double R0, int N, int M) {
351 8 : bool fast = N%2==0 && M%2==0 && D0==0.0;
352 : //if (M==0) M will automatically became 2*NN
353 : constexpr int highestPrecompiledPower=12;
354 : //precompiled rational
355 8 : if(((2*N)==M || M == 0) && fast && N<=highestPrecompiledPower) {
356 : #define FIXEDRATIONALENUM(x) case x: \
357 : return std::pair <switchType,Data> {switchType::rationalfix##x,Data::init(D0,DMAX,R0)};
358 2 : switch (N) {
359 0 : FIXEDRATIONALENUM(12)
360 0 : FIXEDRATIONALENUM(10)
361 0 : FIXEDRATIONALENUM(8)
362 2 : FIXEDRATIONALENUM(6)
363 0 : FIXEDRATIONALENUM(4)
364 0 : FIXEDRATIONALENUM(2)
365 : default:
366 : break;
367 : }
368 : #undef FIXEDRATIONALENUM
369 : }
370 : #undef FIXEDRATIONALENUM
371 : //continue with the 'at runtime implementation'
372 6 : return rational<rationalPow::standard,rationalForm::standard>::init(D0,DMAX,R0,N,M);
373 : }
374 :
375 : struct exponentialSwitch: public baseSwitch<exponentialSwitch> {
376 85 : static std::pair <switchType,Data> init(
377 : const double D0,
378 : const double DMAX,
379 : const double R0) {
380 : return{switchType::exponential,
381 85 : Data::init(D0,DMAX,R0)};
382 : }
383 2997089 : static inline ValueDerivative function(const double rdist) {
384 2997089 : double result = std::exp(-rdist);
385 2997089 : return {result,-result};
386 : }
387 : };
388 :
389 : struct gaussianSwitch: public baseSwitch<gaussianSwitch> {
390 82 : static std::pair <switchType,Data> init(
391 : const double D0,
392 : const double DMAX,
393 : const double R0) {
394 : return{switchType::gaussian,
395 82 : Data::init(D0,DMAX,R0)};
396 : }
397 1313645 : static inline ValueDerivative function(const double rdist) {
398 1313645 : double result = std::exp(-0.5*rdist*rdist);
399 1313645 : return {result,-rdist*result};
400 : }
401 : };
402 :
403 : struct fastgaussianSwitch: public baseSwitch<fastgaussianSwitch> {
404 77 : static std::pair <switchType,Data> init(
405 : const double DMAX) {
406 : return{switchType::fastgaussian,
407 77 : Data::init(0.0,DMAX,1.0)};
408 : }
409 41 : static inline ValueDerivative function(const double rdist) {
410 41 : double result = std::exp(-0.5*rdist*rdist);
411 41 : return {result,-rdist*result};
412 : }
413 38635198 : static inline ValueDerivative calculateSqr(const Data& data,const double distance2) {
414 : double result = 0.0;
415 : double dfunc = 0.0;
416 38635198 : if(distance2<data.dmax_2) {
417 : result=1.0;
418 38522934 : if(distance2 >0.0) {
419 38522868 : result = exp(-0.5*distance2);
420 : dfunc = -result;
421 : //I have to multiply and then divide dfunc by sqrt(distance2)
422 : //by omitting that I get a funny wrong value of the derivative in 0
423 : //by not omitting that I get a funny nan because I divide by 0
424 : //hence the extra if
425 : // stretch:
426 38522868 : result=result*data.stretch+data.shift;
427 : dfunc*=data.stretch;
428 : }
429 : }
430 38635198 : return {result,dfunc};
431 : }
432 : };
433 :
434 : struct smapSwitch: public baseSwitch<smapSwitch> {
435 19 : static std::pair <switchType,Data> init(
436 : const double D0,
437 : const double DMAX,
438 : const double R0,
439 : const int A,
440 : const int B) {
441 19 : auto data =Data::init(D0,DMAX,R0);
442 : data.a= A;
443 : data.b= B;
444 19 : data.c= std::pow(2., static_cast<double>(data.a)/static_cast<double>(data.b) ) - 1.0;
445 19 : data.d= -static_cast<double>(data.b) / static_cast<double>(data.a);
446 19 : return{switchType::smap,data};
447 : }
448 :
449 21911376 : static inline ValueDerivative function(const Data& data,const double rdist) {
450 21911376 : const double sx=data.c*Tools::fastpow( rdist, data.a );
451 21911376 : double result=std::pow( 1.0 + sx, data.d );
452 21911376 : double dfunc=-data.b*sx/rdist*result/(1.0+sx);
453 21911376 : return {result,dfunc};
454 : }
455 : };
456 :
457 : struct cubicSwitch: public baseSwitch<cubicSwitch> {
458 20 : static std::pair <switchType,Data> init(const double D0, const double DMAX) {
459 : return{switchType::cubic,
460 20 : Data::init(D0,DMAX,DMAX-D0)};
461 : }
462 127299 : static inline ValueDerivative function(const double rdist) {
463 127299 : const double tmp1 = rdist - 1.0;
464 127299 : const double tmp2 = 1.0+2.0*rdist;
465 127299 : double dfunc = 2*tmp1*tmp2 + 2*tmp1*tmp1;
466 127299 : return {tmp1*tmp1*tmp2,dfunc};
467 : }
468 : };
469 :
470 : struct tanhSwitch: public baseSwitch<tanhSwitch> {
471 9 : static std::pair <switchType,Data> init(
472 : const double D0,
473 : const double DMAX,
474 : const double R0) {
475 : return{switchType::tanh,
476 9 : Data::init(D0,DMAX,R0)};
477 : }
478 16372 : static inline ValueDerivative function(const double rdist) {
479 16372 : const double tmp1 = std::tanh(rdist);
480 16372 : double dfunc = tmp1 * tmp1 - 1.0;
481 16372 : return {1.0 - tmp1,dfunc };
482 : }
483 : };
484 :
485 : struct cosinusSwitch: public baseSwitch<cosinusSwitch> {
486 9 : static std::pair <switchType,Data> init(
487 : const double D0,
488 : const double DMAX,
489 : const double R0) {
490 : return{switchType::cosinus,
491 9 : Data::init(D0,DMAX,R0)};
492 : }
493 627159 : static inline ValueDerivative function(const Data& data,const double rdist) {
494 : double result = 0.0;
495 : double dfunc=0.0;
496 627159 : if(rdist<=1.0) {
497 : // rdist = (r-r1)/(r2-r1) ; 0.0<=rdist<=1.0 if r1 <= r <=r2; (r2-r1)/(r2-r1)=1
498 297044 : double rdistPI = rdist * PLMD::pi;
499 297044 : result = 0.5 * (std::cos ( rdistPI ) + 1.0);
500 297044 : dfunc = -0.5 * PLMD::pi * std::sin ( rdistPI );
501 : }
502 627159 : return {result,dfunc};
503 : }
504 : };
505 :
506 : struct nativeqSwitch: public baseSwitch<nativeqSwitch> {
507 1160 : static std::pair <switchType,Data> init(
508 : const double D0,
509 : const double DMAX,
510 : const double R0,
511 : const double BETA, // nm-1
512 : const double LAMBDA,// unitless
513 : const double REF) {
514 1160 : auto data=Data::init(D0,DMAX,R0);
515 : data.beta=BETA;
516 : data.lambda=LAMBDA;
517 : data.ref=REF;
518 1160 : return{switchType::nativeq,data};
519 : }
520 439300 : static inline ValueDerivative calculate(const Data& data,const double distance) {
521 : double res = 0.0;
522 : double dfunc = 0.0;
523 439300 : if(distance<=data.dmax) {
524 : res = 1.0;
525 439268 : if(distance > data.d0) {
526 439251 : const double rdist = data.beta*(distance - data.lambda * data.ref);
527 439251 : double exprdist=std::exp(rdist);
528 439251 : res=1.0/(1.0+exprdist);
529 : /*2.9
530 : //need to see if this (5op+assign)
531 : //double exprmdist=1.0 + exprdist;
532 : //dfunc = - (beta *exprdist)/(exprmdist*exprmdist);
533 : //or this (5op but 2 divisions) is faster
534 : dfunc = - beta /(exprdist+ 2.0 +1.0/exprdist);
535 : //this cames from - beta * exprdist/(exprdist*exprdist+ 2.0 *exprdist +1.0)
536 : //dfunc *= invr0;
537 : dfunc /= distance;
538 : */
539 : //2.10
540 439251 : dfunc = - data.beta /(exprdist+ 2.0 +1.0/exprdist) /distance;
541 :
542 439251 : dfunc*=data.stretch;
543 : }
544 439268 : res=res*data.stretch+data.shift;
545 : }
546 439300 : return {res,dfunc};
547 : }
548 : };
549 :
550 : class leptonSwitch {
551 : /// Lepton expression.
552 182 : class funcAndDeriv {
553 : lepton::CompiledExpression expression;
554 : lepton::CompiledExpression deriv;
555 : double* varRef=nullptr;
556 : double* varDevRef=nullptr;
557 : public:
558 60 : funcAndDeriv(const std::string &func) {
559 60 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
560 60 : expression=pe.createCompiledExpression();
561 62 : std::string arg="x";
562 :
563 : {
564 60 : auto vars=expression.getVariables();
565 60 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
566 60 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
567 :
568 60 : if(found_x2) {
569 : arg="x2";
570 : }
571 60 : if (vars.size()==0) {
572 : // this is necessary since in some cases lepton thinks a variable is not present even though it is present
573 : // e.g. func=0*x
574 0 : varRef=nullptr;
575 60 : } else if(vars.size()==1 && (found_x || found_x2)) {
576 58 : varRef=&expression.getVariableReference(arg);
577 : } else {
578 4 : plumed_error()
579 : <<"Please declare a function with only ONE argument that can only be x or x2. Your function is: "
580 4 : << func;
581 : }
582 : }
583 :
584 60 : lepton::ParsedExpression ped=lepton::Parser::parse(func).differentiate(arg).optimize(lepton::Constants());
585 58 : deriv=ped.createCompiledExpression();
586 : {
587 58 : auto vars=expression.getVariables();
588 58 : if (vars.size()==0) {
589 0 : varDevRef=nullptr;
590 : } else {
591 58 : varDevRef=&deriv.getVariableReference(arg);
592 : }
593 : }
594 :
595 62 : }
596 124 : funcAndDeriv (const funcAndDeriv& other):
597 124 : expression(other.expression),
598 124 : deriv(other.deriv) {
599 124 : std::string arg="x";
600 :
601 : {
602 124 : auto vars=expression.getVariables();
603 124 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
604 124 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
605 :
606 124 : if(found_x2) {
607 : arg="x2";
608 : }
609 124 : if (vars.size()==0) {
610 0 : varRef=nullptr;
611 124 : } else if(vars.size()==1 && (found_x || found_x2)) {
612 124 : varRef=&expression.getVariableReference(arg);
613 : }// UB: I assume that the function is already correct
614 : }
615 :
616 : {
617 124 : auto vars=expression.getVariables();
618 124 : if (vars.size()==0) {
619 0 : varDevRef=nullptr;
620 : } else {
621 124 : varDevRef=&deriv.getVariableReference(arg);
622 : }
623 : }
624 124 : }
625 :
626 : funcAndDeriv& operator= (const funcAndDeriv& other) {
627 : if(this != &other) {
628 : expression = other.expression;
629 : deriv = other.deriv;
630 : std::string arg="x";
631 :
632 : {
633 : auto vars=expression.getVariables();
634 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
635 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
636 :
637 : if(found_x2) {
638 : arg="x2";
639 : }
640 : if (vars.size()==0) {
641 : varRef=nullptr;
642 : } else if(vars.size()==1 && (found_x || found_x2)) {
643 : varRef=&expression.getVariableReference(arg);
644 : }// UB: I assume that the function is already correct
645 : }
646 :
647 : {
648 : auto vars=expression.getVariables();
649 : if (vars.size()==0) {
650 : varDevRef=nullptr;
651 : } else {
652 : varDevRef=&deriv.getVariableReference(arg);
653 : }
654 : }
655 : }
656 : return *this;
657 : }
658 :
659 6753859 : std::pair<double,double> operator()(double const x) const {
660 : //FAQ: why this works? this thing is const and you are modifying things!
661 : //Actually I am modifying something that is pointed at, not my pointers,
662 : //so I am not mutating the state of this!
663 6753859 : if(varRef) {
664 6753859 : *varRef=x;
665 : }
666 6753859 : if(varDevRef) {
667 6753859 : *varDevRef=x;
668 : }
669 : return std::make_pair(
670 6753859 : expression.evaluate(),
671 6753859 : deriv.evaluate());
672 : }
673 :
674 : auto& getVariables() const {
675 58 : return expression.getVariables();
676 : }
677 : auto& getVariables_derivative() const {
678 : return deriv.getVariables();
679 : }
680 : };
681 : /// Function for lepton
682 : std::string lepton_func;
683 : /// \warning Since lepton::CompiledExpression is mutable, a vector is necessary for multithreading!
684 : std::vector <funcAndDeriv> expressions{};
685 : /// Set to true if lepton only uses x2
686 : bool leptonx2=false;
687 : Data data;
688 : public:
689 49 : leptonSwitch(double D0, double DMAX, double R0, const std::string & func)
690 49 : :lepton_func(func),
691 49 : expressions(OpenMP::getNumThreads(), lepton_func),
692 47 : data(Data::init(D0,DMAX,R0)) {
693 : //this is a bit odd, but it works
694 : auto vars=expressions[0].getVariables();
695 47 : leptonx2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
696 47 : }
697 :
698 11 : leptonSwitch(const leptonSwitch& other)
699 11 : :lepton_func(other.lepton_func),
700 11 : expressions(OpenMP::getNumThreads(), lepton_func),
701 11 : data(other.data) {
702 : //this is a bit odd, but it works
703 : auto vars=expressions[0].getVariables();
704 11 : leptonx2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
705 11 : }
706 :
707 : // the implicit assignment is what is needed (the copy constructor above only rebuilds the expressions)
708 : leptonSwitch& operator=(const leptonSwitch& other) = default;
709 :
710 6106994 : inline double calculate(const double distance,double&dfunc) const {
711 6106994 : double res = 0.0;
712 6106994 : dfunc = 0.0;
713 6106994 : if(leptonx2) {
714 2 : res= calculateSqr(distance*distance,dfunc);
715 : } else {
716 6106992 : if(distance<=data.dmax) {
717 5811749 : res = 1.0;
718 5811749 : const double rdist = (distance-data.d0)*data.invr0;
719 5811749 : if(rdist > 0.0) {
720 5505757 : const unsigned t=OpenMP::getThreadNum();
721 5505757 : plumed_assert(t<expressions.size());
722 5505757 : std::tie(res,dfunc) = expressions[t](rdist);
723 5505757 : dfunc *= data.invr0;
724 5505757 : dfunc /= distance;
725 : }
726 5811749 : res=res*data.stretch+data.shift;
727 5811749 : dfunc*=data.stretch;
728 : }
729 : }
730 6106994 : return res;
731 : }
732 :
733 7354820 : double calculateSqr(const double distance2,double&dfunc) const {
734 7354820 : double result =0.0;
735 7354820 : dfunc=0.0;
736 7354820 : if(leptonx2) {
737 1248110 : if(distance2<=data.dmax_2) {
738 1248102 : const unsigned t=OpenMP::getThreadNum();
739 1248102 : const double rdist_2 = distance2*data.invr0_2;
740 1248102 : plumed_assert(t<expressions.size());
741 1248102 : std::tie(result,dfunc) = expressions[t](rdist_2);
742 : // chain rule:
743 1248102 : dfunc*=2*data.invr0_2;
744 : // stretch:
745 1248102 : result=result*data.stretch+data.shift;
746 1248102 : dfunc*=data.stretch;
747 : }
748 : } else {
749 6106710 : result = calculate(std::sqrt(distance2),dfunc);
750 : }
751 7354820 : return result;
752 : }
753 : const Data & getData() const {
754 22 : return data;
755 : }
756 22 : void setupStretch() {
757 22 : if(data.dmax!=std::numeric_limits<double>::max()) {
758 22 : data.stretch=1.0;
759 22 : data.shift=0.0;
760 : double dummy;
761 22 : double s0=calculate(0.0,dummy);
762 22 : double sd=calculate(data.dmax,dummy);
763 22 : data.stretch=1.0/(s0-sd);
764 22 : data.shift=-sd*data.stretch;
765 : }
766 22 : }
767 : void removeStretch() {
768 : data.stretch=1.0;
769 : data.shift=0.0;
770 : }
771 23 : std::string description() const {
772 : using namespace switchContainers;
773 23 : std::ostringstream ostr;
774 23 : ostr<<1.0/data.invr0
775 : <<". Using "
776 46 : << typeToString(switchType::lepton)
777 46 : <<" switching function with parameters d0="<< data.d0
778 23 : <<" func=" << lepton_func;
779 23 : return ostr.str();
780 23 : }
781 : };
782 :
783 : //call to calculate with no inheritance
784 262 : ValueDerivative calculate(const switchType type,
785 : const Data& data,
786 : const double rdist) {
787 : #define SWITCHCALL(x) case switchType::x: return x##Switch::calculate(data, rdist);
788 : #define RATCALL(x) case switchType::rationalfix##x:return fixedRational<x>::calculate(data, rdist);
789 262 : switch (type) {
790 0 : RATCALL(12)
791 0 : RATCALL(10)
792 0 : RATCALL(8)
793 22 : RATCALL(6)
794 0 : RATCALL(4)
795 0 : RATCALL(2)
796 22 : case switchType::rational:
797 22 : return rational<rationalPow::standard,rationalForm::standard>::calculate(data, rdist);
798 22 : case switchType::rationalFast:
799 22 : return rational<rationalPow::fast,rationalForm::standard>::calculate(data, rdist);
800 22 : case switchType::rationalSimple:
801 22 : return rational<rationalPow::standard,rationalForm::simplified>::calculate(data, rdist);
802 0 : case switchType::rationalSimpleFast:
803 0 : return rational<rationalPow::fast,rationalForm::simplified>::calculate(data, rdist);
804 22 : SWITCHCALL(exponential)
805 22 : SWITCHCALL(gaussian)
806 22 : SWITCHCALL(fastgaussian)
807 22 : SWITCHCALL(smap)
808 22 : SWITCHCALL(cubic)
809 22 : SWITCHCALL(tanh)
810 20 : SWITCHCALL(cosinus)
811 22 : SWITCHCALL(nativeq)
812 : default:
813 : break;
814 : }
815 : #undef SWITCHCALL
816 : #undef RATCALL
817 0 : return {0.0,0.0};
818 : }
819 :
820 : //call to calculateSqr with no inheritance
821 240 : ValueDerivative calculateSqr(const switchType type,
822 : const Data& data,
823 : const double rdist2) {
824 : #define SWITCHCALL(x) case switchType::x: return x##Switch::calculateSqr(data, rdist2);
825 : #define RATCALL(x) case switchType::rationalfix##x:return fixedRational<x>::calculateSqr(data, rdist2);
826 240 : switch (type) {
827 0 : RATCALL(12)
828 0 : RATCALL(10)
829 0 : RATCALL(8)
830 20 : RATCALL(6)
831 0 : RATCALL(4)
832 0 : RATCALL(2)
833 20 : case switchType::rational:
834 20 : return rational<rationalPow::standard,rationalForm::standard>::calculateSqr(data, rdist2);
835 20 : case switchType::rationalFast:
836 20 : return rational<rationalPow::fast,rationalForm::standard>::calculateSqr(data, rdist2);
837 20 : case switchType::rationalSimple:
838 20 : return rational<rationalPow::standard,rationalForm::simplified>::calculateSqr(data, rdist2);
839 0 : case switchType::rationalSimpleFast:
840 0 : return rational<rationalPow::fast,rationalForm::simplified>::calculateSqr(data, rdist2);
841 20 : SWITCHCALL(exponential)
842 20 : SWITCHCALL(gaussian)
843 20 : SWITCHCALL(fastgaussian)
844 20 : SWITCHCALL(smap)
845 20 : SWITCHCALL(cubic)
846 20 : SWITCHCALL(tanh)
847 20 : SWITCHCALL(cosinus)
848 20 : SWITCHCALL(nativeq)
849 : default:
850 : break;
851 : }
852 : #undef SWITCHCALL
853 : #undef RATCALL
854 0 : return {0.0,0.0};
855 : }
856 :
857 : //call to setupStretch with no inheritance
858 11 : void setupStretch(switchType type, Data& data) {
859 11 : if(data.dmax!=std::numeric_limits<double>::max()) {
860 11 : data.stretch=1.0;
861 11 : data.shift=0.0;
862 11 : double s0=calculate(type,data,0.0).first;
863 11 : double sd=calculate(type,data,data.dmax).first;
864 11 : data.stretch=1.0/(s0-sd);
865 11 : data.shift=-sd*data.stretch;
866 : }
867 11 : }
868 :
869 0 : void removeStretch(Data& data) {
870 0 : data.stretch=1.0;
871 0 : data.shift=0.0;
872 0 : }
873 :
874 1960 : std::string typeToString(switchType type) {
875 : #define DEFAULTPRINT(x) case switchType::x: return #x;
876 1960 : switch (type) {
877 : case switchType::rationalfix12:
878 : case switchType::rationalfix10:
879 : case switchType::rationalfix8:
880 : case switchType::rationalfix6:
881 : case switchType::rationalfix4:
882 : case switchType::rationalfix2:
883 : case switchType::rationalFast:
884 : case switchType::rationalSimple:
885 : case switchType::rationalSimpleFast:
886 564 : DEFAULTPRINT(rational)
887 81 : DEFAULTPRINT(exponential)
888 78 : DEFAULTPRINT(gaussian)
889 17 : DEFAULTPRINT(fastgaussian)
890 15 : DEFAULTPRINT(smap)
891 16 : DEFAULTPRINT(cubic)
892 5 : DEFAULTPRINT(tanh)
893 5 : DEFAULTPRINT(cosinus)
894 1156 : DEFAULTPRINT(nativeq)
895 23 : DEFAULTPRINT(lepton)
896 : case switchType::not_initialized:
897 0 : return "not initialized!";
898 : }
899 : #undef DEFAULTPRINT
900 0 : return "";
901 : }
902 :
903 1937 : std::string description(switchType type, const Data& data) {
904 1937 : std::ostringstream ostr;
905 1937 : ostr<<1.0/data.invr0
906 : <<". Using "
907 1937 : << typeToString(type)
908 3874 : <<" switching function with parameters d0="<< data.d0;
909 1937 : switch (type) {
910 : #define RATFIX(N) case switchType::rationalfix##N: \
911 : ostr<< " " << " nn=" << N << " mm=" <<N*2; \
912 : break;
913 16 : RATFIX(12)
914 0 : RATFIX(10)
915 2 : RATFIX(8)
916 324 : RATFIX(6)
917 1 : RATFIX(4)
918 0 : RATFIX(2)
919 : #undef RATFIX
920 6 : case switchType::rational:
921 6 : ostr << " nn=" << data.nn << " mm=" << data.mm;
922 : break;
923 15 : case switchType::smap:
924 15 : ostr<<" a="<<data.a<<" b="<<data.b;
925 : break;
926 16 : case switchType::cubic:
927 16 : ostr<<" dmax="<<data.dmax;
928 : break;
929 1156 : case switchType::nativeq:
930 1156 : ostr<<" beta="<<data.beta<<" lambda="<<data.lambda<<" ref="<<data.ref;
931 : break;
932 : default:
933 :
934 : break;
935 : }
936 1937 : return ostr.str();
937 1937 : }
938 :
939 : class SwitchInterface_lepton: public Switch {
940 : leptonSwitch function;
941 : public:
942 49 : SwitchInterface_lepton(double D0, double DMAX, double R0, const std::string & func):
943 49 : function(D0,DMAX,R0,func) {}
944 11 : SwitchInterface_lepton(const SwitchInterface_lepton& other):
945 11 : function(other.function) {}
946 : SwitchInterface_lepton& operator=(const SwitchInterface_lepton& other) = default;
947 240 : double calculate(double distance, double& dfunc) const override {
948 240 : return function.calculate(distance,dfunc);
949 : }
950 7354818 : double calculateSqr(double distance2, double& dfunc) const override {
951 7354818 : return function.calculateSqr(distance2,dfunc);
952 : }
953 22 : void setupStretch() override {
954 22 : function.setupStretch();
955 22 : }
956 22 : const Data& getData() const override {
957 22 : return function.getData();
958 : }
959 11 : switchType getType() const override {
960 11 : return switchType::lepton;
961 : }
962 23 : std::string description() const override {
963 23 : return function.description();
964 : }
965 : };
966 :
967 : } // namespace switchContainers
968 :
969 2659 : SwitchingFunction::SwitchingFunction()=default;
970 2668 : SwitchingFunction::~SwitchingFunction()=default;
971 :
972 1 : SwitchingFunction::SwitchingFunction(const SwitchingFunction& other) {
973 1 : copyFunction(other);
974 1 : }
975 :
976 8 : SwitchingFunction::SwitchingFunction(SwitchingFunction&& other):
977 8 : function(std::move(other.function)) {}
978 :
979 332 : SwitchingFunction& SwitchingFunction::operator=(const SwitchingFunction& other) {
980 332 : if (this != &other) {
981 332 : copyFunction(other);
982 : }
983 332 : return *this;
984 : }
985 :
986 0 : SwitchingFunction& SwitchingFunction::operator=(SwitchingFunction&& other) {
987 : function.reset();
988 : function=std::move(other.function);
989 0 : return *this;
990 : }
991 :
992 333 : void SwitchingFunction::copyFunction(const SwitchingFunction& other) {
993 333 : if (this == &other) {
994 0 : return; // nothing to do
995 : }
996 : function.reset();
997 333 : if (!other.function) {
998 : //now both have the function uninitialized
999 : return; // nothing to copy
1000 : }
1001 : using namespace switchContainers;
1002 333 : const auto settings = std::make_pair(other.function->getType(),other.function->getData());
1003 : #define SWITCHCALL(x) case switchType::x: \
1004 : function = std::make_unique<SwitchInterface<x##Switch>>(settings); \
1005 : break;
1006 : #define RATCALL(x) case switchType::rationalfix##x:\
1007 : function = std::make_unique<SwitchInterface<fixedRational<x>>>(settings); \
1008 : break;
1009 333 : switch (settings.first) {
1010 0 : RATCALL(12)
1011 0 : RATCALL(10)
1012 2 : RATCALL(8)
1013 138 : RATCALL(6)
1014 1 : RATCALL(4)
1015 0 : RATCALL(2)
1016 1 : case switchType::rational:
1017 1 : function = std::make_unique<SwitchInterface<rational<rationalPow::standard,rationalForm::standard>>>(settings);
1018 1 : break;
1019 0 : case switchType::rationalFast:
1020 0 : function = std::make_unique<SwitchInterface<rational<rationalPow::fast,rationalForm::standard>>>(settings);
1021 0 : break;
1022 104 : case switchType::rationalSimple:
1023 104 : function = std::make_unique<SwitchInterface<rational<rationalPow::standard,rationalForm::simplified>>>(settings);
1024 104 : break;
1025 0 : case switchType::rationalSimpleFast:
1026 0 : function = std::make_unique<SwitchInterface<rational<rationalPow::fast,rationalForm::simplified>>>(settings);
1027 0 : break;
1028 28 : SWITCHCALL(exponential)
1029 32 : SWITCHCALL(gaussian)
1030 0 : SWITCHCALL(fastgaussian)
1031 1 : SWITCHCALL(smap)
1032 12 : SWITCHCALL(cubic)
1033 3 : SWITCHCALL(tanh)
1034 0 : SWITCHCALL(cosinus)
1035 0 : SWITCHCALL(nativeq)
1036 11 : case switchType::lepton:
1037 11 : function = std::make_unique<SwitchInterface_lepton>(*dynamic_cast<SwitchInterface_lepton*>(other.function.get()));
1038 : case switchType::not_initialized:
1039 : break;
1040 : }
1041 : #undef SWITCHCALL
1042 : #undef RATCALL
1043 : }
1044 :
1045 :
1046 0 : void SwitchingFunction::registerKeywords( Keywords& keys ) {
1047 0 : keys.add("compulsory","R_0","the value of R_0 in the switching function");
1048 0 : keys.add("compulsory","D_0","0.0","the value of D_0 in the switching function");
1049 0 : keys.add("optional","D_MAX","the value at which the switching function can be assumed equal to zero");
1050 0 : keys.add("compulsory","NN","6","the value of n in the switching function (only needed for TYPE=RATIONAL)");
1051 0 : keys.add("compulsory","MM","0","the value of m in the switching function (only needed for TYPE=RATIONAL); 0 implies 2*NN");
1052 0 : keys.add("compulsory","A","the value of a in the switching function (only needed for TYPE=SMAP)");
1053 0 : keys.add("compulsory","B","the value of b in the switching function (only needed for TYPE=SMAP)");
1054 0 : }
1055 :
1056 2025 : void SwitchingFunction::set(const std::string & definition,std::string& errormsg) {
1057 : function.reset();
1058 2025 : std::vector<std::string> data=Tools::getWords(definition);
1059 : #define CHECKandPARSE(datastring,keyword,variable,errormsg) \
1060 : if(Tools::findKeyword(datastring,keyword) && !Tools::parse(datastring,keyword,variable))\
1061 : errormsg="could not parse " keyword; //adjacent strings are automagically concatenated
1062 : #define REQUIREDPARSE(datastring,keyword,variable,errormsg) \
1063 : if(!Tools::parse(datastring,keyword,variable))\
1064 : errormsg=keyword " is required for " + name ; //adjacent strings are automagically concatenated
1065 :
1066 2025 : if( data.size()<1 ) {
1067 : errormsg="missing all input for switching function";
1068 : return;
1069 : }
1070 2025 : std::string name=data[0];
1071 : data.erase(data.begin());
1072 2025 : double d0=0.0;
1073 2025 : double dmax=std::numeric_limits<double>::max();
1074 2593 : CHECKandPARSE(data,"D_0",d0,errormsg);
1075 2573 : CHECKandPARSE(data,"D_MAX",dmax,errormsg);
1076 :
1077 2025 : bool dostretch=false;
1078 2025 : Tools::parseFlag(data,"STRETCH",dostretch); // this is ignored now
1079 2025 : dostretch=true;
1080 2025 : bool dontstretch=false;
1081 2025 : Tools::parseFlag(data,"NOSTRETCH",dontstretch); // this is ignored now
1082 2025 : if(dontstretch) {
1083 100 : dostretch=false;
1084 : }
1085 : using namespace switchContainers;
1086 2025 : if(name=="CUBIC") {
1087 : //cubic is the only switch type that only uses d0 and dmax
1088 36 : function= std::make_unique<SwitchInterface<cubicSwitch>>(
1089 54 : cubicSwitch::init(d0,dmax));
1090 : } else {
1091 2007 : double r0=0.0;
1092 4014 : REQUIREDPARSE(data,"R_0",r0,errormsg);
1093 2007 : if(name=="RATIONAL") {
1094 530 : int nn=6;
1095 530 : int mm=0;
1096 880 : CHECKandPARSE(data,"NN",nn,errormsg);
1097 864 : CHECKandPARSE(data,"MM",mm,errormsg);
1098 1060 : function = switchContainers::rationalFactory(d0,dmax,r0,nn,mm);
1099 :
1100 1477 : } else if(name=="SMAP") {
1101 17 : int a=0;
1102 17 : int b=0;
1103 : //in the original a and b are "default=0",
1104 : //but you divide by a and b during the initialization!
1105 : //better an error message than an UB, so no default
1106 34 : REQUIREDPARSE(data,"A",a,errormsg);
1107 34 : REQUIREDPARSE(data,"B",b,errormsg);
1108 34 : function= std::make_unique<SwitchInterface<smapSwitch>>(
1109 51 : smapSwitch::init(d0,dmax,r0,a,b));
1110 1460 : } else if(name=="Q") {
1111 1158 : double beta = 50.0; // nm-1
1112 1158 : double lambda = 1.8; // unitless
1113 : double ref;
1114 3474 : CHECKandPARSE(data,"BETA",beta,errormsg);
1115 3474 : CHECKandPARSE(data,"LAMBDA",lambda,errormsg);
1116 2316 : REQUIREDPARSE(data,"REF",ref,errormsg);
1117 : //the original error message was not standard
1118 : // if(!Tools::parse(data,"REF",ref))
1119 : // errormsg="REF (reference distaance) is required for native Q";
1120 2316 : function= std::make_unique<SwitchInterface<nativeqSwitch>>(
1121 3474 : nativeqSwitch::init(d0,dmax,r0,beta,lambda,ref));
1122 302 : } else if(name=="EXP") {
1123 166 : function= std::make_unique<SwitchInterface<exponentialSwitch>>(
1124 249 : exponentialSwitch::init(d0,dmax,r0));
1125 219 : } else if(name=="GAUSSIAN") {
1126 155 : if ( r0==1.0 && d0==0.0 ) {
1127 150 : function= std::make_unique<SwitchInterface<fastgaussianSwitch>>(
1128 225 : fastgaussianSwitch::init(dmax));
1129 : } else {
1130 160 : function= std::make_unique<SwitchInterface<gaussianSwitch>>(
1131 240 : gaussianSwitch::init(d0,dmax,r0));
1132 : }
1133 64 : } else if(name=="TANH") {
1134 14 : function= std::make_unique<SwitchInterface<tanhSwitch>>(
1135 21 : tanhSwitch::init(d0,dmax,r0));
1136 57 : } else if(name=="COSINUS") {
1137 14 : function= std::make_unique<SwitchInterface<cosinusSwitch>>(
1138 21 : cosinusSwitch::init(d0,dmax,r0));
1139 50 : } else if((name=="MATHEVAL" || name=="CUSTOM")) {
1140 : std::string func;
1141 49 : Tools::parse(data,"FUNC",func);
1142 49 : function= std::make_unique<SwitchInterface_lepton>(d0,dmax,r0,func);
1143 : } else {
1144 3 : errormsg="cannot understand switching function type '"+name+"'";
1145 : }
1146 : }
1147 : #undef CHECKandPARSE
1148 : #undef REQUIREDPARSE
1149 :
1150 2023 : if( !data.empty() ) {
1151 : errormsg="found the following rogue keywords in switching function input : ";
1152 0 : for(unsigned i=0; i<data.size(); ++i) {
1153 2 : errormsg = errormsg + data[i] + " ";
1154 : }
1155 : }
1156 :
1157 2023 : if(dostretch && dmax!=std::numeric_limits<double>::max()) {
1158 225 : function->setupStretch();
1159 : }
1160 2025 : }
1161 :
1162 1960 : std::string SwitchingFunction::description() const {
1163 : // if this error is necessary, something went wrong in the constructor
1164 : // plumed_merror("Unknown switching function type");
1165 1960 : return function->description();
1166 : }
1167 :
1168 92383774 : double SwitchingFunction::calculateSqr(double distance2,double&dfunc)const {
1169 92383774 : return function->calculateSqr( distance2, dfunc);
1170 : }
1171 :
1172 129800219 : double SwitchingFunction::calculate(double distance,double&dfunc)const {
1173 129800219 : plumed_massert(function,"you are trying to use an unset SwitchingFunction");
1174 129800219 : return function->calculate( distance, dfunc);
1175 : }
1176 :
1177 86 : void SwitchingFunction::set(const int nn,int mm, const double r0, const double d0) {
1178 : function.reset();
1179 86 : if(mm == 0) {
1180 81 : mm = 2*nn;
1181 : }
1182 86 : double dmax=d0+r0*std::pow(0.00001,1./(nn-mm));
1183 172 : function = switchContainers::rationalFactory(d0,dmax,r0,nn,mm);
1184 86 : function->setupStretch();
1185 86 : }
1186 :
1187 41 : double SwitchingFunction::get_r0() const {
1188 41 : return 1.0/function->getData().invr0;
1189 : }
1190 :
1191 8 : double SwitchingFunction::get_d0() const {
1192 8 : return function->getData().d0;
1193 : }
1194 :
1195 539215961 : double SwitchingFunction::get_dmax() const {
1196 539215961 : return function->getData().dmax;
1197 : }
1198 :
1199 49800138 : double SwitchingFunction::get_dmax2() const {
1200 49800138 : return function->getData().dmax_2;
1201 : }
1202 :
1203 0 : const switchContainers::Data& SwitchingFunction::get_data() const {
1204 0 : return function->getData();
1205 : }
1206 :
1207 0 : void SwitchingFunctionAccelerable::registerKeywords( Keywords& keys ) {
1208 0 : SwitchingFunction::registerKeywords(keys);
1209 0 : }
1210 :
1211 24 : void SwitchingFunctionAccelerable::set(const std::string & definition,std::string& errormsg) {
1212 24 : std::vector<std::string> data=Tools::getWords(definition);
1213 : #define CHECKandPARSE(datastring,keyword,variable,errormsg) \
1214 : if(Tools::findKeyword(datastring,keyword) && !Tools::parse(datastring,keyword,variable))\
1215 : errormsg="could not parse " keyword; //adjacent strings are automagically concatenated
1216 : #define REQUIREDPARSE(datastring,keyword,variable,errormsg) \
1217 : if(!Tools::parse(datastring,keyword,variable))\
1218 : errormsg=keyword " is required for " + name ; //adjacent strings are automagically concatenated
1219 :
1220 24 : if( data.size()<1 ) {
1221 : errormsg="missing all input for switching function";
1222 : return;
1223 : }
1224 24 : std::string name=data[0];
1225 : data.erase(data.begin());
1226 24 : double d0=0.0;
1227 24 : double dmax=std::numeric_limits<double>::max();
1228 44 : CHECKandPARSE(data,"D_0",d0,errormsg);
1229 68 : CHECKandPARSE(data,"D_MAX",dmax,errormsg);
1230 :
1231 24 : bool dostretch=false;
1232 24 : Tools::parseFlag(data,"STRETCH",dostretch); // this is ignored now
1233 24 : dostretch=true;
1234 24 : bool dontstretch=false;
1235 24 : Tools::parseFlag(data,"NOSTRETCH",dontstretch); // this is ignored now
1236 24 : if(dontstretch) {
1237 12 : dostretch=false;
1238 : }
1239 : using namespace switchContainers;
1240 24 : if(name=="CUBIC") {
1241 : //cubic is the only switch type that only uses d0 and dmax
1242 2 : std::tie(type,switchData) = cubicSwitch::init(d0,dmax);
1243 : } else {
1244 22 : double r0=0.0;
1245 44 : REQUIREDPARSE(data,"R_0",r0,errormsg);
1246 22 : if(name=="RATIONAL") {
1247 8 : int nn=6;
1248 8 : int mm=0;
1249 20 : CHECKandPARSE(data,"NN",nn,errormsg);
1250 16 : CHECKandPARSE(data,"MM",mm,errormsg);
1251 8 : std::tie(type,switchData) = switchContainers::initRational(d0,dmax,r0,nn,mm);
1252 :
1253 14 : } else if(name=="SMAP") {
1254 2 : int a=0;
1255 2 : int b=0;
1256 : //in the original a and b are "default=0",
1257 : //but you divide by a and b during the initialization!
1258 : //better an error message than an UB, so no default
1259 4 : REQUIREDPARSE(data,"A",a,errormsg);
1260 4 : REQUIREDPARSE(data,"B",b,errormsg);
1261 2 : std::tie(type,switchData) = smapSwitch::init(d0,dmax,r0,a,b);
1262 12 : } else if(name=="Q") {
1263 2 : double beta = 50.0; // nm-1
1264 2 : double lambda = 1.8; // unitless
1265 : double ref;
1266 6 : CHECKandPARSE(data,"BETA",beta,errormsg);
1267 6 : CHECKandPARSE(data,"LAMBDA",lambda,errormsg);
1268 4 : REQUIREDPARSE(data,"REF",ref,errormsg);
1269 : //the original error message was not standard
1270 : // if(!Tools::parse(data,"REF",ref))
1271 : // errormsg="REF (reference distaance) is required for native Q";
1272 2 : std::tie(type,switchData) = nativeqSwitch::init(d0,dmax,r0,beta,lambda,ref);
1273 10 : } else if(name=="EXP") {
1274 2 : std::tie(type,switchData) = exponentialSwitch::init(d0,dmax,r0);
1275 8 : } else if(name=="GAUSSIAN") {
1276 4 : if ( r0==1.0 && d0==0.0 ) {
1277 2 : std::tie(type,switchData) = fastgaussianSwitch::init(dmax);
1278 : } else {
1279 2 : std::tie(type,switchData) = gaussianSwitch::init(d0,dmax,r0);
1280 : }
1281 4 : } else if(name=="TANH") {
1282 2 : std::tie(type,switchData) = tanhSwitch::init(d0,dmax,r0);
1283 2 : } else if(name=="COSINUS") {
1284 2 : std::tie(type,switchData) = cosinusSwitch::init(d0,dmax,r0);
1285 0 : } else if((name=="MATHEVAL" || name=="CUSTOM")) {
1286 : std::string func;
1287 0 : Tools::parse(data,"FUNC",func);
1288 : errormsg="the custom switching function is not comatible with acceleration";
1289 : } else {
1290 0 : errormsg="cannot understand switching function type '"+name+"'";
1291 : }
1292 : }
1293 : #undef CHECKandPARSE
1294 : #undef REQUIREDPARSE
1295 :
1296 24 : if( !data.empty() ) {
1297 : errormsg="found the following rogue keywords in switching function input : ";
1298 0 : for(unsigned i=0; i<data.size(); ++i) {
1299 0 : errormsg = errormsg + data[i] + " ";
1300 : }
1301 : }
1302 :
1303 24 : if(dostretch && dmax!=std::numeric_limits<double>::max()) {
1304 11 : setupStretch(type,switchData);
1305 : }
1306 24 : }
1307 :
1308 0 : std::string SwitchingFunctionAccelerable::description() const {
1309 : // if this error is necessary, something went wrong in the constructor
1310 : // plumed_merror("Unknown switching function type");
1311 0 : return switchContainers::description(type,switchData);
1312 : }
1313 :
1314 240 : double SwitchingFunctionAccelerable::calculateSqr(double distance2,double&dfunc)const {
1315 : double result;
1316 240 : std::tie(result, dfunc) = switchContainers::calculateSqr(type,switchData,distance2);
1317 240 : return result;
1318 : }
1319 :
1320 240 : double SwitchingFunctionAccelerable::calculate(double distance,double&dfunc)const {
1321 : //massert do not go with openacc
1322 : // plumed_massert(init,"you are trying to use an unset SwitchingFunction");
1323 : double result;
1324 240 : std::tie(result, dfunc) = switchContainers::calculate(type,switchData,distance);
1325 240 : return result;
1326 : }
1327 :
1328 0 : void SwitchingFunctionAccelerable::set(const int nn,int mm, const double r0, const double d0) {
1329 0 : if(mm == 0) {
1330 0 : mm = 2*nn;
1331 : }
1332 0 : double dmax=d0+r0*std::pow(0.00001,1./(nn-mm));
1333 0 : std::tie(type,switchData) = switchContainers::initRational(d0,dmax,r0,nn,mm);
1334 0 : setupStretch(type,switchData);
1335 0 : }
1336 :
1337 0 : double SwitchingFunctionAccelerable::get_r0() const {
1338 0 : return 1.0/switchData.invr0;
1339 : }
1340 :
1341 0 : double SwitchingFunctionAccelerable::get_d0() const {
1342 0 : return switchData.d0;
1343 : }
1344 :
1345 0 : double SwitchingFunctionAccelerable::get_dmax() const {
1346 0 : return switchData.dmax;
1347 : }
1348 :
1349 0 : double SwitchingFunctionAccelerable::get_dmax2() const {
1350 0 : return switchData.dmax_2;
1351 : }
1352 :
1353 0 : const switchContainers::Data& SwitchingFunctionAccelerable::get_data() const {
1354 0 : return switchData;
1355 : }
1356 0 : void SwitchingFunctionAccelerable::toACCDevice() const {
1357 : #pragma acc enter data copyin(this[0:1],type)
1358 0 : switchData.toACCDevice();
1359 0 : }
1360 :
1361 0 : void SwitchingFunctionAccelerable::removeFromACCDevice() const {
1362 0 : switchData.removeFromACCDevice();
1363 : #pragma acc exit data delete(type,this[0:1])
1364 0 : }
1365 : }// Namespace PLMD
|
