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
 |
 |
# Calculate the contact matrix between atom 1 and all the other atoms cmat:CONTACT_MATRIXGROUPA=1 GROUPB=2-100Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Create a vector containing 99 ones ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=99 # Calculate the coordination number for atom 1 by multiplying the contact matrix above # by a vector of ones coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 # Now evaluate the numerator in the expression above. This will output 6 scalars. sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones # Now evaluate the q_lm(i) values by dividing by the coordination number of the atom q1-rmn-n1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-n1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # And output these six scalars to a fileif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=q1-rmn-n1,q1-imn-n1,q1-rmn-0,q1-imn-0,q1-rmn-p1,q1-imn-p1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# Calculate the contact matrix between atom 1 and all the other atoms cmat:CONTACT_MATRIXGROUPA=1 GROUPB=2-100Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Create a vector containing 99 ones ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=99 # Calculate the coordination number for atom 1 by multiplying the contact matrix above # by a vector of ones coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 # Now evaluate the numerator in the expression above. This will output 6 scalars. sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones # Now evaluate the q_lm(i) values by dividing by the coordination number of the atom q1-rmn-n1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-n1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # Calculate the sum of the squares of the q_lm(i) values q1_2:if the output of your function is periodic then you should specify the periodicity of the functionCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=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:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO # And output this single scalar that results to a fileif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=q1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
Alternatively, we could compute the following quantity:
Click on the labels of the actions for more information on what each action computes
|
|
# Calculate the contact matrix between atom 1 and all the other atoms cmat:CONTACT_MATRIXGROUPA=1 GROUPB=2-100Adjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Create a vector containing 99 ones ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=99 # Calculate the coordination number for atom 1 by multiplying the contact matrix above # by a vector of ones coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 # Now evaluate the numerators in the expression above. This will output 6 scalars. sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones # Now take the sum of the squares of the numerators q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 # Square root q1:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO # And divide by the number of neighbours anorm:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the function
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
|
|
q1:Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More detailsSPECIESA=1this keyword is used for colvars such as the coordination numberSPECIESB=2-100this keyword is used for colvars such as the coordination numberSWITCH={RATIONAL D_0=2.0 R_0=1.0}the switching function that it used in the construction of the contact matrixPrint quantities to a file. More detailsARG=q1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
q4:Q4Calculate fourth order Steinhardt parameters. This action is a shortcut. More detailsSPECIES=1-100this keyword is used for colvars such as coordination numberSWITCH={RATIONAL D_0=2.0 R_0=1.0}the switching function that it used in the construction of the contact matrixPrint quantities to a file. More detailsARG=q4the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Create a vector containing 100 ones ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 # Calculate the coordination numbers for the atoms by multiplying the contact matrix above # by a vector of ones coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 # Now evaluate the numerator in the expression above. This will output 6 vectors each of which has 100 components. sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-n1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # Calculate the sum of the squares of the q_lm(i) values q1_2:if the output of your function is periodic then you should specify the periodicity of the functionCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=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:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO # Calculate the average Q1 values. q1_mean:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1the input to this functionPERIODIC=NO # And output this single scalar that results to a fileif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=q1_meanthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
q1:Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More detailsSPECIES=1-100this keyword is used for colvars such as coordination numberSWITCH={RATIONAL D_0=2.0 R_0=1.0}the switching function that it used in the construction of the contact matrixMEANcalculate the mean of all the quantitiesPrint quantities to a file. More detailsARG=q1_meanthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Create a vector containing 100 ones ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 # Calculate the coordination numbers for the atoms by multiplying the contact matrix above # by a vector of ones coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 # Now evaluate the numerator in the expression above. This will output 6 vectors each of which has 100 components. sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-n1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # Now calculate the mean for each m value. The output from these six actions are scalars q1-rms-n1:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1-rmn-n1the input to this functionPERIODIC=NO q1-ims-n1:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1-imn-n1the input to this functionPERIODIC=NO q1-rms-0:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1-rmn-0the input to this functionPERIODIC=NO q1-ims-0:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1-imn-0the input to this functionPERIODIC=NO q1-rms-p1:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1-rmn-p1the input to this functionPERIODIC=NO q1-ims-p1:if the output of your function is periodic then you should specify the periodicity of the functionMEANCalculate the arithmetic mean of the elements in a vector More detailsARG=q1-imn-p1the input to this functionPERIODIC=NO # Now calculate the sum of the squares q1_mean2:if the output of your function is periodic then you should specify the periodicity of the functionCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=q1-rms-n1,q1-ims-n1,q1-rms-0,q1-ims-0,q1-rms-p1,q1-ims-p1 # Take the square root q1_mean:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_mean2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO # Print the final valueif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=q1_meanthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
q1:Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More detailsSPECIES=1-100this keyword is used for colvars such as coordination numberSWITCH={RATIONAL D_0=2.0 R_0=1.0}the switching function that it used in the construction of the contact matrixVMEANcalculate the norm of the mean vectorPrint quantities to a file. More detailsARG=q1_vmeanthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# Calculate the Steinhardt parameters q6:Q6Calculate sixth order Steinhardt parameters. This action is a shortcut. More detailsSPECIES=1-100this keyword is used for colvars such as coordination numberSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Take a local average of the Steinhardt parameters lq6:the switching function that it used in the construction of the contact matrixLOCAL_AVERAGECalculate averages over spherical regions centered on atoms This action is a shortcut. More detailsSPECIES=q6this keyword is used for colvars such as coordination numberSWITCH={RATIONAL D_0=3.0 R_0=1.5} # This prints out the vector of 100 locally averaged order paraemtersthe switching function that it used in the construction of the contact matrixPrint quantities to a file. More detailsARG=lq6the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat,ones # Now calculate a second contact matrix, which we will use for our local averaging. cmat2:the input for this action is the scalar output from one or more other actionsCONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={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:specify the switching function to use between two sets of indistinguishable atomsMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=prod,coord,coord2the input to this functionFUNC=(x+y)/(1+zthe function you wish to evaluatePERIODIC=NO # This will print out the 100 average values for the local averages of the coordination numbersif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=lavthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones q1-rmn-n1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-n1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-n1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-0:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-0,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-rmn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.rm-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO q1-imn-p1:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=sp.im-p1,coordthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # Now lets calculate the contact matrix that we will use for the local averaging cmat2:if the output of your function is periodic then you should specify the periodicity of the functionCONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={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:specify the switching function to use between two sets of indistinguishable atomsMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat2,ones # This action constructs a 100 x 6 matrix that contains the q_lm(i) values for all the atoms stack:the input for this action is the scalar output from one or more other actionsVSTACKCreate a matrix by stacking vectors together More detailsARG=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:the input for this action is the scalar output from one or more other actionsMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=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:the input for this action is the scalar output from one or more other actionsVSTACKCreate a matrix by stacking vectors together More detailsARG=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:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=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:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=vpstack,vpstack,unormthe input to this functionFUNC=(x+y)/(1+zthe function you wish to evaluatePERIODIC=NO # Do an element-wise square of the matrix above. This is not doing matrix multiplication. av2:if the output of your function is periodic then you should specify the periodicity of the functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=avthe input to this functionFUNC=x*xthe function you wish to evaluatePERIODIC=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:if the output of your function is periodic then you should specify the periodicity of the functionMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=av2,lones # And these are the norms of the vectors of locally-averaged q_lm(i) values. aq1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=aq1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO # We now print out our 100 locally averaged q_l values.if the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=aq1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones # This calculates the sum of the squares of the unormalized \overline{q}_lm(i) values. q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=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:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO # Create a 100 x 6 matrix to hold all the unofrmalised \overline{q}_lm(i) values. stack:if the output of your function is periodic then you should specify the periodicity of the functionVSTACKCreate a matrix by stacking vectors together More detailsARG=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:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=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:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=stack,unormthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # And finally print out the 100 x 6 matrix of overline{q}_lm(i) values.if the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=oq1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 q1:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO stack:if the output of your function is periodic then you should specify the periodicity of the functionVSTACKCreate a matrix by stacking vectors together More detailsARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 lones:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=q1,lones oq1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=stack,unormthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # Now transpose the matrix of \overline{q}_lm(i) values oq1T:if the output of your function is periodic then you should specify the periodicity of the functionTRANSPOSECalculate the transpose of a matrix More detailsARG=oq1 # Calculate the matrix product s:the input for this action is the scalar output from one or more other actionsMATRIX_PRODUCTCalculate the product of two matrices More detailsARG=oq1,oq1T # And output the 100 x 100 matrix of inner productsthe input for this action is the scalar output from one or more other actionsPrint quantities to a file. More detailsARG=sthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# 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 detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 q1:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO stack:if the output of your function is periodic then you should specify the periodicity of the functionVSTACKCreate a matrix by stacking vectors together More detailsARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 lones:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=q1,lones oq1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=stack,unormthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO # Now transpose the matrix of \overline{q}_lm(i) values oq1T:if the output of your function is periodic then you should specify the periodicity of the functionTRANSPOSECalculate the transpose of a matrix More detailsARG=oq1 # Calculate a second contact matrix. This will contain the \sigma(r_ij) values in the second equation above. cmat2:the input for this action is the scalar output from one or more other actionsCONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Calculate the matrix product, S, described above s:specify the switching function to use between two sets of indistinguishable atomsMATRIX_PRODUCTCalculate the product of two matrices More detailsARG=oq1,oq1T # Now take the element-wise product of the contact matrix and the matrix of inner products. f:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=cmat2,sthe input to this functionFUNC=x*ythe function you wish to evaluatePERIODIC=NO # And output the 100 x 100 matrixif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=fthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# These lines compute the matrix F cmat:CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 q1:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO stack:if the output of your function is periodic then you should specify the periodicity of the functionVSTACKCreate a matrix by stacking vectors together More detailsARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 lones:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=q1,lones oq1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=stack,unormthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO oq1T:if the output of your function is periodic then you should specify the periodicity of the functionTRANSPOSECalculate the transpose of a matrix More detailsARG=oq1 cmat2:the input for this action is the scalar output from one or more other actionsCONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={RATIONAL D_0=2.0 R_0=1.0} s:specify the switching function to use between two sets of indistinguishable atomsMATRIX_PRODUCTCalculate the product of two matrices More detailsARG=oq1,oq1T f:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=cmat2,sthe input to this functionFUNC=x*ythe function you wish to evaluatePERIODIC=NO # Add all the elemnets of the matrix f together fs:if the output of your function is periodic then you should specify the periodicity of the functionSUMCalculate the sum of the arguments More detailsARG=fthe input to this functionPERIODIC=NO # This outputs a single scalarif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=fsthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# These lines compute the matrix F cmat:CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 q1:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO stack:if the output of your function is periodic then you should specify the periodicity of the functionVSTACKCreate a matrix by stacking vectors together More detailsARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 lones:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=q1,lones oq1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=stack,unormthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO oq1T:if the output of your function is periodic then you should specify the periodicity of the functionTRANSPOSECalculate the transpose of a matrix More detailsARG=oq1 cmat2:the input for this action is the scalar output from one or more other actionsCONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={RATIONAL D_0=2.0 R_0=1.0} s:specify the switching function to use between two sets of indistinguishable atomsMATRIX_PRODUCTCalculate the product of two matrices More detailsARG=oq1,oq1T f:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=cmat2,sthe input to this functionFUNC=x*ythe function you wish to evaluatePERIODIC=NO # Tranform all the elements of f with a switching function that is one when f_ij>0.5 fmt:if the output of your function is periodic then you should specify the periodicity of the functionMORE_THANUse a switching function to determine how many of the input variables are more than a certain cutoff. More detailsARG=fthe input to this functionSWITCH={RATIONAL R_0=0.5} # Add all the elemnets of the matrix fmt together fs:This keyword is used if you want to employ an alternative to the continuous swiching function defined aboveSUMCalculate the sum of the arguments More detailsARG=fmtthe input to this functionPERIODIC=NO # This outputs a single scalarif the output of your function is periodic then you should specify the periodicity of the functionPrint quantities to a file. More detailsARG=fsthe input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
cmat:CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableCOMPONENTSalso calculate the components of the vector connecting the atoms in the contact matrixSWITCH={RATIONAL D_0=2.0 R_0=1.0} ones:specify the switching function to use between two sets of indistinguishable atomsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=100 coord:the number of ones that you would like to createMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=cmat.w,ones sh:the input for this action is the scalar output from one or more other actionsSPHERICAL_HARMONICCalculate the values of all the spherical harmonic funtions for a particular value of l. More detailsARG=cmat.x,cmat.y,cmat.z,cmat.wthe input to this functionL=1 sp:the value of the angular momentumMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=sh.*,ones q1_2:the input for this action is the scalar output from one or more other actionsCOMBINECalculate a polynomial combination of a set of other variables. More detailsPERIODIC=NOif the output of your function is periodic then you should specify the periodicity of the functionPOWERS=2,2,2,2,2,2the powers to which you are raising each of the arguments in your functionARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 q1:the input to this functionCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=q1_2the input to this functionFUNC=sqrt(xthe function you wish to evaluatePERIODIC=NO stack:if the output of your function is periodic then you should specify the periodicity of the functionVSTACKCreate a matrix by stacking vectors together More detailsARG=sp.rm-n1,sp.im-n1,sp.rm-0,sp.im-0,sp.rm-p1,sp.im-p1 lones:the input for this action is the scalar output from one or more other actionsONESCreate a constant vector with all elements equal to one This action is a shortcut. More detailsSIZE=6 unorm:the number of ones that you would like to createOUTER_PRODUCTCalculate the outer product matrix of two vectors This action has hidden defaults. More detailsARG=q1,lones oq1:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=stack,unormthe input to this functionFUNC=x/ythe function you wish to evaluatePERIODIC=NO oq1T:if the output of your function is periodic then you should specify the periodicity of the functionTRANSPOSECalculate the transpose of a matrix More detailsARG=oq1 cmat2:the input for this action is the scalar output from one or more other actionsCONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More detailsGROUP=1-100specifies the list of atoms that should be assumed indistinguishableSWITCH={RATIONAL D_0=2.0 R_0=1.0} s:specify the switching function to use between two sets of indistinguishable atomsMATRIX_PRODUCTCalculate the product of two matrices More detailsARG=oq1,oq1T f:the input for this action is the scalar output from one or more other actionsCUSTOMCalculate a combination of variables using a custom expression. More detailsARG=cmat2,sthe input to this functionFUNC=x*ythe function you wish to evaluatePERIODIC=NO # This outputs a vector with 100 elements lq1:if the output of your function is periodic then you should specify the periodicity of the functionMATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More detailsARG=f,ones # This prints the 100 symmetry function values to the file colvarthe input for this action is the scalar output from one or more other actionsPrint quantities to a file. More detailsARG=lq1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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
|
|
# Calculate the steinhardt parameters q1:Q1Calculate 1st order Steinhardt parameters This action is a shortcut. More detailsSPECIES=1-100this keyword is used for colvars such as coordination numberSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Calculate the local steinhard parameters lq1:the switching function that it used in the construction of the contact matrixLOCAL_Q1SPECIES=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 detailsSWITCH={RATIONAL D_0=2.0 R_0=1.0} # Print the local steinhardt parametersThis keyword is used if you want to employ an alternative to the continuous swiching function defined abovePrint quantities to a file. More detailsARG=lq1the input for this action is the scalar output from one or more other actionsFILE=colvarthe name of the file on which to output these quantities
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.