Calculating symmetry functions

In the other blog articles on these pages I have written a great deal about how you can calculate coordination numbers using inputs such as the one below:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix.  This is a square 7x7 matrix
c1: 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-7 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=2.6 NN=6 MM=12}
# Calculate the coordination numbers for the 7 atoms by a vector that contains all ones
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=7
cc: 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=c1,ones
# And output the 7 coordination numbers for the atoms
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=cc FILEthe name of the file on which to output these quantities=colvar

In this article I want to explain how we can use PLUMED input files that are similar to the one above to calculate other input files.

CONTACT_MATRIX and the COMPONENTS flag

In the input at the top of this file I used a CONTACT_MATRIX. This command outputs a matrix, $\mathbf{A}$, in which element $i,j$ is given by:

\[A_{ij} = \sigma(r_{ij})\]

where $\sigma$ is a switching function and $r_{ij}$ is the distance between atom $i$ and atom $j$.

Switching functions like $\sigma(r_{ij})$ are used in the definitions symmetry functions such as those used in the FCCUBIC, TETRAHEDRAL and SIMPLECUBIC actions. These functions have the following general form:

\[s_i = \sum_{j\ne i } f(x_{ij},y_{ij},z_{ij})\sigma(r_{ij})\]

where $f$ is some function of the $x_{ij}$, $y_{ij}$ and $z_{ij}$ components of the vector that connects atom $i$ to atom $j$. To make it easy to implement such functions I have introduced a COMPONENTS flag on the CONTACT_MATRIX command that can be used as shown below:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix and the three matrices that contain the components of the 
# vectors connecting atom i and atom j.  In other words, calculate four 7x7 matrices
c1: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix GROUPspecifies the list of atoms that should be assumed indistinguishable=1-7 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=2.6 NN=6 MM=12}
# Output the four 7x7 matrices calculated by the command above to a file.
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=c1.w,c1.x,c1.y,c1.z FILEthe name of the file on which to output these quantities=colvar

The COMPONENTS flag for the CONTACT_MATRIX action is specfically designed to make calculating symmetry functions such as $s_i$ striaghforward. The elements of the matrix c1.x are thus calculated as:

\[X_{ij} = \begin{cases} x_{ij} & \textrm{if} \quad \sigma(r_{ij})>0 \\ 0 & \textrm{otherwise} \end{cases}\]

The elements of c1.y and c1.z are also calculated in a similar fashion. In other words, the $i,j$ element of these three matrices are only non-zero if atom $j$ is within the first coordination sphere of atom $i$.

Combining CONTACT_MATRIX and CUSTOM

We can use the CONTACT_MATRIX command with the COMPONENTS keyword to calculate a function similar to $s_i$ above as follows:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix and the three matrices that contain the components of the 
# vectors connecting atom i and atom j.  In other words, calculate four 7x7 matrices
c1: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix GROUPspecifies the list of atoms that should be assumed indistinguishable=1-7 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=2.6 NN=6 MM=12}
# This applies the function to the three input matrices element-wise and thus outputs 
# a 7x7 matrix
r: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=c1.x,c1.y,c1.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
# Again we are applying the function element wise to the five input matrices and thus 
# outputting a 7x7 matrix
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=c1.w,c1.x,c1.y,c1.z,r FUNCthe function you wish to evaluate=w*(x^4+y^4+z^4)/(r^4) 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
# We now multiply the 7x7 by a vector of ones to calculate the symmetry function
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=7
s: 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=f,ones 
# And print the sum of the symmetry function to a file
sums: SUMCalculate the sum of the arguments More details ARGthe vector/matrix/grid whose elements shuld be added together=s 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=sums FILEthe name of the file on which to output these quantities=colvar

The graph for this input is shown below:

flowchart TB 
MD(positions from MD)
Box("label=Box 
 PBC 
")
Box -- Box --> c1
linkStyle 0 stroke:red,color:red;
MD --> c1
linkStyle 1 stroke:violet,color:violet;
c1(["label=c1 
 CONTACT_MATRIX 
"])
c1 -- c1.x --> r
linkStyle 2 stroke:red,color:red;
c1 -- c1.y --> r
linkStyle 3 stroke:red,color:red;
c1 -- c1.z --> r
linkStyle 4 stroke:red,color:red;
r(["label=r 
 CUSTOM
FUNC=sqrt(x\*x+y\*y+z\*z) 
"])
c1 -- c1.w --> f
linkStyle 5 stroke:red,color:red;
c1 -- c1.x --> f
linkStyle 6 stroke:red,color:red;
c1 -- c1.y --> f
linkStyle 7 stroke:red,color:red;
c1 -- c1.z --> f
linkStyle 8 stroke:red,color:red;
r -- r --> f
linkStyle 9 stroke:red,color:red;
f(["label=f 
 CUSTOM
FUNC=w\*(x^4+y^4+z^4)/(r^4) 
"])
ones(["label=ones 
 CONSTANT 
"])
f -- f --> s
linkStyle 10 stroke:red,color:red;
ones -- ones --> s
linkStyle 11 stroke:blue,color:blue;
s(["label=s 
 MATRIX_VECTOR_PRODUCT 
"])
s -- s --> sums
linkStyle 12 stroke:blue,color:blue;
sums(["label=sums 
 SUM 
"])
sums -- sums --> 11
11("label=#64;11 
 PRINT
FILE=colvar 
")

The way forces are passed through the code is illustrated below.

flowchart TB 
MD(positions from MD)
Box("label=Box 
 PBC 
")
Box -- Box --> c1
linkStyle 0 stroke:red,color:red;
MD --> c1
linkStyle 1 stroke:violet,color:violet;
c1(["label=c1 
 CONTACT_MATRIX 
"])
c1 -- c1.x --> r
linkStyle 2 stroke:red,color:red;
c1 -- c1.y --> r
linkStyle 3 stroke:red,color:red;
c1 -- c1.z --> r
linkStyle 4 stroke:red,color:red;
r(["label=r 
 CUSTOM
FUNC=sqrt(x\*x+y\*y+z\*z) 
"])
c1 -- c1.w --> f
linkStyle 5 stroke:red,color:red;
c1 -- c1.x --> f
linkStyle 6 stroke:red,color:red;
c1 -- c1.y --> f
linkStyle 7 stroke:red,color:red;
c1 -- c1.z --> f
linkStyle 8 stroke:red,color:red;
r -- r --> f
linkStyle 9 stroke:red,color:red;
f(["label=f 
 CUSTOM
FUNC=w\*(x^4+y^4+z^4)/(r^4) 
"])
ones(["label=ones 
 CONSTANT 
"])
f -- f --> s
linkStyle 10 stroke:red,color:red;
ones -- ones --> s
linkStyle 11 stroke:blue,color:blue;
s(["label=s 
 MATRIX_VECTOR_PRODUCT 
"])
s -- s --> sums
linkStyle 12 stroke:blue,color:blue;
sums(["label=sums 
 SUM 
"])
sums -- sums --> 11
11("label=#64;11 
 PRINT
FILE=colvar 
")

Once again the elements of any matrices and vectors are recomputed during the backward loop when we apply the chain rule.

Extending these symmetry functions

There are ways of defining the neighbourhood of atom $i$ that do not involve using the CONTACT_MATRIX keyword. For example, you can define a matrix in which element $i,j$ is only non-zero if there is a hydrogen bond between atom $i$ and atom $j$ (HBOND_MATRIX or HBPAMM_MATRIX) or, if you have molecules, you can say that molecule $i$ and $j$ are connected if they are within a cutoff of each other and if the two molecules have some favourable orientation with respect to each other. In all these cases you can also use a COMPONENTS keywords to get the direction of the bonds in the first coordination sphere as I have described in the article above. You are thus not forced to use a switching function to define the $\sigma(r_{ij})$ part of the symmetry function that was defined above. You can use other methods to determine whether atom $i$ and atom $j$ are adjacent or not. Consequently, the implementation of these symmetry function in PLUMED is quite flexible and allows users to try many things without modifying the code.