Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "ActionRegister.h"
23 : #include "Function.h"
24 :
25 : #ifdef __PLUMED_HAS_MATHEVAL
26 : #include <matheval.h>
27 : #endif
28 :
29 : #include "lepton/Lepton.h"
30 :
31 : using namespace std;
32 :
33 : namespace PLMD {
34 : namespace function {
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 :
53 : //+PLUMEDOC FUNCTION MATHEVAL
54 : /*
55 : Calculate a combination of variables using a matheval expression.
56 :
57 : This action computes an arbitrary function of one or more precomputed
58 : collective variables. Arguments are chosen with the ARG keyword,
59 : and the function is provided with the FUNC string. Notice that this
60 : string should contain no space. Within FUNC, one can refer to the
61 : arguments as x,y,z, and t (up to four variables provided as ARG).
62 : This names can be customized using the VAR keyword (see examples below).
63 :
64 : If you want a function that depends not only on collective variables
65 : but also on time you can use the \subpage TIME action.
66 :
67 : \attention
68 : The MATHEVAL object only works if one of these conditions is satisfied:
69 : (a) libmatheval is installed on the system and PLUMED has been linked to it or
70 : (b) the environment variable `PLUMED_USE_LEPTON` is set equal to `yes` at runtime
71 : (in this case, the internal lepton library will be used).
72 : Notice that in version v2.5 the internal lepton library will be used by default.
73 :
74 : \par Examples
75 :
76 : The following input tells plumed to perform a metadynamics
77 : using as a CV the difference between two distances.
78 : \plumedfile
79 : dAB: DISTANCE ATOMS=10,12
80 : dAC: DISTANCE ATOMS=10,15
81 : diff: MATHEVAL ARG=dAB,dAC FUNC=y-x PERIODIC=NO
82 : # notice: the previous line could be replaced with the following
83 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
84 : METAD ARG=diff WIDTH=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
85 : \endplumedfile
86 : (see also \ref DISTANCE, \ref COMBINE, and \ref METAD).
87 : Notice that forces applied to diff will be correctly propagated
88 : to atoms 10, 12, and 15.
89 : Also notice that since MATHEVAL is used without the VAR option
90 : the two arguments should be referred to as x and y in the expression FUNC.
91 : For simple functions
92 : such as this one it is possible to use \ref COMBINE, which does
93 : not require libmatheval to be installed on your system.
94 :
95 : The following input tells plumed to print the angle between vectors
96 : identified by atoms 1,2 and atoms 2,3
97 : its square (as computed from the x,y,z components) and the distance
98 : again as computed from the square root of the square.
99 : \plumedfile
100 : DISTANCE LABEL=d1 ATOMS=1,2 COMPONENTS
101 : DISTANCE LABEL=d2 ATOMS=2,3 COMPONENTS
102 : MATHEVAL ...
103 : LABEL=theta
104 : ARG=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
105 : VAR=ax,ay,az,bx,by,bz
106 : FUNC=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz))
107 : PERIODIC=NO
108 : ... MATHEVAL
109 : PRINT ARG=theta
110 : \endplumedfile
111 : (See also \ref PRINT and \ref DISTANCE).
112 :
113 : Notice that the matheval library implements a large number of functions (trigonometric, exp, log, etc).
114 : Among the useful functions, have a look at the step function (that is the Heaviside function).
115 : `step(x)` is defined as 1 when `x` is positive and `0` when x is negative. This allows for
116 : a straightforward implementation of if clauses.
117 :
118 : For example, imagine that you want to implement a restraint that only acts when a
119 : distance is larger than 0.5. You can do it with
120 : \plumedfile
121 : d: DISTANCE ATOMS=10,15
122 : m: MATHEVAL ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
123 : # check the function you are applying:
124 : PRINT ARG=d,n FILE=checkme
125 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
126 : \endplumedfile
127 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
128 :
129 : The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` is:
130 : - If x<0.5 (step(0.5-x)!=0) use 0.5
131 : - If x>0.5 (step(x-0.5)!=0) use x
132 : Notice that the same could have been obtained using an \ref UPPER_WALLS
133 : However, with MATHEVAL you can create way more complex definitions.
134 :
135 : \warning If you apply forces on the variable (as in the previous example) you should
136 : make sure that the variable is continuous!
137 : Conversely, if you are just analyzing a trajectory you can safely use
138 : discontinuous variables.
139 :
140 : A possible continuity check with gnuplot is
141 : \verbatim
142 : # this allow to step function to be used in gnuplot:
143 : gnuplot> step(x)=0.5*(erf(x*10000000)+1)
144 : # here you can test your function
145 : gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)
146 : \endverbatim
147 :
148 : Also notice that you can easily make logical operations on the conditions that you
149 : create. The equivalent of the AND operator is the product: `step(1.0-x)*step(x-0.5)` is
150 : only equal to 1 when x is between 0.5 and 1.0. By combining negation and AND you can obtain an OR. That is,
151 : `1-step(1.0-x)*step(x-0.5)` is only equal to 1 when x is outside the 0.5-1.0 interval.
152 :
153 : MATHEVAL can be used in combination with \ref DISTANCE to implement variants of the
154 : DISTANCE keyword that were present in PLUMED 1.3 and that allowed to compute
155 : the distance of a point from a line defined by two other points, or the progression
156 : along that line.
157 : \plumedfile
158 : # take center of atoms 1 to 10 as reference point 1
159 : p1: CENTER ATOMS=1-10
160 : # take center of atoms 11 to 20 as reference point 2
161 : p2: CENTER ATOMS=11-20
162 : # take center of atoms 21 to 30 as reference point 3
163 : p3: CENTER ATOMS=21-30
164 :
165 : # compute distances
166 : d12: DISTANCE ATOMS=p1,p2
167 : d13: DISTANCE ATOMS=p1,p3
168 : d23: DISTANCE ATOMS=p2,p3
169 :
170 : # compute progress variable of the projection of point p3
171 : # along the vector joining p1 and p2
172 : # notice that progress is measured from the middle point
173 : onaxis: MATHEVAL ARG=d13,d23,d12 FUNC=(0.5*(y^2-x^2)/z) PERIODIC=NO
174 :
175 : # compute between point p3 and the vector joining p1 and p2
176 : fromaxis: MATHEVAL ARG=d13,d23,d12,onaxis VAR=x,y,z,o FUNC=(0.5*(y^2+x^2)-o^2-0.25*z^2) PERIODIC=NO
177 :
178 : PRINT ARG=onaxis,fromaxis
179 :
180 : \endplumedfile
181 :
182 : Notice that these equations have been used to combine \ref RMSD
183 : from different snapshots of a protein so as to define
184 : progression (S) and distance (Z) variables \cite perez2015atp.
185 :
186 :
187 : */
188 : //+ENDPLUMEDOC
189 :
190 :
191 : class Matheval :
192 : public Function
193 : {
194 : const bool use_lepton;
195 : lepton::CompiledExpression expression;
196 : std::vector<lepton::CompiledExpression> expression_deriv;
197 : void* evaluator;
198 : vector<void*> evaluator_deriv;
199 : vector<string> var;
200 : string func;
201 : vector<double> values;
202 : vector<char*> names;
203 : public:
204 : explicit Matheval(const ActionOptions&);
205 : ~Matheval();
206 : void calculate();
207 : static void registerKeywords(Keywords& keys);
208 : };
209 :
210 6714 : PLUMED_REGISTER_ACTION(Matheval,"MATHEVAL")
211 :
212 : //+PLUMEDOC FUNCTION CUSTOM
213 : /*
214 : An alias to the \ref MATHEVAL function.
215 :
216 : \par Examples
217 :
218 : Just replace \ref MATHEVAL with \ref CUSTOM.
219 :
220 : \plumedfile
221 : d: DISTANCE ATOMS=10,15
222 : m: CUSTOM ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
223 : # check the function you are applying:
224 : PRINT ARG=d,n FILE=checkme
225 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
226 : \endplumedfile
227 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
228 :
229 : */
230 : //+ENDPLUMEDOC
231 :
232 : class Custom :
233 : public Matheval {
234 : };
235 :
236 6453 : PLUMED_REGISTER_ACTION(Matheval,"CUSTOM")
237 :
238 265 : void Matheval::registerKeywords(Keywords& keys) {
239 265 : Function::registerKeywords(keys);
240 795 : keys.use("ARG"); keys.use("PERIODIC");
241 1060 : keys.add("compulsory","FUNC","the function you wish to evaluate");
242 1060 : keys.add("optional","VAR","the names to give each of the arguments in the function. If you have up to three arguments in your function you can use x, y and z to refer to them. Otherwise you must use this flag to give your variables names.");
243 265 : }
244 :
245 263 : Matheval::Matheval(const ActionOptions&ao):
246 : Action(ao),
247 : Function(ao),
248 263 : use_lepton(std::getenv("PLUMED_USE_LEPTON")),
249 : expression_deriv(getNumberOfArguments()),
250 : evaluator(NULL),
251 : evaluator_deriv(getNumberOfArguments(),NULL),
252 : values(getNumberOfArguments()),
253 2104 : names(getNumberOfArguments())
254 : {
255 526 : parseVector("VAR",var);
256 263 : if(var.size()==0) {
257 249 : var.resize(getNumberOfArguments());
258 249 : if(getNumberOfArguments()>3)
259 0 : error("Using more than 3 arguments you should explicitly write their names with VAR");
260 249 : if(var.size()>0) var[0]="x";
261 249 : if(var.size()>1) var[1]="y";
262 249 : if(var.size()>2) var[2]="z";
263 : }
264 263 : if(var.size()!=getNumberOfArguments())
265 0 : error("Size of VAR array should be the same as number of arguments");
266 526 : parse("FUNC",func);
267 263 : addValueWithDerivatives();
268 263 : checkRead();
269 :
270 526 : log.printf(" with function : %s\n",func.c_str());
271 263 : log.printf(" with variables :");
272 1645 : for(unsigned i=0; i<var.size(); i++) log.printf(" %s",var[i].c_str());
273 263 : log.printf("\n");
274 :
275 263 : if(use_lepton) {
276 37 : log<<" WARNING: you are using lepton as a replacement for libmatheval\n";
277 74 : lepton::ParsedExpression pe=lepton::Parser::parse(func).optimize(leptonConstants);
278 37 : log<<" function as parsed by lepton: "<<pe<<"\n";
279 37 : expression=pe.createCompiledExpression();
280 117 : for(auto &p: expression.getVariables()) {
281 86 : if(std::find(var.begin(),var.end(),p)==var.end()) {
282 0 : error("variable " + p + " is not defined");
283 : }
284 : }
285 37 : log<<" derivatives as computed by lepton:\n";
286 131 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
287 188 : lepton::ParsedExpression pe=lepton::Parser::parse(func).differentiate(var[i]).optimize(leptonConstants);
288 47 : log<<" "<<pe<<"\n";
289 94 : expression_deriv[i]=pe.createCompiledExpression();
290 : }
291 : } else {
292 : #ifdef __PLUMED_HAS_MATHEVAL
293 226 : evaluator=evaluator_create(const_cast<char*>(func.c_str()));
294 226 : if(!evaluator) error("There was some problem in parsing matheval formula "+func);
295 : char **check_names;
296 : int check_count;
297 226 : evaluator_get_variables(evaluator,&check_names,&check_count);
298 226 : if(check_count!=int(getNumberOfArguments())) {
299 : string sc;
300 0 : Tools::convert(check_count,sc);
301 0 : error("Your function string contains "+sc+" arguments. This should be equal to the number of ARGs");
302 : }
303 878 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
304 : bool found=false;
305 1894 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
306 1568 : if(var[i]==check_names[j])found=true;
307 : }
308 326 : if(!found)
309 0 : error("Variable "+var[i]+" cannot be found in your function string");
310 : }
311 878 : for(unsigned i=0; i<getNumberOfArguments(); i++)
312 978 : evaluator_deriv[i]=evaluator_derivative(evaluator,const_cast<char*>(var[i].c_str()));
313 226 : log.printf(" function as parsed by matheval: %s\n", evaluator_get_string(evaluator));
314 226 : log.printf(" derivatives as computed by matheval:\n");
315 1430 : for(unsigned i=0; i<var.size(); i++) log.printf(" %s\n",evaluator_get_string(evaluator_deriv[i]));
316 : #else
317 : error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
318 : #endif
319 : }
320 263 : }
321 :
322 19747 : void Matheval::calculate() {
323 19747 : if(use_lepton) {
324 457 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
325 : try {
326 487 : expression.getVariableReference(var[i])=getArgument(i);
327 20 : } catch(PLMD::lepton::Exception& exc) {
328 : // this is necessary since in some cases lepton things a variable is not present even though it is present
329 : // e.g. func=0*x
330 : }
331 : }
332 119 : setValue(expression.evaluate());
333 457 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
334 947 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
335 : try {
336 1239 : expression_deriv[i].getVariableReference(var[j])=getArgument(j);
337 317 : } catch(PLMD::lepton::Exception& exc) {
338 : // this is necessary since in some cases lepton things a variable is not present even though it is present
339 : // e.g. func=0*x
340 : }
341 : }
342 338 : setDerivative(i,expression_deriv[i].evaluate());
343 : }
344 : } else {
345 : #ifdef __PLUMED_HAS_MATHEVAL
346 60166 : for(unsigned i=0; i<getNumberOfArguments(); i++) values[i]=getArgument(i);
347 60166 : for(unsigned i=0; i<getNumberOfArguments(); i++) names[i]=const_cast<char*>(var[i].c_str());
348 39256 : setValue(evaluator_evaluate(evaluator,names.size(),&names[0],&values[0]));
349 60166 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
350 40538 : setDerivative(i,evaluator_evaluate(evaluator_deriv[i],names.size(),&names[0],&values[0]));
351 : }
352 : #else
353 : error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
354 : #endif
355 : }
356 19747 : }
357 :
358 1052 : Matheval::~Matheval() {
359 263 : if(evaluator) {
360 : #ifdef __PLUMED_HAS_MATHEVAL
361 226 : evaluator_destroy(evaluator);
362 1430 : for(unsigned i=0; i<evaluator_deriv.size(); i++)evaluator_destroy(evaluator_deriv[i]);
363 : #else
364 : error("MATHEVAL not available, please export PLUMED_USE_LEPTON=yes");
365 : #endif
366 : }
367 526 : }
368 :
369 : }
370 4839 : }
371 :
372 :
|
