Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2020 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 <vector>
27 : #include <limits>
28 :
29 : #define PI 3.14159265358979323846
30 :
31 : using namespace std;
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$
115 : s(r) &= 1 & if r<=d0
116 : s(r) &= 0.5 \left( \cos ( \frac{ r - d_0 }{ r_0 } * PI ) + 1 \right) & if d0<r<=d0+r0
117 : s(r) &= 0 & if r> d0+r0
118 : \f$
119 : </td> <td>
120 : {COSINUS R_0=\f$r_0\f$ D_0=\f$d_0\f$}
121 : </td> <td> </td>
122 : </tr> <tr>
123 : <td> CUSTOM </td> <td>
124 : \f$
125 : s(r) = FUNC
126 : \f$
127 : </td> <td>
128 : {CUSTOM FUNC=1/(1+x^6) R_0=\f$r_0\f$ D_0=\f$d_0\f$}
129 : </td> <td> </td>
130 : </tr>
131 : </table>
132 :
133 : Notice that for backward compatibility we allow using `MATHEVAL` instead of `CUSTOM`.
134 : Also notice that if the a `CUSTOM` switching function only depends on even powers of `x` it can be
135 : made faster by using `x2` as a variable. For instance
136 : \verbatim
137 : {CUSTOM FUNC=1/(1+x2^3) R_0=0.3}
138 : \endverbatim
139 : is equivalent to
140 : \verbatim
141 : {CUSTOM FUNC=1/(1+x^6) R_0=0.3}
142 : \endverbatim
143 : but runs faster. The reason is that there is an expensive square root calculation that can be optimized out.
144 :
145 :
146 : \attention
147 : With the default implementation CUSTOM is slower than other functions
148 : (e.g., it is slower than an equivalent RATIONAL function by approximately a factor 2).
149 : Checkout page \ref Lepton to see how to improve its performance.
150 :
151 : For all the switching functions in the above table one can also specify a further (optional) parameter using the parameter
152 : keyword D_MAX to assert that for \f$r>d_{\textrm{max}}\f$ the switching function can be assumed equal to zero.
153 : In this case the function is brought smoothly to zero by stretching and shifting it.
154 : \verbatim
155 : KEYWORD={RATIONAL R_0=1 D_MAX=3}
156 : \endverbatim
157 : the resulting switching function will be
158 : \f$
159 : s(r) = \frac{s'(r)-s'(d_{max})}{s'(0)-s'(d_{max})}
160 : \f$
161 : where
162 : \f$
163 : s'(r)=\frac{1-r^6}{1-r^{12}}
164 : \f$
165 : Since PLUMED 2.2 this is the default. The old behavior (no stretching) can be obtained with the
166 : NOSTRETCH flag. The NOSTRETCH keyword is only provided for backward compatibility and might be
167 : removed in the future. Similarly, the STRETCH keyword is still allowed but has no effect.
168 :
169 : Notice that switching functions defined with the simplified syntax are never stretched
170 : for backward compatibility. This might change in the future.
171 :
172 : */
173 : //+ENDPLUMEDOC
174 :
175 84 : void SwitchingFunction::registerKeywords( Keywords& keys ) {
176 336 : keys.add("compulsory","R_0","the value of R_0 in the switching function");
177 420 : keys.add("compulsory","D_0","0.0","the value of D_0 in the switching function");
178 336 : keys.add("optional","D_MAX","the value at which the switching function can be assumed equal to zero");
179 420 : keys.add("compulsory","NN","6","the value of n in the switching function (only needed for TYPE=RATIONAL)");
180 420 : keys.add("compulsory","MM","0","the value of m in the switching function (only needed for TYPE=RATIONAL); 0 implies 2*NN");
181 336 : keys.add("compulsory","A","the value of a in the switching funciton (only needed for TYPE=SMAP)");
182 336 : keys.add("compulsory","B","the value of b in the switching funciton (only needed for TYPE=SMAP)");
183 84 : }
184 :
185 1006 : void SwitchingFunction::set(const std::string & definition,std::string& errormsg) {
186 2012 : vector<string> data=Tools::getWords(definition);
187 1006 : if( data.size()<1 ) {
188 : errormsg="missing all input for switching function";
189 0 : return;
190 : }
191 : string name=data[0];
192 : data.erase(data.begin());
193 1006 : invr0=0.0;
194 1006 : invr0_2=0.0;
195 1006 : d0=0.0;
196 1006 : dmax=std::numeric_limits<double>::max();
197 1006 : dmax_2=std::numeric_limits<double>::max();
198 1006 : stretch=1.0;
199 1006 : shift=0.0;
200 1006 : init=true;
201 :
202 : bool present;
203 :
204 2012 : present=Tools::findKeyword(data,"D_0");
205 2012 : if(present && !Tools::parse(data,"D_0",d0)) errormsg="could not parse D_0";
206 :
207 2012 : present=Tools::findKeyword(data,"D_MAX");
208 2012 : if(present && !Tools::parse(data,"D_MAX",dmax)) errormsg="could not parse D_MAX";
209 1006 : if(dmax<std::sqrt(std::numeric_limits<double>::max())) dmax_2=dmax*dmax;
210 1006 : bool dostretch=false;
211 2012 : Tools::parseFlag(data,"STRETCH",dostretch); // this is ignored now
212 1006 : dostretch=true;
213 1006 : bool dontstretch=false;
214 2012 : Tools::parseFlag(data,"NOSTRETCH",dontstretch); // this is ignored now
215 1006 : if(dontstretch) dostretch=false;
216 : double r0;
217 1006 : if(name=="CUBIC") {
218 18 : r0 = dmax - d0;
219 : } else {
220 1976 : bool found_r0=Tools::parse(data,"R_0",r0);
221 988 : if(!found_r0) errormsg="R_0 is required";
222 : }
223 1006 : invr0=1.0/r0;
224 1006 : invr0_2=invr0*invr0;
225 :
226 1006 : if(name=="RATIONAL") {
227 281 : type=rational;
228 281 : nn=6;
229 281 : mm=0;
230 562 : present=Tools::findKeyword(data,"NN");
231 562 : if(present && !Tools::parse(data,"NN",nn)) errormsg="could not parse NN";
232 562 : present=Tools::findKeyword(data,"MM");
233 562 : if(present && !Tools::parse(data,"MM",mm)) errormsg="could not parse MM";
234 281 : if(mm==0) mm=2*nn;
235 281 : fastrational=(nn%2==0 && mm%2==0 && d0==0.0);
236 725 : } else if(name=="SMAP") {
237 10 : type=smap;
238 20 : present=Tools::findKeyword(data,"A");
239 20 : if(present && !Tools::parse(data,"A",a)) errormsg="could not parse A";
240 20 : present=Tools::findKeyword(data,"B");
241 20 : if(present && !Tools::parse(data,"B",b)) errormsg="could not parse B";
242 10 : c=pow(2., static_cast<double>(a)/static_cast<double>(b) ) - 1;
243 10 : d = -static_cast<double>(b) / static_cast<double>(a);
244 : }
245 715 : else if(name=="Q") {
246 570 : type=nativeq;
247 570 : beta = 50.0; // nm-1
248 570 : lambda = 1.8; // unitless
249 1140 : present=Tools::findKeyword(data,"BETA");
250 1140 : if(present && !Tools::parse(data, "BETA", beta)) errormsg="could not parse BETA";
251 1140 : present=Tools::findKeyword(data,"LAMBDA");
252 1140 : if(present && !Tools::parse(data, "LAMBDA", lambda)) errormsg="could not parse LAMBDA";
253 1140 : bool found_ref=Tools::parse(data,"REF",ref); // nm
254 570 : if(!found_ref) errormsg="REF (reference disatance) is required for native Q";
255 :
256 : }
257 145 : else if(name=="EXP") type=exponential;
258 80 : else if(name=="GAUSSIAN") type=gaussian;
259 32 : else if(name=="CUBIC") type=cubic;
260 14 : else if(name=="TANH") type=tanh;
261 12 : else if(name=="COSINUS") type=cosinus;
262 19 : else if((name=="MATHEVAL" || name=="CUSTOM")) {
263 10 : type=leptontype;
264 : std::string func;
265 30 : Tools::parse(data,"FUNC",func);
266 20 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
267 10 : lepton_func=func;
268 10 : expression.resize(OpenMP::getNumThreads());
269 30 : for(auto & e : expression) e=pe.createCompiledExpression();
270 20 : lepton_ref.resize(expression.size());
271 80 : for(unsigned t=0; t<lepton_ref.size(); t++) {
272 : try {
273 60 : lepton_ref[t]=&const_cast<lepton::CompiledExpression*>(&expression[t])->getVariableReference("x");
274 12 : } catch(const PLMD::lepton::Exception& exc) {
275 : try {
276 18 : lepton_ref[t]=&const_cast<lepton::CompiledExpression*>(&expression[t])->getVariableReference("x2");
277 6 : leptonx2=true;
278 0 : } catch(const PLMD::lepton::Exception& exc) {
279 : // this is necessary since in some cases lepton things a variable is not present even though it is present
280 : // e.g. func=0*x
281 0 : lepton_ref[t]=nullptr;
282 : }
283 : }
284 : }
285 10 : std::string arg="x";
286 10 : if(leptonx2) arg="x2";
287 30 : lepton::ParsedExpression ped=lepton::Parser::parse(func).differentiate(arg).optimize(lepton::Constants());
288 10 : expression_deriv.resize(OpenMP::getNumThreads());
289 30 : for(auto & e : expression_deriv) e=ped.createCompiledExpression();
290 20 : lepton_ref_deriv.resize(expression_deriv.size());
291 80 : for(unsigned t=0; t<lepton_ref_deriv.size(); t++) {
292 : try {
293 20 : lepton_ref_deriv[t]=&const_cast<lepton::CompiledExpression*>(&expression_deriv[t])->getVariableReference(arg);
294 0 : } catch(const PLMD::lepton::Exception& exc) {
295 : // this is necessary since in some cases lepton things a variable is not present even though it is present
296 : // e.g. func=3*x
297 0 : lepton_ref_deriv[t]=nullptr;
298 : }
299 : }
300 :
301 : }
302 4 : else errormsg="cannot understand switching function type '"+name+"'";
303 1006 : if( !data.empty() ) {
304 : errormsg="found the following rogue keywords in switching function input : ";
305 0 : for(unsigned i=0; i<data.size(); ++i) errormsg = errormsg + data[i] + " ";
306 : }
307 :
308 1006 : if(dostretch && dmax!=std::numeric_limits<double>::max()) {
309 : double dummy;
310 93 : double s0=calculate(0.0,dummy);
311 93 : double sd=calculate(dmax,dummy);
312 93 : stretch=1.0/(s0-sd);
313 93 : shift=-sd*stretch;
314 : }
315 1006 : plumed_assert(!(leptonx2 && d0!=0.0)) << "You cannot use lepton x2 optimization with d0!=0.0 (d0=" << d0 <<")\n"
316 0 : << "Please rewrite your function using x as a variable";
317 : }
318 :
319 1058 : std::string SwitchingFunction::description() const {
320 2116 : std::ostringstream ostr;
321 2116 : ostr<<1./invr0<<". Using ";
322 1058 : if(type==rational) {
323 338 : ostr<<"rational";
324 720 : } else if(type==exponential) {
325 61 : ostr<<"exponential";
326 659 : } else if(type==nativeq) {
327 570 : ostr<<"nativeq";
328 89 : } else if(type==gaussian) {
329 48 : ostr<<"gaussian";
330 41 : } else if(type==smap) {
331 10 : ostr<<"smap";
332 31 : } else if(type==cubic) {
333 18 : ostr<<"cubic";
334 13 : } else if(type==tanh) {
335 2 : ostr<<"tanh";
336 11 : } else if(type==cosinus) {
337 1 : ostr<<"cosinus";
338 10 : } else if(type==leptontype) {
339 10 : ostr<<"lepton";
340 : } else {
341 0 : plumed_merror("Unknown switching function type");
342 : }
343 1058 : ostr<<" switching function with parameters d0="<<d0;
344 1058 : if(type==rational) {
345 676 : ostr<<" nn="<<nn<<" mm="<<mm;
346 720 : } else if(type==nativeq) {
347 1710 : ostr<<" beta="<<beta<<" lambda="<<lambda<<" ref="<<ref;
348 150 : } else if(type==smap) {
349 20 : ostr<<" a="<<a<<" b="<<b;
350 140 : } else if(type==cubic) {
351 18 : ostr<<" dmax="<<dmax;
352 122 : } else if(type==leptontype) {
353 : ostr<<" func="<<lepton_func;
354 :
355 : }
356 1058 : return ostr.str();
357 : }
358 :
359 41284833 : double SwitchingFunction::do_rational(double rdist,double&dfunc,int nn,int mm)const {
360 : double result;
361 41284833 : if(2*nn==mm) {
362 : // if 2*N==M, then (1.0-rdist^N)/(1.0-rdist^M) = 1.0/(1.0+rdist^N)
363 25189547 : double rNdist=Tools::fastpow(rdist,nn-1);
364 25189547 : double iden=1.0/(1+rNdist*rdist);
365 25189547 : dfunc = -nn*rNdist*iden*iden;
366 : result = iden;
367 : } else {
368 16095286 : if(rdist>(1.-100.0*epsilon) && rdist<(1+100.0*epsilon)) {
369 0 : result=nn/mm;
370 0 : dfunc=0.5*nn*(nn-mm)/mm;
371 : } else {
372 16095286 : double rNdist=Tools::fastpow(rdist,nn-1);
373 16095286 : double rMdist=Tools::fastpow(rdist,mm-1);
374 16095286 : double num = 1.-rNdist*rdist;
375 16095286 : double iden = 1./(1.-rMdist*rdist);
376 16095286 : double func = num*iden;
377 : result = func;
378 16095286 : dfunc = ((-nn*rNdist*iden)+(func*(iden*mm)*rMdist));
379 : }
380 : }
381 41284833 : return result;
382 : }
383 :
384 19709549 : double SwitchingFunction::calculateSqr(double distance2,double&dfunc)const {
385 19709549 : if(fastrational) {
386 7657959 : if(distance2>dmax_2) {
387 144482 : dfunc=0.0;
388 144482 : return 0.0;
389 : }
390 7513477 : const double rdist_2 = distance2*invr0_2;
391 7513477 : double result=do_rational(rdist_2,dfunc,nn/2,mm/2);
392 : // chain rule:
393 7513477 : dfunc*=2*invr0_2;
394 : // stretch:
395 7513477 : result=result*stretch+shift;
396 7513477 : dfunc*=stretch;
397 7513477 : return result;
398 12051590 : } else if(leptonx2) {
399 1248110 : if(distance2>dmax_2) {
400 8 : dfunc=0.0;
401 8 : return 0.0;
402 : }
403 1248102 : const unsigned t=OpenMP::getThreadNum();
404 1248102 : const double rdist_2 = distance2*invr0_2;
405 2496204 : plumed_assert(t<expression.size());
406 1248102 : if(lepton_ref[t]) *lepton_ref[t]=rdist_2;
407 1248102 : if(lepton_ref_deriv[t]) *lepton_ref_deriv[t]=rdist_2;
408 1248102 : double result=expression[t].evaluate();
409 1248102 : dfunc=expression_deriv[t].evaluate();
410 : // chain rule:
411 1248102 : dfunc*=2*invr0_2;
412 : // stretch:
413 1248102 : result=result*stretch+shift;
414 1248102 : dfunc*=stretch;
415 1248102 : return result;
416 : } else {
417 10803480 : double distance=std::sqrt(distance2);
418 10803480 : return calculate(distance,dfunc);
419 : }
420 : }
421 :
422 89759393 : double SwitchingFunction::calculate(double distance,double&dfunc)const {
423 89759393 : plumed_massert(init,"you are trying to use an unset SwitchingFunction");
424 89759393 : if(distance>dmax) {
425 450533 : dfunc=0.0;
426 450533 : return 0.0;
427 : }
428 : // in this case, the lepton object stores only the calculateSqr function
429 : // so we have to implement calculate in terms of calculateSqr
430 89308860 : if(leptonx2) {
431 2 : return calculateSqr(distance*distance,dfunc);
432 : }
433 89308858 : const double rdist = (distance-d0)*invr0;
434 : double result;
435 :
436 89308858 : if(rdist<=0.) {
437 : result=1.;
438 24420664 : dfunc=0.0;
439 : } else {
440 64888194 : if(type==smap) {
441 45232932 : double sx=c*Tools::fastpow( rdist, a );
442 22616466 : result=pow( 1.0 + sx, d );
443 22616466 : dfunc=-b*sx/rdist*result/(1.0+sx);
444 42271728 : } else if(type==rational) {
445 33771356 : result=do_rational(rdist,dfunc,nn,mm);
446 8500372 : } else if(type==exponential) {
447 2486473 : result=exp(-rdist);
448 2486473 : dfunc=-result;
449 6013899 : } else if(type==nativeq) {
450 146570 : double rdist2 = beta*(distance - lambda * ref);
451 146570 : double exprdist=exp(rdist2);
452 146570 : double exprmdist=1.0/exprdist;
453 146570 : result=1./(1.+exprdist);
454 146570 : dfunc=-1.0/(exprmdist+1.0)/(1.+exprdist);
455 5867329 : } else if(type==gaussian) {
456 194880 : result=exp(-0.5*rdist*rdist);
457 194880 : dfunc=-rdist*result;
458 5672449 : } else if(type==cubic) {
459 127132 : double tmp1=rdist-1, tmp2=(1+2*rdist);
460 127132 : result=tmp1*tmp1*tmp2;
461 127132 : dfunc=2*tmp1*tmp2 + 2*tmp1*tmp1;
462 5545317 : } else if(type==tanh) {
463 8000 : double tmp1=std::tanh(rdist);
464 8000 : result = 1.0 - tmp1;
465 8000 : dfunc=-(1-tmp1*tmp1);
466 5537317 : } else if(type==cosinus) {
467 : if(rdist<=0.0) {
468 : // rdist = (r-r1)/(r2-r1) ; rdist<=0.0 if r <=r1
469 : result=1.;
470 : dfunc=0.0;
471 522053 : } else if(rdist<=1.0) {
472 : // rdist = (r-r1)/(r2-r1) ; 0.0<=rdist<=1.0 if r1 <= r <=r2; (r2-r1)/(r2-r1)=1
473 226962 : double tmpcos = cos ( rdist * PI );
474 226962 : double tmpsin = sin ( rdist * PI );
475 226962 : result = 0.5 * (tmpcos + 1.0);
476 226962 : dfunc=-0.5 * PI * tmpsin * invr0;
477 : } else {
478 : result=0.;
479 295091 : dfunc=0.0;
480 : }
481 5015264 : } else if(type==leptontype) {
482 5015264 : const unsigned t=OpenMP::getThreadNum();
483 10030528 : plumed_assert(t<expression.size());
484 5015264 : if(lepton_ref[t]) *lepton_ref[t]=rdist;
485 5015264 : if(lepton_ref_deriv[t]) *lepton_ref_deriv[t]=rdist;
486 5015264 : result=expression[t].evaluate();
487 5015264 : dfunc=expression_deriv[t].evaluate();
488 0 : } else plumed_merror("Unknown switching function type");
489 : // this is for the chain rule:
490 64888194 : dfunc*=invr0;
491 : // this is because calculate() sets dfunc to the derivative divided times the distance.
492 : // (I think this is misleading and I would like to modify it - GB)
493 64888194 : dfunc/=distance;
494 : }
495 :
496 89308858 : result=result*stretch+shift;
497 89308858 : dfunc*=stretch;
498 :
499 89308858 : return result;
500 : }
501 :
502 61 : void SwitchingFunction::set(int nn,int mm,double r0,double d0) {
503 61 : init=true;
504 61 : type=rational;
505 61 : if(mm==0) mm=2*nn;
506 61 : this->nn=nn;
507 61 : this->mm=mm;
508 61 : this->invr0=1.0/r0;
509 61 : this->invr0_2=this->invr0*this->invr0;
510 61 : this->d0=d0;
511 61 : this->dmax=d0+r0*pow(0.00001,1./(nn-mm));
512 61 : this->dmax_2=this->dmax*this->dmax;
513 61 : this->leptonx2=false;
514 61 : this->fastrational=(nn%2==0 && mm%2==0 && d0==0.0);
515 :
516 : double dummy;
517 61 : double s0=calculate(0.0,dummy);
518 61 : double sd=calculate(dmax,dummy);
519 61 : stretch=1.0/(s0-sd);
520 61 : shift=-sd*stretch;
521 61 : }
522 :
523 30 : double SwitchingFunction::get_r0() const {
524 30 : return 1./invr0;
525 : }
526 :
527 6 : double SwitchingFunction::get_d0() const {
528 6 : return d0;
529 : }
530 :
531 117918073 : double SwitchingFunction::get_dmax() const {
532 117918073 : return dmax;
533 : }
534 :
535 23893569 : double SwitchingFunction::get_dmax2() const {
536 23893569 : return dmax_2;
537 : }
538 :
539 : }
540 :
541 :
542 :
|
