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 109 : void SwitchingFunction::registerKeywords( Keywords& keys ) {
177 218 : keys.add("compulsory","R_0","the value of R_0 in the switching function");
178 218 : keys.add("compulsory","D_0","0.0","the value of D_0 in the switching function");
179 218 : keys.add("optional","D_MAX","the value at which the switching function can be assumed equal to zero");
180 218 : keys.add("compulsory","NN","6","the value of n in the switching function (only needed for TYPE=RATIONAL)");
181 218 : keys.add("compulsory","MM","0","the value of m in the switching function (only needed for TYPE=RATIONAL); 0 implies 2*NN");
182 218 : keys.add("compulsory","A","the value of a in the switching function (only needed for TYPE=SMAP)");
183 218 : keys.add("compulsory","B","the value of b in the switching function (only needed for TYPE=SMAP)");
184 109 : }
185 :
186 1377 : void SwitchingFunction::set(const std::string & definition,std::string& errormsg) {
187 1377 : std::vector<std::string> data=Tools::getWords(definition);
188 1377 : if( data.size()<1 ) {
189 : errormsg="missing all input for switching function";
190 : return;
191 : }
192 1377 : std::string name=data[0];
193 : data.erase(data.begin());
194 1377 : invr0=0.0;
195 1377 : invr0_2=0.0;
196 1377 : d0=0.0;
197 1377 : dmax=std::numeric_limits<double>::max();
198 1377 : dmax_2=std::numeric_limits<double>::max();
199 1377 : stretch=1.0;
200 1377 : shift=0.0;
201 1377 : init=true;
202 :
203 : bool present;
204 :
205 1377 : present=Tools::findKeyword(data,"D_0");
206 1711 : if(present && !Tools::parse(data,"D_0",d0)) {
207 : errormsg="could not parse D_0";
208 : }
209 :
210 1377 : present=Tools::findKeyword(data,"D_MAX");
211 1697 : if(present && !Tools::parse(data,"D_MAX",dmax)) {
212 : errormsg="could not parse D_MAX";
213 : }
214 1377 : if(dmax<std::sqrt(std::numeric_limits<double>::max())) {
215 160 : dmax_2=dmax*dmax;
216 : }
217 1377 : bool dostretch=false;
218 1377 : Tools::parseFlag(data,"STRETCH",dostretch); // this is ignored now
219 1377 : dostretch=true;
220 1377 : bool dontstretch=false;
221 1377 : Tools::parseFlag(data,"NOSTRETCH",dontstretch); // this is ignored now
222 1377 : if(dontstretch) {
223 84 : dostretch=false;
224 : }
225 : double r0;
226 1377 : if(name=="CUBIC") {
227 20 : r0 = dmax - d0;
228 : } else {
229 1357 : bool found_r0=Tools::parse(data,"R_0",r0);
230 1357 : if(!found_r0) {
231 : errormsg="R_0 is required";
232 : }
233 : }
234 1377 : invr0=1.0/r0;
235 1377 : invr0_2=invr0*invr0;
236 :
237 1377 : if(name=="RATIONAL") {
238 317 : type=rational;
239 317 : nn=6;
240 317 : mm=0;
241 317 : present=Tools::findKeyword(data,"NN");
242 515 : if(present && !Tools::parse(data,"NN",nn)) {
243 : errormsg="could not parse NN";
244 : }
245 317 : present=Tools::findKeyword(data,"MM");
246 511 : if(present && !Tools::parse(data,"MM",mm)) {
247 : errormsg="could not parse MM";
248 : }
249 317 : if(mm==0) {
250 220 : mm=2*nn;
251 : }
252 387 : fastrational=(nn%2==0 && mm%2==0 && d0==0.0);
253 1060 : } else if(name=="SMAP") {
254 12 : type=smap;
255 12 : present=Tools::findKeyword(data,"A");
256 36 : if(present && !Tools::parse(data,"A",a)) {
257 : errormsg="could not parse A";
258 : }
259 12 : present=Tools::findKeyword(data,"B");
260 36 : if(present && !Tools::parse(data,"B",b)) {
261 : errormsg="could not parse B";
262 : }
263 12 : c=std::pow(2., static_cast<double>(a)/static_cast<double>(b) ) - 1;
264 12 : d = -static_cast<double>(b) / static_cast<double>(a);
265 1048 : } else if(name=="Q") {
266 864 : type=nativeq;
267 864 : beta = 50.0; // nm-1
268 864 : lambda = 1.8; // unitless
269 864 : present=Tools::findKeyword(data,"BETA");
270 2592 : if(present && !Tools::parse(data, "BETA", beta)) {
271 : errormsg="could not parse BETA";
272 : }
273 864 : present=Tools::findKeyword(data,"LAMBDA");
274 2592 : if(present && !Tools::parse(data, "LAMBDA", lambda)) {
275 : errormsg="could not parse LAMBDA";
276 : }
277 864 : bool found_ref=Tools::parse(data,"REF",ref); // nm
278 864 : if(!found_ref) {
279 : errormsg="REF (reference disatance) is required for native Q";
280 : }
281 :
282 184 : } else if(name=="EXP") {
283 67 : type=exponential;
284 117 : } else if(name=="GAUSSIAN") {
285 52 : type=gaussian;
286 65 : } else if(name=="CUBIC") {
287 20 : type=cubic;
288 45 : } else if(name=="TANH") {
289 4 : type=tanh;
290 41 : } else if(name=="COSINUS") {
291 3 : type=cosinus;
292 73 : } else if((name=="MATHEVAL" || name=="CUSTOM")) {
293 37 : type=leptontype;
294 : std::string func;
295 37 : Tools::parse(data,"FUNC",func);
296 39 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(lepton::Constants());
297 37 : lepton_func=func;
298 37 : expression.resize(OpenMP::getNumThreads());
299 111 : for(auto & e : expression) {
300 74 : e=pe.createCompiledExpression();
301 : }
302 37 : lepton_ref.resize(expression.size());
303 107 : for(unsigned t=0; t<lepton_ref.size(); t++) {
304 72 : auto vars=expression[t].getVariables();
305 72 : bool found_x=std::find(vars.begin(),vars.end(),"x")!=vars.end();
306 72 : bool found_x2=std::find(vars.begin(),vars.end(),"x2")!=vars.end();
307 72 : if (vars.size()==0) {
308 : // this is necessary since in some cases lepton thinks a variable is not present even though it is present
309 : // e.g. func=0*x
310 0 : lepton_ref[t]=nullptr;
311 72 : } else if(vars.size()==1 && found_x) {
312 64 : lepton_ref[t]=&const_cast<lepton::CompiledExpression*>(&expression[t])->getVariableReference("x");
313 8 : } else if(vars.size()==1 && found_x2) {
314 8 : lepton_ref[t]=&const_cast<lepton::CompiledExpression*>(&expression[t])->getVariableReference("x2");
315 6 : leptonx2=true;
316 2 : } else if(vars.size()==2 && found_x && found_x2) {
317 2 : plumed_error() << "Cannot use simultaneously x and x2 argument in switching function: "<<func;
318 : } else {
319 2 : plumed_error() << "Something wrong in the arguments for switching function: "<<func;
320 : }
321 : }
322 37 : std::string arg="x";
323 35 : if(leptonx2) {
324 : arg="x2";
325 : }
326 70 : lepton::ParsedExpression ped=lepton::Parser::parse(func).differentiate(arg).optimize(lepton::Constants());
327 35 : expression_deriv.resize(OpenMP::getNumThreads());
328 105 : for(auto & e : expression_deriv) {
329 70 : e=ped.createCompiledExpression();
330 : }
331 35 : lepton_ref_deriv.resize(expression_deriv.size());
332 105 : for(unsigned t=0; t<lepton_ref_deriv.size(); t++) {
333 : try {
334 70 : lepton_ref_deriv[t]=&const_cast<lepton::CompiledExpression*>(&expression_deriv[t])->getVariableReference(arg);
335 0 : } catch(const PLMD::lepton::Exception& exc) {
336 : // this is necessary since in some cases lepton things a variable is not present even though it is present
337 : // e.g. func=3*x
338 0 : lepton_ref_deriv[t]=nullptr;
339 0 : }
340 : }
341 :
342 : } else {
343 2 : errormsg="cannot understand switching function type '"+name+"'";
344 : }
345 1375 : if( !data.empty() ) {
346 : errormsg="found the following rogue keywords in switching function input : ";
347 0 : for(unsigned i=0; i<data.size(); ++i) {
348 2 : errormsg = errormsg + data[i] + " ";
349 : }
350 : }
351 :
352 1375 : if(dostretch && dmax!=std::numeric_limits<double>::max()) {
353 : double dummy;
354 118 : double s0=calculate(0.0,dummy);
355 118 : double sd=calculate(dmax,dummy);
356 118 : stretch=1.0/(s0-sd);
357 118 : shift=-sd*stretch;
358 : }
359 1375 : plumed_assert(!(leptonx2 && d0!=0.0)) << "You cannot use lepton x2 optimization with d0!=0.0 (d0=" << d0 <<")\n"
360 0 : << "Please rewrite your function using x as a variable";
361 1377 : }
362 :
363 1379 : std::string SwitchingFunction::description() const {
364 1379 : std::ostringstream ostr;
365 1379 : ostr<<1./invr0<<". Using ";
366 1379 : if(type==rational) {
367 366 : ostr<<"rational";
368 : } else if(type==exponential) {
369 61 : ostr<<"exponential";
370 : } else if(type==nativeq) {
371 862 : ostr<<"nativeq";
372 : } else if(type==gaussian) {
373 48 : ostr<<"gaussian";
374 : } else if(type==smap) {
375 10 : ostr<<"smap";
376 : } else if(type==cubic) {
377 18 : ostr<<"cubic";
378 : } else if(type==tanh) {
379 2 : ostr<<"tanh";
380 : } else if(type==cosinus) {
381 1 : ostr<<"cosinus";
382 : } else if(type==leptontype) {
383 11 : ostr<<"lepton";
384 : } else {
385 0 : plumed_merror("Unknown switching function type");
386 : }
387 1379 : ostr<<" switching function with parameters d0="<<d0;
388 1379 : if(type==rational) {
389 366 : ostr<<" nn="<<nn<<" mm="<<mm;
390 1013 : } else if(type==nativeq) {
391 862 : ostr<<" beta="<<beta<<" lambda="<<lambda<<" ref="<<ref;
392 151 : } else if(type==smap) {
393 10 : ostr<<" a="<<a<<" b="<<b;
394 141 : } else if(type==cubic) {
395 18 : ostr<<" dmax="<<dmax;
396 123 : } else if(type==leptontype) {
397 11 : ostr<<" func="<<lepton_func;
398 :
399 : }
400 1379 : return ostr.str();
401 1379 : }
402 :
403 41337464 : double SwitchingFunction::do_rational(double rdist,double&dfunc,int nn,int mm)const {
404 : double result;
405 41337464 : if(2*nn==mm) {
406 : // if 2*N==M, then (1.0-rdist^N)/(1.0-rdist^M) = 1.0/(1.0+rdist^N)
407 25193522 : double rNdist=Tools::fastpow(rdist,nn-1);
408 25193522 : double iden=1.0/(1+rNdist*rdist);
409 25193522 : dfunc = -nn*rNdist*iden*iden;
410 : result = iden;
411 : } else {
412 16143942 : if(rdist>(1.-5.0e10*epsilon) && rdist<(1+5.0e10*epsilon)) {
413 10 : const double secDev = ((nn * (mm * mm - 3.0* mm * (-1 + nn ) + nn *(-3 + 2* nn )))/(6.0* mm ));
414 10 : const double x =(rdist-1.0);
415 10 : dfunc=0.5*nn*double(nn-mm)/mm;
416 10 : result = double(nn)/mm+ x * ( dfunc + 0.5 * x * secDev);
417 10 : dfunc = dfunc + x * secDev;
418 10 : } else {
419 16143932 : double rNdist=Tools::fastpow(rdist,nn-1);
420 16143932 : double rMdist=Tools::fastpow(rdist,mm-1);
421 16143932 : double num = 1.-rNdist*rdist;
422 16143932 : double iden = 1./(1.-rMdist*rdist);
423 16143932 : double func = num*iden;
424 : result = func;
425 16143932 : dfunc = ((-nn*rNdist*iden)+(func*(iden*mm)*rMdist));
426 : }
427 : }
428 41337464 : return result;
429 : }
430 :
431 19723321 : double SwitchingFunction::calculateSqr(double distance2,double&dfunc)const {
432 19723321 : if(fastrational) {
433 7658019 : if(distance2>dmax_2) {
434 144494 : dfunc=0.0;
435 144494 : return 0.0;
436 : }
437 7513525 : const double rdist_2 = distance2*invr0_2;
438 7513525 : double result=do_rational(rdist_2,dfunc,nn/2,mm/2);
439 : // chain rule:
440 7513525 : dfunc*=2*invr0_2;
441 : // stretch:
442 7513525 : result=result*stretch+shift;
443 7513525 : dfunc*=stretch;
444 7513525 : return result;
445 12065302 : } else if(leptonx2) {
446 1248110 : if(distance2>dmax_2) {
447 8 : dfunc=0.0;
448 8 : return 0.0;
449 : }
450 1248102 : const unsigned t=OpenMP::getThreadNum();
451 1248102 : const double rdist_2 = distance2*invr0_2;
452 1248102 : plumed_assert(t<expression.size());
453 1248102 : if(lepton_ref[t]) {
454 1248102 : *lepton_ref[t]=rdist_2;
455 : }
456 1248102 : if(lepton_ref_deriv[t]) {
457 1248102 : *lepton_ref_deriv[t]=rdist_2;
458 : }
459 1248102 : double result=expression[t].evaluate();
460 1248102 : dfunc=expression_deriv[t].evaluate();
461 : // chain rule:
462 1248102 : dfunc*=2*invr0_2;
463 : // stretch:
464 1248102 : result=result*stretch+shift;
465 1248102 : dfunc*=stretch;
466 1248102 : return result;
467 : } else {
468 10817192 : double distance=std::sqrt(distance2);
469 10817192 : return calculate(distance,dfunc);
470 : }
471 : }
472 :
473 89172907 : double SwitchingFunction::calculate(double distance,double&dfunc)const {
474 89172907 : plumed_massert(init,"you are trying to use an unset SwitchingFunction");
475 89172907 : if(distance>dmax) {
476 450805 : dfunc=0.0;
477 450805 : return 0.0;
478 : }
479 : // in this case, the lepton object stores only the calculateSqr function
480 : // so we have to implement calculate in terms of calculateSqr
481 88722102 : if(leptonx2) {
482 2 : return calculateSqr(distance*distance,dfunc);
483 : }
484 88722100 : const double rdist = (distance-d0)*invr0;
485 : double result;
486 :
487 88722100 : if(rdist<=0.) {
488 : result=1.;
489 24460073 : dfunc=0.0;
490 : } else {
491 64262027 : if(type==smap) {
492 21789996 : double sx=c*Tools::fastpow( rdist, a );
493 21789996 : result=std::pow( 1.0 + sx, d );
494 21789996 : dfunc=-b*sx/rdist*result/(1.0+sx);
495 : } else if(type==rational) {
496 33823939 : result=do_rational(rdist,dfunc,nn,mm);
497 : } else if(type==exponential) {
498 2486499 : result=std::exp(-rdist);
499 2486499 : dfunc=-result;
500 : } else if(type==nativeq) {
501 292887 : double rdist2 = beta*(distance - lambda * ref);
502 292887 : double exprdist=std::exp(rdist2);
503 292887 : double exprmdist=1.0/exprdist;
504 292887 : result=1./(1.+exprdist);
505 292887 : dfunc=-beta/(exprmdist+1.0)/(1.+exprdist)/invr0;
506 : } else if(type==gaussian) {
507 195733 : result=std::exp(-0.5*rdist*rdist);
508 195733 : dfunc=-rdist*result;
509 : } else if(type==cubic) {
510 127153 : double tmp1=rdist-1, tmp2=(1+2*rdist);
511 127153 : result=tmp1*tmp1*tmp2;
512 127153 : dfunc=2*tmp1*tmp2 + 2*tmp1*tmp1;
513 : } else if(type==tanh) {
514 8025 : double tmp1=std::tanh(rdist);
515 8025 : result = 1.0 - tmp1;
516 8025 : dfunc=-(1-tmp1*tmp1);
517 : } else if(type==cosinus) {
518 : if(rdist<=0.0) {
519 : // rdist = (r-r1)/(r2-r1) ; rdist<=0.0 if r <=r1
520 : result=1.;
521 : dfunc=0.0;
522 522089 : } else if(rdist<=1.0) {
523 : // rdist = (r-r1)/(r2-r1) ; 0.0<=rdist<=1.0 if r1 <= r <=r2; (r2-r1)/(r2-r1)=1
524 226986 : double tmpcos = std::cos ( rdist * PI );
525 226986 : double tmpsin = std::sin ( rdist * PI );
526 226986 : result = 0.5 * (tmpcos + 1.0);
527 226986 : dfunc=-0.5 * PI * tmpsin;
528 : } else {
529 : result=0.;
530 295103 : dfunc=0.0;
531 : }
532 : } else if(type==leptontype) {
533 5015706 : const unsigned t=OpenMP::getThreadNum();
534 5015706 : plumed_assert(t<expression.size());
535 5015706 : if(lepton_ref[t]) {
536 5015706 : *lepton_ref[t]=rdist;
537 : }
538 5015706 : if(lepton_ref_deriv[t]) {
539 5015706 : *lepton_ref_deriv[t]=rdist;
540 : }
541 5015706 : result=expression[t].evaluate();
542 5015706 : dfunc=expression_deriv[t].evaluate();
543 : } else {
544 0 : plumed_merror("Unknown switching function type");
545 : }
546 : // this is for the chain rule (derivative of rdist):
547 64262027 : dfunc*=invr0;
548 : // for any future switching functions, be aware that multiplying invr0 is only correct for functions of rdist = (r-d0)/r0.
549 :
550 : // this is because calculate() sets dfunc to the derivative divided times the distance.
551 : // (I think this is misleading and I would like to modify it - GB)
552 64262027 : dfunc/=distance;
553 : }
554 :
555 88722100 : result=result*stretch+shift;
556 88722100 : dfunc*=stretch;
557 :
558 88722100 : return result;
559 : }
560 :
561 61 : void SwitchingFunction::set(int nn,int mm,double r0,double d0) {
562 61 : init=true;
563 61 : type=rational;
564 61 : if(mm==0) {
565 61 : mm=2*nn;
566 : }
567 61 : this->nn=nn;
568 61 : this->mm=mm;
569 61 : this->invr0=1.0/r0;
570 61 : this->invr0_2=this->invr0*this->invr0;
571 61 : this->d0=d0;
572 61 : this->dmax=d0+r0*std::pow(0.00001,1./(nn-mm));
573 61 : this->dmax_2=this->dmax*this->dmax;
574 61 : this->leptonx2=false;
575 61 : this->fastrational=(nn%2==0 && mm%2==0 && d0==0.0);
576 :
577 : double dummy;
578 61 : double s0=calculate(0.0,dummy);
579 61 : double sd=calculate(dmax,dummy);
580 61 : stretch=1.0/(s0-sd);
581 61 : shift=-sd*stretch;
582 61 : }
583 :
584 30 : double SwitchingFunction::get_r0() const {
585 30 : return 1./invr0;
586 : }
587 :
588 6 : double SwitchingFunction::get_d0() const {
589 6 : return d0;
590 : }
591 :
592 117918073 : double SwitchingFunction::get_dmax() const {
593 117918073 : return dmax;
594 : }
595 :
596 23893569 : double SwitchingFunction::get_dmax2() const {
597 23893569 : return dmax_2;
598 : }
599 :
600 : }
601 :
602 :
603 :
|
