Skip to content

Shortcut: CUSTOM

Module function
Description Usage
Calculate a combination of variables using a custom expression. used in 11 tutorialsused in 54 eggs
output value type
an arbitrary function scalar/vector/matrix/grid

Details and examples

Calculate a combination of variables using a custom expression.

CUSTOM is one of the most useful actions in PLUMED. This action takes in a list of arguments and then uses the lepton mathematical expression parser to evaluate a user defined function of these input arguments. We can thus use this action in the input below to perform a metadynamics simulation using the difference between two distances as a CV.

Click on the labels of the actions for more information on what each action computes
tested on2.11
dAB: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=10,12
dAC: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=10,15
diff: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=dAB,dAC FUNCthe function you wish to evaluate=y-x PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# notice: the previous line could be replaced with the following
# diff: COMBINE ARG=dAB,dAC COEFFICIENTS=-1,1
METADUsed to performed metadynamics on one or more collective variables. More details ARGthe labels of the scalars on which the bias will act=diff SIGMAthe widths of the Gaussian hills=0.1 HEIGHTthe heights of the Gaussian hills=0.5 BIASFACTORuse well tempered metadynamics and use this bias factor=10 PACEthe frequency for hill addition=100

The particular function that should be evaluated from the input arguments is specified using the FUNC keyword. The function provided to the FUNC keyword is written in terms of x and y in the input above. x here deonetes the first argument provided to the ARG keyword, dAB, while y is the second argument provided to the ARG keyword, dAC.

The VAR keyword

If you wish, you can rewrite the example input above more transparantly as:

Click on the labels of the actions for more information on what each action computes
tested on2.11
dAB: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=10,12
dAC: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=10,15
diff: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=dAB,dAC VARthe names to give each of the arguments in the function=dAB,dAC FUNCthe function you wish to evaluate=dAC-dAB PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
METADUsed to performed metadynamics on one or more collective variables. More details ARGthe labels of the scalars on which the bias will act=diff SIGMAthe widths of the Gaussian hills=0.1 HEIGHTthe heights of the Gaussian hills=0.5 BIASFACTORuse well tempered metadynamics and use this bias factor=10 PACEthe frequency for hill addition=100

By using the VAR keyword here we ensure that we can use the labels of the arguments in the mathematical expression that we provide to the FUNC argument. Notice that, if you have four or more arguments you must use the VAR keyword. With three or less arguments the VAR keyword can be ommitted and the symbols x, y and z can be used to denote the first, second and third arguments respectively.

The following input illustrates a case where using the VAR keyword is essential. This input tells PLUMED to print the angle between the vector connecting atoms 1,2 and the vector connecting atoms 2,3.

Click on the labels of the actions for more information on what each action computes
tested on2.11
d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2 COMPONENTS calculate the x, y and z components of the distance separately and store them as label
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=2,3 COMPONENTS calculate the x, y and z components of the distance separately and store them as label
theta: CUSTOMCalculate a combination of variables using a custom expression. More details ...
  ARGthe values input to this function=d1.x,d1.y,d1.z,d2.x,d2.y,d2.z
  VARthe names to give each of the arguments in the function=ax,ay,az,bx,by,bz
  FUNCthe function you wish to evaluate=acos((ax*bx+ay*by+az*bz)/sqrt((ax*ax+ay*ay+az*az)*(bx*bx+by*by+bz*bz)))
  PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
...

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=theta

Functions and constants

The previous examples demonstrates how CUSTOM can be used to evaluate mathematical expressions that involve taking powers (^), adding (+), multiplying (*), dividing (/) and subtracting (-) variables and can control the order in which operations are performed by using parenthesis. 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 that you can use in in the input to the FUNC keyword for CUSTOM.

Function Description
sqrt(x) The square root of x
exp(x) The exponential of x i.e.
log(x) The natural logarithm of x
sin(x) The sine of x
cos(x) The cosine of x
sec(x) The reciprocal of the cosine of .
csc(x) The reciprocal of the sine of .
tan(x) The tangent of x i.e.
cot(x) The reciprocal of the tangent of .
asin(x) The principal arc sine of x. Returns which gives
acos(x) The principal arc cosine of x. Returns which gives
atan(x) The principal arc tangent of x. Returns which gives
atan2(x,y) The principal arg tangent of . Returns which gives
sinh(x) The hyperbolic sine of
cosh(x) the hyperbolic cosine of
tanh(x) The hyperbolic tangent of
erf(x) The error function
erfc(x) The complementary error function
step(x) 1 if and 0 otherwise
delta(x) inf if and 0 otherwise
nandelta(x) nan if and 0 otherwise
square(x) The square of x i.e.
cube(x) The cube of x i.e.
recip(x) The reciprocal of x i.e.
min(x,y) If this function returns . If this function returns
max(x,y) If this function returns . If this function returns
abs(x) The absolute value of x
floor(x) The largest integer that is less than
ceil(x) The smallest integer that is greater than
select(x,y,z) If this returns otherwise this returns
acot(x) Returns the value of which gives
asec(x) Returns the value of which gives
acsc(x) Returns the value of which gives
coth(x) The recipricoal of the hyperbolic tangent of .
sech(x) The recipricoal of the hyperbolic cosine of .
csch(x) The recipricoal of the hyperbolic sine of .
asinh(x) The nonnegative area hyperbolic sine of . Returns which gives
acosh(x) The nonnegative area hyperbolic cosine of . Returns which gives
atanh(x) The nonnegative area hyperbolic tangent of . Returns which gives .
acoth(x) The inverse hyperbolic tangent of is calculated as
asech(x) The inverse hyperbolic secant of is calculated as
acsch(x) The inverse hyperbolic cosecant of is calculated as

Notice, that you can also incorporate the following constants in the expressions that are used in the input to FUNC:

Symbol Description
e Euler's number - the base for the natural logarithm
log2e
log10e
ln2
ln10
pi the circle constant
pi_2
pi_4
`sqrt2
sqrt1_2

The way one of these constants can be used in an expression is illustrated in the following example:

Click on the labels of the actions for more information on what each action computes
tested on2.11
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
# This function is evaluating the area of the circle whose radius is
# given by the distance between atoms 1 and 2.
c: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d FUNCthe function you wish to evaluate=pi*x*x PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=c FILEthe name of the file on which to output these quantities=colvar

The PERIODIC keyword

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. To see how to use this keyword consider the following input:

Click on the labels of the actions for more information on what each action computes
tested on2.11
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2 COMPONENTS calculate the x, y and z components of the distance separately and store them as label
theta: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d.x,d.y FUNCthe function you wish to evaluate=atan2(y,x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=-pi,+pi
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=theta FILEthe name of the file on which to output these quantities=colvar

The above input calculates the torsional angle between the vector connecting atoms 1 and 2 and positive direction around the axis. This quantity must lie between and so we use the PERIODIC keyword to tell PLUMED this fact about the output quantity.

Notice that PLUMED cannot automatically detect if the codomain of the function is periodic. You thus always have to specifiy the periodicity of the codomain whenever you use the CUSTOM command. However, the codomain of most functions that people use are not periodic so you can normally just write PERIODIC=NO

The step function

The step operation (that is the Heaviside function) from the table above allow you to use if clauses in CUSTOM actions. 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 1 if x is less than one and zero otherwise.

Using the step operation multiple times in a function allows you to perform logical operations. For example, the equivalent of the AND operator 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. 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 greater than 1.0.

The example below illustrtes how we can put these ideas of using the step operation into practise.

Click on the labels of the actions for more information on what each action computes
tested on2.11
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=10,15
m: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d FUNCthe function you wish to evaluate=0.5*step(0.5-x)+x*step(x-0.5) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# check the function you are applying:
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d,m FILEthe name of the file on which to output these quantities=checkme
RESTRAINTAdds harmonic and/or linear restraints on one or more variables. More details ARGthe values the harmonic restraint acts upon=d ATthe position of the restraint=0.5 KAPPA specifies that the restraint is harmonic and what the values of the force constants on each of the variables are=10.0

The meaning of the function 0.5*step(0.5-x)+x*step(x-0.5) in this example is: - If x<0.5 (step(0.5-x)!=0) use 0.5 - If x>0.5 (step(x-0.5)!=0) use x Notice that the same result can be achieved by using UPPER_WALLS However, with CUSTOM you can create much more complex definitions.

Notice that we can apply a force on the value m by using the RESTRAINT command as the function we defined in the expression that was passed to the FUNC keyword is continuous. In general, however, you must be careful when using the step, delta, nandelta and select functions as you can easily write expression for discontinuous functions by using these operations. If you want to check if a function you have created using step is continuous you can easily plot the function in gnuplot by using a commands like those shown below

# this allow to step function to be used in gnuplot:
gnuplot> step(x)=0.5*(erf(x*10000000)+1)
# here you can test your function
gnuplot> p 0.5*step(0.5-x)+x*step(x-0.5)

Using TIME as a variable

Notice that you can use CUSTOM to implement a MOVINGRESTRAINT as shown below.

Click on the labels of the actions for more information on what each action computes
tested on2.11
t: TIMEretrieve the time of the simulation to be used elsewhere More details
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d,t FUNCthe function you wish to evaluate=100*(x-((0.2-0.1)*y/100))^2 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
BIASVALUETakes the value of one variable and use it as a bias More details ARGthe labels of the scalar/vector arguments whose values will be used as a bias on the system=f

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.

Using CUSTOM in shortcuts

The CUSTOM action is used in many of the shortcut actions that are implemented in PLUMED. We think that using this action in these places is beneficial as it ensures that the mathematical expressions that are used in the method are visible in the log. We have found that there are many actions that 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 fast.

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 two other points or the progression along that line. In PLUMED 2.10 we can implement these variables using the following input file:

Click on the labels of the actions for more information on what each action computes
tested on2.11
# take center of atoms 1 to 10 as reference point 1
p1: CENTERCalculate the center for a group of atoms, with arbitrary weights. More details ATOMSthe group of atoms that appear in the definition of this center=1-10
# take center of atoms 11 to 20 as reference point 2
p2: CENTERCalculate the center for a group of atoms, with arbitrary weights. More details ATOMSthe group of atoms that appear in the definition of this center=11-20
# take center of atoms 21 to 30 as reference point 3
p3: CENTERCalculate the center for a group of atoms, with arbitrary weights. More details ATOMSthe group of atoms that appear in the definition of this center=21-30

# compute distances
d12: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=p1,p2
d13: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=p1,p3
d23: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=p2,p3

# compute progress variable of the projection of point p3
# along the vector joining p1 and p2
# notice that progress is measured from the middle point
onaxis: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d13,d23,d12 FUNCthe function you wish to evaluate=(0.5*(y^2-x^2)/z) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO

# compute between point p3 and the vector joining p1 and p2
fromaxis: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d13,d23,d12,onaxis VARthe names to give each of the arguments in the function=x,y,z,o FUNCthe function you wish to evaluate=(0.5*(y^2+x^2)-o^2-0.25*z^2) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=onaxis,fromaxis

The equations in this input were also used to combine RMSD values from different snapshots of a protein so as to define progression (S) and distance (Z) variables in the paper that is cited in the bibliography. We can understand how these expressions are derived by noting that , and are the distances between atoms 1 and 3, 2 and 3 and 1 and 2 respectively. The projection of the vector connecting atom 1 to atom 3 onto the vector connecting atom 1 to atom 2 is thus , where theta is the angle between the vector connecting atoms 1 and 3 and the vector connecting atoms 1 and 2. We can arrive at the following expression for by rearranging the cosine rule:

Notice that the value called onaxis in the above input is thus . Adding the factor of ensures that the origin is at the center of the bond connecting atom 1 to atom 2.

The value fromaxis measures the square of the distance from the the line. It is calculated using pythagoras theorem as follows:

Inserting into this expression gives:

Inserting the fact that , which comes from the expression for that was used to calculate onaxis, gets us to the expression that is used to calculate fromaxis.

CUSTOM with vector arguments

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 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 cosine of these three angles and outputs them in a three dimensional vector called c that is printed to the colvar file.

Click on the labels of the actions for more information on what each action computes
tested on2.11
a: ANGLECalculate one or multiple angle/s. More details ATOMS1the list of atoms involved in this collective variable (either 3 or 4 atoms)=1,2,3 ATOMS2the list of atoms involved in this collective variable (either 3 or 4 atoms)=4,5,6 ATOMS3the list of atoms involved in this collective variable (either 3 or 4 atoms)=7,8,9
c: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=a FUNCthe function you wish to evaluate=cos(x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=c FILEthe name of the file on which to output these quantities=colvar

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 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 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 8 and 9.

Click on the labels of the actions for more information on what each action computes
tested on2.11
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMS1the pair of atom that we are calculating the distance between=1,2 ATOMS2the pair of atom that we are calculating the distance between=4,5 ATOMS3the pair of atom that we are calculating the distance between=7,8
a: ANGLECalculate one or multiple angle/s. More details ATOMS1the list of atoms involved in this collective variable (either 3 or 4 atoms)=1,2,3 ATOMS2the list of atoms involved in this collective variable (either 3 or 4 atoms)=4,5,6 ATOMS3the list of atoms involved in this collective variable (either 3 or 4 atoms)=7,8,9
c: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d,a FUNCthe function you wish to evaluate=x*cos(y) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=c FILEthe name of the file on which to output these quantities=colvar

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.

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 must appear after the labels of vector vectors as shown below:

Click on the labels of the actions for more information on what each action computes
tested on2.11
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
a: ANGLECalculate one or multiple angle/s. More details ATOMS1the list of atoms involved in this collective variable (either 3 or 4 atoms)=1,2,3 ATOMS2the list of atoms involved in this collective variable (either 3 or 4 atoms)=1,2,4 ATOMS3the list of atoms involved in this collective variable (either 3 or 4 atoms)=1,2,5 ATOMS4the list of atoms involved in this collective variable (either 3 or 4 atoms)=1,2,6
c: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=a,d FUNCthe function you wish to evaluate=y*cos(x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=c FILEthe name of the file on which to output these quantities=colvar

The multiplication of a scalar by a vector in the above input is done in the usual way. 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 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 result of this operation is a vector f that contains the three input distances from d with 23 added to each of them.

Click on the labels of the actions for more information on what each action computes
tested on2.11
c: CONSTANTCreate a constant value that can be passed to actions This action has hidden defaults. More details VALUEthe single number that you would like to store=23
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMS1the pair of atom that we are calculating the distance between=1,2 ATOMS2the pair of atom that we are calculating the distance between=4,5 ATOMS3the pair of atom that we are calculating the distance between=7,8
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d,c FUNCthe function you wish to evaluate=x+y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=f FILEthe name of the file on which to output these quantities=colvar

CUSTOM, MASK and vector arguments

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 are vectors:

Click on the labels of the actions for more information on what each action computes
tested on2.11
# Fixed virtual atom which serves as the probe volume's center (pos. in nm)
center: FIXEDATOMAdd a virtual atom in a fixed position. This action has hidden defaults. More details ATcoordinates of the virtual atom=2.5,2.5,2.5
ow: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=1-400
hw: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=401-1200
# Vector in which element i is one if atom i is in sphere of interest and zero otherwise
sphere: INSPHEREThis quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a particular, user-specified part of of the cell. More details ATOMSthe group of atoms that you would like to investigate=ow CENTERthe atom whose vicinity we are interested in examining=center RADIUSthe switching function that tells us the extent of the sphereical region of interest. Options for this keyword are explained in the documentation for LESS_THAN.={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
# Calculate coordination with oxygen atoms
ow_ones: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=400
ow_mat: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details GROUPspecifies the list of atoms that should be assumed indistinguishable=ow MASKa vector that is used to used to determine which rows of the adjancency matrix to compute=sphere SWITCHthe input for the switching function that acts upon the distance between each pair of atoms. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.3 D_MAX=0.8}
ow_coord: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=ow_mat,ow_ones
# And transform the coordination numbers by a switching function
ow_mt: MORE_THANUse a switching function to determine how many of the input variables are more than a certain cutoff. More details ARGthe values input to this function=ow_coord MASKthe label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed=sphere SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL R_0=4}
# Calculate the coordination with hydrogen atoms
hw_ones: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=800
hw_mat: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details GROUPAwhen you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB=ow GROUPBwhen you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA=hw MASKa vector that is used to used to determine which rows of the adjancency matrix to compute=sphere SWITCHthe input for the switching function that acts upon the distance between each pair of atoms. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
hw_coord: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=hw_mat,hw_ones
# And transform the coordination numbers by a switching function
hw_mt: MORE_THANUse a switching function to determine how many of the input variables are more than a certain cutoff. More details ARGthe values input to this function=hw_coord MASKthe label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed=sphere SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL R_0=2}
# And multiply the two vectors of transformed coordination numbers together
prod: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=hw_mt,ow_mt FUNCthe function you wish to evaluate=x*y MASKthe label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed=sphere PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Sum of coordination numbers for atoms that are in the sphere of interest
cv: SUMCalculate the sum of the arguments More details ARGthe vector/matrix/grid whose elements shuld be added together=prod PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And print out final CV to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=cv FILEthe name of the file on which to output these quantities=colvar STRIDE the frequency with which the quantities of interest should be output=1

The variable being computed here measures the number of oxygen atoms in a spherical region centered on the point that have a coordination 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 in the CONTACT_MATRIX, [MORE_THAN](MORE_THAN.md and CUSTOM actions we avoid calculating the order parameters that are not in the spherical region of interest in the way that is described in more detail in the documentation for CONTACT_MATRIX.

CUSTOM with matrix arguments

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 is applied to each element of the matrix in turn so if the input matrix is the output matrix is also . The following example illustrates how you can use this functionality to calculate all the angles between a set of bond vectors.

Click on the labels of the actions for more information on what each action computes
tested on2.11
# Calculate the vectors connecting four atoms
d: DISTANCECalculate the distance/s between pairs of atoms. More details COMPONENTS calculate the x, y and z components of the distance separately and store them as label ATOMS1the pair of atom that we are calculating the distance between=1,2 ATOMS2the pair of atom that we are calculating the distance between=3,4 ATOMS3the pair of atom that we are calculating the distance between=5,6 ATOMS4the pair of atom that we are calculating the distance between=7,8
# Calculate the norm of these four vectors
norm: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d.x,d.y,d.z FUNCthe function you wish to evaluate=sqrt(x*x+y*y+z*z) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Now calculate the directors of the vectors
norm_x: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d.x,norm FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
norm_y: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d.y,norm FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
norm_z: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d.z,norm FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And combine all these directors in a matrix
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=norm_x,norm_y,norm_z
# Now calculate the matrix of dot products between directors (these are the cosines of the angles)
stackT: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=stack
cosa: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=stack,stackT
# And finally get the 4x4 matrix of angles and print it to a file
angles: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cosa FUNCthe function you wish to evaluate=acos(x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=angles FILEthe name of the file on which to output these quantities=colvar

Notice that you can pass multiple matrices in the input to a CUSTOM action as illustrated in the example below:

Click on the labels of the actions for more information on what each action computes
tested on2.11
c1: CONSTANTCreate a constant value that can be passed to actions More details VALUESthe numbers that are in your constant value=2,3,4,5 NROWS the number of rows in your input matrix=2 NCOLS the number of columns in your matrix=2
c2: CONSTANTCreate a constant value that can be passed to actions More details VALUESthe numbers that are in your constant value=1,0,0,1 NROWS the number of rows in your input matrix=2 NCOLS the number of columns in your matrix=2
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=c1,c2 FUNCthe function you wish to evaluate=x*y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=f FILEthe name of the file on which to output these quantities=colvar

The four by four matrix that is calculated by the custom action in this input is given by:

which is the element-wise Hadamard product and not the matrix product. 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 x and a matrix that contains the reciprocals of each of the elements in y is computed. If you wish to calcalculate the product of two matrices you should use the MATRIX_PRODUCT comamnd. Similarly, if you want to calculate the product of a matrix and a vector you should use the MATRIX_VECTOR_PRODUCT command.

Lastly, note that you can pass a mixture of scalars and matrices in the input to a CUSTOM command. As with vectors, you can think of any scalars you pass as being converted into matrix in which every element is equal to the input scalar. Furthermore, the labels of the input scalars must appear after the labels fo the input matrices in the input for the ARG keyword.

CUSTOM, MASK and matrix arguments

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 are matrices:

Click on the labels of the actions for more information on what each action computes
tested on2.11
# Calculate the distances between atoms
d_mat: DISTANCE_MATRIXCalculate a matrix of distances between atoms. This action has hidden defaults. More details GROUPthe atoms for which you would like to calculate the adjacency matrix=1-64 CUTOFF ignore distances that have a value larger than this cutoff=4.5 COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix
# Find the six nearest atom to each of the coordinates
nn: NEIGHBORSBuild a matrix with ones in for the N nearest neighbours of an atom This action has hidden defaults. More details ARGthe label of an adjacency/distance matrix that will be used to find the nearest neighbors=d_mat.w NLOWEST in each row of the output matrix set the elements that correspond to the n lowest elements in each row of the input matrix equal to one=6
# Evalulate the lengths of the bonds connecting each atom to its nearest neighbors
r: CUSTOMCalculate a combination of variables using a custom expression. More details ...
   ARGthe values input to this function=d_mat.x,d_mat.y,d_mat.z MASKthe label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed=nn
   FUNCthe function you wish to evaluate=sqrt(x*x+y*y+z*z)
   PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
...
# Evaluate a function of the nearest neighbors
d_wvfunc: CUSTOMCalculate a combination of variables using a custom expression. More details ...
   ARGthe values input to this function=nn,d_mat.x,d_mat.y,d_mat.z,r MASKthe label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed=nn
   FUNCthe function you wish to evaluate=w*(((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3)
   VARthe names to give each of the arguments in the function=w,x,y,z,r
   PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
...
# Calculate the sum of fcc cubic function values for each atom
d_ones: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=64
d: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=d_wvfunc,d_ones
# Calculate the number of neighbours
d_denom: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=nn,d_ones
# Calculate the average value of the fcc cubic function per bonds
d_n: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d,d_denom FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
d_mean: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe vector/matrix/grid whose elements shuld be added together=d_n PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d_mean FILEthe name of the file on which to output these quantities=colv

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 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 (more conventional version) of the above CV the MASK is unecessary as the sparsity pattern of the input matrices are all identical.

Click on the labels of the actions for more information on what each action computes
tested on2.11
cmat: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details GROUPspecifies the list of atoms that should be assumed indistinguishable=1-64 SWITCHthe input for the switching function that acts upon the distance between each pair of atoms. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=1.0 D_MAX=4.5} COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix
# Evalulate the lengths of the bonds connecting each atom to its nearest neighbors
r: CUSTOMCalculate a combination of variables using a custom expression. More details ...
   ARGthe values input to this function=cmat.x,cmat.y,cmat.z
   FUNCthe function you wish to evaluate=sqrt(x*x+y*y+z*z)
   PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
...
# Evaluate a function of the nearest neighbors
d_wvfunc: CUSTOMCalculate a combination of variables using a custom expression. More details ...
   ARGthe values input to this function=cmat.w,cmat.x,cmat.y,cmat.z,r
   FUNCthe function you wish to evaluate=w*(((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3)
   VARthe names to give each of the arguments in the function=w,x,y,z,r
   PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
...
# Calculate the sum of fcc cubic function values for each atom
d_ones: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=64
d: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=d_wvfunc,d_ones
# Calculate the number of neighbours
d_denom: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=cmat.w,d_ones
# Calculate the average value of the fcc cubic function per bonds
d_n: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d,d_denom FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
d_mean: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe vector/matrix/grid whose elements shuld be added together=d_n PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d_mean FILEthe name of the file on which to output these quantities=colv

CUSTOM with grid arguments

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 using as illustrated in the input below:

Click on the labels of the actions for more information on what each action computes
tested on2.11
x: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
hA: HISTOGRAMAccumulate the average probability density along a few CVs from a trajectory. This action is a shortcut and it has hidden defaults. More details ARGthe quantities that are being used to construct the histogram=x GRID_MIN the lower bounds for the grid=0.0 GRID_MAX the upper bounds for the grid=3.0 GRID_BINthe number of bins for the grid=100 BANDWIDTHthe bandwidths for kernel density esimtation=0.1
fE: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=hA FUNCthe function you wish to evaluate=-2.5*log(x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
DUMPGRIDOutput the function on the grid to a file with the PLUMED grid format. More details ARGthe label for the grid that you would like to output=fE FILE the file on which to write the grid=fes.dat

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 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.

Click on the labels of the actions for more information on what each action computes
tested on2.11
# Evaluate the FCC order parameter for all of the atoms
fcc: FCCUBICMeasure how similar the environment around atoms is to that found in a FCC structure. This action is a shortcut and it has hidden defaults. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-5184 SWITCHthe switching function that it used in the construction of the contact matrix. Options for this keyword are explained in the documentation for LESS_THAN.={CUBIC D_0=1.2 D_MAX=1.5} ALPHA The alpha parameter of the angular function that is used for FCCUBIC=27
# Calculate the distance between each atom and the origin on atom 1
dens_dist: DISTANCESCalculate the distances between multiple piars of atoms This action is a shortcut. More details ORIGINcalculate the distance of all the atoms specified using the ATOMS keyword from this point=1 ATOMSthe pairs of atoms that you would like to calculate the angles for=fcc COMPONENTS calculate the x, y and z components of the distance separately and store them as label
# Do a KDE using the FCC order parameters for the weights of each gaussian and the positions of the atoms as the centers
dens_numer: KDECreate a histogram from the input scalar/vector/matrix using KDE This action has hidden defaults. More details HEIGHTSthis keyword takes the label of an action that calculates a vector of values=fcc_n ARGthe label for the value that should be used to construct the histogram=dens_dist.x,dens_dist.y,dens_dist.z GRID_BINthe number of bins for the grid=14,14,28 BANDWIDTHthe bandwidths for kernel density esimtation=1.0,1.0,1.0
# Estimate the density of atoms at each point in the box
ones: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=5184
dens_denom: KDECreate a histogram from the input scalar/vector/matrix using KDE This action has hidden defaults. More details ARGthe label for the value that should be used to construct the histogram=dens_dist.x,dens_dist.y,dens_dist.z GRID_BINthe number of bins for the grid=14,14,28 HEIGHTSthis keyword takes the label of an action that calculates a vector of values=ones BANDWIDTHthe bandwidths for kernel density esimtation=1.0,1.0,1.0
# Now calculate the average value of the order parameter at each point in the box
dens: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=dens_numer,dens_denom FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
DUMPCUBEOutput a three dimensional grid using the Gaussian cube file format. More details ARGthe label for the grid that you would like to output=dens FILE the file on which to write the grid=dens.cube FMTthe format that should be used to output real numbers=%8.4f

Making lepton library faster

If you use a lot of CUSTOM actions or CUSTOM switching functions it can slow down PLUMED. These commands use the lepton library that is included in PLUMED. This library replaces libmatheval since PLUMED 2.5, and by itself it is significantly faster than libmatheval. However, you can make it even faster using a just-in-time compiler. As of PLUMED 2.6, the correct version of ASMJIT is embedded in PLUMED. As of PLUMED 2.8, ASMJIT is enabled by default on supported architectures (X86/X64). You can disable it at runtime setting the environment variable PLUMED_USE_ASMJIT:

export PLUMED_USE_ASMJIT=no

In some case using a custom expression is almost as fast as using a hard-coded function. For instance, with an input that contained the following lines:

Click on the labels of the actions for more information on what each action computes
tested on2.11
c: COORDINATIONCalculate coordination numbers. This action has hidden defaults. More details GROUPAFirst list of atoms=1-108 GROUPBSecond list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)=1-108 R_0The r_0 parameter of the switching function=1
d_fast: COORDINATIONCalculate coordination numbers. More details GROUPAFirst list of atoms=1-108 GROUPBSecond list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)=1-108 SWITCHThis keyword is used if you want to employ an alternative to the continuous switching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={CUSTOM FUNC=1/(1+x2^3) R_0=1}

I (GB) obtained the following timings (on a Macbook laptop):

...
PLUMED: 4A  1 c                                          108     0.126592     0.001172     0.000701     0.002532
PLUMED: 4A  2 d_fast                                      108     0.135210     0.001252     0.000755     0.002623
...

Notice the usage of x2 as a variable for the switching function, which avoids an unnecessary square root calculation (this is done automatically by the hard-coded switching functions when you use only even powers). The asmjit calculation (d_fast) takes less than 10% more than the hard-coded one (c).

Input

The arguments that serve as the input for this action are specified using one or more of the keywords in the following table.

Keyword Type Description
ARG scalar/vector/matrix/grid the values input to this function
MASK vector/matrix the label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed

Full list of keywords

The following table describes the keywords and options that can be used with this action

Keyword Type Default Description
ARG input none the values input to this function
MASK input none the label for a sparse vector/matrix that should be used to determine which elements of the vector/matrix should be computed
PERIODIC compulsory none if the output of your function is periodic then you should specify the periodicity of the function
FUNC compulsory none the function you wish to evaluate
VAR optional not used the names to give each of the arguments in the function

References

More information about how this action can be used is available in the following articles: