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