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 : CUSTOM is one of the most useful actions in PLUMED. This action takes in a list of arguments and then uses
39 : the [lepton mathematical expression parser](https://simtk.org/projects/lepton) to evaluate a user defined function
40 : of these input arguments. We can thus use this action in the input below to perform a metadynamics
41 : simulation using the difference between two distances as a CV.
42 :
43 : ```plumed
44 : dAB: DISTANCE ATOMS=10,12
45 : dAC: DISTANCE ATOMS=10,15
46 : diff: CUSTOM ARG=dAB,dAC FUNC=y-x PERIODIC=NO
47 : # notice: the previous line could be replaced with the following
48 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
49 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
50 : ```
51 :
52 : The particular function that should be evaluated from the input arguments is specified using the `FUNC` keyword.
53 : The function provided to the `FUNC` keyword is written in terms of `x` and `y` in the input above. `x` here deonetes the
54 : first argument provided to the ARG keyword, `dAB`, while `y` is the second argument provided to the `ARG` keyword, `dAC`.
55 :
56 : ## The VAR keyword
57 :
58 : If you wish, you can rewrite the example input above more transparantly as:
59 :
60 : ```plumed
61 : dAB: DISTANCE ATOMS=10,12
62 : dAC: DISTANCE ATOMS=10,15
63 : diff: CUSTOM ARG=dAB,dAC VAR=dAB,dAC FUNC=dAC-dAB PERIODIC=NO
64 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
65 : ```
66 :
67 : By using the `VAR` keyword here we ensure that we can use the labels of the arguments in the mathematical expression that
68 : we provide to the `FUNC` argument. Notice that, if you have four or more arguments you must use the `VAR` keyword. With
69 : three or less arguments the `VAR` keyword can be ommitted and the symbols `x`, `y` and `z` can be used to denote the first,
70 : second and third arguments respectively.
71 :
72 : The following input illustrates a case where using the `VAR` keyword is essential. This input tells PLUMED to
73 : print the angle between the vector connecting atoms 1,2 and the vector connecting atoms 2,3.
74 :
75 : ```plumed
76 : d1: DISTANCE ATOMS=1,2 COMPONENTS
77 : d2: DISTANCE ATOMS=2,3 COMPONENTS
78 : theta: CUSTOM ...
79 : ARG=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
80 : VAR=ax,ay,az,bx,by,bz
81 : FUNC=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz)))
82 : PERIODIC=NO
83 : ...
84 :
85 : PRINT ARG=theta
86 : ```
87 :
88 : ## Functions and constants
89 :
90 : The previous examples demonstrates how CUSTOM can be used to evaluate mathematical expressions that involve taking powers (^), adding (+),
91 : multiplying (*), dividing (/) and subtracting (-) variables and can control the order in which operations are performed by using parenthesis.
92 : In addition to these basic binary operations you can also use a range of mathematical functions in your expressions. The following table lists all the mathematical functions
93 : that you can use in in the input to the `FUNC` keyword for CUSTOM.
94 :
95 : | Function | Description |
96 : |:--------:|:------------|
97 : | sqrt(x) | The square root of x |
98 : | exp(x) | The exponential of x i.e. $e^{x}$ |
99 : | log(x) | The natural logarithm of x |
100 : | sin(x) | The sine of x |
101 : | cos(x) | The cosine of x |
102 : | sec(x) | The reciprocal of the cosine of $x$. $\frac{1}{\cos(x)}$ |
103 : | csc(x) | The reciprocal of the sine of $x$. $\frac{1}{\sin(x)}$ |
104 : | tan(x) | The tangent of x i.e. $\frac{\sin(x)}{\cos(x)}$ |
105 : | cot(x) | The reciprocal of the tangent of $x$. $\frac{1}{\tan(x)}$ |
106 : | asin(x) | The principal arc sine of x. Returns $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $x = \sin(y)$ |
107 : | acos(x) | The principal arc cosine of x. Returns $0 \le y \le \pi$ which gives $x=\cos(y)$ |
108 : | atan(x) | The principal arc tangent of x. Returns $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $x = \tan(y)$ |
109 : | atan2(x,y) | The principal arg tangent of $\frac{x}{y}$. Returns $-\pi \le z \le \pi$ which gives $\frac{x}{y} = \tan(z)$ |
110 : | sinh(x) | The hyperbolic sine of $x$ |
111 : | cosh(x) | the hyperbolic cosine of $x$ |
112 : | tanh(x) | The hyperbolic tangent of $x$ |
113 : | erf(x) | The [error function](https://en.wikipedia.org/wiki/Error_function) $\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2} \textrm{d}t$ |
114 : | erfc(x) | The complementary error function $1-\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2} \textrm{d}t$ |
115 : | step(x) | 1 if $x \ge 0$ and 0 otherwise |
116 : | delta(x) | inf if $x=0$ and 0 otherwise |
117 : | nandelta(x) | nan if $x=0$ and 0 otherwise |
118 : | square(x) | The square of x i.e. $x^2$ |
119 : | cube(x) | The cube of x i.e. $x^3$ |
120 : | recip(x) | The reciprocal of x i.e. $\frac{1}{x}$ |
121 : | min(x,y) | If $x<y$ this function returns $x$. If $y\le x$ this function returns $y$ |
122 : | max(x,y) | If $x>y$ this function returns $x$. If $y\ge x$ this function returns $y$ |
123 : | abs(x) | The absolute value of x |
124 : | floor(x) | The largest integer that is less than $x$ |
125 : | ceil(x) | The smallest integer that is greater than $x$ |
126 : | select(x,y,z) | If $x==0$ this returns $z$ otherwise this returns $y$ |
127 : | acot(x) | Returns the value of $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $\frac{1}{x}=\tan(y)$ |
128 : | asec(x) | Returns the value of $0 \le y \le \pi$ which gives $\frac{1}{x}=\cos(y)$ |
129 : | acsc(x) | Returns the value of $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $\frac{1}{x}=\sin(y)$ |
130 : | coth(x) | The recipricoal of the hyperbolic tangent of $x$. $\frac{1}{\tanh(x)}$ |
131 : | sech(x) | The recipricoal of the hyperbolic cosine of $x$. $\frac{1}{\cosh(x)}$ |
132 : | csch(x) | The recipricoal of the hyperbolic sine of $x$. $\frac{1}{\sinh(x)}$ |
133 : | asinh(x) | The nonnegative area hyperbolic sine of $x$. Returns $y$ which gives $x=\sinh(y)$ |
134 : | acosh(x) | The nonnegative area hyperbolic cosine of $x$. Returns $0\le y \le \infty$ which gives $x=\cosh(y)$ |
135 : | atanh(x) | The nonnegative area hyperbolic tangent of $-1 \le x \le 1$. Returns $y$ which gives $x=\tanh(y)$. |
136 : | acoth(x) | The inverse hyperbolic tangent of $x$ is calculated as $\frac{1}{2}\ln\left( \frac{x+1}{x-1} \right)$ |
137 : | asech(x) | The inverse hyperbolic secant of $x$ is calculated as $\log\left( \sqrt{\frac{1}{x}-1}\sqrt{\frac{1}{x}+1} + \frac{1}{x}\right)$ |
138 : | acsch(x) | The inverse hyperbolic cosecant of $x$ is calculated as $\log\left(\frac{1}{x}+\sqrt{\frac{1}{x^2}+1}\right)$ |
139 :
140 : Notice, that you can also incorporate the following constants in the expressions that are used in the input to FUNC:
141 :
142 : | Symbol | Description |
143 : | :----- |:------------|
144 : | `e` | Euler's number - the base for the natural logarithm |
145 : | `log2e` | $1 / \log(2)$ |
146 : | `log10e` | $1 / \log(10)$ |
147 : | `ln2` | $\log(2)$ |
148 : | `ln10` | $\log(10)$ |
149 : | `pi` | the circle constant $\pi$ |
150 : | `pi_2` | $\pi / 2$ |
151 : | `pi_4` | $\pi / 4$ |
152 : | `sqrt2 | $\sqrt(2)$ |
153 : | `sqrt1_2 ` | $\sqrt(0.5)$ |
154 :
155 : The way one of these constants can be used in an expression is illustrated in the following example:
156 :
157 : ```plumed
158 : d: DISTANCE ATOMS=1,2
159 : # This function is evaluating the area of the circle whose radius is
160 : # given by the distance between atoms 1 and 2.
161 : c: CUSTOM ARG=d FUNC=pi*x*x PERIODIC=NO
162 : PRINT ARG=c FILE=colvar
163 : ```
164 :
165 : ## The step function
166 :
167 : The `step` operation (that is the Heaviside function) from the table above allow you to use if clauses in CUSTOM actions.
168 : As discussed in the table `step(x)` is 1 when `x` is positive and `0` when `x` is negative. So the function `step(1-x)` is
169 : 1 if x is less than one and zero otherwise.
170 :
171 : Using the `step` operation multiple times in a function allows you to perform logical operations. For example, the equivalent of the AND operator
172 : is the product so, for example, `step(1.0-x)*step(x-0.5)` is only equal to 1 when x is greater than 0.5 AND less than 1.0.
173 : By a similar logic we can use the function `1-step(1.0-x)*step(x-0.5)` to create a value that is 1 if x is less than 0.5 OR
174 : greater than 1.0.
175 :
176 : The example below illustrtes how we can put these ideas of using the `step` operation into practise.
177 :
178 : ```plumed
179 : d: DISTANCE ATOMS=10,15
180 : m: CUSTOM ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
181 : # check the function you are applying:
182 : PRINT ARG=d,m FILE=checkme
183 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
184 : ```
185 :
186 : The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` in this example is:
187 : - If x<0.5 (step(0.5-x)!=0) use 0.5
188 : - If x>0.5 (step(x-0.5)!=0) use x
189 : Notice that the same result can be achieved by using [UPPER_WALLS](UPPER_WALLS.md)
190 : However, with CUSTOM you can create much more complex definitions.
191 :
192 : Notice that we can apply a force on the value `m` by using the [RESTRAINT](RESTRAINT.md) command as the function we defined in the expression that was passed to the `FUNC` keyword is continuous.
193 : In general, however, you must be careful when using the `step`, `delta`, `nandelta` and `select` functions as you can easily write expression for
194 : discontinuous functions by using these operations. If you want to check if a function you have created using `step`
195 : is continuous you can easily plot the function in gnuplot by using a commands like those shown below
196 :
197 : ````
198 : # this allow to step function to be used in gnuplot:
199 : gnuplot> step(x)=0.5*(erf(x*10000000)+1)
200 : # here you can test your function
201 : gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)
202 : ````
203 :
204 : ## Using TIME as a variable
205 :
206 : Notice that you can use CUSTOM to implement a [MOVINGRESTRAINT](MOVINGRESTRAINT.md) as shown below.
207 :
208 : ```plumed
209 : t: TIME
210 : d: DISTANCE ATOMS=1,2
211 : f: CUSTOM ARG=d,t FUNC=100*(x-((0.2-0.1)*y/100))^2 PERIODIC=NO
212 : BIASVALUE ARG=f
213 : ```
214 :
215 : In a 100~ps simulation that was run with this input the distance beetween atom 1 and 2 would be forced to increase from 0.1 to 0.2 nm.
216 :
217 : ## Using CUSTOM in shortcuts
218 :
219 : The CUSTOM action is used in many of the [shortcut actions](shortcuts.md) that are implemented in PLUMED. We think that using this action in these places is beneficial as it ensures that the mathematical
220 : expressions that are used in the method are visible in the log. We have found that there are many actions that
221 : are used in relatively few papers. When implementing these actions we think that sharing implementations of these methods that are comprehensible is more important than sharing methods that are
222 : fast.
223 :
224 : As an example of why this is useful consider some of the variants of the DISTANCE keyword that were present in PLUMED 1.3. These variants allowed you to compute the distance between a point and a line defined by
225 : two other points or the progression along that line. In PLUMED 2.10 we can implement these variables using the following input file:
226 :
227 : ```plumed
228 : # take center of atoms 1 to 10 as reference point 1
229 : p1: CENTER ATOMS=1-10
230 : # take center of atoms 11 to 20 as reference point 2
231 : p2: CENTER ATOMS=11-20
232 : # take center of atoms 21 to 30 as reference point 3
233 : p3: CENTER ATOMS=21-30
234 :
235 : # compute distances
236 : d12: DISTANCE ATOMS=p1,p2
237 : d13: DISTANCE ATOMS=p1,p3
238 : d23: DISTANCE ATOMS=p2,p3
239 :
240 : # compute progress variable of the projection of point p3
241 : # along the vector joining p1 and p2
242 : # notice that progress is measured from the middle point
243 : onaxis: CUSTOM ARG=d13,d23,d12 FUNC=(0.5*(y^2-x^2)/z) PERIODIC=NO
244 :
245 : # compute between point p3 and the vector joining p1 and p2
246 : 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
247 :
248 : PRINT ARG=onaxis,fromaxis
249 : ```
250 :
251 : The equations in this input were also used to combine [RMSD](RMSD.md) values from different snapshots of a protein so as to define
252 : progression (S) and distance (Z) variables in the paper that is cited in the bibliography. We can understand how these expressions
253 : are derived by noting that $x$, $y$ and $z$ are the distances between atoms 1 and 3, 2 and 3 and 1 and 2 respectively. The projection
254 : of the vector connecting atom 1 to atom 3 onto the vector connecting atom 1 to atom 2 is thus $x\cos(\theta)$, where theta is the angle
255 : between the vector connecting atoms 1 and 3 and the vector connecting atoms 1 and 2. We can arrive at the following expression for $x\cos(\theta)$
256 : by rearranging the cosine rule:
257 :
258 : $$
259 : x\cos(\theta) = \frac{y^2 - x^2}{z} - \frac{z}{2}
260 : $$
261 :
262 : Notice that the value called `onaxis` in the above input is thus $o=x\cos(\theta) + \frac{z}{2}$. Adding the factor of $\frac{z}{2}$ ensures that the origin
263 : is at the center of the bond connecting atom 1 to atom 2.
264 :
265 : The value `fromaxis` measures the square of the distance from the the line. It is calculated using pythagoras theorem as follows:
266 :
267 : $$
268 : f^2 = y^2 - x^2\cos^2(\theta)
269 : $$
270 :
271 : Inserting $x\cos(\theta) = o - \frac{z}{2}$ into this expression gives:
272 :
273 : $$
274 : f^2 = y^2 - (o -\frac{z}{2})^2 = y^2 - o^2 + oz - \frac{z^2}{4}
275 : $$
276 :
277 : Inserting the fact that $oz = \frac{y^2 - x^2}{2}$, which comes from the expression for $o$ that was used to calculate `onaxis`, gets us to the expression that is used to calculate `fromaxis`.
278 :
279 : ## CUSTOM with vector arguments
280 :
281 : The examples above have shown how CUSTOM operates when the input arguments are scalars. You can also pass vectors in the argument to a CUSTOM action. The function you define will then
282 : be applied to each element of the vector in turn. So, for example in the input below a vector that contains three angles is passed to the CUSTOM action. The CUSTOM action calculates the
283 : cosine of these three angles and outputs them in a three dimensional vector called `c` that is printed to the colvar file.
284 :
285 : ```plumed
286 : a: ANGLE ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
287 : c: CUSTOM ARG=a FUNC=cos(x) PERIODIC=NO
288 : PRINT ARG=c FILE=colvar
289 : ```
290 :
291 : You are not confined to passing a single vector to CUSTOM. The input below shows what you can do by pass two vectors with the same numbers of elements. The first element of the vector
292 : output by the CUSTOM action here contains the projection of the vector connecting atom 1 and 2 on the vector connecting atom 2 and 3, the second element contains the projection of the
293 : vector connecting atoms 4 and 5 on the vector connecting atoms 5 and 6 and the final element contains the projection of the vector connecting atoms 7 and 8 on the vector connecting atoms
294 : 8 and 9.
295 :
296 : ```plumed
297 : d: DISTANCE ATOMS1=1,2 ATOMS2=4,5 ATOMS3=7,8
298 : a: ANGLE ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
299 : c: CUSTOM ARG=d,a FUNC=x*cos(y) PERIODIC=NO
300 : PRINT ARG=c FILE=colvar
301 : ```
302 :
303 : Notice, when we multiply two vectors in CUSTOM the output is a vector. This product that emerges from using a CUSTOM action is __not__ the scalar or cross product of the input vectors.
304 :
305 : Lastly, notice that you can pass a mixture of scalars and vectors in the input to a CUSTOM action as shown below.
306 :
307 : ```plumed
308 : d: DISTANCE ATOMS=1,2
309 : a: ANGLE ATOMS1=1,2,3 ATOMS2=1,2,4 ATOMS3=1,2,5 ATOMS4=1,2,6
310 : c: CUSTOM ARG=d,a FUNC=x*cos(y) PERIODIC=NO
311 : PRINT ARG=c FILE=colvar
312 : ```
313 :
314 : The multiplication of a scalar by a vector in the above input is done in [the usual way](https://en.wikipedia.org/wiki/Scalar_multiplication).
315 : Similarly, dividing a vector by a scalar is equivalent to multiplying the vector by the reciprocal of the scalar. If you write an expression that adds or subtract a
316 : scalar from a vector addition is performed. To understand why consider the following example that adds the constant `c` to the input vector of distances `d`. The
317 : result of this operation is a vector `f` that contains the three input distances from `d` with 23 added to each of them.
318 :
319 : ```plumed
320 : c: CONSTANT VALUE=23
321 : d: DISTANCE ATOMS1=1,2 ATOMS2=4,5 ATOMS3=7,8
322 : f: CUSTOM ARG=d,c FUNC=x+y PERIODIC=NO
323 : PRINT ARG=f FILE=colvar
324 : ```
325 :
326 : ## CUSTOM with matrix arguments
327 :
328 : You can also pass matrices in the argument to a CUSTOM action. These input matrices are treated similarly to input vectors. In other words, any function you define
329 : is applied to each element of the matrix in turn so if the input matrix is $N \times $M$ the output matrix is also $N \times M$. The following example illustrates how you
330 : can use this functionality to calculate all the angles between a set of bond vectors.
331 :
332 : ```plumed
333 : # Calculate the vectors connecting four atoms
334 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
335 : # Calculate the norm of these four vectors
336 : norm: CUSTOM ARG=d.x,d.y,d.z FUNC=sqrt(x*x+y*y+z*z) PERIODIC=NO
337 : # Now calculate the directors of the vectors
338 : norm_x: CUSTOM ARG=d.x,norm FUNC=x/y PERIODIC=NO
339 : norm_y: CUSTOM ARG=d.y,norm FUNC=x/y PERIODIC=NO
340 : norm_z: CUSTOM ARG=d.z,norm FUNC=x/y PERIODIC=NO
341 : # And combine all these directors in a matrix
342 : stack: VSTACK ARG=norm_x,norm_y,norm_z
343 : # Now calculate the matrix of dot products between directors (these are the cosines of the angles)
344 : stackT: TRANSPOSE ARG=stack
345 : cosa: MATRIX_PRODUCT ARG=stack,stackT
346 : # And finally get the 4x4 matrix of angles and print it to a file
347 : angles: CUSTOM ARG=cosa FUNC=acos(x) PERIODIC=NO
348 : PRINT ARG=angles FILE=colvar
349 : ```
350 :
351 : Notice that you can pass multiple $N\times M$ matrices in the input to a CUSTOM action as illustrated in the example below:
352 :
353 : ```plumed
354 : c1: CONSTANT VALUES=2,3,4,5 NROWS=2 NCOLS=2
355 : c2: CONSTANT VALUES=1,0,0,1 NROWS=2 NCOLS=2
356 : f: CUSTOM ARG=c1,c2 FUNC=x*y PERIODIC=NO
357 : PRINT ARG=f FILE=colvar
358 : ```
359 :
360 : The four by four matrix that is calculated by the custom action in this input is given by:
361 :
362 : $$
363 : f = \left(
364 : \begin{matrix}
365 : 2 & 0 \\
366 : 0 & 5
367 : \end{matrix}
368 : \right)
369 : $$
370 :
371 : which is the element-wise [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)) and __not__ the matrix product.
372 : By a similar token if you have a custom command that takes two matrices in input and use `FUNC=x/y` the Hadamard product between the matrix
373 : x and a matrix that contains the reciprocals of each of the elements in y is computed. If you wish to calcalculate the
374 : [product of two matrices](https://en.wikipedia.org/wiki/Matrix_multiplication) you should use the [MATRIX_PRODUCT](MATRIX_PRODUCT.md) comamnd.
375 : Similarly, if you want to calculate the product of a matrix and a vector you should use the [MATRIX_VECTOR_PRODUCT](MATRIX_VECTOR_PRODUCT.md)
376 : command.
377 :
378 : Lastly, note that you can pass a mixture of scalars and $N\times M$ matrices in the input to a CUSTOM command. As with vectors, you can think of
379 : any scalars you pass as being converted into $N\times M$ matrix in which every element is equal to the input scalar.
380 :
381 : ## CUSTOM with grid arguments
382 :
383 : CUSTOM will also accept a function on a grid in input. You might use this feature if you want to calculate a free energy from a histogram
384 : using $F(s) = -k_BT\log[H(s)]$ as illustrated in the input below:
385 :
386 : ```plumed
387 : x: DISTANCE ATOMS=1,2
388 : hA: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
389 : fE: CUSTOM ARG=hA FUNC=-2.5*log(x) PERIODIC=NO
390 : DUMPGRID ARG=fE FILE=fes.dat
391 : ```
392 :
393 : Another way that this functonality is routinely used in PLUMED is in calculating the density field for a symmetry function. The example below shows how you
394 : would do this in practise. The final output here is a function on a grid that tells you the average value of the FCC order parameter at each point in the cell.
395 :
396 : ```plumed
397 : # Evaluate the FCC order parameter for all of the atoms
398 : fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
399 : # Calculate the distance between each atom and the origin on atom 1
400 : dens_dist: DISTANCES ORIGIN=1 ATOMS=fcc COMPONENTS
401 : # Do a KDE using the FCC order parameters for the weights of each gaussian and the positions of the atoms as the centers
402 : dens_numer: KDE HEIGHTS=fcc_n ARG=dens_dist.x,dens_dist.y,dens_dist.z GRID_BIN=14,14,28 BANDWIDTH=1.0,1.0,1.0
403 : # Estimate the density of atoms at each point in the box
404 : ones: ONES SIZE=5184
405 : dens_denom: KDE ARG=dens_dist.x,dens_dist.y,dens_dist.z GRID_BIN=14,14,28 HEIGHTS=ones BANDWIDTH=1.0,1.0,1.0
406 : # Now calculate the average value of the order parameter at each point in the box
407 : dens: CUSTOM ARG=dens_numer,dens_denom FUNC=x/y PERIODIC=NO
408 : DUMPCUBE ARG=dens FILE=dens.cube FMT=%8.4f
409 : ```
410 :
411 : */
412 : //+ENDPLUMEDOC
413 :
414 : //+PLUMEDOC FUNCTION MATHEVAL_SCALAR
415 : /*
416 : Calculate a function of a set of input scalars
417 :
418 : See \ref MATHEVAL
419 :
420 : \par Examples
421 :
422 : */
423 : //+ENDPLUMEDOC
424 :
425 : //+PLUMEDOC FUNCTION CUSTOM_SCALAR
426 : /*
427 : Calculate a function of a set of input scalars
428 :
429 : See \ref CUSTOM
430 :
431 : \par Examples
432 :
433 : */
434 : //+ENDPLUMEDOC
435 :
436 : //+PLUMEDOC FUNCTION MATHEVAL_VECTOR
437 : /*
438 : Calculate a function of a set of input vectors elementwise
439 :
440 : See \ref MATHEVAL
441 :
442 : \par Examples
443 :
444 : */
445 : //+ENDPLUMEDOC
446 :
447 : //+PLUMEDOC FUNCTION CUSTOM_VECTOR
448 : /*
449 : Calculate a function of a set of input vectors elementwise
450 :
451 : See \ref CUSTOM
452 :
453 : \par Examples
454 :
455 : */
456 : //+ENDPLUMEDOC
457 :
458 : //+PLUMEDOC COLVAR CUSTOM_MATRIX
459 : /*
460 : Calculate an arbitrary function piecewise for one or multiple input matrices.
461 :
462 : \par Examples
463 :
464 : */
465 : //+ENDPLUMEDOC
466 :
467 : //+PLUMEDOC COLVAR MATHEVAL_MATRIX
468 : /*
469 : Calculate an arbitrary function piecewise for one or multiple input matrices.
470 :
471 : \par Examples
472 :
473 : */
474 : //+ENDPLUMEDOC
475 :
476 : typedef FunctionShortcut<Custom> CustomShortcut;
477 : PLUMED_REGISTER_ACTION(CustomShortcut,"CUSTOM")
478 : PLUMED_REGISTER_ACTION(CustomShortcut,"MATHEVAL")
479 : typedef FunctionOfScalar<Custom> ScalarCustom;
480 : PLUMED_REGISTER_ACTION(ScalarCustom,"CUSTOM_SCALAR")
481 : PLUMED_REGISTER_ACTION(ScalarCustom,"MATHEVAL_SCALAR")
482 : typedef FunctionOfVector<Custom> VectorCustom;
483 : PLUMED_REGISTER_ACTION(VectorCustom,"CUSTOM_VECTOR")
484 : PLUMED_REGISTER_ACTION(VectorCustom,"MATHEVAL_VECTOR")
485 : typedef FunctionOfMatrix<Custom> MatrixCustom;
486 : PLUMED_REGISTER_ACTION(MatrixCustom,"CUSTOM_MATRIX")
487 : PLUMED_REGISTER_ACTION(MatrixCustom,"MATHEVAL_MATRIX")
488 :
489 : //+PLUMEDOC FUNCTION MATHEVAL
490 : /*
491 : An alias to the CUSTOM function that can also be used to calaculate combinations of variables using a custom expression.
492 :
493 : The documentation for this action is identical to that for [CUSTOM](CUSTOM.md). You can thus use it to evaluate a evaluate an arbitrary function as in the following example input:
494 :
495 : ```plumed
496 : dAB: DISTANCE ATOMS=10,12
497 : dAC: DISTANCE ATOMS=10,15
498 : diff: MATHEVAL ARG=dAB,dAC FUNC=y-x PERIODIC=NO
499 : # notice: the previous line could be replaced with the following
500 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
501 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
502 : ```
503 :
504 : This alias is kept in order to maintain compatibility with previous PLUMED versions.
505 : However, notice that as of PLUMED 2.5 the libmatheval library is not linked anymore,
506 : and that the MATHEVAL action evaluates functions [the Lepton library](https://simtk.org/projects/lepton).
507 :
508 : \par Examples
509 :
510 : Just replace \ref CUSTOM with \ref MATHEVAL.
511 :
512 : \plumedfile
513 : d: DISTANCE ATOMS=10,15
514 : m: MATHEVAL ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
515 : # check the function you are applying:
516 : PRINT ARG=d,m FILE=checkme
517 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
518 : \endplumedfile
519 : (see also \ref DISTANCE, \ref PRINT, and \ref RESTRAINT)
520 :
521 : */
522 : //+ENDPLUMEDOC
523 :
524 10770 : void Custom::registerKeywords(Keywords& keys) {
525 10770 : keys.use("PERIODIC");
526 10770 : keys.add("compulsory","FUNC","the function you wish to evaluate");
527 10770 : 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.");
528 21540 : keys.setValueDescription("scalar/vector/matrix/grid","an arbitrary function");
529 10770 : keys.addDOI("10.1093/nar/gkv872");
530 10770 : }
531 :
532 2996 : void Custom::read( ActionWithArguments* action ) {
533 : // Read in the variables
534 : std::vector<std::string> var;
535 2996 : parseVector(action,"VAR",var);
536 5992 : parse(action,"FUNC",func);
537 2996 : if(var.size()==0) {
538 2937 : var.resize(action->getNumberOfArguments());
539 2937 : if(var.size()>3) {
540 0 : action->error("Using more than 3 arguments you should explicitly write their names with VAR");
541 : }
542 2937 : if(var.size()>0) {
543 : var[0]="x";
544 : }
545 2937 : if(var.size()>1) {
546 : var[1]="y";
547 : }
548 2937 : if(var.size()>2) {
549 : var[2]="z";
550 : }
551 : }
552 2996 : if(var.size()!=action->getNumberOfArguments()) {
553 0 : action->error("Size of VAR array should be the same as number of arguments");
554 : }
555 : // Check for operations that are not multiplication (this can probably be done much more cleverly)
556 2996 : bool onlymultiplication = func.find("*")!=std::string::npos;
557 : // Find first bracket in expression
558 2996 : if( func.find("(")!=std::string::npos ) {
559 1596 : std::size_t br = func.find_first_of("(");
560 1596 : std::string subexpr=func.substr(0,br);
561 1596 : onlymultiplication = func.find("*")!=std::string::npos;
562 1596 : if( subexpr.find("/")!=std::string::npos ) {
563 159 : std::size_t sl = func.find_first_of("/");
564 159 : std::string aa = subexpr.substr(0,sl);
565 : subexpr=aa;
566 : }
567 1596 : if( subexpr.find("+")!=std::string::npos || subexpr.find("-")!=std::string::npos ) {
568 : onlymultiplication=false;
569 : }
570 : // Now work out which vars are in multiplication
571 1489 : if( onlymultiplication ) {
572 2321 : for(unsigned i=0; i<var.size(); ++i) {
573 1465 : if( subexpr.find(var[i])!=std::string::npos &&
574 : action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) {
575 316 : check_multiplication_vars.push_back(i);
576 : }
577 : }
578 : }
579 1400 : } else if( func.find("/")!=std::string::npos ) {
580 : onlymultiplication=true;
581 821 : if( func.find("+")!=std::string::npos || func.find("-")!=std::string::npos ) {
582 : onlymultiplication=false;
583 : }
584 : if( onlymultiplication ) {
585 820 : std::size_t br = func.find_first_of("/");
586 820 : std::string subexpr=func.substr(0,br);
587 2280 : for(unsigned i=0; i<var.size(); ++i) {
588 1460 : if( subexpr.find(var[i])!=std::string::npos &&
589 : action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) {
590 31 : check_multiplication_vars.push_back(i);
591 : }
592 : }
593 : }
594 579 : } else if( func.find("+")!=std::string::npos || func.find("-")!=std::string::npos ) {
595 : onlymultiplication=false;
596 : } else {
597 1028 : for(unsigned i=0; i<var.size(); ++i) {
598 624 : if( action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) {
599 197 : check_multiplication_vars.push_back(i);
600 : }
601 : }
602 : }
603 2996 : if( check_multiplication_vars.size()>0 ) {
604 451 : action->log.printf(" optimizing implementation as function only involves multiplication \n");
605 : }
606 :
607 2996 : action->log.printf(" with function : %s\n",func.c_str());
608 2996 : action->log.printf(" with variables :");
609 7809 : for(unsigned i=0; i<var.size(); i++) {
610 4813 : action->log.printf(" %s",var[i].c_str());
611 : }
612 2996 : action->log.printf("\n");
613 2996 : function.set( func, var, action );
614 2996 : std::vector<double> zeros( action->getNumberOfArguments(), 0 );
615 2996 : double fval = abs(function.evaluate(zeros));
616 2996 : zerowhenallzero=(fval<epsilon );
617 2996 : if( zerowhenallzero ) {
618 1201 : action->log.printf(" not calculating when all arguments are zero \n");
619 : }
620 2996 : }
621 :
622 5 : std::string Custom::getGraphInfo( const std::string& name ) const {
623 10 : return FunctionTemplateBase::getGraphInfo( name ) + + "\n" + "FUNC=" + func;
624 : }
625 :
626 1819 : bool Custom::getDerivativeZeroIfValueIsZero() const {
627 1819 : return check_multiplication_vars.size()>0;
628 : }
629 :
630 966 : std::vector<Value*> Custom::getArgumentsToCheck( const std::vector<Value*>& args ) {
631 966 : std::vector<Value*> fargs( check_multiplication_vars.size() );
632 2093 : for(unsigned i=0; i<check_multiplication_vars.size(); ++i) {
633 1127 : fargs[i] = args[check_multiplication_vars[i]];
634 : }
635 966 : return fargs;
636 : }
637 :
638 20478004 : void Custom::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
639 20478004 : if( args.size()>1 ) {
640 : bool allzero=false;
641 15780962 : if( check_multiplication_vars.size()>0 ) {
642 17313345 : for(unsigned i=0; i<check_multiplication_vars.size(); ++i) {
643 10885579 : if( fabs(args[check_multiplication_vars[i]])<epsilon ) {
644 : allzero=true;
645 : break;
646 : }
647 : }
648 5142171 : } else if( zerowhenallzero ) {
649 2860794 : allzero=(fabs(args[0])<epsilon);
650 3296324 : for(unsigned i=1; i<args.size(); ++i) {
651 3151815 : if( fabs(args[i])>epsilon ) {
652 : allzero=false;
653 : break;
654 : }
655 : }
656 : }
657 13499585 : if( allzero ) {
658 4351594 : vals[0]=0;
659 13352933 : for(unsigned i=0; i<args.size(); i++) {
660 9001339 : derivatives(0,i) = 0.0;
661 : }
662 : return;
663 : }
664 : }
665 16126410 : vals[0] = function.evaluate( args );
666 16126410 : if( !noderiv ) {
667 47278265 : for(unsigned i=0; i<args.size(); i++) {
668 31160198 : derivatives(0,i) = function.evaluateDeriv( i, args );
669 : }
670 : }
671 : }
672 :
673 : }
674 : }
675 :
676 :
|
