Steinhardt parameters in PLUMED

The Steinhardt parameters are arguably some of the most complicated collective variables to use in PLUMED. These variables are difficult to use because researchers use a variety of subtly different variations on the same basic concept. When you are trying to reproduce the values quoted in a particular paper you often find that you need to spend time comparing the choices that were made in the paper with the choices that have been used in the code you are employing. To resolve this problem I have tried to use the new shortcut features that I have implemented in PLUMED to make the implementation of the Steinhardt parameters in PLUMED a little more transparent. Hopefully this more transparent and piecemeal implementation allows users to implement and use all the variations in the literature by modifying the template input files that are generated by the shortcut actions that can be used to run the most commonly used Steinhardt variants.

In the following article I will thus explain how the Steinhardt parameters have been implemented and the variety of things that you can do within them.

Calculating the Steinhardt parameter for a single atom

The Steinhardt parameter is an example of a symmetry function. In other words, it is a function of the position of the atoms in the first coordination sphere around a central atom. If we have a system of $N$ indistinguishable atoms we can thus calculate $N$ separate instances of the Steinhardt parameter.

Steinhardt parameters are complicated as the spherical harmonics are used to calculate them. Consequently, when we calculate the Steinhardt parameter for atom $i$ we evaluate $(2l + 1)$ complex quntities as follows:

\[q_{lm}(i) = \frac{\sum_{j \ne i} \sigma(r_{ij})Y_l^m(\theta_{ij},\phi_{ij}) }{ \sum_{j \ne i} \sigma(r_{ij})}\]

In this expression the sums run over all the specified atoms. $\sigma(r_{ij})$ is a switching function that acts on the distance, $r_{ij}$, between atom $i$ and atom $j$. $Y_l^m(\theta_{ij},\phi_{ij})$ is the spherical harmonic, which is a function of the angles $\theta_{ij}$ and $\phi_{ij}$ the vector that connects atom $i$ to atom $j$ and the $x$ and $z$ axis of the lab frame respectively.

We can calculate $q_{lm}(i)$ values for atom 1 in PLUMED using an input like the one shown below. Here I am using $l=1$ so as to keep the input short.

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix between atom 1 and all the other atoms
cmat: 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=1 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=2-100 COMPONENTS also calculate the components of the vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
# Create a vector containing 99 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=99
# Calculate the coordination number for atom 1 by multiplying the contact matrix above 
# by a vector of ones
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=cmat.w,ones
# Evaluate the Y_l^m values for each vector in the first coordination sphere and multiply these 
# values by \sigma(r_ij).  This action outputs 6 matrices that correspond to the real and imaginary 
# components of Y^l_m for m values of -1, 0 and +1.
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
# Now evaluate the numerator in the expression above.  This will output 6 scalars. 
sp: 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=sh.*,ones
# Now evaluate the q_lm(i) values by dividing by  the coordination number of the atom
q1-rmn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-n1,coord 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
q1-imn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-n1,coord 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
q1-rmn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-0,coord 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
q1-imn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-0,coord 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
q1-rmn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-p1,coord 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
q1-imn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-p1,coord 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 output these six scalars to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1-rmn-n1,q1-imn-n1,q1-rmn-0,q1-imn-0,q1-rmn-p1,q1-imn-p1 FILEthe name of the file on which to output these quantities=colvar

The $q_{lm}(i)$ values that are output by the PLUMED input above are not rotationally invariant. We thus introduce $q_l(i)$ as follows:

\[q_{l}(i) = \sqrt{ \sum_{m=-l}^l | q_{lm}(i) |^2 }\]

This quantity can be calculated using the following PLUMED input:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix between atom 1 and all the other atoms
cmat: 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=1 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=2-100 COMPONENTS also calculate the components of the vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
# Create a vector containing 99 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=99
# Calculate the coordination number for atom 1 by multiplying the contact matrix above 
# by a vector of ones
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=cmat.w,ones
# Evaluate the Y_l^m values for each vector in the first coordination sphere and multiply these 
# values by \sigma(r_ij).  This action outputs 6 matrices that correspond to the real and imaginary 
# components of Y^l_m for m values of -1, 0 and +1. 
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
# Now evaluate the numerator in the expression above.  This will output 6 scalars.  
sp: 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=sh.*,ones
# Now evaluate the q_lm(i) values by dividing by  the coordination number of the atom
q1-rmn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-n1,coord 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
q1-imn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-n1,coord 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
q1-rmn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-0,coord 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
q1-imn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-0,coord 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
q1-rmn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-p1,coord 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
q1-imn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-p1,coord 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
# Calculate the sum of the squares of the q_lm(i) values
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=q1-rmn-n1,q1-imn-n1,q1-rmn-0,q1-imn-0,q1-rmn-p1,q1-imn-p1
# Take the square root of the sum of the squares
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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 output this single scalar that results to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1 FILEthe name of the file on which to output these quantities=colvar

Alternatively, we could compute the following quantity:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix between atom 1 and all the other atoms
cmat: 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=1 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=2-100 COMPONENTS also calculate the components of the vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
# Create a vector containing 99 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=99
# Calculate the coordination number for atom 1 by multiplying the contact matrix above
# by a vector of ones
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=cmat.w,ones
# Evaluate the Y_l^m values for each vector in the first coordination sphere and multiply these
# values by \sigma(r_ij).  This action outputs 6 matrices that correspond to the real and imaginary
# components of Y^l_m for m values of -1, 0 and +1.
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
# Now evaluate the numerators in the expression above.  This will output 6 scalars.
sp: 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=sh.*,ones
# Now take the sum of the squares of the numerators
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
# Square root 
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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 divide by the number of neighbours
anorm: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1,coord 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

This input computes:

\[q_l(i) = \frac{1}{N(i)} \sqrt{\sum_{m=-1}^l N(i)^2 | q_{lm}(i) |^2} \qquad \textrm{where} \qquad N(i) = \sum_{j \ne i} \sigma(r_{ij})\]

This should be equivalent to what is computed by the first input but there may be some small differences in practise because of numerical error. I only mention this alternative claculation method as this is done when the shortcut below is employed:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
q1: Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More details SPECIESAthe list of atoms for which the symmetry function is being calculated=1 SPECIESBthe list of atoms that can be in the environments of each of the atoms for which the symmetry function is being calculated=2-100 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.={RATIONAL D_0=2.0 R_0=1.0}
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1 FILEthe name of the file on which to output these quantities=colvar

As you can see the Q1 action in this second input expands to the input above.

Calculating the Steinhardt parameters for multiple atoms

Consider the following input for calculating Steinhardt parameters:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
q4: Q4Calculate fourth order Steinhardt parameters. This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0}
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q4 FILEthe name of the file on which to output these quantities=colvar

The colvar file that is output by this action will contain 100 $q_l(i)$ values for each time step. If you look at the expanded version of the input you can see that it has the same structure as the inputs that we looked at for calculating the $q_l(i)$ for a single atom in the previous section. The difference is in the contact matrix that is calculated in the first line of the input. When we are calculating a single $q_l(i)$ value this is a $1 \times 99$ matrix. For the input above the contact matrix is a $100 \times 100$ matrix. When we do the MATRIX_VECTOR_PRODUCT we thus have a series of vectors with 100 elements (i.e. one element for each atom). The vectors output by this action are this size as there are 100 rows in the constact matrix and not just one.

The expanded version of the input file is clear because when we apply the SPHERICAL_HARMONIC operation we do it element-wise. The user can use this action as if it is generating $2(2l+1)$ matrices with the same shape as the CONTACT_MATRIX by applying computing $\sigma(r_{ij})Y_l^m(\theta_{ij},\phi_{ij})$ for each bond vector. (PLUMED is (obviously) not computing the Steinhardt parameter in this way as that would be horribly inneficient.) They are thus free to do many different things with the vectors that are computed.

Calculating functions of the distribution of Steinhardt Parameters

When you use the shortcut actions Q1, Q3, Q4 and Q6 you can use keywords such as MEAN, LOWEST, LESS_THAN, MORE_THAN to calculate functions of the distribution of $q_l(i)$ values in the same way that you would use these keywords with other symmetry functions such as COORDINATIONNUMBERS. When you use these keywords with Steinhardt parameters there are subtleties in the literature that you should be know about. These subtleties are easily understood if you look at the full inputs that PLUMED generates when you use these shortcuts. Lets start by considering the following input:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix for the first 100 atoms in the input file
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
# Create a vector containing 100 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=100
# Calculate the coordination numbers for the atoms by multiplying the contact matrix above
# by a vector of ones
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=cmat.w,ones
# Evaluate the Y_l^m values for each vector in the first coordination sphere and multiply these
# values by \sigma(r_ij).  This action outputs 6 matrices that correspond to the real and imaginary
# components of Y^l_m for m values of -1, 0 and +1.
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
# Now evaluate the numerator in the expression above.  This will output 6 vectors each of which has 100 components.
sp: 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=sh.*,ones
# Now evaluate the q_lm(i) values for each atom by dividing by  the coordination number of the corresponding atom
q1-rmn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-n1,coord 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
q1-imn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-n1,coord 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
q1-rmn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-0,coord 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
q1-imn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-0,coord 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
q1-rmn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-p1,coord 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
q1-imn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-p1,coord 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
# Calculate the sum of the squares of the q_lm(i) values
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=q1-rmn-n1,q1-imn-n1,q1-rmn-0,q1-imn-0,q1-rmn-p1,q1-imn-p1
# Take the square root of the sum of the squares.  This outputs a vector with 100 components
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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
# Calculate the average Q1 values.
q1_mean: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And output this single scalar that results to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1_mean FILEthe name of the file on which to output these quantities=colvar

The average that is being calculated here is:

\[\langle q_l \rangle = \frac{1}{N} \sum_{i=1}^N q_l(i)\]

which is also what is computed if we use the following shortcut:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
q1: Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0} MEAN calculate the mean of all the quantities
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1_mean FILEthe name of the file on which to output these quantities=colvar

We do not need to average the $q_l(i)$ values. We could instead calculate the average of the $q_{lm}(i)$ values and the magnitude of the resulting vector like this:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the contact matrix for the first 100 atoms in the input file
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
# Create a vector containing 100 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=100 
# Calculate the coordination numbers for the atoms by multiplying the contact matrix above
# by a vector of ones
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=cmat.w,ones
# Evaluate the Y_l^m values for each vector in the first coordination sphere and multiply these
# values by \sigma(r_ij).  This action outputs 6 matrices that correspond to the real and imaginary
# components of Y^l_m for m values of -1, 0 and +1.
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
# Now evaluate the numerator in the expression above.  This will output 6 vectors each of which has 100 components.
sp: 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=sh.*,ones
# Now evaluate the q_lm(i) values for each atom by dividing by  the coordination number of the corresponding atom
q1-rmn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-n1,coord 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
q1-imn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-n1,coord 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
q1-rmn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-0,coord 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
q1-imn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-0,coord 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
q1-rmn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-p1,coord 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
q1-imn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-p1,coord 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 calculate the mean for each m value.  The output from these six actions are scalars
q1-rms-n1: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1-rmn-n1 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
q1-ims-n1: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1-imn-n1 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
q1-rms-0: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1-rmn-0 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
q1-ims-0: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1-imn-0 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
q1-rms-p1: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1-rmn-p1 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
q1-ims-p1: MEANCalculate the arithmetic mean of the elements in a vector More details ARGthe values input to this function=q1-imn-p1 PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Now calculate the sum of the squares
q1_mean2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=q1-rms-n1,q1-ims-n1,q1-rms-0,q1-ims-0,q1-rms-p1,q1-ims-p1
# Take the square root 
q1_mean: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_mean2 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
# Print the final value
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1_mean FILEthe name of the file on which to output these quantities=colvar

The quantity that has been computed by this input is:

\[\langle q_l \rangle = \sqrt{ \sum_{m=-l}^l | \langle q_{lm} \rangle |^2 } \qquad \textrm{where} \qquad \langle q_{lm} \rangle = \frac{1}{N} \sum_{i=1}^N q_{lm}(i)\]

This quantity can also be calculated using the following shortcut input:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
q1: Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0} VMEAN calculate the norm of the mean vector
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=q1_vmean FILEthe name of the file on which to output these quantities=colvar

These are probably not the only two ways of averaging over all the Steinhardt parameters. I hope, however, that from these two examples you can see just how much flexibility can be achieved by only modifying the input file.

Calculating Leichner-Delago-style Steinhardt order parameters

In this paper Leichner and Delago proposed calculating local averages of the Steinhardt parameters. In PLUMED you can calculate these order parameters using the following shortcut action:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the Steinhardt parameters
q6: Q6Calculate sixth order Steinhardt parameters. This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0}
# Take a local average of the Steinhardt parameters
lq6: LOCAL_AVERAGECalculate averages over spherical regions centered on atoms This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=q6 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.={RATIONAL D_0=3.0 R_0=1.5}
# This prints out the vector of 100 locally averaged order paraemters
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lq6 FILEthe name of the file on which to output these quantities=colvar

The expanded version of this shortcut looks rather complicated so lets consider a simpler calculation first. In the following input I am calculating a locally averaged value for the coordination number:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the coordination numbers of the first 100 atoms in the system in the usual way.
# In other words, compute the contact matrix and multiply this matrix by a vector containing all ones.
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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat,ones
# Now calculate a second contact matrix, which we will use for our local averaging.
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
# Get a new set of coordination numbers from cmat2 by multiplying above matrix by vector of ones.
coord2: 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=cmat2,ones
# Multiply the second contact matrix above by the vector containing the first set of coordination numbers
# this gives us a vector that contains the sum of the coordination numbers for all the atoms in the first 
# coordination spheres
prod: 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=cmat2,coord
# And compute the local average of the coordination number by taking the vector (prod) that contains the sum
# of the coordination numbers in the first coordination sphere.  Add to this quantity the value of the coordination 
# number of the central atom and divide this sum by 1 + the number of atoms in the first coordination sphere.
lav: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=prod,coord,coord2 FUNCthe function you wish to evaluate=(x+y)/(1+z) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# This will print out the 100 average values for the local averages of the coordination numbers
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lav FILEthe name of the file on which to output these quantities=colvar

The point I want to make with this input is that the heart of this process of taking a local average of a symmetry function involves multiplying the vector of symmetry function values by a contact matrix. In other words, the locally averaged coordination number $\widehat{c}(i)$ can be computed as:

\[\widehat{c}(i) = \frac{c(i) + \sum_{j\ne i} \sigma(r_{ij})c(j)}{1 + \sum_{j\ne i} \sigma(r_{ij})}\]

where $\sigma(r_{ij})$ is a switching function that acts on the distance $r_{ij}$ between atom $i$ and atom $j$ and $c(i)$ is the coordination number of atom $i$.

With that in mind lets now consder the full input for the Leichner and Delago order parameters, which is as follows:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# This is the usual input that we have for computing the Steinhardt parameters that is explained above.
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
q1-rmn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-n1,coord 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
q1-imn-n1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-n1,coord 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
q1-rmn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-0,coord 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
q1-imn-0: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-0,coord 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 
q1-rmn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.rm-p1,coord 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
q1-imn-p1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=sp.im-p1,coord 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 lets calculate the contact matrix that we will use for the local averaging
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
# Get a new set of coordination numbers from cmat2 by multiplying above matrix by vector of ones.
coord2: 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=cmat2,ones
# This action constructs a 100 x 6 matrix that contains the q_lm(i) values for all the atoms
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=q1-rmn-n1,q1-imn-n1,q1-rmn-0,q1-imn-0,q1-rmn-p1,q1-imn-p1 
# We now multiply each 100 element vector of q_lm(i) values by the contact matrix.  The result is
# 6 100 element vectors that contain.  Element i in these vector contains the sum of the q_lm(j) values
# for those atoms that are in the first coordination sphere of atom i.
prod: 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=cmat2,q1-rmn-n1,q1-imn-n1,q1-rmn-0,q1-imn-0,q1-rmn-p1,q1-imn-p1
# The 6 vectors output by the command above are combined into a 100 x 6 matrix that contains the sums of 
# the q_lm(i) values for the neighbours of atom i.
vpstack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=prod.q1-rmn-n1,prod.q1-imn-n1,prod.q1-rmn-0,prod.q1-imn-0,prod.q1-rmn-p1,prod.q1-imn-p1
# These two commands create a 100 x 6 element matrix.  Every column of this matrix is identical and contains 
# the 100 coordination numbers for the atoms.
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=coord2,lones
# We now compute our local averages for q_lm(i).  The output here is a 100 x 6 matrix of local average values for q_lm(i)
av: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=vpstack,vpstack,unorm FUNCthe function you wish to evaluate=(x+y)/(1+z) PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Do an element-wise square of the matrix above.  This is not doing matrix multiplication.
av2: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=av FUNCthe function you wish to evaluate=x*x PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# We now multiply our 100x6 matrix by a vector with 6 elements.  The result is a vector with 100 elements - one element
# for each atom in the system.  This is essentially the square of the norm of the vectors of locally-averaged q_lm(i) values.
aq1_2: 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=av2,lones 
# And these are the norms of the vectors of locally-averaged q_lm(i) values.
aq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=aq1_2 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
# We now print out our 100 locally averaged q_l values.
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=aq1 FILEthe name of the file on which to output these quantities=colvar

Notice that we do the same here as we did with the coordination number input. In other words, we calculate locally-averaged versions of the $q_{lm}(i)$ parameters using:

\[\widehat{q}{lm}(i) = \frac{q_{lm}(i) + \sum_{j\ne i} \sigma(r_{ij})q_{lm}(j)}{1 + \sum_{j\ne i} \sigma(r_{ij})}\]

The final, locally-averaged version of $\widehat{q}_l(i)$ is then defined as:

\[\widehat{q}_{l}(i) = \sqrt{ \sum_{m=-l}^l | \widehat{q}_{lm}(i) |^2 }\]

which is the familiar summing over the $l$ values that we have seen in earlier sectioms to make the final symmetry function rotationally invariant.

Calculating ten-Wolde-Frenkel-style Steinhard order parameters

Before discussing how the ten-Wolde-Frenkel order parameters are computed it is useful to introduce yet another variant on the Steinhardt parameter. This variant is computed as follows

\[\overline{q}_{lm}(i) = \frac{1}{\overline{q}_l} \sum_{j \ne i } \sigma(r_{ij})Y_l^m(\theta_{ij},\phi_{ij}) \qquad \textrm{where} \qquad \overline{q}_l(i) = \sqrt{ \sum_{m=-l}^l |q_{lm}(i)|^2 }\]

As always $\sigma(r_{ij})$ is a switching function that acts on the distance, $r_{ij}$, between atom $i$ and atom $j$. $Y_l^m(\theta_{ij},\phi_{ij})$ is then the spherical harmonic function which depends on the direction the bond between atom $i$ and atom $j$ point in. Notice that this version of the Steinhardt parameter has the followign property:

\[\sum_{m=-l}^l \overline{q}_{lm}(i)^{*} \overline{q}_{lm}(i) = 1\]

We can compute the $\overline{q}_{lm}(i)$ values in PLUMED using the following input:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# This is the usual input that we have for computing the Steinhardt parameters that is explained above.
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones 
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
# This calculates the sum of the squares of the unormalized \overline{q}_lm(i) values.
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
# This calculates \overline{q}_l(i)
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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
# Create a 100 x 6 matrix to hold all the unofrmalised \overline{q}_lm(i) values.
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
# Create a second 100 x 6 matrix.  All the columns of this matrix are identical.  They just contain the 100 \overline{q}_l(i) values.
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones
# We can now calculate the \overline{q}_lm(i) values by dividing the unormalised values by the \overline{q}_l(i).
# The output here is a 100 x 6 matrix
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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 print out the 100 x 6 matrix of overline{q}_lm(i) values.
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=oq1 FILEthe name of the file on which to output these quantities=colvar

We introduce $\overline{q}_{lm}(i)$ as we can use the following inner product to measure the similarity between the environment around atom $i$ and the environment around atom $j$:

\[S_{ij} = \sum_{m=-l}^l \overline{q}_{lm}(i)^{*} \overline{q}_{lm}(j)\]

Notice furthermore that if we have $N$ atoms we can store all the real and imaginary parts of the $\overline{q}_{lm}(i)$ values in an $N \times 2l$ matrix $\mathbf{Q}$. The inner product above can be computed for every pair of atoms in our system by doing the following piece of matrix algebra:

\[\mathbf{S} = \mathbf{Q} \mathbf{Q}^T\]

$S$ here is thus an $N\times N$ square matrix. To calculate this matrix in PLUMED you would do the following:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# These lines compute the matrix of \overline{q}_lm(i) values in the way described above
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones  
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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 
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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 transpose the matrix of \overline{q}_lm(i) values
oq1T: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=oq1
# Calculate the matrix product
s: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=oq1,oq1T
# And output the 100 x 100 matrix of inner products
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=s FILEthe name of the file on which to output these quantities=colvar

The matrix $\mathbf{S}$ defined above is not the one used when you are computing the ten-Wolde and Frenkel style order parametres. Element $i,j$ of the matrix, $\mathbf{F}$, that is used when computing the ten-Wolde and Frenkel style order parameters is large when atoms $i$ and $j$ are within a certain cutoff of each other AND when the environment around the two atoms are similar. Element $i,j$ of this matrix is thus computed as:

\[F_{ij} = \sigma(r_{ij}) \sum_{m=-l}^l \overline{q}_{lm}(i)^{*} \overline{q}_{lm}(j)\]

where $\sigma$ is a switching function that acts on the distance $r_{ij}$ between atom $i$ and atom $j$. The full matrix can be computed in PLUMED using the following input:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# These lines compute the matrix of \overline{q}_lm(i) values in the way described above
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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 
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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 transpose the matrix of \overline{q}_lm(i) values
oq1T: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=oq1
# Calculate a second contact matrix.  This will contain the \sigma(r_ij) values in the second equation above.
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
# Calculate the matrix product, S, described above
s: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=oq1,oq1T
# Now take the element-wise product of the contact matrix and the matrix of inner products.
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat2,s 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 output the 100 x 100 matrix 
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

Once you ahve computed the matrix $F$ there are many things that you can do with it. You can for instance sum all the elements of the matrix using an input like this:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# These lines compute the matrix F
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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
oq1T: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=oq1
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
s: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=oq1,oq1T
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat2,s 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
# Add all the elemnets of the matrix f together
fs: SUMCalculate the sum of the arguments More details ARGthe values input to this function=f PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# This outputs a single scalar
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=fs FILEthe name of the file on which to output these quantities=colvar

Some people prefer to apply a switching function to the elements of $F$ before summing like this:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# These lines compute the matrix F
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 vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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
oq1T: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=oq1
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
s: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=oq1,oq1T
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat2,s 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
# Tranform all the elements of f with a switching function that is one when f_ij>0.5
fmt: 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=f SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL R_0=0.5}
# Add all the elemnets of the matrix fmt together
fs: SUMCalculate the sum of the arguments More details ARGthe values input to this function=fmt PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# This outputs a single scalar 
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=fs FILEthe name of the file on which to output these quantities=colvar

You could even apply the switching function to $S$ and only then multiply the transformed version of $S$ by the second contact matrix. Alternatively, you can calculate a new symmetry function for each of the atoms in your system by multiplying $F$ by a vector of all ones as following:

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 COMPONENTS also calculate the components of the vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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
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=cmat.w,ones
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones    
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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
oq1T: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=oq1
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
s: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=oq1,oq1T
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat2,s 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  
# This outputs a vector with 100 elements 
lq1: 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 
# This prints the 100 symmetry function values to the file colvar
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lq1 FILEthe name of the file on which to output these quantities=colvar

There is a shortcut that does the same as this last input which works as follows:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the steinhardt parameters
q1: Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0}
# Calculate the local steinhard parameters
lq1: LOCAL_Q1Calculate the local degree of order around an atoms by taking the average dot product between the q_1 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere. This action is a shortcut. More details SPECIESthe label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters=q1 SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL D_0=2.0 R_0=1.0}
# Print the local steinhardt parameters
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lq1 FILEthe name of the file on which to output these quantities=colvar

Please note that if you use the following input:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the steinhardt parameters
q1: Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0}
# Calculate the local steinhard parameters
lq1: LOCAL_Q1Calculate the local degree of order around an atoms by taking the average dot product between the q_1 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere. This action is a shortcut. More details SPECIESthe label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters=q1 SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL D_0=2.0 R_0=1.0} MEAN calculate the mean of all the quantities
# Print the local steinhardt parameters
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lq1.mean FILEthe name of the file on which to output these quantities=colvar

the quantity in the vector that is averaged is:

\[s_i = \frac{ \sum_j F_{ij} }{ \sum_j \sigma(r_{ij}) }\]

This is not the same vector that is output in the last but one example input above. If you would like to output a vector that contains all the values for the quantity above you can calculate it using the following long 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 COMPONENTS also calculate the components of the vectors 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.={RATIONAL D_0=2.0 R_0=1.0}
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 
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=cmat.w,ones
sh: SPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More details ARGthe values input to this function=cmat.x,cmat.y,cmat.z,cmat.w Lthe value of the angular momentum=1
sp: 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=sh.*,ones
q1_2: COMBINECalculate a polynomial combination of a set of other variables. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO POWERS the powers to which you are raising each of the arguments in your function=2,2,2,2,2,2 ARGthe values input to this function=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
q1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=q1_2 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
stack: VSTACKCreate a matrix by stacking vectors together More details ARGthe values that you would like to stack together to construct the output matrix=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1
lones: 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=6
unorm: OUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More details ARGthe labels of the two vectors from which the outer product is being computed=q1,lones 
oq1: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=stack,unorm 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
oq1T: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=oq1
cmat2: 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.={RATIONAL D_0=2.0 R_0=1.0}
s: MATRIX_PRODUCTCalculate the product of two matrices More details ARGthe label of the two matrices from which the product is calculated=oq1,oq1T
f: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cmat2,s 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  
# This outputs a vector with 100 elements 
lq1: 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 
# This calculates the coordination number again, which is used as the demominator in the expression above
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=cmat2,ones
# And now we normalise by dividing by the coordination number
lq1_av: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=lq1,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
# This prints the 100 symmetry function values to the file colvar
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lq1_av FILEthe name of the file on which to output these quantities=colvar

Alternatively, if you want to use the shortcut you can use:

Click on the labels of the actions for more information on what each action computes
tested onv2.9
tested onmaster
# Calculate the steinhardt parameters
q1: Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More details SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-100 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.={RATIONAL D_0=2.0 R_0=1.0}
# Calculate the local steinhard parameters
lq1: LOCAL_Q1Calculate the local degree of order around an atoms by taking the average dot product between the q_1 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere. This action is a shortcut. More details SPECIESthe label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters=q1 SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above={RATIONAL D_0=2.0 R_0=1.0}
# Print the local steinhardt parameters
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=lq1_av FILEthe name of the file on which to output these quantities=colvar

Conclusion

The implementation of the Steinhardt parameters in previous versions of PLUMED grew rather organically, which made the code confusing and difficult to use. I hope this article convinces readers that the changes that have been made are worthwhile. This new structure exposes the linear algebra that is being used to compute all the variants on these particular order parameters in the input file. It is thus easier to determine what PLUMED is actually computing.

All the old inputs should still work as I have included shortcuts that reproduce new-style input files that reproduce what PLUMED used to do in the past. However, I hope that folks will not use these shortcuts in the future as when these new inputs are used it is much easier for folks to understand what has been done.