Action: HBOND_MATRIX

Module	adjmat
Description	Usage
Adjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them.
output value	type
a matrix containing the weights for the bonds between each pair of atoms	matrix

Details and examples

Adjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them.

A useful tool for developing complex collective variables is the notion of the so called adjacency matrix. An adjacency matrix is an $N \times N$ matrix in which the $i$ th, $j$ th element tells you whether or not the $i$ th and $j$ th atoms/molecules from a set of $N$ atoms/molecules are adjacent or not. As detailed in the documentation for CONTACT_MATRIX there are then a range of further analyses that you can perform on these matrices.

For this action the elements of the adjacency matrix are calculated using:

$a_{ij} = \sigma_{oo}( |\mathbf{r}_{ij}| ) \sum_{k=1}^N \sigma_{oh}( |\mathbf{r}_{ik}| ) \sigma_{\theta}( \theta_{kij} )$

This expression was derived by thinking about how to detect if there is a hydrogen bond between atoms $i$ and $j$ . The notion is that if the hydrogen bond is present atoms $i$ and $j$ should be within a certain cutoff distance. In addition, there should be a hydrogen within a certain cutoff distance of atom $i$ and this hydrogen should lie on or close to the vector connecting atoms $i$ and $j$ . As such $\sigma_{oo}(r_{ij})$ is a switching function that acts on the modulus of the vector connecting atom $i$ to atom $j$ . The sum over $k$ then runs over all the hydrogen atoms that are specified using using HYDROGEN keyword. $\sigma_{oh}(r_{ik})$ is a switching function that acts on the modulus of the vector connecting atom $i$ to atom $k$ and $\sigma_{\theta}(\theta_{kij})$ is a switching function that acts on the angle between the vector connecting atoms $i$ and $j$ and the vector connecting atoms $i$ and $k$ .

It is important to note that hydrogen bonds, unlike regular bonds, are asymmetric. In other words, the hydrogen atom does not sit at the mid point between the two other atoms in this three-center bond. As a result of this adjacency matrices calculated using HBOND_MATRIX are not symmetric like those calculated by CONTACT_MATRIX.

Each water molecule can participate in a hydrogen bond in one of two ways. It can either donate one of its hydrogen atom to the neighboring oxygen or it can accept a bond between the hydrogen of a neighboring water molecule and its own oxygen.

The following input can be used to analyze the number of hydrogen bonds each of the oxygen atoms in a box of water donates. This information is output in an xyz files which contains five columns of data. The first four of these columns are a label for the atom and the x, y and z position of the oxygen. The last column is then the number of hydrogen bond that water molecule donates.

Click on the labels of the actions for more information on what each action computes

matThe HBOND_MATRIX action with label mat calculates the following quantities: Quantity    Type    Description  
mat matrix a matrix containing the weights for the bonds between each pair of atoms: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action has hidden defaults. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
...
mat: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action uses the defaults shown here. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
 NL_CUTOFF The cutoff for the neighbor list=0.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=1
...
onesThe ONES action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector a vector of ones with the required number of elements: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=64
# PLUMED interprets the command:
# ones: ONES SIZE=64
# as follows (Click the red comment above to revert to the short version of the input):
onesThe CONSTANT action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector the constant value that was read from the plumed input: CONSTANTCreate a constant value that can be passed to actions More details NOLOG do not report all the read in scalars in the log VALUESthe numbers that are in your constant value=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
# --- End of included input --- rsumsThe MATRIX_VECTOR_PRODUCT action with label rsums calculates the following quantities: Quantity    Type    Description  
rsums vector the vector that is obtained by taking the product between the matrix and the vector that were input: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=mat,ones
DUMPATOMSDump selected atoms on a file. More details ATOMSthe atom indices whose positions you would like to print out=1-192:3 ARGthe labels of vectors that should be output in the xyz file=rsums FILEfile on which to output coordinates; extension is automatically detected=donors.xyz

If you want to calculate the number of hydorgen bonds each of the oxygen atoms accepts you would use an input like the one below:

Click on the labels of the actions for more information on what each action computes

matThe HBOND_MATRIX action with label mat calculates the following quantities: Quantity    Type    Description  
mat matrix a matrix containing the weights for the bonds between each pair of atoms: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action has hidden defaults. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
...
mat: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action uses the defaults shown here. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
 NL_CUTOFF The cutoff for the neighbor list=0.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=1
...
matTThe TRANSPOSE action with label matT calculates the following quantities: Quantity    Type    Description  
matT matrix the transpose of the input matrix: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=mat
onesThe ONES action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector a vector of ones with the required number of elements: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=64
# PLUMED interprets the command:
# ones: ONES SIZE=64
# as follows (Click the red comment above to revert to the short version of the input):
onesThe CONSTANT action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector the constant value that was read from the plumed input: CONSTANTCreate a constant value that can be passed to actions More details NOLOG do not report all the read in scalars in the log VALUESthe numbers that are in your constant value=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
# --- End of included input --- rsumsThe MATRIX_VECTOR_PRODUCT action with label rsums calculates the following quantities: Quantity    Type    Description  
rsums vector the vector that is obtained by taking the product between the matrix and the vector that were input: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=matT,ones
DUMPATOMSDump selected atoms on a file. More details ATOMSthe atom indices whose positions you would like to print out=1-192:3 ARGthe labels of vectors that should be output in the xyz file=rsums FILEfile on which to output coordinates; extension is automatically detected=acceptors.xyz

Consequently, if you want the total number of hydrogen bonds each oxygen atom participates in you need to use the following input:

Click on the labels of the actions for more information on what each action computes

matThe HBOND_MATRIX action with label mat calculates the following quantities: Quantity    Type    Description  
mat matrix a matrix containing the weights for the bonds between each pair of atoms: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action has hidden defaults. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
...
mat: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action uses the defaults shown here. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
 NL_CUTOFF The cutoff for the neighbor list=0.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=1
...
matTThe TRANSPOSE action with label matT calculates the following quantities: Quantity    Type    Description  
matT matrix the transpose of the input matrix: TRANSPOSECalculate the transpose of a matrix More details ARGthe label of the vector or matrix that should be transposed=mat
hbmatThe CUSTOM action with label hbmat calculates the following quantities: Quantity    Type    Description  
hbmat matrix the matrix obtained by doing an element-wise application of an arbitrary function to the input matrix: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=mat,matT FUNCthe function you wish to evaluate=x+y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
onesThe ONES action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector a vector of ones with the required number of elements: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=64
# PLUMED interprets the command:
# ones: ONES SIZE=64
# as follows (Click the red comment above to revert to the short version of the input):
onesThe CONSTANT action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector the constant value that was read from the plumed input: CONSTANTCreate a constant value that can be passed to actions More details NOLOG do not report all the read in scalars in the log VALUESthe numbers that are in your constant value=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
# --- End of included input --- rsumsThe MATRIX_VECTOR_PRODUCT action with label rsums calculates the following quantities: Quantity    Type    Description  
rsums vector the vector that is obtained by taking the product between the matrix and the vector that were input: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=hbmat,ones
DUMPATOMSDump selected atoms on a file. More details ATOMSthe atom indices whose positions you would like to print out=1-192:3 ARGthe labels of vectors that should be output in the xyz file=rsums FILEfile on which to output coordinates; extension is automatically detected=hbonds.xyz

Notice that in all the inputs above the $r_{ij}$ and $r_{ik}$ values that enter the formula above are calculated in a way that takes the periodic boundary conditions into account. If you want to ignore the periodic boundary conditions you can use the NOPBC flag as shown below.

Click on the labels of the actions for more information on what each action computes

matThe HBOND_MATRIX action with label mat calculates the following quantities: Quantity    Type    Description  
mat matrix a matrix containing the weights for the bonds between each pair of atoms: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action has hidden defaults. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
   NOPBC don't use pbc
...
mat: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action uses the defaults shown here. More details ...
   DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
   SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
   ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
   NOPBC don't use pbc
 NL_CUTOFF The cutoff for the neighbor list=0.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=1
...

COMPONENTS flag

If you add the flag COMPONENTS to the input as shown below:

Click on the labels of the actions for more information on what each action computes

c4The HBOND_MATRIX action with label c4 calculates the following quantities: Quantity    Type    Description  
c4.w matrix a matrix containing the weights for the bonds between each pair of atoms
c4.x matrix the projection of the bond on the x axis
c4.y matrix the projection of the bond on the y axis
c4.z matrix the projection of the bond on the z axis: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action has hidden defaults. More details ...
  DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
  SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
  ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
  COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix
...
c4: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action uses the defaults shown here. More details ...
  DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
  SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
  ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
  COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix
 NL_CUTOFF The cutoff for the neighbor list=0.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=1
...

then four matrices with the labels c4.w, c4.x, c4.y and c4.z are output by the action. The matrix with the label c4.w is the adjacency matrix that would be output if you had not added the COMPONENTS flag. The $i,j$ component of the matrices c4.x, c4.y and c4.z contain the $x$ , $y$ and $z$ components of the vector connecting atoms $j$ and $k$ . Importantly, however, the components of these vectors are only stored in c4.x, c4.y and c4.z if the elements of c4.w are non-zero. Using the COMPONENTS flag in this way ensures that you can use HBOND_MATRIX in tandem with many of the functionalities that are part of the symfunc module. Remember, however, that the $i,j$ element of the HBOND_MATRIX is only non-zero if atom $i$ donates a hydrogen bond to atom $j$ . You cannot use HBOND_MATRIX to identify the set of atoms that each atom is hydrogen bonded to.

The MASK keyword

You use the MASK keyword with HBOND_MATRIX in the same way that is used in CONTACT_MATRIX. This keyword thus expects a vector in input, which tells HBOND_MATRIX that it is safe to not calculate certain rows of the output matrix. An example where this keyword is used is shown below:

Click on the labels of the actions for more information on what each action computes

# Fixed virtual atom which serves as the probe volume's center (pos. in nm)
centerThe FIXEDATOM action with label center calculates the following quantities: Quantity    Type    Description  
center atoms virtual atom calculated by FIXEDATOM action: FIXEDATOMAdd a virtual atom in a fixed position. This action has hidden defaults. More details ATcoordinates of the virtual atom=2.5,2.5,2.5
center: FIXEDATOMAdd a virtual atom in a fixed position. This action uses the defaults shown here. More details ATcoordinates of the virtual atom=2.5,2.5,2.5  SET_MASS mass of the virtual atom=1 SET_CHARGE charge of the virtual atom=0
# Vector in which element i is one if atom i is in sphere of interest and zero otherwise
sphereThe INSPHERE action with label sphere calculates the following quantities: Quantity    Type    Description  
sphere vector vector of numbers between 0 and 1 that measure the degree to which each atom is within the volume of interest: INSPHEREThis quantity can be used to calculate functions of the distribution of collective variables for the atoms that lie in a particular, user-specified part of of the cell. More details ATOMSthe group of atoms that you would like to investigate=1-192:3 CENTERthe atom whose vicinity we are interested in examining=center RADIUSthe switching function that tells us the extent of the sphereical region of interest. Options for this keyword are explained in the documentation for LESS_THAN.={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
# Calculates cooordination numbers
cmapThe HBOND_MATRIX action with label cmap calculates the following quantities: Quantity    Type    Description  
cmap matrix a matrix containing the weights for the bonds between each pair of atoms: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action has hidden defaults. More details ...
  DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
  SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
  ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi} MASKa vector that is used to used to determine which rows of the adjancency matrix to compute=sphere
...
cmap: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. This action uses the defaults shown here. More details ...
  DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
  SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30}
  ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi} MASKa vector that is used to used to determine which rows of the adjancency matrix to compute=sphere
 NL_CUTOFF The cutoff for the neighbor list=0.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=1
...
onesThe ONES action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector a vector of ones with the required number of elements: ONESCreate a constant vector with all elements equal to one This action is a shortcut. More details SIZEthe number of ones that you would like to create=64
# PLUMED interprets the command:
# ones: ONES SIZE=64
# as follows (Click the red comment above to revert to the short version of the input):
onesThe CONSTANT action with label ones calculates the following quantities: Quantity    Type    Description  
ones vector the constant value that was read from the plumed input: CONSTANTCreate a constant value that can be passed to actions More details NOLOG do not report all the read in scalars in the log VALUESthe numbers that are in your constant value=1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
# --- End of included input --- ccThe MATRIX_VECTOR_PRODUCT action with label cc calculates the following quantities: Quantity    Type    Description  
cc vector the vector that is obtained by taking the product between the matrix and the vector that were input: MATRIX_VECTOR_PRODUCTCalculate the product of the matrix and the vector More details ARGthe label for the matrix and the vector/scalar that are being multiplied=cmap,ones
# Multiply coordination numbers by sphere vector
prodThe CUSTOM action with label prod calculates the following quantities: Quantity    Type    Description  
prod vector the vector obtained by doing an element-wise application of an arbitrary function to the input vectors: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=cc,sphere FUNCthe function you wish to evaluate=x*y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Sum of coordination numbers for atoms that are in the sphere of interest
numerThe SUM action with label numer calculates the following quantities: Quantity    Type    Description  
numer scalar the SUM of the elements in the input value: SUMCalculate the sum of the arguments More details ARGthe vector/matrix/grid whose elements shuld be added together=prod PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Number of atoms that are in sphere of interest
denomThe SUM action with label denom calculates the following quantities: Quantity    Type    Description  
denom scalar the SUM of the elements in the input value: SUMCalculate the sum of the arguments More details ARGthe vector/matrix/grid whose elements shuld be added together=sphere PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# Average coordination number for atoms in sphere of interest
avThe CUSTOM action with label av calculates the following quantities: Quantity    Type    Description  
av vector the vector obtained by doing an element-wise application of an arbitrary function to the input vectors: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=prod,sphere FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# And print out final CV to a file
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=av FILEthe name of the file on which to output these quantities=colvar STRIDE the frequency with which the quantities of interest should be output=1

This input calculates the average number of hydrogen bonds each of the aatoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$ donate.

Optimisation details

Adjacency matrices are sparse. Each atom is only be connected to a small number of neighbours and the vast majority of the elements of the contact matrix are thus zero. To reduce the amount of memory that PLUMED requires PLUMED uses sparse matrix storage. Consequently, whenever you calculate and store a contact matrix only the elements of the matrix that are non-zero are stored. The same thing holds for the additional values that are created when you use the COMPONENTS flag. The components of the vectors connecting atoms are only stored when the elements of c4.w are non-zero.

We can also use the sparsity of the adjacency matrix to make the time required to compute a contact matrix scale linearly rather than quadratically with the number of atoms. Element $i,j$ of the contact matrix is only non-zero if two atoms are within a cutoff, $r_c$ . We can determine that many pairs of atoms are further appart than $r_c$ without computing the distance between these atoms by using divide and conquer strategies such as linked lists and neighbour lists. To turn on these features you need to set the D_MAX parameter in the switching functions. The value you pass to the D_MAX keyword is used as the cutoff in the link cell algorithm.

In theory we could further optimize the implementation of the HBOND_MATRIX action by exploiting neighbor lists. If we were to do this we would likely add two further keywords as shown below:

Click on the labels of the actions for more information on what each action computes

cmapThe HBOND_MATRIX action with label cmap calculates the following quantities: Quantity    Type    Description  
cmap matrix a matrix containing the weights for the bonds between each pair of atoms: HBOND_MATRIXAdjacency matrix in which two atoms are adjacent if there is a hydrogen bond between them. More details ...
  DONORSThe list of atoms which can donate a hydrogen bond=1-192:3 ACCEPTORSThe list of atoms which can accept a hydrogen bond=1-192:3 HYDROGENSThe list of atoms that can form the bridge between the two interesting parts of the structure=2-192:3,3-192:3
  SWITCHThe switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=3.20 D_MAX=4.0} HSWITCHThe switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=2.30 D_MAX=3.5}
  ASWITCHA switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.167pi}
  NL_CUTOFF The cutoff for the neighbor list=5.0 NL_STRIDE The frequency with which we are updating the atoms in the neighbor list=5
...

The NL_CUTOFF keyword would be used to specify the cutoff (in nm) to use when constructing neighbor lists. This value would need to be slightly larger than the D_MAX parameter for the switching function that acts on the acceptor donor distance. The NL_STRIDE keyword would then be used to specify how frequently the neighbour list should be updated. Thus far we have not found it necessary to implement this algorithm. We have been happy with the performance even if we use the linked list algorithm to update the neighbors on every step. If you feel that you need this CV to perform better please get in touch as adding a neighbor list for this action should be relatively straightforward.

Input

The arguments and atoms that serve as the input for this action are specified using one or more of the keywords in the following table.

Keyword	Type	Description
MASK	vector	a vector that is used to used to determine which rows of the adjancency matrix to compute
DONORS	atoms	The list of atoms which can donate a hydrogen bond
ACCEPTORS	atoms	The list of atoms which can accept a hydrogen bond
HYDROGENS	atoms	The list of atoms that can form the bridge between the two interesting parts of the structure

Output components

This action can calculate the values in the following table when the associated keyword is included in the input for the action. These values can be referenced elsewhere in the input by using this Action's label followed by a dot and the name of the value required from the list below.

Name	Type	Keyword	Description
w	matrix	COMPONENTS	a matrix containing the weights for the bonds between each pair of atoms
x	matrix	COMPONENTS	the projection of the bond on the x axis
y	matrix	COMPONENTS	the projection of the bond on the y axis
z	matrix	COMPONENTS	the projection of the bond on the z axis

Full list of keywords

The following table describes the keywords and options that can be used with this action

Keyword	Type	Default	Description
MASK	input	none	a vector that is used to used to determine which rows of the adjancency matrix to compute
DONORS	input	none	The list of atoms which can donate a hydrogen bond
ACCEPTORS	input	none	The list of atoms which can accept a hydrogen bond
HYDROGENS	input	none	The list of atoms that can form the bridge between the two interesting parts of the structure
NL_CUTOFF	compulsory	0.0	The cutoff for the neighbor list
NL_STRIDE	compulsory	1	The frequency with which we are updating the atoms in the neighbor list
COMPONENTS	optional	false	also calculate the components of the vector connecting the atoms in the contact matrix
NOPBC	optional	false	don't use pbc
SWITCH	optional	not used	The switchingfunction that specifies how close a pair of atoms must be together for there to be a hydrogen bond between them. Options for this keyword are explained in the documentation for LESS_THAN.
HSWITCH	optional	not used	The switchingfunction that specifies how close the hydrogen must be to the donor atom of the hydrogen bond for it to be considered a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.
ASWITCH	optional	not used	A switchingfunction that is used to specify what the angle between the vector connecting the donor atom to the acceptor atom and the vector connecting the donor atom to the hydrogen must be in order for it considered to be a hydrogen bond. Options for this keyword are explained in the documentation for LESS_THAN.
SERIAL	optional	false	do the calculation in serial. Further information about this flag can be found here.
USEGPU	optional	false	run this calculation on the GPU. Further information about this flag can be found here.

deprecated keywords

The keywords in the following table can still be used with this action but have been deprecated

Keyword	Description
GROUP	the atoms for which you would like to calculate the adjacency matrix

References

More information about how this action can be used is available in the following articles:

G. A. Tribello, F. Giberti, G. C. Sosso, M. Salvalaglio, M. Parrinello, Analyzing and Driving Cluster Formation in Atomistic Simulations. Journal of Chemical Theory and Computation. 13, 1317–1327 (2017)

Quantity	Type	Description
mat	matrix	a matrix containing the weights for the bonds between each pair of atoms