Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 : using namespace std;
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 :
35 : //+PLUMEDOC FUNCTION MATHEVAL
36 : /*
37 : Calculate a combination of variables using a matheval expression.
38 :
39 : This action computes an arbitrary function of one or more precomputed
40 : collective variables. Arguments are chosen with the ARG keyword,
41 : and the function is provided with the FUNC string. Notice that this
42 : string should contain no space. Within FUNC, one can refer to the
43 : arguments as x,y,z, and t (up to four variables provided as ARG).
44 : This names can be customized using the VAR keyword (see examples below).
45 :
46 : If you want a function that depends not only on collective variables
47 : but also on time you can use the \subpage TIME action.
48 :
49 : \attention
50 : The MATHEVAL object only works if libmatheval is installed on the system and
51 : PLUMED has been linked to it
52 :
53 : \par Examples
54 :
55 : The following input tells plumed to perform a metadynamics
56 : using as a CV the difference between two distances.
57 : \verbatim
58 : dAB: DISTANCE ATOMS=10,12
59 : dAC: DISTANCE ATOMS=10,15
60 : diff: MATHEVAL ARG=dAB,dAC FUNC=y-x PERIODIC=NO
61 : # notice: the previous line could be replaced with the following
62 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
63 : METAD ARG=diff WIDTH=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
64 : \endverbatim
65 : (see also \ref DISTANCE, \ref COMBINE, and \ref METAD).
66 : Notice that forces applied to diff will be correctly propagated
67 : to atoms 10, 12, and 15.
68 : Also notice that since MATHEVAL is used without the VAR option
69 : the two arguments should be referred to as x and y in the expression FUNC.
70 : For simple functions
71 : such as this one it is possible to use \ref COMBINE, which does
72 : not require libmatheval to be installed on your system.
73 :
74 : The following input tells plumed to print the angle between vectors
75 : identified by atoms 1,2 and atoms 2,3
76 : its square (as computed from the x,y,z components) and the distance
77 : again as computed from the square root of the square.
78 : \verbatim
79 : DISTANCE LABEL=d1 ATOMS=1,2 COMPONENTS
80 : DISTANCE LABEL=d2 ATOMS=2,3 COMPONENTS
81 : MATHEVAL ...
82 : LABEL=theta
83 : ARG=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
84 : VAR=ax,ay,az,bx,by,bz
85 : FUNC=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz))
86 : PERIODIC=NO
87 : ... MATHEVAL
88 : PRINT ARG=theta
89 : \endverbatim
90 : (See also \ref PRINT and \ref DISTANCE).
91 :
92 : Notice that the matheval library implements a large number of functions (trigonometric, exp, log, etc).
93 : Among the useful functions, have a look at the step function (that is the Heaviside function).
94 : `step(x)` is defined as 1 when `x` is positive and `0` when x is negative. This allows for
95 : a straightforward implementation of if clauses.
96 :
97 : For example, imagine that you want to implement a restraint that only acts when a
98 : distance is larger than 0.5. You can do it with
99 : \verbatim
100 : d: DISTANCE ATOMS=10,15
101 : m: MATHEVAL ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
102 : # check the function you are applying:
103 : PRINT ARG=d,n FILE=checkme
104 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
105 : \endverbatim
106 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
107 :
108 : The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` is:
109 : - If x<0.5 (step(0.5-x)!=0) use 0.5
110 : - If x>0.5 (step(x-0.5)!=0) use x
111 : Notice that the same could have been obtained using an \ref UPPER_WALLS
112 : However, with MATHEVAL you can create way more complex definitions.
113 :
114 : \warning If you apply forces on the variable (as in the previous example) you should
115 : make sure that the variable is continuous!
116 : Conversely, if you are just analyzing a trajectory you can safely use
117 : discontinuous variables.
118 :
119 : A possible continuity check with gnuplot is
120 : \verbatim
121 : # this allow to step function to be used in gnuplot:
122 : gnuplot> step(x)=0.5*(erf(x*10000000)+1)
123 : # here you can test your function
124 : gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)
125 : \endverbatim
126 :
127 : Also notice that you can easily make logical operations on the conditions that you
128 : create. The equivalent of the AND operator is the product: `step(1.0-x)*step(x-0.5)` is
129 : 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,
130 : `1-step(1.0-x)*step(x-0.5)` is only equal to 1 when x is outside the 0.5-1.0 interval.
131 :
132 : MATHEVAL can be used in combination with \ref DISTANCE to implement variants of the
133 : DISTANCE keyword that were present in PLUMED 1.3 and that allowed to compute
134 : the distance of a point from a line defined by two other points, or the progression
135 : along that line.
136 : \verbatim
137 : # take center of atoms 1 to 10 as reference point 1
138 : p1: CENTER ATOMS=1-10
139 : # take center of atoms 11 to 20 as reference point 2
140 : p2: CENTER ATOMS=11-20
141 : # take center of atoms 21 to 30 as reference point 3
142 : p3: CENTER ATOMS=21-30
143 :
144 : # compute distances
145 : d12: DISTANCE ATOMS=p1,p2
146 : d13: DISTANCE ATOMS=p1,p3
147 : d23: DISTANCE ATOMS=p2,p3
148 :
149 : # compute progress variable of the projection of point p3
150 : # along the vector joining p1 and p2
151 : # notice that progress is measured from the middle point
152 : onaxis: MATHEVAL ARG=d13,d23,d12 FUNC=(0.5*(y^2-x^2)/z) PERIODIC=NO
153 :
154 : # compute between point p3 and the vector joining p1 and p2
155 : 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
156 :
157 : PRINT ARG=onaxis,fromaxis
158 :
159 : \endverbatim
160 :
161 : Notice that these equations have been used to combine \ref RMSD
162 : from different snapshots of a protein so as to define
163 : progression (S) and distance (Z) variables \cite perez2015atp.
164 :
165 :
166 : */
167 : //+ENDPLUMEDOC
168 :
169 :
170 : class Matheval :
171 : public Function
172 : {
173 : void* evaluator;
174 : vector<void*> evaluator_deriv;
175 : vector<string> var;
176 : string func;
177 : vector<double> values;
178 : vector<char*> names;
179 : public:
180 : explicit Matheval(const ActionOptions&);
181 : ~Matheval();
182 : void calculate();
183 : static void registerKeywords(Keywords& keys);
184 : };
185 :
186 : #ifdef __PLUMED_HAS_MATHEVAL
187 2557 : PLUMED_REGISTER_ACTION(Matheval,"MATHEVAL")
188 :
189 35 : void Matheval::registerKeywords(Keywords& keys) {
190 35 : Function::registerKeywords(keys);
191 35 : keys.use("ARG"); keys.use("PERIODIC");
192 35 : keys.add("compulsory","FUNC","the function you wish to evaluate");
193 35 : 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.");
194 35 : }
195 :
196 34 : Matheval::Matheval(const ActionOptions&ao):
197 : Action(ao),
198 : Function(ao),
199 34 : evaluator_deriv(getNumberOfArguments()),
200 34 : values(getNumberOfArguments()),
201 102 : names(getNumberOfArguments())
202 : {
203 34 : parseVector("VAR",var);
204 34 : if(var.size()==0) {
205 21 : var.resize(getNumberOfArguments());
206 21 : if(getNumberOfArguments()>3)
207 0 : error("Using more than 3 arguments you should explicitly write their names with VAR");
208 21 : if(var.size()>0) var[0]="x";
209 21 : if(var.size()>1) var[1]="y";
210 21 : if(var.size()>2) var[2]="z";
211 : }
212 34 : if(var.size()!=getNumberOfArguments())
213 0 : error("Size of VAR array should be the same as number of arguments");
214 34 : parse("FUNC",func);
215 34 : addValueWithDerivatives();
216 34 : checkRead();
217 :
218 34 : evaluator=evaluator_create(const_cast<char*>(func.c_str()));
219 :
220 34 : if(!evaluator) error("There was some problem in parsing matheval formula "+func);
221 :
222 : char **check_names;
223 : int check_count;
224 34 : evaluator_get_variables(evaluator,&check_names,&check_count);
225 34 : if(check_count!=int(getNumberOfArguments())) {
226 0 : string sc;
227 0 : Tools::convert(check_count,sc);
228 0 : error("Your function string contains "+sc+" arguments. This should be equal to the number of ARGs");
229 : }
230 168 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
231 134 : bool found=false;
232 726 : for(unsigned j=0; j<getNumberOfArguments(); j++) {
233 592 : if(var[i]==check_names[j])found=true;
234 : }
235 134 : if(!found)
236 0 : error("Variable "+var[i]+" cannot be found in your function string");
237 : }
238 :
239 168 : for(unsigned i=0; i<getNumberOfArguments(); i++)
240 134 : evaluator_deriv[i]=evaluator_derivative(evaluator,const_cast<char*>(var[i].c_str()));
241 :
242 :
243 34 : log.printf(" with function : %s\n",func.c_str());
244 34 : log.printf(" with variables :");
245 34 : for(unsigned i=0; i<var.size(); i++) log.printf(" %s",var[i].c_str());
246 34 : log.printf("\n");
247 34 : log.printf(" function as parsed by matheval: %s\n", evaluator_get_string(evaluator));
248 34 : log.printf(" derivatives as computed by matheval:\n");
249 34 : for(unsigned i=0; i<var.size(); i++) log.printf(" %s\n",evaluator_get_string(evaluator_deriv[i]));
250 34 : }
251 :
252 236 : void Matheval::calculate() {
253 236 : for(unsigned i=0; i<getNumberOfArguments(); i++) values[i]=getArgument(i);
254 236 : for(unsigned i=0; i<getNumberOfArguments(); i++) names[i]=const_cast<char*>(var[i].c_str());
255 236 : setValue(evaluator_evaluate(evaluator,names.size(),&names[0],&values[0]));
256 :
257 1113 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
258 877 : setDerivative(i,evaluator_evaluate(evaluator_deriv[i],names.size(),&names[0],&values[0]));
259 : }
260 236 : }
261 :
262 136 : Matheval::~Matheval() {
263 34 : evaluator_destroy(evaluator);
264 34 : for(unsigned i=0; i<evaluator_deriv.size(); i++)evaluator_destroy(evaluator_deriv[i]);
265 102 : }
266 :
267 : #endif
268 :
269 : }
270 2523 : }
271 :
272 :
|
