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 : #include <cmath>
32 : #include <set>
33 :
34 : namespace PLMD {
35 : namespace function {
36 :
37 : //+PLUMEDOC FUNCTION CUSTOM
38 : /*
39 : Calculate a combination of variables using a custom expression.
40 :
41 : CUSTOM is one of the most useful actions in PLUMED. This action takes in a list of arguments and then uses
42 : the [lepton mathematical expression parser](https://simtk.org/projects/lepton) to evaluate a user defined function
43 : of these input arguments. We can thus use this action in the input below to perform a metadynamics
44 : simulation using the difference between two distances as a CV.
45 :
46 : ```plumed
47 : dAB: DISTANCE ATOMS=10,12
48 : dAC: DISTANCE ATOMS=10,15
49 : diff: CUSTOM ARG=dAB,dAC FUNC=y-x PERIODIC=NO
50 : # notice: the previous line could be replaced with the following
51 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
52 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
53 : ```
54 :
55 : The particular function that should be evaluated from the input arguments is specified using the `FUNC` keyword.
56 : The function provided to the `FUNC` keyword is written in terms of `x` and `y` in the input above. `x` here deonetes the
57 : first argument provided to the ARG keyword, `dAB`, while `y` is the second argument provided to the `ARG` keyword, `dAC`.
58 :
59 : ## The VAR keyword
60 :
61 : If you wish, you can rewrite the example input above more transparantly as:
62 :
63 : ```plumed
64 : dAB: DISTANCE ATOMS=10,12
65 : dAC: DISTANCE ATOMS=10,15
66 : diff: CUSTOM ARG=dAB,dAC VAR=dAB,dAC FUNC=dAC-dAB PERIODIC=NO
67 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
68 : ```
69 :
70 : By using the `VAR` keyword here we ensure that we can use the labels of the arguments in the mathematical expression that
71 : we provide to the `FUNC` argument. Notice that, if you have four or more arguments you must use the `VAR` keyword. With
72 : three or less arguments the `VAR` keyword can be ommitted and the symbols `x`, `y` and `z` can be used to denote the first,
73 : second and third arguments respectively.
74 :
75 : The following input illustrates a case where using the `VAR` keyword is essential. This input tells PLUMED to
76 : print the angle between the vector connecting atoms 1,2 and the vector connecting atoms 2,3.
77 :
78 : ```plumed
79 : d1: DISTANCE ATOMS=1,2 COMPONENTS
80 : d2: DISTANCE ATOMS=2,3 COMPONENTS
81 : theta: CUSTOM ...
82 : ARG=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
83 : VAR=ax,ay,az,bx,by,bz
84 : FUNC=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz)))
85 : PERIODIC=NO
86 : ...
87 :
88 : PRINT ARG=theta
89 : ```
90 :
91 : ## Functions and constants
92 :
93 : The previous examples demonstrates how CUSTOM can be used to evaluate mathematical expressions that involve taking powers (^), adding (+),
94 : multiplying (*), dividing (/) and subtracting (-) variables and can control the order in which operations are performed by using parenthesis.
95 : 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
96 : that you can use in the input to the `FUNC` keyword for CUSTOM.
97 :
98 : | Function | Description |
99 : |:--------:|:------------|
100 : | sqrt(x) | The square root of x |
101 : | exp(x) | The exponential of x i.e. $e^{x}$ |
102 : | log(x) | The natural logarithm of x |
103 : | sin(x) | The sine of x |
104 : | cos(x) | The cosine of x |
105 : | sec(x) | The reciprocal of the cosine of $x$. $\frac{1}{\cos(x)}$ |
106 : | csc(x) | The reciprocal of the sine of $x$. $\frac{1}{\sin(x)}$ |
107 : | tan(x) | The tangent of x i.e. $\frac{\sin(x)}{\cos(x)}$ |
108 : | cot(x) | The reciprocal of the tangent of $x$. $\frac{1}{\tan(x)}$ |
109 : | asin(x) | The principal arc sine of x. Returns $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $x = \sin(y)$ |
110 : | acos(x) | The principal arc cosine of x. Returns $0 \le y \le \pi$ which gives $x=\cos(y)$ |
111 : | atan(x) | The principal arc tangent of x. Returns $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $x = \tan(y)$ |
112 : | atan2(x,y) | The principal arg tangent of $\frac{x}{y}$. Returns $-\pi \le z \le \pi$ which gives $\frac{x}{y} = \tan(z)$ |
113 : | sinh(x) | The hyperbolic sine of $x$ |
114 : | cosh(x) | the hyperbolic cosine of $x$ |
115 : | tanh(x) | The hyperbolic tangent of $x$ |
116 : | 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$ |
117 : | erfc(x) | The complementary error function $1-\frac{2}{\sqrt{\pi}}\int_0^x e^{-t^2} \textrm{d}t$ |
118 : | step(x) | 1 if $x \ge 0$ and 0 otherwise |
119 : | delta(x) | inf if $x=0$ and 0 otherwise |
120 : | nandelta(x) | nan if $x=0$ and 0 otherwise |
121 : | square(x) | The square of x i.e. $x^2$ |
122 : | cube(x) | The cube of x i.e. $x^3$ |
123 : | recip(x) | The reciprocal of x i.e. $\frac{1}{x}$ |
124 : | min(x,y) | If $x<y$ this function returns $x$. If $y\le x$ this function returns $y$ |
125 : | max(x,y) | If $x>y$ this function returns $x$. If $y\ge x$ this function returns $y$ |
126 : | abs(x) | The absolute value of x |
127 : | floor(x) | The largest integer that is less than $x$ |
128 : | ceil(x) | The smallest integer that is greater than $x$ |
129 : | select(x,y,z) | If $x==0$ this returns $z$ otherwise this returns $y$ |
130 : | acot(x) | Returns the value of $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $\frac{1}{x}=\tan(y)$ |
131 : | asec(x) | Returns the value of $0 \le y \le \pi$ which gives $\frac{1}{x}=\cos(y)$ |
132 : | acsc(x) | Returns the value of $-\frac{\pi}{2} \le y \le \frac{\pi}{2}$ which gives $\frac{1}{x}=\sin(y)$ |
133 : | coth(x) | The recipricoal of the hyperbolic tangent of $x$. $\frac{1}{\tanh(x)}$ |
134 : | sech(x) | The recipricoal of the hyperbolic cosine of $x$. $\frac{1}{\cosh(x)}$ |
135 : | csch(x) | The recipricoal of the hyperbolic sine of $x$. $\frac{1}{\sinh(x)}$ |
136 : | asinh(x) | The nonnegative area hyperbolic sine of $x$. Returns $y$ which gives $x=\sinh(y)$ |
137 : | acosh(x) | The nonnegative area hyperbolic cosine of $x$. Returns $0\le y \le \infty$ which gives $x=\cosh(y)$ |
138 : | atanh(x) | The nonnegative area hyperbolic tangent of $-1 \le x \le 1$. Returns $y$ which gives $x=\tanh(y)$. |
139 : | acoth(x) | The inverse hyperbolic tangent of $x$ is calculated as $\frac{1}{2}\ln\left( \frac{x+1}{x-1} \right)$ |
140 : | 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)$ |
141 : | acsch(x) | The inverse hyperbolic cosecant of $x$ is calculated as $\log\left(\frac{1}{x}+\sqrt{\frac{1}{x^2}+1}\right)$ |
142 :
143 : Notice, that you can also incorporate the following constants in the expressions that are used in the input to FUNC:
144 :
145 : | Symbol | Description |
146 : | :----- |:------------|
147 : | `e` | Euler's number - the base for the natural logarithm |
148 : | `log2e` | $1 / \log(2)$ |
149 : | `log10e` | $1 / \log(10)$ |
150 : | `ln2` | $\log(2)$ |
151 : | `ln10` | $\log(10)$ |
152 : | `pi` | the circle constant $\pi$ |
153 : | `pi_2` | $\pi / 2$ |
154 : | `pi_4` | $\pi / 4$ |
155 : | `sqrt2 | $\sqrt(2)$ |
156 : | `sqrt1_2 ` | $\sqrt(0.5)$ |
157 :
158 : The way one of these constants can be used in an expression is illustrated in the following example:
159 :
160 : ```plumed
161 : d: DISTANCE ATOMS=1,2
162 : # This function is evaluating the area of the circle whose radius is
163 : # given by the distance between atoms 1 and 2.
164 : c: CUSTOM ARG=d FUNC=pi*x*x PERIODIC=NO
165 : PRINT ARG=c FILE=colvar
166 : ```
167 :
168 : ## The PERIODIC keyword
169 :
170 : If the function you have evaluated has a periodic codomain then you have to use the PERIODIC keyword to specify the codomain of the function.
171 : To see how to use this keyword consider the following input:
172 :
173 : ```plumed
174 : d: DISTANCE ATOMS=1,2 COMPONENTS
175 : theta: CUSTOM ARG=d.x,d.y FUNC=atan2(y,x) PERIODIC=-pi,+pi
176 : PRINT ARG=theta FILE=colvar
177 : ```
178 :
179 : The above input calculates the torsional angle between the vector connecting atoms 1 and 2 and positive $x$ direction around the $z$ axis.
180 : This quantity must lie between $-\pi$ and $\pi$ so we use the PERIODIC keyword to tell PLUMED this fact about the output quantity.
181 :
182 : Notice that PLUMED cannot automatically detect if the codomain of the function is periodic. You thus always have to specifiy the
183 : periodicity of the codomain whenever you use the CUSTOM command. However, the codomain of most functions that people use are not periodic
184 : so you can normally just write `PERIODIC=NO`
185 :
186 : ## The step function
187 :
188 : The `step` operation (that is the Heaviside function) from the table above allow you to use if clauses in CUSTOM actions.
189 : 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
190 : 1 if x is less than one and zero otherwise.
191 :
192 : Using the `step` operation multiple times in a function allows you to perform logical operations. For example, the equivalent of the AND operator
193 : 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.
194 : 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
195 : greater than 1.0.
196 :
197 : The example below illustrtes how we can put these ideas of using the `step` operation into practise.
198 :
199 : ```plumed
200 : d: DISTANCE ATOMS=10,15
201 : m: CUSTOM ARG=d FUNC=0.5*step(0.5-x)+x*step(x-0.5) PERIODIC=NO
202 : # check the function you are applying:
203 : PRINT ARG=d,m FILE=checkme
204 : RESTRAINT ARG=d AT=0.5 KAPPA=10.0
205 : ```
206 :
207 : The meaning of the function `0.5*step(0.5-x)+x*step(x-0.5)` in this example is:
208 : - If x<0.5 (step(0.5-x)!=0) use 0.5
209 : - If x>0.5 (step(x-0.5)!=0) use x
210 : Notice that the same result can be achieved by using [UPPER_WALLS](UPPER_WALLS.md)
211 : However, with CUSTOM you can create much more complex definitions.
212 :
213 : 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.
214 : In general, however, you must be careful when using the `step`, `delta`, `nandelta` and `select` functions as you can easily write expression for
215 : discontinuous functions by using these operations. If you want to check if a function you have created using `step`
216 : is continuous you can easily plot the function in gnuplot by using a commands like those shown below
217 :
218 : ````
219 : # this allow to step function to be used in gnuplot:
220 : gnuplot> step(x)=0.5*(erf(x*10000000)+1)
221 : # here you can test your function
222 : gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)
223 : ````
224 :
225 : ## Using TIME as a variable
226 :
227 : Notice that you can use CUSTOM to implement a [MOVINGRESTRAINT](MOVINGRESTRAINT.md) as shown below.
228 :
229 : ```plumed
230 : t: TIME
231 : d: DISTANCE ATOMS=1,2
232 : f: CUSTOM ARG=d,t FUNC=100*(x-((0.2-0.1)*y/100))^2 PERIODIC=NO
233 : BIASVALUE ARG=f
234 : ```
235 :
236 : 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.
237 :
238 : ## Using CUSTOM in shortcuts
239 :
240 : 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
241 : expressions that are used in the method are visible in the log. We have found that there are many actions that
242 : 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
243 : fast.
244 :
245 : 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
246 : two other points or the progression along that line. In PLUMED 2.10 we can implement these variables using the following input file:
247 :
248 : ```plumed
249 : # take center of atoms 1 to 10 as reference point 1
250 : p1: CENTER ATOMS=1-10
251 : # take center of atoms 11 to 20 as reference point 2
252 : p2: CENTER ATOMS=11-20
253 : # take center of atoms 21 to 30 as reference point 3
254 : p3: CENTER ATOMS=21-30
255 :
256 : # compute distances
257 : d12: DISTANCE ATOMS=p1,p2
258 : d13: DISTANCE ATOMS=p1,p3
259 : d23: DISTANCE ATOMS=p2,p3
260 :
261 : # compute progress variable of the projection of point p3
262 : # along the vector joining p1 and p2
263 : # notice that progress is measured from the middle point
264 : onaxis: CUSTOM ARG=d13,d23,d12 FUNC=(0.5*(y^2-x^2)/z) PERIODIC=NO
265 :
266 : # compute between point p3 and the vector joining p1 and p2
267 : 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
268 :
269 : PRINT ARG=onaxis,fromaxis
270 : ```
271 :
272 : The equations in this input were also used to combine [RMSD](RMSD.md) values from different snapshots of a protein so as to define
273 : progression (S) and distance (Z) variables in the paper that is cited in the bibliography. We can understand how these expressions
274 : 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
275 : 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
276 : 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)$
277 : by rearranging the cosine rule:
278 :
279 : $$
280 : x\cos(\theta) = \frac{y^2 - x^2}{z} - \frac{z}{2}
281 : $$
282 :
283 : 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
284 : is at the center of the bond connecting atom 1 to atom 2.
285 :
286 : The value `fromaxis` measures the square of the distance from the line. It is calculated using pythagoras theorem as follows:
287 :
288 : $$
289 : f^2 = y^2 - x^2\cos^2(\theta)
290 : $$
291 :
292 : Inserting $x\cos(\theta) = o - \frac{z}{2}$ into this expression gives:
293 :
294 : $$
295 : f^2 = y^2 - (o -\frac{z}{2})^2 = y^2 - o^2 + oz - \frac{z^2}{4}
296 : $$
297 :
298 : 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`.
299 :
300 : ## CUSTOM with vector arguments
301 :
302 : 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
303 : 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
304 : cosine of these three angles and outputs them in a three dimensional vector called `c` that is printed to the colvar file.
305 :
306 : ```plumed
307 : a: ANGLE ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
308 : c: CUSTOM ARG=a FUNC=cos(x) PERIODIC=NO
309 : PRINT ARG=c FILE=colvar
310 : ```
311 :
312 : 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
313 : 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
314 : 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
315 : 8 and 9.
316 :
317 : ```plumed
318 : d: DISTANCE ATOMS1=1,2 ATOMS2=4,5 ATOMS3=7,8
319 : a: ANGLE ATOMS1=1,2,3 ATOMS2=4,5,6 ATOMS3=7,8,9
320 : c: CUSTOM ARG=d,a FUNC=x*cos(y) PERIODIC=NO
321 : PRINT ARG=c FILE=colvar
322 : ```
323 :
324 : 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.
325 :
326 : Lastly, notice that you can pass a mixture of vectors and scalars in the input to a CUSTOM action. However, the labels of the scalar values that appear in the input for the ARG keyword
327 : must appear after the labels of vector vectors as shown below:
328 :
329 : ```plumed
330 : d: DISTANCE ATOMS=1,2
331 : a: ANGLE ATOMS1=1,2,3 ATOMS2=1,2,4 ATOMS3=1,2,5 ATOMS4=1,2,6
332 : c: CUSTOM ARG=a,d FUNC=y*cos(x) PERIODIC=NO
333 : PRINT ARG=c FILE=colvar
334 : ```
335 :
336 : 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).
337 : 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
338 : 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
339 : result of this operation is a vector `f` that contains the three input distances from `d` with 23 added to each of them.
340 :
341 : ```plumed
342 : c: CONSTANT VALUE=23
343 : d: DISTANCE ATOMS1=1,2 ATOMS2=4,5 ATOMS3=7,8
344 : f: CUSTOM ARG=d,c FUNC=x+y PERIODIC=NO
345 : PRINT ARG=f FILE=colvar
346 : ```
347 :
348 : ## CUSTOM, MASK and vector arguments
349 :
350 : The example input below illustrates a way you can use the MASK keyword to improve PLUMED performance for inputs that use CUSTOM actions that take arguments that
351 : are vectors:
352 :
353 : ```plumed
354 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
355 : center: FIXEDATOM AT=2.5,2.5,2.5
356 : ow: GROUP ATOMS=1-400
357 : hw: GROUP ATOMS=401-1200
358 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
359 : sphere: INSPHERE ATOMS=ow CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
360 : # Calculate coordination with oxygen atoms
361 : ow_ones: ONES SIZE=400
362 : ow_mat: CONTACT_MATRIX GROUP=ow MASK=sphere SWITCH={RATIONAL R_0=0.3 D_MAX=0.8}
363 : ow_coord: MATRIX_VECTOR_PRODUCT ARG=ow_mat,ow_ones
364 : # And transform the coordination numbers by a switching function
365 : ow_mt: MORE_THAN ARG=ow_coord MASK=sphere SWITCH={RATIONAL R_0=4}
366 : # Calculate the coordination with hydrogen atoms
367 : hw_ones: ONES SIZE=800
368 : hw_mat: CONTACT_MATRIX GROUPA=ow GROUPB=hw MASK=sphere SWITCH={RATIONAL R_0=0.1 D_MAX=0.5}
369 : hw_coord: MATRIX_VECTOR_PRODUCT ARG=hw_mat,hw_ones
370 : # And transform the coordination numbers by a switching function
371 : hw_mt: MORE_THAN ARG=hw_coord MASK=sphere SWITCH={RATIONAL R_0=2}
372 : # And multiply the two vectors of transformed coordination numbers together
373 : prod: CUSTOM ARG=hw_mt,ow_mt FUNC=x*y MASK=sphere PERIODIC=NO
374 : # Sum of coordination numbers for atoms that are in the sphere of interest
375 : cv: SUM ARG=prod PERIODIC=NO
376 : # And print out final CV to a file
377 : PRINT ARG=cv FILE=colvar STRIDE=1
378 : ```
379 :
380 : The variable being computed here measures the number of oxygen atoms in a spherical region centered on the point $(2.5,2.5,2.5)$ that have a coordination
381 : that is greater than 4 with other oxygen atoms and a coordination that is greater than 2 with other hydrogen atoms. Notice how by using the MASK keyword
382 : in the [CONTACT_MATRIX](CONTACT_MATRIX.md), [MORE_THAN](MORE_THAN.md and CUSTOM actions we avoid calculating the order parameters that are not in the spherical
383 : region of interest in the way that is described in more detail in the documentation for [CONTACT_MATRIX](CONTACT_MATRIX.md).
384 :
385 :
386 : ## CUSTOM with matrix arguments
387 :
388 : 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
389 : 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
390 : can use this functionality to calculate all the angles between a set of bond vectors.
391 :
392 : ```plumed
393 : # Calculate the vectors connecting four atoms
394 : d: DISTANCE COMPONENTS ATOMS1=1,2 ATOMS2=3,4 ATOMS3=5,6 ATOMS4=7,8
395 : # Calculate the norm of these four vectors
396 : norm: CUSTOM ARG=d.x,d.y,d.z FUNC=sqrt(x*x+y*y+z*z) PERIODIC=NO
397 : # Now calculate the directors of the vectors
398 : norm_x: CUSTOM ARG=d.x,norm FUNC=x/y PERIODIC=NO
399 : norm_y: CUSTOM ARG=d.y,norm FUNC=x/y PERIODIC=NO
400 : norm_z: CUSTOM ARG=d.z,norm FUNC=x/y PERIODIC=NO
401 : # And combine all these directors in a matrix
402 : stack: VSTACK ARG=norm_x,norm_y,norm_z
403 : # Now calculate the matrix of dot products between directors (these are the cosines of the angles)
404 : stackT: TRANSPOSE ARG=stack
405 : cosa: MATRIX_PRODUCT ARG=stack,stackT
406 : # And finally get the 4x4 matrix of angles and print it to a file
407 : angles: CUSTOM ARG=cosa FUNC=acos(x) PERIODIC=NO
408 : PRINT ARG=angles FILE=colvar
409 : ```
410 :
411 : Notice that you can pass multiple $N\times M$ matrices in the input to a CUSTOM action as illustrated in the example below:
412 :
413 : ```plumed
414 : c1: CONSTANT VALUES=2,3,4,5 NROWS=2 NCOLS=2
415 : c2: CONSTANT VALUES=1,0,0,1 NROWS=2 NCOLS=2
416 : f: CUSTOM ARG=c1,c2 FUNC=x*y PERIODIC=NO
417 : PRINT ARG=f FILE=colvar
418 : ```
419 :
420 : The four by four matrix that is calculated by the custom action in this input is given by:
421 :
422 : $$
423 : f = \left(
424 : \begin{matrix}
425 : 2 & 0 \\
426 : 0 & 5
427 : \end{matrix}
428 : \right)
429 : $$
430 :
431 : which is the element-wise [Hadamard product](https://en.wikipedia.org/wiki/Hadamard_product_(matrices)) and __not__ the matrix product.
432 : 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
433 : x and a matrix that contains the reciprocals of each of the elements in y is computed. If you wish to calcalculate the
434 : [product of two matrices](https://en.wikipedia.org/wiki/Matrix_multiplication) you should use the [MATRIX_PRODUCT](MATRIX_PRODUCT.md) comamnd.
435 : 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)
436 : command.
437 :
438 : 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
439 : any scalars you pass as being converted into $N\times M$ matrix in which every element is equal to the input scalar. Furthermore, the labels of
440 : the input scalars must appear __after__ the labels fo the input matrices in the input for the ARG keyword.
441 :
442 : ## CUSTOM, MASK and matrix arguments
443 :
444 : The example input below illustrates a way you can use the MASK keyword to improve PLUMED performance for inputs that use CUSTOM actions that take arguments that
445 : are matrices:
446 :
447 : ```plumed
448 : # Calculate the distances between atoms
449 : d_mat: DISTANCE_MATRIX GROUP=1-64 LINKCELL_CUTOFF=4.5 COMPONENTS
450 : # Find the six nearest atom to each of the coordinates
451 : nn: NEIGHBORS ARG=d_mat.w NLOWEST=6
452 : # Evaluate the lengths of the bonds connecting each atom to its nearest neighbors
453 : r: CUSTOM ...
454 : ARG=d_mat.x,d_mat.y,d_mat.z MASK=nn
455 : FUNC=sqrt(x*x+y*y+z*z)
456 : PERIODIC=NO
457 : ...
458 : # Evaluate a function of the nearest neighbors
459 : d_wvfunc: CUSTOM ...
460 : ARG=nn,d_mat.x,d_mat.y,d_mat.z,r MASK=nn
461 : FUNC=w*(((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3)
462 : VAR=w,x,y,z,r
463 : PERIODIC=NO
464 : ...
465 : # Calculate the sum of fcc cubic function values for each atom
466 : d_ones: ONES SIZE=64
467 : d: MATRIX_VECTOR_PRODUCT ARG=d_wvfunc,d_ones
468 : # Calculate the number of neighbours
469 : d_denom: MATRIX_VECTOR_PRODUCT ARG=nn,d_ones
470 : # Calculate the average value of the fcc cubic function per bonds
471 : d_n: CUSTOM ARG=d,d_denom FUNC=x/y PERIODIC=NO
472 : d_mean: MEAN ARG=d_n PERIODIC=NO
473 : PRINT ARG=d_mean FILE=colv
474 : ```
475 :
476 : Using the MASK here ensures that we only calculate the CUSTOM functions for the six nearest neighbors of each of the input atoms. However, it is important to note that
477 : this input is rather unsual. In general you rarely need to use the MASK keyword for CUSTOM actions that take matrices in input. For example, if we use the following
478 : (more conventional version) of the above CV the MASK is unnecessary as the sparsity pattern of the input matrices are all identical.
479 :
480 : ```plumed
481 : cmat: CONTACT_MATRIX GROUP=1-64 SWITCH={RATIONAL R_0=1.0 D_MAX=4.5} COMPONENTS
482 : # Evaluate the lengths of the bonds connecting each atom to its nearest neighbors
483 : r: CUSTOM ...
484 : ARG=cmat.x,cmat.y,cmat.z
485 : FUNC=sqrt(x*x+y*y+z*z)
486 : PERIODIC=NO
487 : ...
488 : # Evaluate a function of the nearest neighbors
489 : d_wvfunc: CUSTOM ...
490 : ARG=cmat.w,cmat.x,cmat.y,cmat.z,r
491 : FUNC=w*(((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3)
492 : VAR=w,x,y,z,r
493 : PERIODIC=NO
494 : ...
495 : # Calculate the sum of fcc cubic function values for each atom
496 : d_ones: ONES SIZE=64
497 : d: MATRIX_VECTOR_PRODUCT ARG=d_wvfunc,d_ones
498 : # Calculate the number of neighbours
499 : d_denom: MATRIX_VECTOR_PRODUCT ARG=cmat.w,d_ones
500 : # Calculate the average value of the fcc cubic function per bonds
501 : d_n: CUSTOM ARG=d,d_denom FUNC=x/y PERIODIC=NO
502 : d_mean: MEAN ARG=d_n PERIODIC=NO
503 : PRINT ARG=d_mean FILE=colv
504 : ```
505 :
506 : ## CUSTOM with grid arguments
507 :
508 : 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
509 : using $F(s) = -k_BT\log[H(s)]$ as illustrated in the input below:
510 :
511 : ```plumed
512 : x: DISTANCE ATOMS=1,2
513 : hA: HISTOGRAM ARG=x GRID_MIN=0.0 GRID_MAX=3.0 GRID_BIN=100 BANDWIDTH=0.1
514 : fE: CUSTOM ARG=hA FUNC=-2.5*log(x) PERIODIC=NO
515 : DUMPGRID ARG=fE FILE=fes.dat
516 : ```
517 :
518 : 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
519 : 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.
520 :
521 : ```plumed
522 : # Evaluate the FCC order parameter for all of the atoms
523 : fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27
524 : # Calculate the distance between each atom and the origin on atom 1
525 : dens_dist: DISTANCES ORIGIN=1 ATOMS=fcc COMPONENTS
526 : # Do a KDE using the FCC order parameters for the weights of each gaussian and the positions of the atoms as the centers
527 : 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
528 : # Estimate the density of atoms at each point in the box
529 : ones: ONES SIZE=5184
530 : 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
531 : # Now calculate the average value of the order parameter at each point in the box
532 : dens: CUSTOM ARG=dens_numer,dens_denom FUNC=x/y PERIODIC=NO
533 : DUMPCUBE ARG=dens FILE=dens.cube FMT=%8.4f STRIDE=1
534 : ```
535 :
536 : ## Making lepton library faster
537 :
538 : If you use a lot of [CUSTOM](CUSTOM.md) actions or CUSTOM [switching functions](LESS_THAN.md) it can slow down PLUMED.
539 : These commands use the lepton library that is included in PLUMED.
540 : This library replaces libmatheval since PLUMED 2.5, and by itself it is significantly faster than libmatheval.
541 : However, you can make it even faster using a [just-in-time compiler](https://github.com/asmjit/asmjit.git).
542 : As of PLUMED 2.6, the correct version of ASMJIT is embedded in PLUMED.
543 : As of PLUMED 2.8, ASMJIT is enabled by default on X86/X64.
544 : Since PLUMED 2.11, it is also enabled by default on AArch64, including Apple Silicon.
545 : It can be disabled at build time with `./configure --disable-asmjit`.
546 : You can disable it at runtime setting the environment variable `PLUMED_USE_ASMJIT`:
547 :
548 : ```bash
549 : export PLUMED_USE_ASMJIT=no
550 : ```
551 :
552 : In some case using a custom expression is almost as fast as using a hard-coded
553 : function. For instance, with an input that contained the following lines:
554 :
555 : ```plumed
556 : c: COORDINATION GROUPA=1-108 GROUPB=1-108 R_0=1
557 : d_fast: COORDINATION GROUPA=1-108 GROUPB=1-108 SWITCH={CUSTOM FUNC=1/(1+x2^3) R_0=1}
558 : ```
559 :
560 : I (GB) obtained the following timings (on a Macbook laptop):
561 :
562 : ````
563 : ...
564 : PLUMED: 4A 1 c 108 0.126592 0.001172 0.000701 0.002532
565 : PLUMED: 4A 2 d_fast 108 0.135210 0.001252 0.000755 0.002623
566 : ...
567 : ````
568 :
569 : Notice the usage of `x2` as a variable for the switching function, which
570 : avoids an unnecessary square root calculation (this is done automatically by the hard-coded switching functions
571 : when you use only even powers). The asmjit calculation (`d_fast`) takes less than 10% more than the hard-coded
572 : one (`c`).
573 :
574 :
575 : */
576 : //+ENDPLUMEDOC
577 :
578 : typedef FunctionShortcut<Custom> CustomShortcut;
579 : PLUMED_REGISTER_ACTION(CustomShortcut,"CUSTOM")
580 : PLUMED_REGISTER_ACTION(CustomShortcut,"MATHEVAL")
581 : typedef FunctionOfScalar<Custom> ScalarCustom;
582 : PLUMED_REGISTER_ACTION(ScalarCustom,"CUSTOM_SCALAR")
583 : PLUMED_REGISTER_ACTION(ScalarCustom,"MATHEVAL_SCALAR")
584 : typedef FunctionOfVector<Custom> VectorCustom;
585 : PLUMED_REGISTER_ACTION(VectorCustom,"CUSTOM_VECTOR")
586 : PLUMED_REGISTER_ACTION(VectorCustom,"MATHEVAL_VECTOR")
587 : typedef FunctionOfMatrix<Custom> MatrixCustom;
588 : PLUMED_REGISTER_ACTION(MatrixCustom,"CUSTOM_MATRIX")
589 : PLUMED_REGISTER_ACTION(MatrixCustom,"MATHEVAL_MATRIX")
590 :
591 : //+PLUMEDOC FUNCTION MATHEVAL
592 : /*
593 : An alias to the CUSTOM function that can also be used to calaculate combinations of variables using a custom expression.
594 :
595 : 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:
596 :
597 : ```plumed
598 : dAB: DISTANCE ATOMS=10,12
599 : dAC: DISTANCE ATOMS=10,15
600 : diff: MATHEVAL ARG=dAB,dAC FUNC=y-x PERIODIC=NO
601 : # notice: the previous line could be replaced with the following
602 : # diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
603 : METAD ARG=diff SIGMA=0.1 HEIGHT=0.5 BIASFACTOR=10 PACE=100
604 : ```
605 :
606 : This alias is kept in order to maintain compatibility with previous PLUMED versions.
607 : However, notice that as of PLUMED 2.5 the libmatheval library is not linked anymore,
608 : and that the MATHEVAL action evaluates functions [the Lepton library](https://simtk.org/projects/lepton).
609 :
610 : */
611 : //+ENDPLUMEDOC
612 :
613 11281 : void Custom::registerKeywords(Keywords& keys) {
614 11281 : if( keys.getDisplayName()=="MATHEVAL") {
615 2692 : keys.setDeprecated("CUSTOM");
616 : }
617 11281 : keys.use("PERIODIC");
618 11281 : keys.add("compulsory","FUNC","the function you wish to evaluate");
619 11281 : 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.");
620 11281 : keys.add("hidden","MASKED_INPUT_ALLOWED","turns on that you are allowed to use masked inputs ");
621 22562 : keys.setValueDescription("scalar/vector/matrix/grid","an arbitrary function");
622 11281 : keys.addDOI("10.1093/nar/gkv872");
623 11281 : }
624 :
625 : /// Collect the names of the variables v such that the expression is of the form v*g(...),
626 : /// i.e. it is zero whenever v is zero (as long as g is finite)
627 6029 : static void findMultiplicativeFactors( const lepton::ExpressionTreeNode& node, std::set<std::string>& factors ) {
628 6029 : const lepton::Operation& op = node.getOperation();
629 6029 : switch( op.getId() ) {
630 2262 : case lepton::Operation::VARIABLE:
631 2262 : factors.insert( op.getName() );
632 2262 : break;
633 603 : case lepton::Operation::MULTIPLY:
634 603 : findMultiplicativeFactors( node.getChildren()[0], factors );
635 603 : findMultiplicativeFactors( node.getChildren()[1], factors );
636 603 : break;
637 987 : case lepton::Operation::DIVIDE:
638 : case lepton::Operation::NEGATE:
639 : case lepton::Operation::SQUARE:
640 : case lepton::Operation::CUBE:
641 987 : findMultiplicativeFactors( node.getChildren()[0], factors );
642 987 : break;
643 562 : case lepton::Operation::MULTIPLY_CONSTANT:
644 562 : if( std::isfinite( dynamic_cast<const lepton::Operation::MultiplyConstant&>(op).getValue() ) ) {
645 562 : findMultiplicativeFactors( node.getChildren()[0], factors );
646 : }
647 : break;
648 0 : case lepton::Operation::POWER_CONSTANT:
649 0 : if( dynamic_cast<const lepton::Operation::PowerConstant&>(op).getValue()>0 ) {
650 0 : findMultiplicativeFactors( node.getChildren()[0], factors );
651 : }
652 : break;
653 : default:
654 : break;
655 : }
656 6029 : }
657 :
658 3274 : void Custom::read( Custom& f, ActionWithArguments* action, FunctionOptions& options ) {
659 : // Read in the variables
660 3274 : unsigned nargs = action->getNumberOfArguments();
661 3274 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(action);
662 3274 : if( av && av->getNumberOfMasks()>0 ) {
663 56 : nargs = nargs - av->getNumberOfMasks();
664 : }
665 3274 : action->parseVector("VAR",f.var);
666 6548 : action->parse("FUNC",f.func);
667 3274 : if(f.var.size()==0) {
668 3152 : f.var.resize(nargs);
669 3152 : if(f.var.size()>3) {
670 0 : action->error("Using more than 3 arguments you should explicitly write their names with VAR");
671 : }
672 3152 : if(f.var.size()>0) {
673 : f.var[0]="x";
674 : }
675 3152 : if(f.var.size()>1) {
676 : f.var[1]="y";
677 : }
678 3152 : if(f.var.size()>2) {
679 : f.var[2]="z";
680 : }
681 : }
682 3274 : if(f.var.size()!=nargs) {
683 0 : action->error("Size of VAR array should be the same as number of arguments");
684 : }
685 3274 : action->log.printf(" with function : %s\n",f.func.c_str());
686 3274 : action->log.printf(" with variables :");
687 8357 : for(unsigned i=0; i<f.var.size(); i++) {
688 5083 : action->log.printf(" %s",f.var[i].c_str());
689 55913 : for(std::map<std::string,double>::const_iterator iter=lepton::Constants().begin(); iter!=lepton::Constants().end(); ++iter ) {
690 50830 : if( f.var[i]==iter->first ) {
691 0 : action->error("variable " + f.var[i] + " clashes with lepton constant " + iter->first );
692 : }
693 : }
694 : }
695 3274 : action->log.printf("\n");
696 3274 : f.function.set( f.func, f.var, action );
697 : // Find the arguments that multiply the whole expression. If one of these
698 : // arguments is zero the function is zero. If, in addition, the derivative of
699 : // that argument is zero whenever its value is zero (e.g. a switching function
700 : // beyond its cutoff) all the derivatives of the function are also zero, so the
701 : // evaluation can be skipped. This is decided on the lepton parse tree: the
702 : // previous string based check got wrong expressions such as x*cos(y)+1 or exp(x)*y.
703 : std::set<std::string> factors;
704 3274 : findMultiplicativeFactors( lepton::Parser::parse(f.func).optimize(lepton::Constants()).getRootNode(), factors );
705 8357 : for(unsigned i=0; i<f.var.size(); ++i) {
706 2190 : if( factors.count(f.var[i])>0 && action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) {
707 1056 : f.check_multiplication_vars.push_back(i);
708 : }
709 : }
710 3274 : if( f.check_multiplication_vars.size()>0 ) {
711 913 : action->log.printf(" optimizing implementation as function only involves multiplication \n");
712 : }
713 :
714 : // If the function is zero when all the arguments are zero and all the arguments have zero
715 : // derivatives when they are zero, value and derivatives are zero and need not be computed.
716 : // If some argument has a non-zero derivative when it is zero the derivatives of the function
717 : // must be computed (e.g. FUNC=x+y), so this shortcut is not used.
718 : bool allderivzero=true;
719 4790 : for(unsigned i=0; i<nargs; ++i) {
720 4122 : if( !action->getPntrToArgument(i)->isDerivativeZeroWhenValueIsZero() ) {
721 : allderivzero=false;
722 : break;
723 : }
724 : }
725 3274 : std::vector<double> zeros( nargs, 0 );
726 3274 : double fval = fabs(f.function.evaluate(zeros));
727 3274 : f.zerowhenallzero=(allderivzero && fval<epsilon );
728 3274 : if( f.zerowhenallzero ) {
729 216 : action->log.printf(" not calculating when all arguments are zero \n");
730 : }
731 3274 : options.derivativeZeroIfValueIsZero = f.check_multiplication_vars.size()>0;
732 3274 : }
733 :
734 4 : std::string Custom::getFunctionString( const std::string& func ) {
735 4 : std::string outstr="", mult="*";
736 20 : for(unsigned i=0; i<func.length(); ++i) {
737 16 : if( func[i]==mult[0] ) {
738 : outstr += "\\*";
739 : } else {
740 : outstr += func[i];
741 : }
742 : }
743 4 : return outstr;
744 : }
745 :
746 34886415 : void Custom::calc( const Custom& func, bool noderiv, const View<const double>& args, FunctionOutput& funcout ) {
747 34886415 : if( args.size()>1 ) {
748 : bool allzero=false;
749 25235932 : if( func.check_multiplication_vars.size()>0 ) {
750 33371481 : for(unsigned i=0; i<func.check_multiplication_vars.size(); ++i) {
751 19004198 : if( fabs(args[func.check_multiplication_vars[i]])<epsilon ) {
752 : allzero=true;
753 : break;
754 : }
755 : }
756 7661877 : } else if( func.zerowhenallzero ) {
757 2519897 : allzero=(fabs(args[0])<epsilon);
758 2703226 : for(unsigned i=1; i<args.size(); ++i) {
759 2616038 : if( fabs(args[i])>epsilon ) {
760 : allzero=false;
761 : break;
762 : }
763 : }
764 : }
765 20093952 : if( allzero ) {
766 3293374 : funcout.values[0]=0;
767 3293374 : if( noderiv ) {
768 : return;
769 : }
770 4545080 : for(unsigned i=0; i<args.size(); i++) {
771 3030056 : funcout.derivs[0][i] = 0.0;
772 : }
773 : return;
774 : }
775 : }
776 31593041 : funcout.values[0] = func.function.evaluate( args );
777 31593041 : if( !noderiv ) {
778 15549522 : for(unsigned i=0; i<args.size(); i++) {
779 10205526 : funcout.derivs[0][i] = func.function.evaluateDeriv( i, args );
780 : }
781 : }
782 : }
783 :
784 : }
785 : }
786 :
787 :
|
