Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2019 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 : #ifdef __PLUMED_HAS_MATHEVAL
30 : #include <matheval.h>
31 : #endif
32 :
33 : using namespace std;
34 : namespace PLMD {
35 :
36 3226 : static std::map<string, double> leptonConstants= {
37 : {"e", std::exp(1.0)},
38 : {"log2e", 1.0/std::log(2.0)},
39 : {"log10e", 1.0/std::log(10.0)},
40 : {"ln2", std::log(2.0)},
41 : {"ln10", std::log(10.0)},
42 : {"pi", pi},
43 : {"pi_2", pi*0.5},
44 : {"pi_4", pi*0.25},
45 : // {"1_pi", 1.0/pi},
46 : // {"2_pi", 2.0/pi},
47 : // {"2_sqrtpi", 2.0/std::sqrt(pi)},
48 : {"sqrt2", std::sqrt(2.0)},
49 : {"sqrt1_2", std::sqrt(0.5)}
50 : };
51 :
52 : //+PLUMEDOC INTERNAL switchingfunction
53 : /*
54 : Functions that measure whether values are less than a certain quantity.
55 :
56 : Switching functions \f$s(r)\f$ take a minimum of one input parameter \f$d_0\f$.
57 : 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.
58 : The various switching functions available in plumed differ in terms of how this decay is performed.
59 :
60 : Where there is an accepted convention in the literature (e.g. \ref COORDINATION) on the form of the
61 : switching function we use the convention as the default. However, the flexibility to use different
62 : switching functions is always present generally through a single keyword. This keyword generally
63 : takes an input with the following form:
64 :
65 : \verbatim
66 : KEYWORD={TYPE <list of parameters>}
67 : \endverbatim
68 :
69 : The following table contains a list of the various switching functions that are available in plumed 2
70 : together with an example input.
71 :
72 : <table align=center frame=void width=95%% cellpadding=5%%>
73 : <tr>
74 : <td> TYPE </td> <td> FUNCTION </td> <td> EXAMPLE INPUT </td> <td> DEFAULT PARAMETERS </td>
75 : </tr> <tr> <td>RATIONAL </td> <td>
76 : \f$
77 : s(r)=\frac{ 1 - \left(\frac{ r - d_0 }{ r_0 }\right)^{n} }{ 1 - \left(\frac{ r - d_0 }{ r_0 }\right)^{m} }
78 : \f$
79 : </td> <td>
80 : {RATIONAL R_0=\f$r_0\f$ D_0=\f$d_0\f$ NN=\f$n\f$ MM=\f$m\f$}
81 : </td> <td> \f$d_0=0.0\f$, \f$n=6\f$, \f$m=2n\f$ </td>
82 : </tr> <tr>
83 : <td> EXP </td> <td>
84 : \f$
85 : s(r)=\exp\left(-\frac{ r - d_0 }{ r_0 }\right)
86 : \f$
87 : </td> <td>
88 : {EXP R_0=\f$r_0\f$ D_0=\f$d_0\f$}
89 : </td> <td> \f$d_0=0.0\f$ </td>
90 : </tr> <tr>
91 : <td> GAUSSIAN </td> <td>
92 : \f$
93 : s(r)=\exp\left(-\frac{ (r - d_0)^2 }{ 2r_0^2 }\right)
94 : \f$
95 : </td> <td>
96 : {GAUSSIAN R_0=\f$r_0\f$ D_0=\f$d_0\f$}
97 : </td> <td> \f$d_0=0.0\f$ </td>
98 : </tr> <tr>
99 : <td> SMAP </td> <td>
100 : \f$
101 : s(r) = \left[ 1 + ( 2^{a/b} -1 )\left( \frac{r-d_0}{r_0} \right)^a \right]^{-b/a}
102 : \f$
103 : </td> <td>
104 : {SMAP R_0=\f$r_0\f$ D_0=\f$d_0\f$ A=\f$a\f$ B=\f$b\f$}
105 : </td> <td> \f$d_0=0.0\f$ </td>
106 : </tr> <tr>
107 : <td> Q </td> <td>
108 : \f$
109 : s(r) = \frac{1}{1 + \exp(\beta(r_{ij} - \lambda r_{ij}^0))}
110 : \f$
111 : </td> <td>
112 : {Q REF=\f$r_{ij}^0\f$ BETA=\f$\beta\f$ LAMBDA=\f$\lambda\f$ }
113 : </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>
114 : </tr> <tr>
115 : <td> CUBIC </td> <td>
116 : \f$
117 : s(r) = (y-1)^2(1+2y) \qquad \textrm{where} \quad y = \frac{r - r_1}{r_0-r_1}
118 : \f$
119 : </td> <td>
120 : {CUBIC D_0=\f$r_1\f$ D_MAX=\f$r_0\f$}
121 : </td> <td> </td>
122 : </tr> <tr>
123 : <td> TANH </td> <td>
124 : \f$
125 : s(r) = 1 - \tanh\left( \frac{ r - d_0 }{ r_0 } \right)
126 : \f$
127 : </td> <td>
128 : {TANH R_0=\f$r_0\f$ D_0=\f$d_0\f$}
129 : </td> <td> </td>
130 : </tr> <tr>
131 : <td> MATHEVAL </td> <td>
132 : \f$
133 : s(r) = FUNC
134 : \f$
135 : </td> <td>
136 : {MATHEVAL FUNC=1/(1+x^6) R_0=\f$r_0\f$ D_0=\f$d_0\f$}
137 : </td> <td> </td>
138 : </tr>
139 : </table>
140 :
141 : \attention
142 : Similarly to the \ref MATHEVAL function, the MATHEVAL switching function
143 : only works if one of these two conditions is satisfied:
144 : (a) libmatheval is installed on the system and PLUMED has been linked to it or
145 : (b) the environment variable `PLUMED_USE_LEPTON` is set equal to `yes` at runtime
146 : (in this case, the internal lepton library will be used).
147 : Notice that in version v2.5 the internal lepton library will be used by default.
148 : Also notice that using MATHEVAL is much slower than using e.g. RATIONAL.
149 : Thus, the MATHEVAL switching function is useful to perform quick
150 : tests on switching functions with arbitrary form before proceeding to their
151 : implementation in C++.
152 :
153 : For all the switching functions in the above table one can also specify a further (optional) parameter using the parameter
154 : keyword D_MAX to assert that for \f$r>d_{\textrm{max}}\f$ the switching function can be assumed equal to zero.
155 : In this case the function is brought smoothly to zero by stretching and shifting it.
156 : \verbatim
157 : KEYWORD={RATIONAL R_0=1 D_MAX=3}
158 : \endverbatim
159 : the resulting switching function will be
160 : \f$
161 : s(r) = \frac{s'(r)-s'(d_{max})}{s'(0)-s'(d_{max})}
162 : \f$
163 : where
164 : \f$
165 : s'(r)=\frac{1-r^6}{1-r^{12}}
166 : \f$
167 : Since PLUMED 2.2 this is the default. The old behavior (no stretching) can be obtained with the
168 : NOSTRETCH flag. The NOSTRETCH keyword is only provided for backward compatibility and might be
169 : removed in the future. Similarly, the STRETCH keyword is still allowed but has no effect.
170 :
171 : Notice that switching functions defined with the simplified syntax are never stretched
172 : for backward compatibility. This might change in the future.
173 :
174 : */
175 : //+ENDPLUMEDOC
176 :
177 84 : void SwitchingFunction::registerKeywords( Keywords& keys ) {
178 336 : keys.add("compulsory","R_0","the value of R_0 in the switching function");
179 420 : keys.add("compulsory","D_0","0.0","the value of D_0 in the switching function");
180 336 : keys.add("optional","D_MAX","the value at which the switching function can be assumed equal to zero");
181 420 : keys.add("compulsory","NN","6","the value of n in the switching function (only needed for TYPE=RATIONAL)");
182 420 : keys.add("compulsory","MM","0","the value of m in the switching function (only needed for TYPE=RATIONAL); 0 implies 2*NN");
183 336 : keys.add("compulsory","A","the value of a in the switching funciton (only needed for TYPE=SMAP)");
184 336 : keys.add("compulsory","B","the value of b in the switching funciton (only needed for TYPE=SMAP)");
185 84 : }
186 :
187 659 : void SwitchingFunction::set(const std::string & definition,std::string& errormsg) {
188 1318 : vector<string> data=Tools::getWords(definition);
189 659 : if( data.size()<1 ) {
190 : errormsg="missing all input for switching function";
191 0 : return;
192 : }
193 : string name=data[0];
194 : data.erase(data.begin());
195 659 : invr0=0.0;
196 659 : invr0_2=0.0;
197 659 : d0=0.0;
198 659 : dmax=std::numeric_limits<double>::max();
199 659 : dmax_2=std::numeric_limits<double>::max();
200 659 : stretch=1.0;
201 659 : shift=0.0;
202 659 : init=true;
203 :
204 : bool present;
205 :
206 1318 : present=Tools::findKeyword(data,"D_0");
207 1318 : if(present && !Tools::parse(data,"D_0",d0)) errormsg="could not parse D_0";
208 :
209 1318 : present=Tools::findKeyword(data,"D_MAX");
210 1318 : if(present && !Tools::parse(data,"D_MAX",dmax)) errormsg="could not parse D_MAX";
211 659 : if(dmax<std::sqrt(std::numeric_limits<double>::max())) dmax_2=dmax*dmax;
212 659 : bool dostretch=false;
213 1318 : Tools::parseFlag(data,"STRETCH",dostretch); // this is ignored now
214 659 : dostretch=true;
215 659 : bool dontstretch=false;
216 1318 : Tools::parseFlag(data,"NOSTRETCH",dontstretch); // this is ignored now
217 659 : if(dontstretch) dostretch=false;
218 : double r0;
219 659 : if(name=="CUBIC") {
220 17 : r0 = dmax - d0;
221 : } else {
222 1284 : bool found_r0=Tools::parse(data,"R_0",r0);
223 642 : if(!found_r0) errormsg="R_0 is required";
224 : }
225 659 : invr0=1.0/r0;
226 659 : invr0_2=invr0*invr0;
227 :
228 659 : if(name=="RATIONAL") {
229 252 : type=rational;
230 252 : nn=6;
231 252 : mm=0;
232 504 : present=Tools::findKeyword(data,"NN");
233 504 : if(present && !Tools::parse(data,"NN",nn)) errormsg="could not parse NN";
234 504 : present=Tools::findKeyword(data,"MM");
235 504 : if(present && !Tools::parse(data,"MM",mm)) errormsg="could not parse MM";
236 252 : if(mm==0) mm=2*nn;
237 407 : } else if(name=="SMAP") {
238 0 : type=smap;
239 0 : present=Tools::findKeyword(data,"A");
240 0 : if(present && !Tools::parse(data,"A",a)) errormsg="could not parse A";
241 0 : present=Tools::findKeyword(data,"B");
242 0 : if(present && !Tools::parse(data,"B",b)) errormsg="could not parse B";
243 0 : c=pow(2., static_cast<double>(a)/static_cast<double>(b) ) - 1;
244 0 : d = -static_cast<double>(b) / static_cast<double>(a);
245 : }
246 407 : else if(name=="Q") {
247 278 : type=nativeq;
248 278 : beta = 50.0; // nm-1
249 278 : lambda = 1.8; // unitless
250 556 : present=Tools::findKeyword(data,"BETA");
251 556 : if(present && !Tools::parse(data, "BETA", beta)) errormsg="could not parse BETA";
252 556 : present=Tools::findKeyword(data,"LAMBDA");
253 556 : if(present && !Tools::parse(data, "LAMBDA", lambda)) errormsg="could not parse LAMBDA";
254 556 : bool found_ref=Tools::parse(data,"REF",ref); // nm
255 278 : if(!found_ref) errormsg="REF (reference disatance) is required for native Q";
256 :
257 : }
258 129 : else if(name=="EXP") type=exponential;
259 69 : else if(name=="GAUSSIAN") type=gaussian;
260 25 : else if(name=="CUBIC") type=cubic;
261 8 : else if(name=="TANH") type=tanh;
262 8 : else if((name=="MATHEVAL" || name=="CUSTOM") && std::getenv("PLUMED_USE_LEPTON")) {
263 2 : type=leptontype;
264 : std::string func;
265 6 : Tools::parse(data,"FUNC",func);
266 4 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(leptonConstants);
267 2 : lepton_func=func;
268 2 : expression.resize(OpenMP::getNumThreads());
269 6 : for(auto & e : expression) e=pe.createCompiledExpression();
270 8 : lepton::ParsedExpression ped=lepton::Parser::parse(func).differentiate("x").optimize(leptonConstants);
271 2 : expression_deriv.resize(OpenMP::getNumThreads());
272 6 : for(auto & e : expression_deriv) e=ped.createCompiledExpression();
273 : }
274 : #ifdef __PLUMED_HAS_MATHEVAL
275 5 : else if(name=="MATHEVAL" || name=="CUSTOM") {
276 3 : type=matheval;
277 : std::string func;
278 6 : Tools::parse(data,"FUNC",func);
279 3 : const unsigned nt=OpenMP::getNumThreads();
280 3 : plumed_assert(nt>0);
281 3 : evaluator.resize(nt);
282 15 : for(unsigned i=0; i<nt; i++) {
283 12 : evaluator[i]=evaluator_create(const_cast<char*>(func.c_str()));
284 : }
285 : char **check_names;
286 : int check_count;
287 3 : evaluator_get_variables(evaluator[0],&check_names,&check_count);
288 3 : if(check_count!=1) {
289 : errormsg="wrong number of arguments in MATHEVAL switching function";
290 : return;
291 : }
292 6 : if(std::string(check_names[0])!="x") {
293 : errormsg ="argument should be named 'x'";
294 : return;
295 : }
296 3 : evaluator_deriv.resize(nt);
297 15 : for(unsigned i=0; i<nt; i++) {
298 18 : evaluator_deriv[i]=evaluator_derivative(evaluator[i],const_cast<char*>("x"));
299 : }
300 : }
301 : #endif
302 4 : else errormsg="cannot understand switching function type '"+name+"'";
303 659 : 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 659 : if(dostretch && dmax!=std::numeric_limits<double>::max()) {
309 : double dummy;
310 91 : double s0=calculate(0.0,dummy);
311 91 : double sd=calculate(dmax,dummy);
312 91 : stretch=1.0/(s0-sd);
313 91 : shift=-sd*stretch;
314 : }
315 : }
316 :
317 704 : std::string SwitchingFunction::description() const {
318 1408 : std::ostringstream ostr;
319 1408 : ostr<<1./invr0<<". Using ";
320 704 : if(type==rational) {
321 302 : ostr<<"rational";
322 402 : } else if(type==exponential) {
323 56 : ostr<<"exponential";
324 346 : } else if(type==nativeq) {
325 278 : ostr<<"nativeq";
326 68 : } else if(type==gaussian) {
327 44 : ostr<<"gaussian";
328 24 : } else if(type==smap) {
329 0 : ostr<<"smap";
330 24 : } else if(type==cubic) {
331 17 : ostr<<"cubic";
332 7 : } else if(type==tanh) {
333 2 : ostr<<"tanh";
334 5 : } else if(type==leptontype) {
335 2 : ostr<<"lepton";
336 : #ifdef __PLUMED_HAS_MATHEVAL
337 3 : } else if(type==matheval) {
338 3 : ostr<<"matheval";
339 : #endif
340 : } else {
341 0 : plumed_merror("Unknown switching function type");
342 : }
343 704 : ostr<<" swiching function with parameters d0="<<d0;
344 704 : if(type==rational) {
345 604 : ostr<<" nn="<<nn<<" mm="<<mm;
346 402 : } else if(type==nativeq) {
347 834 : ostr<<" beta="<<beta<<" lambda="<<lambda<<" ref="<<ref;
348 124 : } else if(type==smap) {
349 0 : ostr<<" a="<<a<<" b="<<b;
350 124 : } else if(type==cubic) {
351 17 : ostr<<" dmax="<<dmax;
352 107 : } else if(type==leptontype) {
353 : ostr<<" func="<<lepton_func;
354 : #ifdef __PLUMED_HAS_MATHEVAL
355 105 : } else if(type==matheval) {
356 3 : ostr<<" func="<<evaluator_get_string(evaluator[0]);
357 : #endif
358 :
359 : }
360 704 : return ostr.str();
361 : }
362 :
363 12031290 : double SwitchingFunction::do_rational(double rdist,double&dfunc,int nn,int mm)const {
364 : double result;
365 12031290 : if(2*nn==mm) {
366 : // if 2*N==M, then (1.0-rdist^N)/(1.0-rdist^M) = 1.0/(1.0+rdist^N)
367 11992172 : double rNdist=Tools::fastpow(rdist,nn-1);
368 11992172 : double iden=1.0/(1+rNdist*rdist);
369 11992172 : dfunc = -nn*rNdist*iden*iden;
370 : result = iden;
371 : } else {
372 39118 : if(rdist>(1.-100.0*epsilon) && rdist<(1+100.0*epsilon)) {
373 0 : result=nn/mm;
374 0 : dfunc=0.5*nn*(nn-mm)/mm;
375 : } else {
376 39118 : double rNdist=Tools::fastpow(rdist,nn-1);
377 39118 : double rMdist=Tools::fastpow(rdist,mm-1);
378 39118 : double num = 1.-rNdist*rdist;
379 39118 : double iden = 1./(1.-rMdist*rdist);
380 39118 : double func = num*iden;
381 : result = func;
382 39118 : dfunc = ((-nn*rNdist*iden)+(func*(iden*mm)*rMdist));
383 : }
384 : }
385 12031290 : return result;
386 : }
387 :
388 11522409 : double SwitchingFunction::calculateSqr(double distance2,double&dfunc)const {
389 11522409 : if(type==rational && nn%2==0 && mm%2==0 && d0==0.0) {
390 8270429 : if(distance2>dmax_2) {
391 118214 : dfunc=0.0;
392 118214 : return 0.0;
393 : }
394 8152215 : const double rdist_2 = distance2*invr0_2;
395 8152215 : double result=do_rational(rdist_2,dfunc,nn/2,mm/2);
396 : // chain rule:
397 8152215 : dfunc*=2*invr0_2;
398 : // stretch:
399 8152215 : result=result*stretch+shift;
400 8152215 : dfunc*=stretch;
401 8152215 : return result;
402 : } else {
403 3251980 : double distance=std::sqrt(distance2);
404 3251980 : return calculate(distance,dfunc);
405 : }
406 : }
407 :
408 30596499 : double SwitchingFunction::calculate(double distance,double&dfunc)const {
409 30596499 : plumed_massert(init,"you are trying to use an unset SwitchingFunction");
410 30596499 : if(distance>dmax) {
411 155450 : dfunc=0.0;
412 155450 : return 0.0;
413 : }
414 30441049 : const double rdist = (distance-d0)*invr0;
415 : double result;
416 :
417 30441049 : if(rdist<=0.) {
418 : result=1.;
419 23808145 : dfunc=0.0;
420 : } else {
421 6632904 : if(type==smap) {
422 0 : double sx=c*pow( rdist, a );
423 0 : result=pow( 1.0 + sx, d );
424 0 : dfunc=-b*sx/rdist*result/(1.0+sx);
425 6632904 : } else if(type==rational) {
426 3879075 : result=do_rational(rdist,dfunc,nn,mm);
427 2753829 : } else if(type==exponential) {
428 1206975 : result=exp(-rdist);
429 1206975 : dfunc=-result;
430 1546854 : } else if(type==nativeq) {
431 278 : double rdist2 = beta*(distance - lambda * ref);
432 278 : double exprdist=exp(rdist2);
433 278 : double exprmdist=1.0/exprdist;
434 278 : result=1./(1.+exprdist);
435 278 : dfunc=-1.0/(exprmdist+1.0)/(1.+exprdist);
436 1546576 : } else if(type==gaussian) {
437 163430 : result=exp(-0.5*rdist*rdist);
438 163430 : dfunc=-rdist*result;
439 1383146 : } else if(type==cubic) {
440 126996 : double tmp1=rdist-1, tmp2=(1+2*rdist);
441 126996 : result=tmp1*tmp1*tmp2;
442 126996 : dfunc=2*tmp1*tmp2 + 2*tmp1*tmp1;
443 1256150 : } else if(type==tanh) {
444 8000 : double tmp1=std::tanh(rdist);
445 8000 : result = 1.0 - tmp1;
446 8000 : dfunc=-(1-tmp1*tmp1);
447 1248150 : } else if(type==leptontype) {
448 51 : const unsigned t=OpenMP::getThreadNum();
449 102 : plumed_assert(t<expression.size());
450 : try {
451 102 : const_cast<lepton::CompiledExpression*>(&expression[t])->getVariableReference("x")=rdist;
452 102 : const_cast<lepton::CompiledExpression*>(&expression_deriv[t])->getVariableReference("x")=rdist;
453 0 : } catch(PLMD::lepton::Exception& exc) {
454 : // this is necessary since in some cases lepton things a variable is not present even though it is present
455 : // e.g. func=0*x
456 : }
457 51 : result=expression[t].evaluate();
458 51 : dfunc=expression_deriv[t].evaluate();
459 : #ifdef __PLUMED_HAS_MATHEVAL
460 1248099 : } else if(type==matheval) {
461 1248099 : const unsigned it=OpenMP::getThreadNum();
462 2496198 : result=evaluator_evaluate_x(evaluator[it],rdist);
463 1248099 : dfunc=evaluator_evaluate_x(evaluator_deriv[it],rdist);
464 : #endif
465 0 : } else plumed_merror("Unknown switching function type");
466 : // this is for the chain rule:
467 6632904 : dfunc*=invr0;
468 : // this is because calculate() sets dfunc to the derivative divided times the distance.
469 : // (I think this is misleading and I would like to modify it - GB)
470 6632904 : dfunc/=distance;
471 : }
472 :
473 30441049 : result=result*stretch+shift;
474 30441049 : dfunc*=stretch;
475 :
476 30441049 : return result;
477 : }
478 :
479 718 : SwitchingFunction::SwitchingFunction():
480 : init(false),
481 : type(rational),
482 : invr0(0.0),
483 : d0(0.0),
484 : dmax(0.0),
485 : nn(6),
486 : mm(0),
487 : a(0.0),
488 : b(0.0),
489 : c(0.0),
490 : d(0.0),
491 : lambda(0.0),
492 : beta(0.0),
493 : ref(0.0),
494 : invr0_2(0.0),
495 : dmax_2(0.0),
496 : stretch(1.0),
497 718 : shift(0.0)
498 : {
499 718 : }
500 :
501 166741966 : SwitchingFunction::SwitchingFunction(const SwitchingFunction&sf):
502 166741966 : init(sf.init),
503 166741966 : type(sf.type),
504 166741966 : invr0(sf.invr0),
505 166741966 : d0(sf.d0),
506 166741966 : dmax(sf.dmax),
507 166741966 : nn(sf.nn),
508 166741966 : mm(sf.mm),
509 166741966 : a(sf.a),
510 166741966 : b(sf.b),
511 166741966 : c(sf.c),
512 166741966 : d(sf.d),
513 166741966 : lambda(sf.lambda),
514 166741966 : beta(sf.beta),
515 166741966 : ref(sf.ref),
516 166741966 : invr0_2(sf.invr0_2),
517 166741966 : dmax_2(sf.dmax_2),
518 166741966 : stretch(sf.stretch),
519 3501581286 : shift(sf.shift)
520 : {
521 : #ifdef __PLUMED_HAS_MATHEVAL
522 166741966 : if(sf.evaluator.size()>0) {
523 0 : const unsigned nt=OpenMP::getNumThreads();
524 0 : plumed_assert(nt>0);
525 0 : evaluator.resize(nt);
526 0 : evaluator_deriv.resize(nt);
527 0 : for(unsigned i=0; i<nt; i++) {
528 0 : evaluator[i]=evaluator_create(evaluator_get_string(sf.evaluator[0]));
529 0 : evaluator_deriv[i]=evaluator_create(evaluator_get_string(sf.evaluator_deriv[0]));
530 : }
531 : }
532 : #endif
533 166741966 : }
534 :
535 0 : SwitchingFunction & SwitchingFunction::operator=(const SwitchingFunction& sf) {
536 0 : if(&sf==this) return *this;
537 0 : init=sf.init;
538 0 : type=sf.type;
539 0 : invr0=sf.invr0;
540 0 : d0=sf.d0;
541 0 : dmax=sf.dmax;
542 0 : nn=sf.nn;
543 0 : mm=sf.mm;
544 0 : a=sf.a;
545 0 : b=sf.b;
546 0 : c=sf.c;
547 0 : d=sf.d;
548 0 : lambda=sf.lambda;
549 0 : beta=sf.beta;
550 0 : ref=sf.ref;
551 0 : invr0_2=sf.invr0_2;
552 0 : dmax_2=sf.dmax_2;
553 0 : stretch=sf.stretch;
554 0 : shift=sf.shift;
555 : #ifdef __PLUMED_HAS_MATHEVAL
556 0 : if(sf.evaluator.size()>0) {
557 0 : const unsigned nt=OpenMP::getNumThreads();
558 0 : plumed_assert(nt>0);
559 0 : evaluator.resize(nt);
560 0 : evaluator_deriv.resize(nt);
561 0 : for(unsigned i=0; i<nt; i++) {
562 0 : evaluator[i]=evaluator_create(evaluator_get_string(sf.evaluator[0]));
563 0 : evaluator_deriv[i]=evaluator_create(evaluator_get_string(sf.evaluator_deriv[0]));
564 : }
565 : }
566 : #endif
567 : return *this;
568 : }
569 :
570 :
571 54 : void SwitchingFunction::set(int nn,int mm,double r0,double d0) {
572 54 : init=true;
573 54 : type=rational;
574 54 : if(mm==0) mm=2*nn;
575 54 : this->nn=nn;
576 54 : this->mm=mm;
577 54 : this->invr0=1.0/r0;
578 54 : this->invr0_2=this->invr0*this->invr0;
579 54 : this->d0=d0;
580 54 : this->dmax=d0+r0*pow(0.00001,1./(nn-mm));
581 54 : this->dmax_2=this->dmax*this->dmax;
582 :
583 : double dummy;
584 54 : double s0=calculate(0.0,dummy);
585 54 : double sd=calculate(dmax,dummy);
586 54 : stretch=1.0/(s0-sd);
587 54 : shift=-sd*stretch;
588 54 : }
589 :
590 6 : double SwitchingFunction::get_r0() const {
591 6 : return 1./invr0;
592 : }
593 :
594 6 : double SwitchingFunction::get_d0() const {
595 6 : return d0;
596 : }
597 :
598 117918069 : double SwitchingFunction::get_dmax() const {
599 117918069 : return dmax;
600 : }
601 :
602 23893569 : double SwitchingFunction::get_dmax2() const {
603 23893569 : return dmax_2;
604 : }
605 :
606 500228052 : SwitchingFunction::~SwitchingFunction() {
607 : #ifdef __PLUMED_HAS_MATHEVAL
608 333485386 : for(unsigned i=0; i<evaluator.size(); i++) evaluator_destroy(evaluator[i]);
609 333485386 : for(unsigned i=0; i<evaluator_deriv.size(); i++) evaluator_destroy(evaluator_deriv[i]);
610 : #endif
611 166742684 : }
612 :
613 :
614 4839 : }
615 :
616 :
617 :
|
