diff --git a/src/symfunc/AngularTetra.cpp b/src/symfunc/AngularTetra.cpp
index 35fcd1dea6..220d935b82 100644
--- a/src/symfunc/AngularTetra.cpp
+++ b/src/symfunc/AngularTetra.cpp
@@ -28,7 +28,27 @@
/*
Calculate the angular tetra CV
-\par Examples
+This shortcut calculates a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html). The particular function that is being
+evaluated for the coordination sphere here is as follows:
+
+$$
+s_i = 1 - \frac{3}{8}\sum_{j=2}^4 \sum_{k=1}^{j-1} \left( \cos(\theta_{jik} + \frac{1}{2} \right)^2
+$$
+
+In this expression the 4 atoms in the sums over $j$ and $k$ are the four atoms that are nearest to atom $i$. $\theta_{jik}$ is the angle between the vectors connecting
+atoms $i$ and $j$ and the vector connecting atoms $i$ and $k$. The CV is large if the four atoms nearest atom $i$ are arranged on the vertices of a regular tetrahedron
+and small otherwise. The following example shows how you can use this action to measure the degree of tetrahedral order in a system.
+
+```plumed
+# Calculate a vector that contains 64 values for the symmetry function.
+# Sum the elements of the vector and calculate the mean value on the atoms from this sum.
+acv: TETRA_ANGULAR SPECIES=1-64 SUM MEAN
+# Print out the positions of the 64 atoms for which the symmetry function was calculated
+# to an xyz file along with the values of the symmetry function
+DUMPATOMS ATOMS=1-64 ARG=acv FILE=mcolv.xyz
+# Print out the average value of the symmetry function and the sum of all the symmetry functions
+PRINT ARG=acv_sum,acv_mean FILE=colvar
+```
*/
//+ENDPLUMEDOC
diff --git a/src/symfunc/AtomicSMAC.cpp b/src/symfunc/AtomicSMAC.cpp
index 6a937c52b3..842066e19f 100644
--- a/src/symfunc/AtomicSMAC.cpp
+++ b/src/symfunc/AtomicSMAC.cpp
@@ -30,7 +30,30 @@
/*
Calculate the atomic smac CV
-\par Examples
+This shortcut offers another example of a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html).
+This action was inspired by [SMAC](SMAC.md) and the distribution of angles in the first coordination sphere around an atom to determine if the
+environment around the atom is ordered. For atom $i$ the symmetry function is calculated using:
+
+$$
+s_i = [ 1 - \gamma(c_i) ] \frac{ \sum_j \sum{k \ne j} \sigma(r_{ij})\sigma(r_{ik}) \sum_n G\left( \frac{ \theta_{jik} - \phi_n}{b_n}\right) }{\sum{k \ne j} \sigma(r_{ij})\sigma(r_{ik})} \qquad \textrm{where} \qquad c_i = \sum_j \sigma(r_{ij})
+$$
+
+In this expression $r_{ij}$ is the distance between atom $i$ and atom $j$ and $\sigma$ is a switching function that acts upon this distance. $c_i$ is thus the number of atoms that are within
+a certain cutoff of atom $i$ and $\gamma$ is another switching function that acts upon this quantity. This switching function ensures that the symmetry function is zero for atoms that are
+regions where the density is low. $\theta_{jik}$ is the angle between the vector connecting atoms $i$ and $j$ and the vector connecting atoms $i$ and $k$. This angle is the argument for the
+set of Gaussian kernel functions, $G$, that are centered on $\phi_n$ and that have bandwidths of $b_n$. The function above is thus determining if the angles between the bonds in the first coordination
+sphere around atom $i$ are similar to the $\phi_n$ values that have been specified by the user or not.
+
+The following example demonstrates how this symmetry function can be used in practise.
+
+```plumed
+smac: ATOMIC_SMAC SPECIES=1-64 KERNEL1={GAUSSIAN CENTER=pi/2 SIGMA=1.0} KERNEL2={GAUSSIAN CENTER=pi/4 SIGMA=1.0} SWITCH_COORD={EXP R_0=4.0} SWITCH={RATIONAL R_0=2.0 D_0=2.0} SUM
+PRINT ARG=smac.* FILE=colvar
+```
+
+The input above would calculate 64 instances of $s_i$ using the formula above. In each of these two Gaussian Kernels are used in the sum over $n$. The parameters for the switching function
+$\sigma$ are specified using the SWITCH keyword, while the parameters for $\gamma$ are specified using SWITCH_COORD. As you can see if you expand the shortcut in the input above, the 64 values
+for $s_i$ are stored in a vector. All the elements of this vector are then added together to produce the single quantity that is output in the colvar file.
*/
//+ENDPLUMEDOC
@@ -48,12 +71,10 @@ PLUMED_REGISTER_ACTION(AtomicSMAC,"ATOMIC_SMAC")
void AtomicSMAC::registerKeywords(Keywords& keys) {
ActionShortcut::registerKeywords( keys );
- keys.add("optional","SPECIES","");
- keys.add("optional","SPECIESA","");
- keys.add("optional","SPECIESB","");
- keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
- "The following provides information on the \\ref switchingfunction that are available. "
- "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
+ keys.add("atoms-3","SPECIES","the list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments");
+ keys.add("atoms-4","SPECIESA","the list of atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment.");
+ keys.add("atoms-4","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which the symmetry function is being calculated.");
+ keys.add("optional","SWITCH","the switching function that it used in the construction of the contact matrix");
keys.add("numbered","KERNEL","The kernels used in the function of the angle");
keys.add("optional","SWITCH_COORD","This keyword is used to define the coordination switching function.");
keys.reset_style("KERNEL","optional");
diff --git a/src/symfunc/CoordShellVectorFunction.cpp b/src/symfunc/CoordShellVectorFunction.cpp
index e8029af0b6..8d5b8b3b10 100644
--- a/src/symfunc/CoordShellVectorFunction.cpp
+++ b/src/symfunc/CoordShellVectorFunction.cpp
@@ -37,7 +37,57 @@ namespace symfunc {
/*
Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom
-\par Examples
+This shortcut allows you to calculate the sum for an arbitrary function of the bond vectors that connect an atom to each of its neighbours.
+In other words, this action allows you to compute the following:
+
+$$
+s_i = \sum_{i \ne j} \sigma(r_{ij}) f(x_{ij}, y_{ij}, z_{ij}, r_{ij}) )
+$$
+
+In this expression, $x_{ij}, y_{ij}, z_{ij}$ are the components of the vector connecting atoms $i$ and $j$ and $r_{ij}$ is the magnitude of this vector.
+$\sigma(r_{ij})$ is then a switching function that ensures that the aritrary function $f$ is only evaluated for if atom $j$ is within a certain cutoff distance
+of atom $i$.
+
+Below you can see a simple example that shows how this action can be used in practise.
+
+```plumed
+cv: COORDINATION_SHELL_FUNCTION SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
+PRINT ARG=cv FILE=colvar
+```
+
+The above input calculates 64 $s_i$ values - one $s_i$ values for each of the atoms specified using the SPECIES keyword. These 64 numbers are then output to a file.
+The function of x, y, z and r to be evaluated is specified using the FUNCTION keyword. Obviously, if your function does not depend on all four of these variables
+they can be excluded from your function.
+
+As discussed in [this paper](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.92.180102) it is sometimes useful to rotate the bond vectors before computing the
+arbitrary function $f$ in the above expression. To perform such rotations you use the PHI, THETA and PSI keywords. The $x_{ij}, y_{ij}$ and $z_{ij}$ values that enter $f$ in the
+expression above are then calculated as:
+
+$$
+\left(
+\begin{matrix}
+x_{ij} \\
+y_{ij} \\
+z_{ij}
+\end{matrix}
+\right) =
+\left(
+\begin{matrix}
+\cos(\psi)\cos(\phi) - \cos(\theta)\sin(\phi)\sin(\psi) & \cos(\psi)*\sin(\phi)+\cos(\theta)*\cos(\phi)*\sin(\psi) & \sin(\psi)*sin(\theta) \\
+-\sin(\psi)*\cos(\phi)-\cos(\theta)*\sin(\phi)*\cos(\psi) & -\sin(\psi)*\sin(\phi)+\cos(\theta)*\cos(\phi)*std::cos(\psi), & \cos(\psi)*\sin(\theta) \\
+\sin(\theta)*\sin(\phi) & \sin(\theta)*\cos(\phi) & \cos(\theta)
+\end{matrix}
+\left(
+\begin{matrix}
+x_{ij}' \\
+y_{ij}' \\
+z_{ij}'
+\end{matrix}
+\right)
+$$
+
+$x_{ij}', y_{ij}'$ and $z_{ij}'$ in this expression are the bond vectors that are calculated in the lab frame. The matrix in the above expression is thus just a
+[rotation matrix](https://en.wikipedia.org/wiki/Rotation_matrix) that converts the lab frame to some frame of interest.
*/
@@ -47,8 +97,30 @@ Calculate an arbitrary function of all the bond vectors in the first coordinatio
/*
Calculate an arbitrary function of all the bond vectors in the first coordination sphere of an atom and take an average
-\par Examples
+This shortcut allows you to calculate the average for an arbitrary function of the bond vectors that connect an atom to each of its neighbours.
+In other words, this action allows you to compute the following:
+
+$$
+s_i = \frac{\sum_{i \ne j} \sigma(r_{ij}) f(x_{ij}, y_{ij}, z_{ij}, r_{ij}) )}{\sum_{i \ne j} \sigma(r_{ij})}
+$$
+In this expression, $x_{ij}, y_{ij}, z_{ij}$ are the components of the vector connecting atoms $i$ and $j$ and $r_{ij}$ is the magnitude of this vector.
+$\sigma(r_{ij})$ is then a switching function that ensures that the aritrary function $f$ is only evaluated for if atom $j$ is within a certain cutoff distance
+of atom $i$.
+
+Below you can see a simple example that shows how this action can be used in practise.
+
+```plumed
+cv: COORDINATION_SHELL_AVERAGE SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} FUNCTION=((x+y+z)/r)^3+((x-y-z)/r)^3+((-x+y-z)/r)^3+((-x-y+z)/r)^3
+PRINT ARG=cv FILE=colvar
+```
+
+The above input calculates 64 $s_i$ values - one $s_i$ values for each of the atoms specified using the SPECIES keyword. These 64 numbers are then output to a file.
+The function of x, y, z and r to be evaluated is specified using the FUNCTION keyword. Obviously, if your function does not depend on all four of these variables
+they can be excluded from your function.
+
+Notice that you can you can rotate the bond vectors before computing the
+arbitrary function $f$ in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
*/
//+ENDPLUMEDOC
@@ -57,35 +129,38 @@ Calculate an arbitrary function of all the bond vectors in the first coordinatio
/*
Calculate whether or not the coordination spheres of atoms are arranged as they would be in a simple cubic structure.
-We can measure how similar the environment around atom \f$i\f$ is to a simple cubic structure is by evaluating
+This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
+which we can use to measure how similar the environment around atom $i$ is to a simple cubic structure. We perform this comparison by evaluating
the following quantity:
-\f[
+$$
s_i = \frac{ \sum_{i \ne j} \sigma(r_{ij}) \left[ \frac{ x_{ij}^4 + y_{ij}^4 + z_{ij}^4 }{r_{ij}^4} \right] }{ \sum_{i \ne j} \sigma(r_{ij}) }
-\f]
+$$
-In this expression \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$ are the \f$x\f$, \f$y\f$ and \f$z\f$ components of the vector connecting atom \f$i\f$ to
-atom \f$j\f$ and \f$r_{ij}\f$ is the magnitude of this vector. \f$\sigma(r_{ij})\f$ is a \ref switchingfunction that acts on the distance between atom \f$i\f$ and atom \f$j\f$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
+In this expression $x_{ij}$, $y_{ij}$ and $z_{ij}$ are the $x$, $y$ and $z$ components of the vector connecting atom $i$ to
+atom $j$ and $r_{ij}$ is the magnitude of this vector. $\sigma(r_{ij})$ is a switching function that acts on the distance between atom $i$ and atom $j$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
over all of the other atoms in the system ensures that we are calculating an average
-of the function of \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$ for the atoms in the first coordination sphere around atom \f$i\f$.
+of the function of $x_{ij}$, $y_{ij}$ and $z_{ij}$ for the atoms in the first coordination sphere around atom $i$.
This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
-the average value for the atoms in your system, the number of atoms that have an \f$s_i\f$ value that is more that some target and
+the average value for the atoms in your system, the number of atoms that have an $s_i$ value that is more that some target and
so on. Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
-
-\par Examples
-
The following input tells plumed to calculate the simple cubic parameter for the atoms 1-100 with themselves.
The mean value is then calculated.
-\plumedfile
+
+```plumed
SIMPLECUBIC SPECIES=1-100 R_0=1.0 MEAN
-\endplumedfile
+```
The following input tells plumed to look at the ways atoms 1-100 are within 3.0 are arranged about atoms
from 101-110. The number of simple cubic parameters that are greater than 0.8 is then output
-\plumedfile
+
+```plumed
SIMPLECUBIC SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=0.8 NN=6 MM=12 D_0=0}
-\endplumedfile
+```
+
+Notice that you can you can rotate the bond vectors before computing the
+function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
*/
//+ENDPLUMEDOC
@@ -94,40 +169,42 @@ SIMPLECUBIC SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=0.8
/*
Calculate the degree to which the environment about ions has a tetrahedral order.
-We can measure the degree to which the atoms in the first coordination shell around any atom, \f$i\f$ is
-is arranged like a tetrahedron using the following function.
+This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
+which we can use to measure the degree to which the atoms in the first coordination shell around any atom, $i$ is
+is arranged like a tetrahedron. We perform this comparison by evaluating the following function.
-\f[
+$$
s(i) = \frac{1}{\sum_j \sigma( r_{ij} )} \sum_j \sigma( r_{ij} )\left[ \frac{(x_{ij} + y_{ij} + z_{ij})^3}{r_{ij}^3} +
\frac{(x_{ij} - y_{ij} - z_{ij})^3}{r_{ij}^3} +
\frac{(-x_{ij} + y_{ij} - z_{ij})^3}{r_{ij}^3} +
\frac{(-x_{ij} - y_{ij} + z_{ij})^3}{r_{ij}^3} \right]
-\f]
+$$
-Here \f$r_{ij}\f$ is the magnitude of the vector connecting atom \f$i\f$ to atom \f$j\f$ and \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$
-are its three components. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
-atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that the function is equal to one
-when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
-
-\par Examples
+Here $r_{ij}$ is the magnitude of the vector connecting atom $i$ to atom $j$ and $x_{ij}$, $y_{ij}$ and $z_{ij}$
+are its three components. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
+atoms $i$ and $j$. The parameters of this function should be set so that the function is equal to one
+when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
The following command calculates the average value of the TETRAHEDRAL parameter for a set of 64 atoms all of the same type
and outputs this quantity to a file called colvar.
-\plumedfile
+```plumed
tt: TETRAHEDRAL SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=tt.mean FILE=colvar
-\endplumedfile
+```
The following command calculates the number of TETRAHEDRAL parameters that are greater than 0.8 in a set of 10 atoms.
In this calculation it is assumed that there are two atom types A and B and that the first coordination sphere of the
10 atoms of type A contains atoms of type B. The formula above is thus calculated for ten different A atoms and within
-it the sum over \f$j\f$ runs over 40 atoms of type B that could be in the first coordination sphere.
+it the sum over $j$ runs over 40 atoms of type B that could be in the first coordination sphere.
-\plumedfile
+```plumed
tt: TETRAHEDRAL SPECIESA=1-10 SPECIESB=11-40 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=0.8}
PRINT ARG=tt.* FILE=colvar
-\endplumedfile
+```
+
+Notice that you can you can rotate the bond vectors before computing the
+function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
*/
//+ENDPLUMEDOC
@@ -155,6 +232,8 @@ void CoordShellVectorFunction::registerKeywords( Keywords& keys ) {
keys.setValueDescription("vector","the symmetry function for each of the specified atoms");
keys.needsAction("CONTACT_MATRIX"); keys.needsAction("FCCUBIC_FUNC"); keys.needsAction("CUSTOM");
keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
+ keys.addDOI("10.1103/PhysRevB.81.125416"); keys.addDOI("10.1103/PhysRevB.92.180102");
+ keys.addDOI("10.1063/1.4997180"); keys.addDOI("10.1063/1.5134461");
}
CoordShellVectorFunction::CoordShellVectorFunction(const ActionOptions& ao):
diff --git a/src/symfunc/CoordinationNumbers.cpp b/src/symfunc/CoordinationNumbers.cpp
index ef664f9603..d8cfd05292 100644
--- a/src/symfunc/CoordinationNumbers.cpp
+++ b/src/symfunc/CoordinationNumbers.cpp
@@ -37,43 +37,50 @@ namespace symfunc {
Calculate the coordination numbers of atoms so that you can then calculate functions of the distribution of
coordination numbers such as the minimum, the number less than a certain quantity and so on.
-So that the calculated coordination numbers have continuous derivatives the following function is used:
+The coordination number of a atom $i$ is the number of atoms that are within a certain cutoff distance of it.
+This quantity can be calculated as follows:
-\f[
-s = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
-\f]
+$$
+s_i = \sum_j \sigma(r_{ij})
+$$
-If R_POWER is set, this will use the product of pairwise distance
-raised to the R_POWER with the coordination number function defined
-above. This was used in White and Voth \cite white2014efficient as a
-way of indirectly biasing radial distribution functions. Note that in
-that reference this function is referred to as moments of coordination
-number, but here we call them powers to distinguish from the existing
-MOMENTS keyword of Multicolvars.
+where $r_{ij}$ is the distance between atoms $i$ and $j$ and $\sigma$ is a switching function. The typical switching
+function that is used in metadynamics is this one:
-\par Examples
+$$
+s(r) = \frac{ 1 - \left(\frac{r-d_0}{r_0}\right)^n } { 1 - \left(\frac{r-d_0}{r_0}\right)^m }
+$$
-The following input tells plumed to calculate the coordination numbers of atoms 1-100 with themselves.
-The minimum coordination number is then calculated.
-\plumedfile
-COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
-\endplumedfile
+The following example shows how you can use this shortcut action to calculate and print the coordination numbers of
+one hundred atoms with themselves:
-The following input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
+```plumed
+c: COORDINATIONNUMBER SPECIES=1-100 R_0=1.0
+DUMPXYZ ATOMS=c ARG=c FILE=coords.xyz
+```
+
+This input will produce an output file called coords that contains the coordination numbers of the 100 input atoms. The cutoff
+that is used to calculate the coordination number in this case is 1.0.
+
+The vectors that are output by the COORDINATIONNUMBER shortcut can be used in the input for many other functions that are within
+PLUMED. In addition, in order to ensure compatibility with older versions of PLUMED you can add additional keywords on the input
+line for COORDINATIONNUMBER in order to calculate various functions of the input. For example, the following input tells plumed
+ato calculate the coordination numbers of atoms 1-100 with themselves. The minimum coordination number is then calculated.
+
+```plumed
+c: COORDINATIONNUMBER SPECIES=1-100 R_0=1.0 MIN={BETA=0.1}
+```
+
+By constrast, this input tells plumed to calculate how many atoms from 1-100 are within 3.0 of each of the atoms
from 101-110. In the first 101 is the central atom, in the second 102 is the central atom and so on. The
-number of coordination numbers more than 6 is then computed.
-\plumedfile
-COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
-\endplumedfile
-
-The following input tells plumed to calculate the mean coordination number of all atoms with themselves
-and its powers. An explicit cutoff is set for each of 8.
-\plumedfile
-cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
-cn1: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
-cn2: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
-PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
-\endplumedfile
+number of coordination numbers that are more than 6 is then computed.
+
+```plumed
+c: COORDINATIONNUMBER SPECIESA=101-110 SPECIESB=1-100 R_0=3.0 MORE_THAN={RATIONAL R_0=6.0 NN=6 MM=12 D_0=0}
+```
+
+Notice that these inputs both use shortcuts. If you expand the inputs above you can determine the set of actions
+that are being used to calculate each of the quantities of interest.
*/
//+ENDPLUMEDOC
@@ -82,7 +89,25 @@ PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
/*
Calculate moments of the distribution of distances in the first coordination sphere
-\par Examples
+This is the CV that was developed by White and Voth and is described in the paper in the bibliograhy below. This action provides a way of indirectly biasing radial distribution functions and computes the following function
+
+$$
+s_i = \sum_j \sigma(r_{ij})r_{ij}^k
+$$
+
+where $k$ is the value that is input using the R_POWER keyword, $r_{ij}$ is the distance between atoms $i$ and $j$ and $\sigma$ is a switching function.
+
+The following example shows how this action can be used.
+
+```plumed
+cn0: COORDINATIONNUMBER SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} MEAN
+cn1: COORDINATION_MOMENTS SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=1 MEAN
+cn2: COORDINATION_MOMENTS SPECIES=1-10 SWITCH={RATIONAL R_0=1.0 D_MAX=8} R_POWER=2 MEAN
+PRINT ARG=cn0.mean,cn1.mean,cn2.mean STRIDE=1 FILE=cn_out
+```
+
+As you can see, it works similarlly to [COORDINATIONNUMBER](COORDINATIONNUMBER.md).
+
*/
//+ENDPLUMEDOC
@@ -93,18 +118,9 @@ PLUMED_REGISTER_ACTION(CoordinationNumbers,"COORDINATION_MOMENTS")
void CoordinationNumbers::shortcutKeywords( Keywords& keys ) {
ActionShortcut::registerKeywords( keys );
- keys.add("atoms-3","SPECIES","this keyword is used for colvars such as coordination number. In that context it specifies that plumed should calculate "
- "one coordination number for each of the atoms specified. Each of these coordination numbers specifies how many of the "
- "other specified atoms are within a certain cutoff of the central atom. You can specify the atoms here as another multicolvar "
- "action or using a MultiColvarFilter or ActionVolume action. When you do so the quantity is calculated for those atoms specified "
- "in the previous multicolvar. This is useful if you would like to calculate the Steinhardt parameter for those atoms that have a "
- "coordination number more than four for example");
- keys.add("atoms-4","SPECIESA","this keyword is used for colvars such as the coordination number. In that context it species that plumed should calculate "
- "one coordination number for each of the atoms specified in SPECIESA. Each of these cooordination numbers specifies how many "
- "of the atoms specifies using SPECIESB is within the specified cutoff. As with the species keyword the input can also be specified "
- "using the label of another multicolvar");
- keys.add("atoms-4","SPECIESB","this keyword is used for colvars such as the coordination number. It must appear with SPECIESA. For a full explanation see "
- "the documentation for that keyword");
+ keys.add("atoms-3","SPECIES","the list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments");
+ keys.add("atoms-4","SPECIESA","the list of atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESB, which specifies the atoms that are in the environment.");
+ keys.add("atoms-4","SPECIESB","the list of atoms that can be in the environments of each of the atoms for which the symmetry function is being calculated. This keyword must be used in conjunction with SPECIESA, which specifies the atoms for which the symmetry function is being calculated.");
keys.add("compulsory","NN","6","The n parameter of the switching function ");
keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
@@ -113,6 +129,7 @@ void CoordinationNumbers::shortcutKeywords( Keywords& keys ) {
keys.linkActionInDocs("SWITCH","LESS_THAN");
multicolvar::MultiColvarShortcuts::shortcutKeywords( keys );
keys.needsAction("CONTACT_MATRIX"); keys.needsAction("GROUP");
+ keys.addDOI("10.1021/ct500320c");
}
void CoordinationNumbers::expandMatrix( const bool& components, const std::string& lab, const std::string& sp_str,
diff --git a/src/symfunc/CylindricalHarmonic.cpp b/src/symfunc/CylindricalHarmonic.cpp
index 737dd9e058..78bf219ced 100644
--- a/src/symfunc/CylindricalHarmonic.cpp
+++ b/src/symfunc/CylindricalHarmonic.cpp
@@ -33,8 +33,41 @@ namespace symfunc {
/*
Calculate the cylindrical harmonic function
-\par Examples
-
+This action allows you to the value of the following complex function. The action outputs
+two components that are the real and imaginary parts of the following function:
+
+$$
+z = w (\frac{x}{r} + \frac{y}{r} i )^n \qquad \textrm{where} \qquad r = \sqrt(x^2 + y^2}
+$$
+
+In this expression $n$ is a parameter that is specified using the DEGREE keyword. $x$ and $y$ are the input arguments and $w$ is an optional input weight, which is set equal to
+one if only two arguments are provided in input. At present, the arguments for this action must be matrices.
+These arguments must all have the same shape as the two output components will also be matrices that are
+calculated by applying the function above to each of the elements of the input matrix in turn.
+
+The following intput provides an example that demonstrates how this function is used:
+
+```plumed
+d: DISTANCE_MATRIX GROUP=1-10 COMPONENTS
+c: CYLINDRICAL_HARMONIC DEGREE=6 ARG=d.x,d.y
+PRINT ARG=c.rm FILE=real_part
+PRINT ARG=c.im FILE=imaginary_part
+```
+
+The DISTANCE_MATRIX command in the above input computes 3 $10\times10$ matrices. Two of these $10\times10$ matrices are used in the input to the cylindrical harmonic command,
+which in turn outputs two $10\times10$ matrices that contain the real and imaginary parts when the function above is applied element-wise to the above input. These two $10\times10$
+matrices are then output to two separate files.
+
+In the above example the weights for every distance is set equal to one. The following example shows how an argument can be used to set the $w$ values to use when computing the function
+above.
+
+```plumed
+s: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0}
+sc: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0} COMPONENTS
+c: CYLINDRICAL_HARMONIC DEGREE=6 ARG=sc.x,sc.y,s
+PRINT ARG=c.rm FILE=real_part
+PRINT ARG=c.im FILE=imaginary_part
+```
*/
//+ENDPLUMEDOC
diff --git a/src/symfunc/Fccubic.cpp b/src/symfunc/Fccubic.cpp
index a3c2b08993..16ecd72e54 100644
--- a/src/symfunc/Fccubic.cpp
+++ b/src/symfunc/Fccubic.cpp
@@ -52,38 +52,40 @@ Measure how similar the environment around atoms is to that found in a FCC struc
/*
Measure how similar the environment around atoms is to that found in a FCC structure.
-This CV was introduced in this article \cite fcc-michele-1 and again in this article \cite fcc-michele-2
-This CV essentially determines whether the environment around any given atom is similar to that found in
-the FCC structure or not. The function that is used to make this determination is as follows:
+This shortcut is an example of a [COORDINATION_SHELL_AVERAGE](COORDINATION_SHELL_AVERAGE.md),
+which we can use to measure how similar the environment around atom $i$ is to an fcc structure.
+The function that is used to make this determination is as follows:
-\f[
+$$
s_i = \frac{ \sum_{i \ne j} \sigma(r_{ij}) \left\{ a\left[ \frac{(x_{ij}y_{ij})^4 + (x_{ij}z_{ij})^4 + (y_{ij}z_{ij})^4}{r_{ij}^8} - \frac{\alpha (x_{ij}y_{ij}z_{ij})^4}{r_{ij}^{12}} \right] + b \right\} }{ \sum_{i \ne j} \sigma(r_{ij}) }
-\f]
+$$
-In this expression \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$ are the \f$x\f$, \f$y\f$ and \f$z\f$ components of the vector connecting atom \f$i\f$ to
-atom \f$j\f$ and \f$r_{ij}\f$ is the magnitude of this vector. \f$\sigma(r_{ij})\f$ is a \ref switchingfunction that acts on the distance between
-atom \f$i\f$ and atom \f$j\f$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
+
+In this expression $x_{ij}$, $y_{ij}$ and $z_{ij}$ are the $x$, $y$ and $z$ components of the vector connecting atom $i$ to
+atom $j$ and $r_{ij}$ is the magnitude of this vector. $\sigma(r_{ij})$ is a switching function that acts on the distance between
+atom $i$ and atom $j$ and its inclusion in the numerator and the denominator of the above expression as well as the fact that we are summing
over all of the other atoms in the system ensures that we are calculating an average
-of the function of \f$x_{ij}\f$, \f$y_{ij}\f$ and \f$z_{ij}\f$ for the atoms in the first coordination sphere around atom \f$i\f$. Lastly, \f$\alpha\f$
-is a parameter that can be set by the user, which by default is equal to three. The values of \f$a\f$ and \f$b\f$ are calculated from \f$\alpha\f$ using:
+of the function of $x_{ij}$, $y_{ij}$ and $z_{ij}$ for the atoms in the first coordination sphere around atom $i$. Lastly, $\alpha$
+is a parameter that can be set by the user, which by default is equal to three. The values of $a$ and $b$ are calculated from $\alpha$ using:
-\f[
+$$
a = \frac{ 80080}{ 2717 + 16 \alpha} \qquad \textrm{and} \qquad b = \frac{ 16(\alpha - 143) }{2717 + 16\alpha}
-\f]
-
-This quantity is once again a multicolvar so you can compute it for multiple atoms using a single PLUMED action and then compute
-the average value for the atoms in your system, the number of atoms that have an \f$s_i\f$ value that is more that some target and
-so on. Notice also that you can rotate the reference frame if you are using a non-standard unit cell.
+$$
-\par Examples
+This action was been used in all the articles in the bibliography. We thus wrote an explict
+action to calculate it in PLUMED instead of using a shortcut as we did for [SIMPLECUBIC](SIMPLECUBIC.md) so that we could get
+good computational performance.
The following input calculates the FCCUBIC parameter for the 64 atoms in the system
and then calculates and prints the average value for this quantity.
-\plumedfile
-FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN LABEL=d
+```plumed
+d: FCCUBIC SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} MEAN
PRINT ARG=d.* FILE=colv
-\endplumedfile
+```
+
+Notice that you can you can rotate the bond vectors before computing the
+function in the above expression as is discussed in the documentation for [COORDINATION_SHELL_FUNCTION](COORDINATION_SHELL_FUNCTION.md)
*/
//+ENDPLUMEDOC
diff --git a/src/symfunc/HexaticParameter.cpp b/src/symfunc/HexaticParameter.cpp
index 32716b1202..50ce9c52d8 100644
--- a/src/symfunc/HexaticParameter.cpp
+++ b/src/symfunc/HexaticParameter.cpp
@@ -36,9 +36,57 @@ namespace symfunc {
/*
Calculate the hexatic order parameter
-\bug Virial is not working currently
+This [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) can be used to understand phase transitions in two dimensional systems.
+The symmetry function for atom $k$ is calculated using:
-\par Examples
+$$
+s_k = \left| \frac{\sum_j \sigma(r_{kj}) e^{6i\theta_j} }{\sum_j \sigma(r_{kj})} \right|
+$$
+
+In this expression, the sum run over all the atoms of interest and $r_{kj}$ is the distance between atom $k$ and atom $j$ and $\sigma$ is
+a switching function that acts upon this distance. $\theta_j$ is the angle between either the $x$, $y$ or $z$ axis and the bond connecting
+atom $k$ and atom $j$. This angle is multiplied by the imaginary number $i$ - the square root of minus one. In the code, we thus calculate
+$e^{i\theta_j}$ as follows:
+
+$$
+e^{i\theta_j) = \frac{x_{kj}}{r_{kj}} + i \frac{y_{kj}}{r_{kj}}
+$$
+
+We then take the 6th power of this complex number directly before compupting the magnitude by multiplying the result by its complex conjugate.
+Notice, furthermore, that we can replace $x_{kj}$ or $y_{kj}$ with $z_{kj}$ by using PLANE=xz or PLANE=yz in place of PLANE=xy.
+
+An example that shows how you can use this shortcut is shown below:
+
+```plumed
+hex: HEXACTIC_PARAMETER SPECIES=1-400 PLANE=xy SWITCH={RATIONAL D_0=1.4 R_0=0.2} MEAN
+PRINT ARG=hex.mean FILE=colvar
+```
+
+As you can see if you expand the shortcut above, this input calculates the quantity defined in the equation above for the 400 atoms in the simulated system and stores them in a vector.
+The elements of this vector are then added together so the mean value can be computed.
+
+In papers where symmetry functions similar to this one have been used a switching function is not employed. The sums over $j$ in the expression above are replaced by sums over the
+six nearest neighbours to each atom. If you would like to calculate this quantity using PLUMED you can use an input like this:
+
+```plumed
+dmat: DISTANCE_MATRIX GROUP=1-400 CUTOFF=3.0 COMPONENTS
+neigh: NEIGHBORS ARG=dmat.w NLOWEST=6
+harm: CYLINDRICAL_HARMONIC DEGREE=6 ARG=dmat.x,dmat.y
+rprod: CUSTOM ARG=neigh,harm.rm FUNC=x*y PERIODIC=NO
+iprod: CUSTOM ARG=neigh,harm.im FUNC=x*y PERIODIC=NO
+hex2_ones: ONES SIZE=400
+hex2_denom: MATRIX_VECTOR_PRODUCT ARG=neigh,hex2_ones
+harm_rm: MATRIX_VECTOR_PRODUCT ARG=rprod,hex2_ones
+harm_im: MATRIX_VECTOR_PRODUCT ARG=iprod,hex2_ones
+hex2_rmn: CUSTOM ARG=harm_rm,hex2_denom FUNC=x/y PERIODIC=NO
+hex2_imn: CUSTOM ARG=harm_im,hex2_denom FUNC=x/y PERIODIC=NO
+DUMPATOMS ATOMS=1-400 ARG=hex2_rmn,hex2_imn,hex2_denom FILE=hexparam.xyz
+```
+
+This input outputs the values of the order parameters for all the atoms to an extended xyz file .
+
+> ![CAUTION]
+> Virial is not working currently
*/
@@ -66,6 +114,7 @@ void HexacticParameter::registerKeywords( Keywords& keys ) {
keys.needsAction("CYLINDRICAL_HARMONIC_MATRIX"); keys.needsAction("ONES");
keys.needsAction("MATRIX_VECTOR_PRODUCT"); keys.needsAction("CUSTOM");
keys.needsAction("MEAN"); keys.needsAction("SUM"); keys.needsAction("COMBINE");
+ keys.addDOI("https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.215701");
}
HexacticParameter::HexacticParameter( const ActionOptions& ao):
diff --git a/src/symfunc/LocalAverage.cpp b/src/symfunc/LocalAverage.cpp
index f97ad199d7..582ac70ac6 100644
--- a/src/symfunc/LocalAverage.cpp
+++ b/src/symfunc/LocalAverage.cpp
@@ -36,11 +36,11 @@ namespace symfunc {
/*
Calculate averages over spherical regions centered on atoms
-As is explained in this video certain multicolvars
+As is explained in this video certain PLUMED actions
calculate one scalar quantity or one vector for each of the atoms in the system. For example
-\ref COORDINATIONNUMBER measures the coordination number of each of the atoms in the system and \ref Q4 measures
+[COORDINATIONNUMBER](COORDINATIONNUMBER.md) measures the coordination number of each of the atoms in the system and [Q4](Q4.md) measures
the 4th order Steinhardt parameter for each of the atoms in the system. These quantities provide tell us something about
-the disposition of the atoms in the first coordination sphere of each of the atoms of interest. Lechner and Dellago \cite dellago-q6
+the disposition of the atoms in the first coordination sphere of each of the atoms of interest. In the paper in the bibliography Lechner and Dellago
have suggested that one can probe local order in a system by taking the average value of such symmetry functions over
the atoms within a spherical cutoff of each of these atoms in the systems. When this is done with Steinhardt parameters
they claim this gives a coordinate that is better able to distinguish solid and liquid configurations of Lennard-Jones atoms.
@@ -48,40 +48,40 @@ they claim this gives a coordinate that is better able to distinguish solid and
You can calculate such locally averaged quantities within plumed by using the LOCAL_AVERAGE command. This command calculates
the following atom-centered quantities:
-\f[
+$$
s_i = \frac{ c_i + \sum_j \sigma(r_{ij})c_j }{ 1 + \sum_j \sigma(r_{ij}) }
-\f]
+$$
-where the \f$c_i\f$ and \f$c_j\f$ values can be for any one of the symmetry functions that can be calculated using plumed
-multicolvars. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
-atoms \f$i\f$ and \f$j\f$. Lechner and Dellago suggest that the parameters of this function should be set so that it the function is equal to one
-when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
+where the $c_i$ and $c_j$ values can be any vector of [symmetry functions](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html)
+ that can be calculated using plumed multicolvars. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
+atoms $i$ and $j$. Lechner and Dellago suggest that the parameters of this function should be set so that it the function is equal to one
+when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
-The \f$s_i\f$ quantities calculated using the above command can be again thought of as atom-centred symmetry functions. They
-thus operate much like multicolvars. You can thus calculate properties of the distribution of \f$s_i\f$ values using MEAN, LESS_THAN, HISTOGRAM
-and so on. You can also probe the value of these averaged variables in regions of the box by using the command in tandem with the
-\ref AROUND command.
+To see how this works in practice consider the following example input.
-\par Examples
-
-This example input calculates the coordination numbers for all the atoms in the system. These coordination numbers are then averaged over
-spherical regions. The number of averaged coordination numbers that are greater than 4 is then output to a file.
-
-\plumedfile
-COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=d1
-LOCAL_AVERAGE ARG=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4} LABEL=la
+```plumed
+d1: COORDINATIONNUMBER SPECIES=1-64 D_0=1.3 R_0=0.2
+la: LOCAL_AVERAGE ARG=d1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MORE_THAN={RATIONAL R_0=4}
PRINT ARG=la.* FILE=colvar
-\endplumedfile
+```
-This example input calculates the \f$q_4\f$ (see \ref Q4) vectors for each of the atoms in the system. These vectors are then averaged
-component by component over a spherical region. The average value for this quantity is then outputeed to a file. This calculates the
-quantities that were used in the paper by Lechner and Dellago \cite dellago-q6
+This input calculates the coordination numbers for all the atoms in the system. These coordination numbers are then averaged over
+spherical regions. The number of averaged coordination numbers that are greater than 4 is then output to a file. Furthermore, if you
+expand the input above you can see how the LOCAL_AVERAGE command is a shortcut action that expands to a longer input that you should be able to
+interpret.
-\plumedfile
+What Lechner and Dellago did in their paper was a little more complicated than this first example. To reproduce what they did you would use
+an input something like this:
+
+```plumed
Q4 SPECIES=1-64 SWITCH={RATIONAL D_0=1.3 R_0=0.2} LABEL=q4
LOCAL_AVERAGE ARG=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=la
PRINT ARG=la.* FILE=colvar
-\endplumedfile
+```
+
+This example input calculates the [Q4](Q4.md) vectors for each of the atoms in the system. These vectors are then averaged
+component by component over a spherical region. The average value for this quantity is then outputeed to a file. If you want
+to understand more about the shortcut that is used here you can read [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html).
*/
//+ENDPLUMEDOC
@@ -100,7 +100,7 @@ void LocalAverage::registerKeywords( Keywords& keys ) {
CoordinationNumbers::shortcutKeywords( keys );
keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
keys.needsAction("VSTACK"); keys.needsAction("CUSTOM"); keys.needsAction("OUTER_PRODUCT");
- keys.setValueDescription("vector","the values of the local averages");
+ keys.setValueDescription("vector","the values of the local averages"); keys.addDOI("10.1063/1.2977970");
}
LocalAverage::LocalAverage(const ActionOptions&ao):
diff --git a/src/symfunc/LocalCrystalinity.cpp b/src/symfunc/LocalCrystalinity.cpp
index e910481c24..852dd6eec7 100644
--- a/src/symfunc/LocalCrystalinity.cpp
+++ b/src/symfunc/LocalCrystalinity.cpp
@@ -36,8 +36,25 @@ namespace symfunc {
/*
Calculate the local crystalinity symmetry function
-\par Examples
+This shortcut provides an implementation of the local crystalinity order parameter that is described in the paper in the bibliography.
+To use this CV you define a series of unit vectors, $g_k$, using multiple instances of the GVECTOR keyword. This allows you to define a
+[symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) for the $i$th atom as:
+$$
+s_i = \sum_k \left| \frac{\sum_j \sigma(|r_{ij}|) e^{ig_k r_{ij}}}{\sum_j \sigma(|r_{ij}|) \right|^2
+$$
+
+In this expression $r_{ij}$ is the vector connecting atom $i$ to atom $j$ and $\sigma$ is a switching function that acts upon the magnidue of this vector, $|r_{ij}|$.
+The following input is an example that shows how this symmetry function can be used in practice.
+
+```plumed
+lc: LOCAL_CRYSTALINITY SPECIES=1-64 SWITCH={RATIONAL D_0=3.0 R_0=1.5} GVECTOR1=1,1,1 GVECTOR2=1,0.5,0.5 GVECTOR3=0.5,1.0,1.0 SUM
+PRINT ARG=lc_sum FILE=colvar
+```
+
+As you can see if you expand the shortcut in this input, the sum over $k$ in the above expression has three terms in this input as 3 GVECTORS are specified.
+Sixty four values for the expression above are computed. These sixty four numbers are then added together in order to give a global mesuare of the crystallinity
+for the simulated system.
*/
//+ENDPLUMEDOC
@@ -57,6 +74,7 @@ void LocalCrystallinity::registerKeywords( Keywords& keys ) {
keys.setValueDescription("vector","the value of the local crystalinity for each of the input atoms");
keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT");
keys.needsAction("COMBINE"); keys.needsAction("CUSTOM");
+ keys.addDOI("10.1063/1.4822877");
}
LocalCrystallinity::LocalCrystallinity( const ActionOptions& ao):
diff --git a/src/symfunc/LocalSteinhardt.cpp b/src/symfunc/LocalSteinhardt.cpp
index 643f07871d..2663c57817 100644
--- a/src/symfunc/LocalSteinhardt.cpp
+++ b/src/symfunc/LocalSteinhardt.cpp
@@ -33,7 +33,76 @@ namespace symfunc {
/*
Calculate 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.
-\par Examples
+The [Q1](Q1.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
+around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
+measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
+examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
+atoms in the surrounding liquid will remain pretty much the same. As such if one models a small nucleus embedded in a very large amount of solution/melt any
+change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
+because the number of atoms is relatively small.
+
+When the average [Q1](Q1.md) parameter is used to bias the dynamics a problems
+can occur. These averaged coordinates cannot distinguish between the correct,
+single-nucleus pathway and a concerted pathway in which all the atoms rearrange
+themselves into their solid-like configuration simultaneously. This second type
+of pathway would be impossible in reality because there is a large entropic
+barrier that prevents concerted processes like this from happening. However,
+in the finite sized systems that are commonly simulated this barrier is reduced
+substantially. As a result in simulations where average Steinhardt parameters
+are biased there are often quite dramatic system size effects
+
+If one wants to simulate nucleation using some form on
+biased dynamics what is really required is an order parameter that measures:
+
+- Whether or not the coordination spheres around atoms are ordered
+- Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
+
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
+introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
+Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
+like this one to help you compute CVs that have been widely used in the literature easily.
+
+This particular shortcut allows you to compute the LOCAL_Q1 parameter for a particular atom, which is a number that measures the extent to
+which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
+quantity for each of the atoms in the system:
+
+$$
+ s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-1}^1 q_{1m}^{*}(i)q_{1m}(j) }{ \sum_j \sigma( r_{ij} ) }
+$$
+
+where $q_{1m}(i)$ and $q_{1m}(j)$ are the 1st order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
+conjugation. The function
+$\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
+so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
+of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
+adjacent atoms are correlated.
+
+The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q1 parameter for the 64 Lennard Jones atoms in the system under study and prints this
+quantity to a file called colvar.
+
+```plumed
+q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq1: LOCAL_Q1 SPECIES=q1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
+PRINT ARG=lq1.mean FILE=colvar
+```
+
+The following input calculates the distribution of LOCAL_Q1 parameters at any given time and outputs this information to a file.
+
+```plumed
+q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq1: LOCAL_Q1 SPECIES=q1 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
+PRINT ARG=lq1.* FILE=colvar
+```
+
+The following calculates the LOCAL_Q1 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
+are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
+
+```plumed
+q1a: Q1 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
+q1b: Q1 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
+w1: LOCAL_Q1 SPECIES=q1a,q1b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
+PRINT ARG=w1.* FILE=colvar
+```
*/
//+ENDPLUMEDOC
@@ -42,7 +111,7 @@ Calculate the local degree of order around an atoms by taking the average dot pr
/*
Calculate the local degree of order around an atoms by taking the average dot product between the q_3 vector on the central atom and the q_3 vector on the atoms in the first coordination sphere.
-The \ref Q3 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
+The [Q3](Q3.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
@@ -50,7 +119,7 @@ atoms in the surrounding liquid will remain pretty much the same. As such if one
change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
because the number of atoms is relatively small.
-When the average \ref Q3 parameter is used to bias the dynamics a problems
+When the average [Q3](Q3.md) parameter is used to bias the dynamics a problems
can occur. These averaged coordinates cannot distinguish between the correct,
single-nucleus pathway and a concerted pathway in which all the atoms rearrange
themselves into their solid-like configuration simultaneously. This second type
@@ -66,51 +135,52 @@ biased dynamics what is really required is an order parameter that measures:
- Whether or not the coordination spheres around atoms are ordered
- Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
-\ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
-LOCAL_Q3 is another variable that can be used in these sorts of calculations. The LOCAL_Q3 parameter for a particular atom is a number that measures the extent to
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
+introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
+Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
+like this one to help you compute CVs that have been widely used in the literature easily.
+
+This particular shortcut allows you to compute the LOCAL_Q3 parameter for a particular atom, which is a number that measures the extent to
which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
quantity for each of the atoms in the system:
-\f[
+$$
s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-3}^3 q_{3m}^{*}(i)q_{3m}(j) }{ \sum_j \sigma( r_{ij} ) }
-\f]
+$$
-where \f$q_{3m}(i)\f$ and \f$q_{3m}(j)\f$ are the 3rd order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes complex
+where $q_{3m}(i)$ and $q_{3m}(j)$ are the 3rd order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
conjugation. The function
-\f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set
-so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise. The sum in the numerator
-of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
-adjacent atoms is correlated.
-
-\par Examples
+$\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
+so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
+of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
+adjacent atoms are correlated.
-The following command calculates the average value of the LOCAL_Q3 parameter for the 64 Lennard Jones atoms in the system under study and prints this
+The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q3 parameter for the 64 Lennard Jones atoms in the system under study and prints this
quantity to a file called colvar.
-\plumedfile
-Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
-LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq3
+```plumed
+q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq3: LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=lq3.mean FILE=colvar
-\endplumedfile
+```
The following input calculates the distribution of LOCAL_Q3 parameters at any given time and outputs this information to a file.
-\plumedfile
-Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q3
-LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq3
+```plumed
+q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq3: LOCAL_Q3 SPECIES=q3 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
PRINT ARG=lq3.* FILE=colvar
-\endplumedfile
+```
The following calculates the LOCAL_Q3 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
-\plumedfile
-Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3a
-Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q3b
-
-LOCAL_Q3 SPECIES=q3a,q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w3
+```plumed
+q3a: Q3 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
+q3b: Q3 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
+w3: LOCAL_Q3 SPECIES=q3a,q3b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=w3.* FILE=colvar
-\endplumedfile
+```
*/
//+ENDPLUMEDOC
@@ -119,7 +189,7 @@ PRINT ARG=w3.* FILE=colvar
/*
Calculate the local degree of order around an atoms by taking the average dot product between the q_4 vector on the central atom and the q_4 vector on the atoms in the first coordination sphere.
-The \ref Q4 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
+The [Q4](Q4.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
@@ -127,7 +197,7 @@ atoms in the surrounding liquid will remain pretty much the same. As such if one
change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
because the number of atoms is relatively small.
-When the average \ref Q4 parameter is used to bias the dynamics a problems
+When the average [Q4](Q4.md) parameter is used to bias the dynamics a problems
can occur. These averaged coordinates cannot distinguish between the correct,
single-nucleus pathway and a concerted pathway in which all the atoms rearrange
themselves into their solid-like configuration simultaneously. This second type
@@ -143,51 +213,52 @@ biased dynamics what is really required is an order parameter that measures:
- Whether or not the coordination spheres around atoms are ordered
- Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
-\ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
-LOCAL_Q4 is another variable that can be used in these sorts of calculations. The LOCAL_Q4 parameter for a particular atom is a number that measures the extent to
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
+introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
+Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
+like this one to help you compute CVs that have been widely used in the literature easily.
+
+This particular shortcut allows you to compute the LOCAL_Q4 parameter for a particular atom, which is a number that measures the extent to
which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
quantity for each of the atoms in the system:
-\f[
+$$
s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-4}^4 q_{4m}^{*}(i)q_{4m}(j) }{ \sum_j \sigma( r_{ij} ) }
-\f]
-
-where \f$q_{4m}(i)\f$ and \f$q_{4m}(j)\f$ are the fourth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
-complex conjugation. The function
-\f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set
-so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise. The sum in the numerator
-of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
-adjacent atoms is correlated.
+$$
-\par Examples
+where $q_{4m}(i)$ and $q_{4m}(j)$ are the 4th order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
+conjugation. The function
+$\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
+so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
+of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
+adjacent atoms are correlated.
-The following command calculates the average value of the LOCAL_Q4 parameter for the 64 Lennard Jones atoms in the system under study and prints this
+The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q4 parameter for the 64 Lennard Jones atoms in the system under study and prints this
quantity to a file called colvar.
-\plumedfile
-Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
-LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq4
+```plumed
+q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq4: LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=lq4.mean FILE=colvar
-\endplumedfile
+```
The following input calculates the distribution of LOCAL_Q4 parameters at any given time and outputs this information to a file.
-\plumedfile
-Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q4
-LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq4
+```plumed
+q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq4: LOCAL_Q4 SPECIES=q4 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
PRINT ARG=lq4.* FILE=colvar
-\endplumedfile
+```
The following calculates the LOCAL_Q4 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
-\plumedfile
-Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4a
-Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q4b
-
-LOCAL_Q4 SPECIES=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w4
+```plumed
+q4a: Q4 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
+q4b: Q4 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
+w4: LOCAL_Q4 SPECIES=q4a,q4b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=w4.* FILE=colvar
-\endplumedfile
+```
*/
//+ENDPLUMEDOC
@@ -196,7 +267,7 @@ PRINT ARG=w4.* FILE=colvar
/*
Calculate the local degree of order around an atoms by taking the average dot product between the q_6 vector on the central atom and the q_6 vector on the atoms in the first coordination sphere.
-The \ref Q6 command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
+The [Q6](Q6.md) command allows one to calculate one complex vectors for each of the atoms in your system that describe the degree of order in the coordination sphere
around a particular atom. The difficulty with these vectors comes when combining the order parameters from all of the individual atoms/molecules so as to get a
measure of the global degree of order for the system. The simplest way of doing this - calculating the average Steinhardt parameter - can be problematic. If one is
examining nucleation say only the order parameters for those atoms in the nucleus will change significantly when the nucleus forms. The order parameters for the
@@ -204,7 +275,7 @@ atoms in the surrounding liquid will remain pretty much the same. As such if one
change in the average order parameter will be negligible. Substantial changes in the value of this average can be observed in simulations of nucleation but only
because the number of atoms is relatively small.
-When the average \ref Q6 parameter is used to bias the dynamics a problems
+When the average [Q6](Q6.md) parameter is used to bias the dynamics a problems
can occur. These averaged coordinates cannot distinguish between the correct,
single-nucleus pathway and a concerted pathway in which all the atoms rearrange
themselves into their solid-like configuration simultaneously. This second type
@@ -220,51 +291,52 @@ biased dynamics what is really required is an order parameter that measures:
- Whether or not the coordination spheres around atoms are ordered
- Whether or not the atoms that are ordered are clustered together in a crystalline nucleus
-\ref LOCAL_AVERAGE and \ref NLINKS are variables that can be combined with the Steinhardt parameters allow to calculate variables that satisfy these requirements.
-LOCAL_Q6 is another variable that can be used in these sorts of calculations. The LOCAL_Q6 parameter for a particular atom is a number that measures the extent to
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html) a variety of variations on the Steinhardt parameters have been
+introduced to better describe nucleation. That page also shows how PLUMED provides you with flexibility that you can use to implement new combinations of the
+Steinhardt parameters. However, the inputs that you would need to write to implement common symmetry functions are rather complex so we also provide shortcuts
+like this one to help you compute CVs that have been widely used in the literature easily.
+
+This particular shortcut allows you to compute the LOCAL_Q6 parameter for a particular atom, which is a number that measures the extent to
which the orientation of the atoms in the first coordination sphere of an atom match the orientation of the central atom. It does this by calculating the following
quantity for each of the atoms in the system:
-\f[
+$$
s_i = \frac{ \sum_j \sigma( r_{ij} ) \sum_{m=-6}^6 q_{6m}^{*}(i)q_{6m}(j) }{ \sum_j \sigma( r_{ij} ) }
-\f]
-
-where \f$q_{6m}(i)\f$ and \f$q_{6m}(j)\f$ are the sixth order Steinhardt vectors calculated for atom \f$i\f$ and atom \f$j\f$ respectively and the asterisk denotes
-complex conjugation. The function
-\f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set
-so that it the function is equal to one when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise. The sum in the numerator
-of this expression is the dot product of the Steinhardt parameters for atoms \f$i\f$ and \f$j\f$ and thus measures the degree to which the orientations of these
-adjacent atoms is correlated.
+$$
-\par Examples
+where $q_{6m}(i)$ and $q_{6m}(j)$ are the 6th order Steinhardt vectors calculated for atom $i$ and atom $j$ respectively and the asterisk denotes complex
+conjugation. The function
+$\sigma( r_{ij} )$ is a switching function that acts on the distance between atoms $i$ and $j$. The parameters of this function should be set
+so that it the function is equal to one when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise. The sum in the numerator
+of this expression is the dot product of the Steinhardt parameters for atoms $i$ and $j$ and thus measures the degree to which the orientations of these
+adjacent atoms are correlated.
-The following command calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
+The following input shows how this works in practice. This input calculates the average value of the LOCAL_Q6 parameter for the 64 Lennard Jones atoms in the system under study and prints this
quantity to a file called colvar.
-\plumedfile
-Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
-LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LABEL=lq6
+```plumed
+q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=lq6.mean FILE=colvar
-\endplumedfile
+```
The following input calculates the distribution of LOCAL_Q6 parameters at any given time and outputs this information to a file.
-\plumedfile
-Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 LABEL=q6
-LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=lq6
+```plumed
+q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2
+lq6: LOCAL_Q6 SPECIES=q6 SWITCH={RATIONAL D_0=1.3 R_0=0.2} HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
PRINT ARG=lq6.* FILE=colvar
-\endplumedfile
+```
The following calculates the LOCAL_Q6 parameters for atoms 1-5 only. For each of these atoms comparisons of the geometry of the coordination sphere
are done with those of all the other atoms in the system. The final quantity is the average and is outputted to a file
-\plumedfile
-Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6a
-Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2 LABEL=q6b
-
-LOCAL_Q6 SPECIES=q6a,q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN LOWMEM LABEL=w6
+```plumed
+q6a: Q6 SPECIESA=1-5 SPECIESB=1-64 D_0=1.3 R_0=0.2
+q6b: Q6 SPECIESA=6-64 SPECIESB=1-64 D_0=1.3 R_0=0.2
+w6: LOCAL_Q6 SPECIES=q6a,q6b SWITCH={RATIONAL D_0=1.3 R_0=0.2} MEAN
PRINT ARG=w6.* FILE=colvar
-\endplumedfile
+```
*/
//+ENDPLUMEDOC
@@ -285,9 +357,9 @@ PLUMED_REGISTER_ACTION(LocalSteinhardt,"LOCAL_Q6")
void LocalSteinhardt::registerKeywords( Keywords& keys ) {
ActionShortcut::registerKeywords( keys );
- keys.add("optional","SPECIES","");
- keys.add("optional","SPECIESA","");
- keys.add("optional","SPECIESB","");
+ keys.add("optional","SPECIES","the label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters");
+ keys.add("optional","SPECIESA","the label of the action that computes the Steinhardt parameters for which you would like to calculate local steinhardt parameters");
+ keys.add("optional","SPECIESB","the label of the action that computes the Steinhardt parameters that you would like to use when calculating the loal steinhardt parameters");
keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
"The following provides information on the \\ref switchingfunction that are available. "
"When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
diff --git a/src/symfunc/RadialTetra.cpp b/src/symfunc/RadialTetra.cpp
index 01798963a9..de2504b81d 100644
--- a/src/symfunc/RadialTetra.cpp
+++ b/src/symfunc/RadialTetra.cpp
@@ -31,7 +31,27 @@
/*
Calculate the radial tetra CV
-\par Examples
+This shortcut calculates a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html). The particular function that is being
+evaluated for the coordination sphere here is as follows:
+
+$$
+s_i = 1 - \frac{\sum_{j=1}^4 r_{ij}^2 - z_i\sum_{j=1}^4 r_{ij}}{12 z_i^2} \qquad \textrm{where} \qquad z_i = \frac{1}{4} \sum_{j=1}^4 r_{ij}
+$$
+
+In this expression the 4 atoms in the sums over $j$ are the four atoms that are nearest to atom $i$ and $r_{ij}$ is the distance between atoms $i$ and $j$.
+The CV is large if the four atoms nearest atom $i$ are arranged on the vertices of a regular tetrahedron
+and small otherwise. The following example shows how you can use this action to measure the degree of tetrahedral order in a system.
+
+```plumed
+# Calculate a vector that contains 64 values for the symmetry function.
+# Sum the elements of the vector and calculate the mean value on the atoms from this sum.
+acv: TETRA_RADIAL SPECIES=1-64 SUM MEAN
+# Print out the positions of the 64 atoms for which the symmetry function was calculated
+# to an xyz file along with the values of the symmetry function
+DUMPATOMS ATOMS=1-64 ARG=acv FILE=mcolv.xyz
+# Print out the average value of the symmetry function and the sum of all the symmetry functions
+PRINT ARG=acv_sum,acv_mean FILE=colvar
+```
*/
//+ENDPLUMEDOC
diff --git a/src/symfunc/SMAC.cpp b/src/symfunc/SMAC.cpp
index 989207333a..ad8cbc174a 100644
--- a/src/symfunc/SMAC.cpp
+++ b/src/symfunc/SMAC.cpp
@@ -30,8 +30,50 @@
/*
Calculate the SMAC order parameter for a set of molecules
-\par Examples
+This shortcut action provides a quantity that can be thought of as a [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) for molecules.
+You can thus use it to detect whether the molecules in a sample are ordered in the way that they would expect to be ordered in the crystal. If you are using this CV you normally use
+the [COM](COM.md) command to define centers of mass for your molecules and the [DISTANCE](DISTANCE.md) or [PLANE](PLANE.md) command to define each molecules orientation. You can thus calculate
+the relative orientation molecules $i$ and $j$ by calculating the torsional angle, $\theta_{ij}$, between the two molecular orientations around the vector $r_{ij}$ connecting the centers of
+mass of the two molecules. The value of the molecular order parameter for molecule $i$ is thus calculated as:
+$$
+s_i = frac{1-\gamma(c_i)}{c_i} \sum_j \sigma(|r_{ij}|) \sum_k K\left( \frac{\theta_{ij} - \phi_k}{b} \right) \qquad \textrm{where} \qquad c_i = \sum_j \sigma(|r_{ij}|)
+$$
+
+In this expression $\sigma$ is a switching function that acts on the distance between the centers of mass of molecules $i$ and $j$. $c_i$ is thus the number of molecules that are within
+a certain cutoff of molecule $i$ and $\gamma$ is another switching function that acts upon this quantity. This switching function ensures that the symmetry function is zero for atoms that are
+regions where the density of molecules are low. $K$ is then a kernel function with bandwidth $b$ that is centered at $\phi_k$. The function above is thus only large if molecule $i$ is surrounded
+by molecules whose relative orientations are as the user has requested by specifying $\phi_k$ parameters.
+
+The following example illustrates how the SMAC order parameter in PLUMED is used:
+
+```plumed
+m3: DISTANCES ...
+ ATOMS1=9,10 LOCATION1=9
+ ATOMS2=89,90 LOCATION2=89
+ ATOMS3=473,474 LOCATION3=473
+ ATOMS4=1161,1162 LOCATION4=1161
+ ATOMS5=1521,1522 LOCATION5=1521
+ ATOMS6=1593,1594 LOCATION6=1593
+ ATOMS7=1601,1602 LOCATION7=1601
+ ATOMS8=2201,2202 LOCATION8=2201
+ COMPONENTS
+...
+
+s2: SMAC SPECIES=m3 KERNEL1={GAUSSIAN CENTER=0 SIGMA=0.480} KERNEL2={GAUSSIAN CENTER=pi SIGMA=0.480} SWITCH={RATIONAL R_0=0.6} MORE_THAN={RATIONAL R_0=0.7} SWITCH_COORD={EXP R_0=4}
+
+PRINT ARG=s2_morethan FILE=colvar
+```
+
+Here the orientations of the molecules are specified by calculating the vectors connecting pairs of atoms in the molecules. The LOCATION keywords in the distance command are used to specify the
+positions of the molecules from which the $r_{ij}$ vectors in the above expression are calculated. The sum over $k$ in the above expression has two terms corresponding to the two Gaussian kernels
+that have been specified using KERNEL keywords. The SWITCH keyword has been used to specify the parameters of the switching function $\sigma$, while the SWITCH_COORD keyword has been used to specify
+the parameters of the switching function $\gamma$.
+
+A vector of 8 values for $s_i$ is calculated. The elements of this vector are transformed by a switching function that is one if the $s_i$ value is larger than 0.7. The eight elements of the resulting vector
+of transformed $s_i$ are then added together to determine the final quantity that is output to the colvar file.
+
+Incidentally, the authors who designed the SMAC symmetry function have forgotten what the letters in this acronym stand for.
*/
//+ENDPLUMEDOC
@@ -65,6 +107,7 @@ void SMAC::registerKeywords(Keywords& keys) {
keys.needsAction("VSTACK"); keys.needsAction("TRANSPOSE"); keys.needsAction("CONTACT_MATRIX");
keys.needsAction("TORSIONS_MATRIX"); keys.needsAction("COMBINE"); keys.needsAction("CUSTOM");
keys.needsAction("ONES"); keys.needsAction("MATRIX_VECTOR_PRODUCT"); keys.needsAction("MORE_THAN");
+ keys.addDOI("10.1016/j.ces.2014.08.032"); keys.addDOI("10.1021/acs.jctc.6b01073");
}
SMAC::SMAC(const ActionOptions& ao):
diff --git a/src/symfunc/SphericalHarmonic.cpp b/src/symfunc/SphericalHarmonic.cpp
index 8e1c60ed7a..5fa007194f 100644
--- a/src/symfunc/SphericalHarmonic.cpp
+++ b/src/symfunc/SphericalHarmonic.cpp
@@ -33,8 +33,58 @@ namespace symfunc {
/*
Calculate the values of all the spherical harmonic funtions for a particular value of l.
-\par Examples
-
+This action allows you to the all the [spherical harmonic](https://en.wikipedia.org/wiki/Spherical_harmonics) functions for a particular input
+$L$ value. As discussed in more detail in the article provided above the spherical harmonics this action calculates have the following general form:
+
+$$
+Y_l^m(\theta,\phi) = wNe^{im\phi} P_l^m(\cos\theta)
+$$
+
+where $N$ is a normalisation constant, $P_l^m$ is tn associated Legendre Polynomial and $w$ is an optional weight that can be passed in an input argument or simply set equal to one.
+$e^{i\phi}$ and $\cos(\theta) are computed from the other input arguments, $x, y$ and $z$ as follows:
+
+$$
+e^{i\phi} = \frac{x}{r} + i\frac{y}{r} \qquad \textrm{and} \qquad \cos(\theta) = \frac{z}{r} \qquad \textrm{where} \qquad r = \sqrt{x^2 + y^2 + z^2}
+$$
+
+At present, the arguments for this action must be matrices. However, it would be easy to add functionality that would allow you to compute this function for scalar or vector input.
+The arguments must all have the same shape as the two output components will also be matrices that are
+calculated by applying the function above to each of the elements of the input matrix in turn. The number of components output will be equal to $2(2L+1)$ and will contain
+the real and imaginary parts of the $Y_l^m$ functions with the the $2l+1$ possible $m$ values.
+
+The following intput provides an example that demonstrates how this function is used:
+
+```plumed
+d: DISTANCE_MATRIX GROUP=1-10 COMPONENTS
+c: SPHERICAL_HARMONIC L=1 ARG=d.x,d.y,d.z
+PRINT ARG=c.rm-n1 FILE=real_part_m-1
+PRINT ARG=c.im-n1 FILE=imaginary_part_m-1
+PRINT ARG=c.rm-0 FILE=real_part_m0
+PRINT ARG=c.im-0 FILE=imaginary_part_m0
+PRINT ARG=c.rm-p1 FILE=real_part_m+1
+PRINT ARG=c.im-p1 FILE=imaginary_part_m+1
+```
+
+The DISTANCE_MATRIX command in the above input computes 3 $10\times10$ matrices. These 3 $10\times10$ matrices are used in the input to the sphierical harmonic command,
+which in turn outputs 6 $10\times10$ matrices that contain the real and imaginary parts when the three spherical harmonic functions with $l=1$ are applied element-wise to the above input. These six $10\times10$
+matrices are then output to six separate files.
+
+In the above example the weights for every distance is set equal to one. The following example shows how an argument can be used to set the $w$ values to use when computing the function
+above.
+
+```plumed
+s: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0}
+sc: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0} COMPONENTS
+c: SPHERICAL_HARMONIC L=1 ARG=sc.x,sc.y,sc.z,s
+PRINT ARG=c.rm-n1 FILE=real_part_m-1
+PRINT ARG=c.im-n1 FILE=imaginary_part_m-1
+PRINT ARG=c.rm-0 FILE=real_part_m0
+PRINT ARG=c.im-0 FILE=imaginary_part_m0
+PRINT ARG=c.rm-p1 FILE=real_part_m+1
+PRINT ARG=c.im-p1 FILE=imaginary_part_m+1
+```
+
+This function is used in the calculation of the Steinhardt order parameters, which are described in detail [here](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html).
*/
//+ENDPLUMEDOC
diff --git a/src/symfunc/Steinhardt.cpp b/src/symfunc/Steinhardt.cpp
index 85b36c5abd..e21e6b4376 100644
--- a/src/symfunc/Steinhardt.cpp
+++ b/src/symfunc/Steinhardt.cpp
@@ -36,7 +36,64 @@ namespace symfunc {
/*
Calculate 1st order Steinhardt parameters
-\par Examples
+The 1st order Steinhardt parameters allow us to measure the degree to which the first coordination shell
+around an atom is ordered with the atoms aranged on a line. The Steinhardt parameter for atom, $i$ is complex vector whose components are
+calculated using the following formula:
+
+$$
+q_{1m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{1m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
+$$
+
+where $Y_{1m}$ is one of the 1st order spherical harmonics so $m$ is a number that runs from $-1$ to
+$+1$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
+atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
+when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
+
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
+be used to measure the degree of order in the system in a variety of different ways. The
+simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
+
+\f[
+Q_1(i) = \sqrt{ \sum_{m=-1}^1 q_{1m}(i)^{*} q_{1m}(i) }
+\f]
+
+This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
+the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
+that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
+
+```plumed
+q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
+PRINT ARG=q1.mean FILE=colvar
+```
+
+Another similar command is provided below. This time the histogram of Q1 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
+to a file called colvar:
+
+```plumed
+q1: Q1 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
+PRINT ARG=q1.* FILE=colvar
+```
+
+The command below could be used to measure the Q1 parameters that describe the arrangement of chlorine ions around the
+sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
+with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q1 parameter is calculated and output to a
+file called colvar
+
+```plumed
+q1: Q1 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
+PRINT ARG=q1.mean FILE=colvar
+```
+
+If you simply want to examine the values of the Q1 parameters for each of the atoms in your system you can do so by exploiting the
+command [DUMPATOMS](DUMPATOMS.xyz) as shown in the example below. The following output file will output a file in an extended xyz format
+called q1.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
+atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q1 parameter, columns
+6-12 will contain the real parts of the director of the $q_{1m}$ vector while columns 12-19 will contain the imaginary parts of this director.
+
+```plumed
+q1: Q1 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
+DUMPATOMS ATOMS=q1 ARG=q1_anorm FILE=q1.xyz
+```
*/
//+ENDPLUMEDOC
@@ -46,71 +103,63 @@ Calculate 1st order Steinhardt parameters
Calculate 3rd order Steinhardt parameters.
The 3rd order Steinhardt parameters allow us to measure the degree to which the first coordination shell
-around an atom is ordered. The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
+around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
calculated using the following formula:
-\f[
+$$
q_{3m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{3m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
-\f]
+$$
-where \f$Y_{3m}\f$ is one of the 3rd order spherical harmonics so \f$m\f$ is a number that runs from \f$-3\f$ to
-\f$+3\f$. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
-atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that it the function is equal to one
-when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
+where $Y_{3m}$ is one of the 3rd order spherical harmonics so $m$ is a number that runs from $-3$ to
+$+3$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
+atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
+when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
-The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways. The
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
+be used to measure the degree of order in the system in a variety of different ways. The
simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
\f[
Q_3(i) = \sqrt{ \sum_{m=-3}^3 q_{3m}(i)^{*} q_{3m}(i) }
\f]
-This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
-the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
-that is investigated.
-
-Other measures of order can be taken by averaging the components of the individual \f$q_3\f$ vectors individually or by taking dot products of
-the \f$q_{3}\f$ vectors on adjacent atoms. More information on these variables can be found in the documentation for \ref LOCAL_Q3,
-\ref LOCAL_AVERAGE and \ref NLINKS.
+This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
+the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
+that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
-\par Examples
-
-The following command calculates the average Q3 parameter for the 64 atoms in a box of Lennard Jones and prints this
-quantity to a file called colvar:
-
-\plumedfile
-Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q3
+```plumed
+q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
PRINT ARG=q3.mean FILE=colvar
-\endplumedfile
+```
-The following command calculates the histogram of Q3 parameters for the 64 atoms in a box of Lennard Jones and prints these
-quantities to a file called colvar:
+Another similar command is provided below. This time the histogram of Q3 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
+to a file called colvar:
-\plumedfile
-Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q3
+```plumed
+q3: Q3 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
PRINT ARG=q3.* FILE=colvar
-\endplumedfile
+```
-The following command could be used to measure the Q3 parameters that describe the arrangement of chlorine ions around the
+The command below could be used to measure the Q3 parameters that describe the arrangement of chlorine ions around the
sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
-with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions. Once again the average Q3 parameter is calculated and output to a
+with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q3 parameter is calculated and output to a
file called colvar
-\plumedfile
-Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q3
+```plumed
+q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
PRINT ARG=q3.mean FILE=colvar
-\endplumedfile
+```
If you simply want to examine the values of the Q3 parameters for each of the atoms in your system you can do so by exploiting the
-command \ref DUMPATOMS as shown in the example below. The following output file will output a file in an extended xyz format
+command [DUMPATOMS](DUMPATOMS.xyz) as shown in the example below. The following output file will output a file in an extended xyz format
called q3.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q3 parameter, columns
-6-12 will contain the real parts of the director of the \f$q_{3m}\f$ vector while columns 12-19 will contain the imaginary parts of this director.
+6-12 will contain the real parts of the director of the $q_{3m}$ vector while columns 12-19 will contain the imaginary parts of this director.
-\plumedfile
+```plumed
q3: Q3 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
DUMPATOMS ATOMS=q3 ARG=q3_anorm FILE=q3.xyz
-\endplumedfile
+```
*/
//+ENDPLUMEDOC
@@ -119,72 +168,64 @@ DUMPATOMS ATOMS=q3 ARG=q3_anorm FILE=q3.xyz
/*
Calculate fourth order Steinhardt parameters.
-The fourth order Steinhardt parameters allow us to measure the degree to which the first coordination shell
-around an atom is ordered. The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
+The 4th order Steinhardt parameters allow us to measure the degree to which the first coordination shell
+around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
calculated using the following formula:
-\f[
+$$
q_{4m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{4m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
-\f]
+$$
-where \f$Y_{4m}\f$ is one of the fourth order spherical harmonics so \f$m\f$ is a number that runs from \f$-4\f$ to
-\f$+4\f$. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
-atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that it the function is equal to one
-when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
+where $Y_{4m}$ is one of the 4th order spherical harmonics so $m$ is a number that runs from $-4$ to
+$+4$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
+atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
+when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
-The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways. The
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
+be used to measure the degree of order in the system in a variety of different ways. The
simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
\f[
Q_4(i) = \sqrt{ \sum_{m=-4}^4 q_{4m}(i)^{*} q_{4m}(i) }
\f]
-This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
-the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
-that is investigated.
+This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
+the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
+that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
-Other measures of order can be taken by averaging the components of the individual \f$q_4\f$ vectors individually or by taking dot products of
-the \f$q_{4}\f$ vectors on adjacent atoms. More information on these variables can be found in the documentation for \ref LOCAL_Q4,
-\ref LOCAL_AVERAGE and \ref NLINKS.
-
-\par Examples
-
-The following command calculates the average Q4 parameter for the 64 atoms in a box of Lennard Jones and prints this
-quantity to a file called colvar:
-
-\plumedfile
-Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q4
+```plumed
+q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
PRINT ARG=q4.mean FILE=colvar
-\endplumedfile
+```
-The following command calculates the histogram of Q4 parameters for the 64 atoms in a box of Lennard Jones and prints these
-quantities to a file called colvar:
+Another similar command is provided below. This time the histogram of Q4 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
+to a file called colvar:
-\plumedfile
-Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q4
+```plumed
+q4: Q4 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
PRINT ARG=q4.* FILE=colvar
-\endplumedfile
+```
-The following command could be used to measure the Q4 parameters that describe the arrangement of chlorine ions around the
+The command below could be used to measure the Q4 parameters that describe the arrangement of chlorine ions around the
sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
-with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions. Once again the average Q4 parameter is calculated and output to a
+with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q4 parameter is calculated and output to a
file called colvar
-\plumedfile
-Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q4
+```plumed
+q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
PRINT ARG=q4.mean FILE=colvar
-\endplumedfile
+```
If you simply want to examine the values of the Q4 parameters for each of the atoms in your system you can do so by exploiting the
-command \ref DUMPATOMS as shown in the example below. The following output file will output a file in an extended xyz format
-called q$.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
+command [DUMPATOMS](DUMPATOMS.xyz) as shown in the example below. The following output file will output a file in an extended xyz format
+called q4.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q4 parameter, columns
-6-15 will contain the real parts of the director of the \f$q_{6m}\f$ vector while columns 15-24 will contain the imaginary parts of this director.
+6-12 will contain the real parts of the director of the $q_{4m}$ vector while columns 12-19 will contain the imaginary parts of this director.
-\plumedfile
+```plumed
q4: Q4 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
DUMPATOMS ATOMS=q4 ARG=q4_anorm FILE=q4.xyz
-\endplumedfile
+```
*/
//+ENDPLUMEDOC
@@ -193,72 +234,64 @@ DUMPATOMS ATOMS=q4 ARG=q4_anorm FILE=q4.xyz
/*
Calculate sixth order Steinhardt parameters.
-The sixth order Steinhardt parameters allow us to measure the degree to which the first coordination shell
-around an atom is ordered. The Steinhardt parameter for atom, \f$i\f$ is complex vector whose components are
+The 6th order Steinhardt parameters allow us to measure the degree to which the first coordination shell
+around an atom is ordered. The Steinhardt parameter for atom, $i$ is complex vector whose components are
calculated using the following formula:
-\f[
+$$
q_{6m}(i) = \frac{\sum_j \sigma( r_{ij} ) Y_{6m}(\mathbf{r}_{ij}) }{\sum_j \sigma( r_{ij} ) }
-\f]
+$$
-where \f$Y_{6m}\f$ is one of the sixth order spherical harmonics so \f$m\f$ is a number that runs from \f$-6\f$ to
-\f$+6\f$. The function \f$\sigma( r_{ij} )\f$ is a \ref switchingfunction that acts on the distance between
-atoms \f$i\f$ and \f$j\f$. The parameters of this function should be set so that it the function is equal to one
-when atom \f$j\f$ is in the first coordination sphere of atom \f$i\f$ and is zero otherwise.
+where $Y_{6m}$ is one of the 6th order spherical harmonics so $m$ is a number that runs from $-6$ to
+$+6$. The function $\sigma( r_{ij} )$ is a switching function that acts on the distance between
+atoms $i$ and $j$. The parameters of this function should be set so that it the function is equal to one
+when atom $j$ is in the first coordination sphere of atom $i$ and is zero otherwise.
-The Steinhardt parameters can be used to measure the degree of order in the system in a variety of different ways. The
+As discussed on [this page](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html), the Steinhardt parameters can
+be used to measure the degree of order in the system in a variety of different ways. The
simplest way of measuring whether or not the coordination sphere is ordered is to simply take the norm of the above vector i.e.
\f[
Q_6(i) = \sqrt{ \sum_{m=-6}^6 q_{6m}(i)^{*} q_{6m}(i) }
\f]
-This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, when
-the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with this colvar it is the distribution of these normed quantities
-that is investigated.
-
-Other measures of order can be taken by averaging the components of the individual \f$q_6\f$ vectors individually or by taking dot products of
-the \f$q_{6}\f$ vectors on adjacent atoms. More information on these variables can be found in the documentation for \ref LOCAL_Q6,
-\ref LOCAL_AVERAGE and \ref NLINKS.
-
-\par Examples
-
-The following command calculates the average Q6 parameter for the 64 atoms in a box of Lennard Jones and prints this
-quantity to a file called colvar:
+This norm is small when the coordination shell is disordered and larger when the coordination shell is ordered. Furthermore, in inputs like
+the one shown below where the keywords LESS_THAN, MIN, MAX, HISTOGRAM, MEAN and so on are used with it is the distribution of these normed quantities
+that is investigated. You can investigate precisely what calculation is performed here by expaning the shortcuts in the input below.
-\plumedfile
-Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN LABEL=q6
+```plumed
+q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 MEAN
PRINT ARG=q6.mean FILE=colvar
-\endplumedfile
+```
-The following command calculates the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones and prints these
-quantities to a file called colvar:
+Another similar command is provided below. This time the histogram of Q6 parameters for the 64 atoms in a box of Lennard Jones is computed and printed to
+to a file called colvar:
-\plumedfile
-Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1} LABEL=q6
+```plumed
+q6: Q6 SPECIES=1-64 D_0=1.3 R_0=0.2 HISTOGRAM={GAUSSIAN LOWER=0.0 UPPER=1.0 NBINS=20 SMEAR=0.1}
PRINT ARG=q6.* FILE=colvar
-\endplumedfile
+```
-The following command could be used to measure the Q6 parameters that describe the arrangement of chlorine ions around the
+The command below could be used to measure the Q6 parameters that describe the arrangement of chlorine ions around the
sodium atoms in sodium chloride. The imagined system here is composed of 64 NaCl formula units and the atoms are arranged in the input
-with the 64 Na\f$^+\f$ ions followed by the 64 Cl\f$-\f$ ions. Once again the average Q6 parameter is calculated and output to a
+with the 64 Na$^+$ ions followed by the 64 Cl$-$ ions. Once again the average Q6 parameter is calculated and output to a
file called colvar
-\plumedfile
-Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN LABEL=q6
+```plumed
+q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
PRINT ARG=q6.mean FILE=colvar
-\endplumedfile
+```
If you simply want to examine the values of the Q6 parameters for each of the atoms in your system you can do so by exploiting the
-command \ref DUMPATOMS as shown in the example below. The following output file will output a file in an extended xyz format
+command [DUMPATOMS](DUMPATOMS.xyz) as shown in the example below. The following output file will output a file in an extended xyz format
called q6.xyz for each frame of the analyzed MD trajectory. The first column in this file will contain a dummy name for each of the
atoms, columns 2-4 will then contain the x, y and z positions of the atoms, column 5 will contain the value of the Q6 parameter, columns
-6-19 will contain the real parts of the director of the \f$q_{6m}\f$ vector while columns 20-33 will contain the imaginary parts of this director.
+6-12 will contain the real parts of the director of the $q_{6m}$ vector while columns 12-19 will contain the imaginary parts of this director.
-\plumedfile
+```plumed
q6: Q6 SPECIESA=1-64 SPECIESB=65-128 D_0=1.3 R_0=0.2 MEAN
-DUMPATOMS ARG=q6_anorm ATOMS=q6 FILE=q6.xyz
-\endplumedfile
+DUMPATOMS ATOMS=q6 ARG=q6_anorm FILE=q6.xyz
+```
*/
//+ENDPLUMEDOC
diff --git a/src/symfunc/ThreeBodyGFunctions.cpp b/src/symfunc/ThreeBodyGFunctions.cpp
index 2517513f11..ed0e455625 100644
--- a/src/symfunc/ThreeBodyGFunctions.cpp
+++ b/src/symfunc/ThreeBodyGFunctions.cpp
@@ -31,7 +31,65 @@ namespace symfunc {
/*
Calculate functions of the coordinates of the coordinates of all pairs of bonds in the first coordination sphere of an atom
-\par Examples
+This shortcut can be used to calculate [symmetry function](https://www.plumed-tutorials.org/lessons/23/001/data/SymmetryFunction.html) that
+are like those defined by Behler in the paper that is cited in the bibliography below. The particular symmetry functions that are computed
+by this shortcut are the angular ones that are functions of the set pairs of atoms in the coordination sphere of the central atom. One of
+the angular symmetry functions that Behler introduces is:
+
+$$
+G^5_i = 2^{1-\zeta} \sum_{j,k\ne i} (1 + \lambda\cos\theta_{ijk})^\zeta e^{-\nu(R_{ij}^2 + R_{ik}^2)} f_c(R_{ij}) f_c(R_{ik})
+$$
+
+In this expression $\zeta$, $\nu$ and $\lambda$ are all parameters. $f_c$ is a switching function which acts upon $R_{ij}$, the distance between atom $i$ and atom $j$, and
+$R_{ik}$, the distance between atom $i$ and atom $k$. $\theta_{ijk}$ is then the angle between the vector that points from atom $i$ to atom $j$ and the vector that points from
+atom $i$ to atom $k$. THe input below can be used to get PLUMED to calculate the 100 values for this symmetry function for the 100 atoms in a system.
+
+```plumed
+# Calculate the contact matrix and the x,y and z components of the bond vectors
+# This action calculates 4 100x100 matrices
+cmat: CONTACT_MATRIX GROUP=1-100 SWITCH={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS
+
+# Compute the symmetry function for the 100 atoms from the 4 100x100 matrices output
+# by cmat. The output from this action is a vector with 100 elements
+beh3: GSYMFUNC_THREEBODY ...
+ WEIGHT=cmat.w ARG=cmat.x,cmat.y,cmat.z
+ FUNCTION1={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3*cos(ajik))^3 LABEL=g5}
+...
+
+# Print the 100 symmetry function values to a file
+PRINT ARG=beh3.g5 FILE=colvar
+```
+
+The GSYMFUNC_THREEBODY action sums over all the distinct triples of atoms that are identified in the contact matrix. This action uses the same functionality as [CUSTOM](CUSTOM.md) and can thus compute any
+function of the following four quantities:
+
+* `rij` - the distance between atom $i$ and atom $j$
+* `rik` - the distance between atom $i$ and atom $k$
+* `rjk` - the distance between atom $j$ and atom $k$
+* `ajik` - the angle between the vector connecting atom $i$ to atom $j$ and the vector connecting atom $i$ to atom $k$.
+
+Furthermore we can calculate more than one function of these four quantities at a time as illustrated by the input below:
+
+```plumed
+# Calculate the contact matrix and the x,y and z components of the bond vectors
+# This action calculates 4 100x100 matrices
+cmat: CONTACT_MATRIX GROUP=1-100 SWITCH={CUSTOM R_0=4.5 D_MAX=4.5 FUNC=0.5*(cos(pi*x)+1)} COMPONENTS
+
+# Compute the 4 symmetry function below for the 100 atoms from the 4 100x100 matrices output
+# by cmat. The output from this action is a vector with 100 elements
+beh3: GSYMFUNC_THREEBODY ...
+ WEIGHT=cmat.w ARG=cmat.x,cmat.y,cmat.z
+ FUNCTION1={FUNC=0.25*(cos(pi*sqrt(rjk)/4.5)+1)*exp(-0.1*(rij+rik+rjk))*(1+2*cos(ajik))^2 LABEL=g4}
+ FUNCTION2={FUNC=0.25*exp(-0.1*(rij+rik))*(1+3.5*cos(ajik))^3 LABEL=g5}
+ FUNCTION3={FUNC=0.125*(1+6.6*cos(ajik))^4 LABEL=g6}
+ FUNCTION4={FUNC=sin(3.0*(ajik-1)) LABEL=g7}
+...
+
+# Print the 4 sets of 100 symmetry function values to a file
+PRINT ARG=beh3.g4,beh3.g5,beh3.g6,beh3.g7 FILE=colvar
+```
+
+You can read more about how to calculate more Behler-type symmetry functions [here](https://www.plumed-tutorials.org/lessons/23/001/data/Behler.html).
*/
//+ENDPLUMEDOC
@@ -55,7 +113,7 @@ void ThreeBodyGFunctions::registerKeywords( Keywords& keys ) {
keys.addInputKeyword("compulsory","ARG","matrix","three matrices containing the bond vectors of interest");
keys.addInputKeyword("compulsory","WEIGHT","matrix","the matrix that contains the weights that should be used for each connection");
keys.add("numbered","FUNCTION","the parameters of the function you would like to compute");
- ActionWithValue::useCustomisableComponents( keys );
+ ActionWithValue::useCustomisableComponents( keys ); keys.addDOI("10.1063/1.3553717");
}
ThreeBodyGFunctions::ThreeBodyGFunctions(const ActionOptions&ao):