Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-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 "Custom.h"
23 : #include "core/ActionRegister.h"
24 : #include "FunctionShortcut.h"
25 : #include "FunctionOfScalar.h"
26 : #include "FunctionOfVector.h"
27 : #include "FunctionOfMatrix.h"
28 : #include "tools/OpenMP.h"
29 : #include "tools/LeptonCall.h"
30 :
31 : namespace PLMD {
32 : namespace function {
33 :
34 : //+PLUMEDOC FUNCTION CUSTOM
35 : /*
36 : Calculate a combination of variables using a custom expression.
37 :
38 : This action computes an arbitrary function of one or more
39 : collective variables. Arguments are chosen with the ARG keyword,
40 : and the function is provided with the FUNC string. Notice that this
41 : string should contain no space. Within FUNC, one can refer to the
42 : arguments as x,y,z, and t (up to four variables provided as ARG).
43 : This names can be customized using the VAR keyword (see examples below).
44 :
45 : This function is implemented using the Lepton library, that allows to evaluate
46 : algebraic expressions and to automatically differentiate them.
47 :
48 : If you want a function that depends not only on collective variables
49 : but also on time you can use the \subpage TIME action.
50 :
51 : \par Examples
52 :
53 : The following input tells plumed to perform a metadynamics
54 : using as a CV the difference between two distances.
55 : \plumedfile
56 : dAB: DISTANCE ATOMS=10,12
57 : dAC: DISTANCE ATOMS=10,15
58 : diff: CUSTOM ARG=dAB,dAC FUNC=y-x PERIODIC=NO
59 : # notice: the previous line could be replaced with the following
60 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
61 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
62 : \endplumedfile
63 : (see also \ref DISTANCE, \ref COMBINE, and \ref METAD).
64 : Notice that forces applied to diff will be correctly propagated
65 : to atoms 10, 12, and 15.
66 : Also notice that since CUSTOM is used without the VAR option
67 : the two arguments should be referred to as x and y in the expression FUNC.
68 : For simple functions
69 : such as this one it is possible to use \ref COMBINE.
70 :
71 : The following input tells plumed to print the angle between vectors
72 : identified by atoms 1,2 and atoms 2,3
73 : its square (as computed from the x,y,z components) and the distance
74 : again as computed from the square root of the square.
75 : \plumedfile
76 : DISTANCE LABEL=d1 ATOMS=1,2 COMPONENTS
77 : DISTANCE LABEL=d2 ATOMS=2,3 COMPONENTS
78 : CUSTOM ...
79 : LABEL=theta
80 : ARG=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
81 : VAR=ax,ay,az,bx,by,bz
82 : FUNC=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz)))
83 : PERIODIC=NO
84 : ... CUSTOM
85 : PRINT ARG=theta
86 : \endplumedfile
87 : (See also \ref PRINT and \ref DISTANCE).
88 :
89 : Notice that this action implements a large number of functions (trigonometric, exp, log, etc).
90 : Among the useful functions, have a look at the step function (that is the Heaviside function).
91 : `step(x)` is defined as 1 when `x` is positive and `0` when x is negative. This allows for
92 : a straightforward implementation of if clauses.
93 :
94 : For example, imagine that you want to implement a restraint that only acts when a
95 : distance is larger than 0.5. You can do it with
96 : \plumedfile
97 : d: DISTANCE ATOMS=10,15
98 : m: CUSTOM ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
99 : # check the function you are applying:
100 : PRINT ARG=d,m FILE=checkme
101 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
102 : \endplumedfile
103 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
104 :
105 : The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` is:
106 : - If x<0.5 (step(0.5-x)!=0) use 0.5
107 : - If x>0.5 (step(x-0.5)!=0) use x
108 : Notice that the same could have been obtained using an \ref UPPER_WALLS
109 : However, with CUSTOM you can create way more complex definitions.
110 :
111 : \warning If you apply forces on the variable (as in the previous example) you should
112 : make sure that the variable is continuous!
113 : Conversely, if you are just analyzing a trajectory you can safely use
114 : discontinuous variables.
115 :
116 : A possible continuity check with gnuplot is
117 : \verbatim
118 : # this allow to step function to be used in gnuplot:
119 : gnuplot> step(x)=0.5*(erf(x*10000000)+1)
120 : # here you can test your function
121 : gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)
122 : \endverbatim
123 :
124 : Also notice that you can easily make logical operations on the conditions that you
125 : create. The equivalent of the AND operator is the product: `step(1.0-x)*step(x-0.5)` is
126 : 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,
127 : `1-step(1.0-x)*step(x-0.5)` is only equal to 1 when x is outside the 0.5-1.0 interval.
128 :
129 : CUSTOM can be used in combination with \ref DISTANCE to implement variants of the
130 : DISTANCE keyword that were present in PLUMED 1.3 and that allowed to compute
131 : the distance of a point from a line defined by two other points, or the progression
132 : along that line.
133 : \plumedfile
134 : # take center of atoms 1 to 10 as reference point 1
135 : p1: CENTER ATOMS=1-10
136 : # take center of atoms 11 to 20 as reference point 2
137 : p2: CENTER ATOMS=11-20
138 : # take center of atoms 21 to 30 as reference point 3
139 : p3: CENTER ATOMS=21-30
140 :
141 : # compute distances
142 : d12: DISTANCE ATOMS=p1,p2
143 : d13: DISTANCE ATOMS=p1,p3
144 : d23: DISTANCE ATOMS=p2,p3
145 :
146 : # compute progress variable of the projection of point p3
147 : # along the vector joining p1 and p2
148 : # notice that progress is measured from the middle point
149 : onaxis: CUSTOM ARG=d13,d23,d12 FUNC=(0.5*(y^2-x^2)/z) PERIODIC=NO
150 :
151 : # compute between point p3 and the vector joining p1 and p2
152 : fromaxis: CUSTOM 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
153 :
154 : PRINT ARG=onaxis,fromaxis
155 :
156 : \endplumedfile
157 :
158 : Notice that these equations have been used to combine \ref RMSD
159 : from different snapshots of a protein so as to define
160 : progression (S) and distance (Z) variables \cite perez2015atp.
161 :
162 :
163 : */
164 : //+ENDPLUMEDOC
165 :
166 : //+PLUMEDOC FUNCTION MATHEVAL_SCALAR
167 : /*
168 : Calculate a function of a set of input scalars
169 :
170 : See \ref MATHEVAL
171 :
172 : \par Examples
173 :
174 : */
175 : //+ENDPLUMEDOC
176 :
177 : //+PLUMEDOC FUNCTION CUSTOM_SCALAR
178 : /*
179 : Calculate a function of a set of input scalars
180 :
181 : See \ref CUSTOM
182 :
183 : \par Examples
184 :
185 : */
186 : //+ENDPLUMEDOC
187 :
188 : //+PLUMEDOC FUNCTION MATHEVAL_VECTOR
189 : /*
190 : Calculate a function of a set of input vectors elementwise
191 :
192 : See \ref MATHEVAL
193 :
194 : \par Examples
195 :
196 : */
197 : //+ENDPLUMEDOC
198 :
199 : //+PLUMEDOC FUNCTION CUSTOM_VECTOR
200 : /*
201 : Calculate a function of a set of input vectors elementwise
202 :
203 : See \ref CUSTOM
204 :
205 : \par Examples
206 :
207 : */
208 : //+ENDPLUMEDOC
209 :
210 : //+PLUMEDOC COLVAR CUSTOM_MATRIX
211 : /*
212 : Calculate an arbitrary function piecewise for one or multiple input matrices.
213 :
214 : \par Examples
215 :
216 : */
217 : //+ENDPLUMEDOC
218 :
219 : //+PLUMEDOC COLVAR MATHEVAL_MATRIX
220 : /*
221 : Calculate an arbitrary function piecewise for one or multiple input matrices.
222 :
223 : \par Examples
224 :
225 : */
226 : //+ENDPLUMEDOC
227 :
228 : typedef FunctionShortcut<Custom> CustomShortcut;
229 : PLUMED_REGISTER_ACTION(CustomShortcut,"CUSTOM")
230 : PLUMED_REGISTER_ACTION(CustomShortcut,"MATHEVAL")
231 : typedef FunctionOfScalar<Custom> ScalarCustom;
232 : PLUMED_REGISTER_ACTION(ScalarCustom,"CUSTOM_SCALAR")
233 : PLUMED_REGISTER_ACTION(ScalarCustom,"MATHEVAL_SCALAR")
234 : typedef FunctionOfVector<Custom> VectorCustom;
235 : PLUMED_REGISTER_ACTION(VectorCustom,"CUSTOM_VECTOR")
236 : PLUMED_REGISTER_ACTION(VectorCustom,"MATHEVAL_VECTOR")
237 : typedef FunctionOfMatrix<Custom> MatrixCustom;
238 : PLUMED_REGISTER_ACTION(MatrixCustom,"CUSTOM_MATRIX")
239 : PLUMED_REGISTER_ACTION(MatrixCustom,"MATHEVAL_MATRIX")
240 :
241 : //+PLUMEDOC FUNCTION MATHEVAL
242 : /*
243 : An alias to the CUSTOM function that can also be used to calaculate combinations of variables using a custom expression.
244 :
245 : Documentation for this action is identical to that for \ref CUSTOM
246 :
247 : This alias is kept in order to maintain compatibility with previous PLUMED versions.
248 : However, notice that as of PLUMED 2.5 the libmatheval library is not linked anymore,
249 : and the \ref MATHEVAL function is implemented using the Lepton library.
250 :
251 : \par Examples
252 :
253 : Just replace \ref CUSTOM with \ref MATHEVAL.
254 :
255 : \plumedfile
256 : d: DISTANCE ATOMS=10,15
257 : m: MATHEVAL ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
258 : # check the function you are applying:
259 : PRINT ARG=d,m FILE=checkme
260 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
261 : \endplumedfile
262 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
263 :
264 : */
265 : //+ENDPLUMEDOC
266 :
267 10806 : void Custom::registerKeywords(Keywords& keys) {
268 10806 : keys.use("PERIODIC");
269 21612 : keys.add("compulsory","FUNC","the function you wish to evaluate");
270 21612 : 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.");
271 10806 : keys.setValueDescription("an arbitrary function");
272 10806 : }
273 :
274 3005 : void Custom::read( ActionWithArguments* action ) {
275 : // Read in the variables
276 9015 : std::vector<std::string> var; parseVector(action,"VAR",var); parse(action,"FUNC",func);
277 3005 : if(var.size()==0) {
278 2946 : var.resize(action->getNumberOfArguments());
279 2946 : if(var.size()>3) action->error("Using more than 3 arguments you should explicitly write their names with VAR");
280 2946 : if(var.size()>0) var[0]="x";
281 2946 : if(var.size()>1) var[1]="y";
282 2946 : if(var.size()>2) var[2]="z";
283 : }
284 3005 : if(var.size()!=action->getNumberOfArguments()) action->error("Size of VAR array should be the same as number of arguments");
285 : // Check for operations that are not multiplication (this can probably be done much more cleverly)
286 3005 : bool onlymultiplication = func.find("*")!=std::string::npos;
287 : // Find first bracket in expression
288 3005 : if( func.find("(")!=std::string::npos ) {
289 1605 : std::size_t br = func.find_first_of("("); std::string subexpr=func.substr(0,br); onlymultiplication = func.find("*")!=std::string::npos;
290 1605 : if( subexpr.find("/")!=std::string::npos ) { std::size_t sl = func.find_first_of("/"); std::string aa = subexpr.substr(0,sl); subexpr=aa; }
291 1605 : if( subexpr.find("+")!=std::string::npos || subexpr.find("-")!=std::string::npos ) onlymultiplication=false;
292 : // Now work out which vars are in multiplication
293 1498 : if( onlymultiplication ) {
294 2321 : for(unsigned i=0; i<var.size(); ++i) {
295 1465 : if( subexpr.find(var[i])!=std::string::npos &&
296 316 : action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) check_multiplication_vars.push_back(i);
297 : }
298 : }
299 1400 : } else if( func.find("/")!=std::string::npos ) {
300 821 : onlymultiplication=true; if( func.find("+")!=std::string::npos || func.find("-")!=std::string::npos ) onlymultiplication=false;
301 : if( onlymultiplication ) {
302 820 : std::size_t br = func.find_first_of("/"); std::string subexpr=func.substr(0,br);
303 2280 : for(unsigned i=0; i<var.size(); ++i) {
304 1460 : if( subexpr.find(var[i])!=std::string::npos &&
305 31 : action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) check_multiplication_vars.push_back(i);
306 : }
307 : }
308 579 : } else if( func.find("+")!=std::string::npos || func.find("-")!=std::string::npos ) {
309 : onlymultiplication=false;
310 : } else {
311 1028 : for(unsigned i=0; i<var.size(); ++i) {
312 624 : if( action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) check_multiplication_vars.push_back(i);
313 : }
314 : }
315 3005 : if( check_multiplication_vars.size()>0 ) {
316 451 : action->log.printf(" optimizing implementation as function only involves multiplication \n");
317 : }
318 :
319 3005 : action->log.printf(" with function : %s\n",func.c_str());
320 3005 : action->log.printf(" with variables :");
321 7827 : for(unsigned i=0; i<var.size(); i++) action->log.printf(" %s",var[i].c_str());
322 3005 : action->log.printf("\n"); function.set( func, var, action );
323 3005 : std::vector<double> zeros( action->getNumberOfArguments(), 0 ); double fval = abs(function.evaluate(zeros));
324 3005 : zerowhenallzero=(fval<epsilon );
325 3005 : if( zerowhenallzero ) action->log.printf(" not calculating when all arguments are zero \n");
326 3005 : }
327 :
328 5 : std::string Custom::getGraphInfo( const std::string& name ) const {
329 10 : return FunctionTemplateBase::getGraphInfo( name ) + + "\n" + "FUNC=" + func;
330 : }
331 :
332 1819 : bool Custom::getDerivativeZeroIfValueIsZero() const {
333 1819 : return check_multiplication_vars.size()>0;
334 : }
335 :
336 966 : std::vector<Value*> Custom::getArgumentsToCheck( const std::vector<Value*>& args ) {
337 966 : std::vector<Value*> fargs( check_multiplication_vars.size() );
338 2093 : for(unsigned i=0; i<check_multiplication_vars.size(); ++i) fargs[i] = args[check_multiplication_vars[i]];
339 966 : return fargs;
340 : }
341 :
342 20478013 : void Custom::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
343 20478013 : if( args.size()>1 ) {
344 : bool allzero=false;
345 15780962 : if( check_multiplication_vars.size()>0 ) {
346 17313345 : for(unsigned i=0; i<check_multiplication_vars.size(); ++i) {
347 10885579 : if( fabs(args[check_multiplication_vars[i]])<epsilon ) { allzero=true; break; }
348 : }
349 5142171 : } else if( zerowhenallzero ) {
350 2860794 : allzero=(fabs(args[0])<epsilon);
351 3296324 : for(unsigned i=1; i<args.size(); ++i) {
352 3151815 : if( fabs(args[i])>epsilon ) { allzero=false; break; }
353 : }
354 : }
355 13499585 : if( allzero ) {
356 13352933 : vals[0]=0; for(unsigned i=0; i<args.size(); i++) derivatives(0,i) = 0.0;
357 : return;
358 : }
359 : }
360 16126419 : vals[0] = function.evaluate( args );
361 16126419 : if( !noderiv ) {
362 47278265 : for(unsigned i=0; i<args.size(); i++) derivatives(0,i) = function.evaluateDeriv( i, args );
363 : }
364 : }
365 :
366 : }
367 : }
368 :
369 :
|
