Implementing Behler-Parinello symmetry functions

There are other codes that implement the Behler-Parinello symmetry functions and to my knowledge no one has used these functions as a CV for a metadynamics simulation. It is, therefore, not unreasonable to ask if an implementation of these CVs in PLUMED is really necessary. The implementation I will describe below was done relatively quickly and I hope it demonstrates how new functionality can be quickly prototyped in the modified version of PLUMED that I have described in these pages.

Implementing $G^1$ symmetry functions

The $G^1$ symmetry function that Behler introduces in the paper that I have linked above is:

\[G^1_i = \sum_{j \ne i} f_c(R_{ij})\]

In this expression, $f_c$ is a switching function and $R_{ij}$ is the distance between atom $i$ and $j$. We have seen this function already, however, it is just the coordination number. To calculate it using PLUMED we can use the following input.

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
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-100 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.={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)}
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=100
g1: 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,ones
# This will print out the 100 coordination number values calculated by the input above.
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=g1 FILEthe name of the file on which to output these quantities=colvar

We have thus already implemented the first of Behler’s symmetry functions.

Implmenting radial symmetry functions

The $G^2$ and $G^3$ symmetry functions that Behler introduces in his paper are:

\[G^2_i = \sum_{j \ne i} e^{-\nu(R_{ij} -R_s)^2} f_c(R_{ij}) \qquad \textrm{and} \qquad G^3_i = \sum_{j \ne i} \cos(\kappa R_{ij}) f_c(R_{ij})\]

$\nu$, $R_s$ and $\kappa$ here are parameters, while $f_c$ is a switching function that acts on the distance $R_{ij}$ between atom $i$ and $j$. In calculating these two functions we need to do a sum of a function of all the $R_{ij}$ values. We thus have everything we need to calculate these variables within PLUMED and need to implement nothing new. We can simply use the following code to compute $G^2$ and $G^3$ with $\nu=1$, $r_s=3$ and $\kappa=1$:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix.  Element i,j of this matrix tells you if atoms i and j are within a certain cutoff of each other
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-100 COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix 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.={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)}
# Now calculate a matrix with all the R_ij values
cmatr: 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
# Compute the quantity in the summation in the expression for G^2.  The output here is a 100x100 matrix.
g2_f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmatr,cmat.w FUNCthe function you wish to evaluate=y*exp(-(x-3)^2) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Compute the quantity in the summation in the expression for G^3.  The output here is a 100x100 matrix.
g3_f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmatr,cmat.w 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
# Now multply the matrices above by a vector of all ones to do the summations
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=100
g2: 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=g2_f,ones
g3: 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=g3_f,ones
# Print out the 100 values for g2
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=g2 FILEthe name of the file on which to output these quantities=g2_file
# Print out the 100 values for g3
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=g3 FILEthe name of the file on which to output these quantities=g3_file

Notice the flexibility of this implementation for the Behler symmetry function you can basically compute any weighted function of the bond vectors in the first coordination sphere using an input like this one. A similar input (see below) could be used to calculate the FCC cubic parameter:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix.  Element i,j of this matrix tells you if atoms i and j are within a certain cutoff of each other
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-100 COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix 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.={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)}
# Now calculate a matrix with all the R_ij values
cmatr: 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
# Now calculate the fcc cubic parameter for each bond.  This outputs a matrix
fcc_f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmatr FUNCthe function you wish to evaluate=((x^4)*(y^4)+(x^4)*(z^4)+(y^4)*(z^4))/(r^8)-27*(x^4)*(y^4)*(z^4)/(r^12) VARthe names to give each of the arguments in the function=x,y,z,r PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And multiply the above by the weights to get another 100 x 100 matrix
wfcc_f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=fcc_f,cmat.w 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
# Now sum the rows of the matrix above to get a vector with 100 elements.  One symmetry function for each atom. 
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=100
fcc_u: 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=wfcc_f,ones 
# This is the coordination number
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,ones
# We devide fcc_u by the coordination number to get the average value of the function above for the atoms in the first coordination sphere.
fcc: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=fcc_u,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
# And finally we print the 100 symmetry function values to a file called colvar.
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=fcc FILEthe name of the file on which to output these quantities=colvar

The code is faster if the calculations that are done with CUSTOM keywords in the input above are implemented directly in C++. Even so the fact that you can so quickly prototype complicated CVs like these is (I hope) useful.

Implementing angular symmetry functionsa

One of the angular symmetry functions that Behler introduces is:

\[G^5_i = 2^{1-\zeta} \sum_{j,k\ne i} (1 + \lambda\cos\theta_{ijk})^\zeta e^{-\nu(R_{ij}^2 + R_{ik}^2)} f_c(R_{ij}) f_c(R_{ik})\]

In this expression $\zeta$, $\nu$ and $\lambda$ are all parameters. $f_c$ is a switching function which acts upon $R_{ij}$, the distance between atom $i$ and atom $j$, and $R_{ik}$, the distance between atom $i$ and atom $k$. $\theta_{ijk}$ is then the angle between the vector that points from atom $i$ to atom $j$ and the vector that points from atom $i$ to atom $k$.

Even this complicated function can be calculated directly in PLUMED. The input below caluclates the $G^5_i$ value for atom 1 with $\zeta=3$, $\lambda=3 and $\nu=0.1$.

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the distances between atom one and all the other atoms
d1: 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=1,3 ATOMS3the pair of atom that we are calculating the distance between=1,4 ATOMS4the pair of atom that we are calculating the distance between=1,5 ATOMS5the pair of atom that we are calculating the distance between=1,7 ATOMS6the pair of atom that we are calculating the distance between=1,8
# Now transform the distances above by a switching function
d1lt: LESS_THANUse a switching function to determine how many of the input variables are less than a certain cutoff. More details ARGthe values input to this function=d1 SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)}
# Calculate the vectors between atom one and all the other atoms
d1c: 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=1,3 ATOMS3the pair of atom that we are calculating the distance between=1,4 ATOMS4the pair of atom that we are calculating the distance between=1,5 ATOMS5the pair of atom that we are calculating the distance between=1,7 ATOMS6the pair of atom that we are calculating the distance between=1,8 
# Now calculate the lengths of all the vectors connecting atom 1 to the other atoms
d1r2: COMBINECalculate a polynomial combination of a set of other variables. More details ARGthe values input to this function=d1c.x,d1c.y,d1c.z POWERS the powers to which you are raising each of the arguments in your function=2,2,2 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
d1r: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d1r2 FUNCthe function you wish to evaluate=sqrt(x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And find unit vectors that connect atom 1 to all its neighbours
d1ux: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d1c.x,d1r 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
d1uy: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d1c.y,d1r 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
d1uz: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=d1c.z,d1r 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 put the three unit vectors into 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=d1ux,d1uz,d1uz
# Now create a matrix in which element i,j = f_c(R_ij) f_c(R_ik)
wmat: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ELEMENTS_ON_DIAGONAL_ARE_ZERO set all diagonal elements to zero ARGthe labels of the two vectors from which the outer product is being computed=d1lt,d1lt
# And a matrix containing the cosines of the angles between the two vectors
stackT: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=stack
cmat: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=stack,stackT
# Now compute a matrix containing (1 + \lambda cos(\theta_ijk))^3
pmat: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO FUNCthe function you wish to evaluate=(1+3*cos(x))^3 
# And a matrix containing in which element j,k is R_ij^2 + R_ik^2 
smat: OUTER_PRODUCTCalculate the outer product matrix of two vectors More details ARGthe labels of the two vectors from which the outer product is being computed=d1r2,d1r2 FUNC the function of the input vectors that should be put in the elements of the outer product=x+y 
# Now take the exponential of this matrix
emat: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=smat FUNCthe function you wish to evaluate=exp(-0.1*x) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And combine everything to get G^5
g5mat: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat,emat,wmat FUNCthe function you wish to evaluate=x*y*z PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Now sum all the elements of g5mat
g5s: SUMCalculate the sum of the arguments More details ARGthe vector/matrix/grid whose elements shuld be added together=g5mat PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And multiply by 0.5 and 2^{1-\zeta}.  We need to multiply by 0.5 here to avoid double counting
g5: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=g5s FUNCthe function you wish to evaluate=0.5*0.25*x PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And print the final scalar to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=g5 FILEthe name of the file on which to output these quantities=colvar

This input computes $G^5$ for a single atom and is rather cumbersome. The calculation here is not particularly well optimised as we have to calculate the distance between atom 1 and all the other atoms. In other words, when an input similar to this one is used PLUMED is not using link cells or neighbour lists to optimse the calculation. This way of implementing $G^5$ is thus not really suited to production applications. I thus wrote a second implementation of the $G^5$ function. The input for this alternative implementation is 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 x,y and z components of the bond vectors
# This action calculates 4 100x100 matrices
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-100 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.={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix

# Compute the symmetry function for the 100 atoms from the 4 100x100 matrices output
# by cmat.  The output from this action is a vector with 100 elements
beh3: GSYMFUNC_THREEBODYCalculate functions of the coordinates of the coordinates of all pairs of bonds in the first coordination sphere of an atom More details ...
    WEIGHTthe matrix that contains the weights that should be used for each connection=cmat.w ARGthree matrices containing the bond vectors of interest=cmat.x,cmat.y,cmat.z
    FUNCTION1the parameters of the function you would like to compute={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3*cos(ajik))^3 LABEL=g5}
...

# Print the 100 symmetry function values to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=beh3.g5 FILEthe name of the file on which to output these quantities=colvar

The GSYMFUNC_THREEBODY action sums over all the distinct triples of atoms that are identified in the contact matrix. This action auses lepton and can thus compute any function of the following four quantities:

Furthermore we can calculate more than one function of these four quantities at a time as illustrated by the input 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 x,y and z components of the bond vectors
# This action calculates 4 100x100 matrices
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-100 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.={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix

# Compute the 4 symmetry function below for the 100 atoms from the 4 100x100 matrices output
# by cmat.  The output from this action is a vector with 100 elements
beh3: GSYMFUNC_THREEBODYCalculate functions of the coordinates of the coordinates of all pairs of bonds in the first coordination sphere of an atom More details ...
    WEIGHTthe matrix that contains the weights that should be used for each connection=cmat.w ARGthree matrices containing the bond vectors of interest=cmat.x,cmat.y,cmat.z
    FUNCTION1the parameters of the function you would like to compute={FUNC=0.25*(cos(pi*sqrt(rjk)/4.5)+1)*exp(-0.1*(rij+rik+rjk))*(1+2*cos(ajik))^2 LABEL=g4}
    FUNCTION2the parameters of the function you would like to compute={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3.5*cos(ajik))^3 LABEL=g5}
    FUNCTION3the parameters of the function you would like to compute={FUNC=0.125*(1+6.6*cos(ajik))^4 LABEL=g6}
    FUNCTION4the parameters of the function you would like to compute={FUNC=sin(3.0*(ajik-1)) LABEL=g7}
...

# Print the 4 sets of 100 symmetry function values to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=beh3.g4,beh3.g5,beh3.g6,beh3.g7 FILEthe name of the file on which to output these quantities=colvar

As we saw for the radial symmetry function, we thus again have an implementation of these angular symmetry functions that can be quickly used to prototype new function types.

Conclusions

I hope the examples above have convinced you that the PLUMED input file now provides the flexibility required to prototype CVs that are really rather complicated. Using CUSTOM actions when implementing these CVs is probably not optimal when it comes to computational performance but I am not convinced that this matters. I think our objective should always be to have less code as this makes PLUMED easier to maintain.

It is also the case that many of the CVs that are implemented in PLUMED are not used by many people. The understanding of what particular actions do is thus lost with time. When CVs are implemented in the way above the series of action that are performed within them are more transparent. This transparency allows allows other to understand what has been done in the past and reduces the requirements to write large amounts of documentation.

Lastly, if users need code that runs faster they can always reimplement expensive CVs that are implemented in PLUMED. Furthermore, in any effort to reimplement CVs they are helped by the fact that PLUMED offers reference values that they can test their new code against.