Action: TOPOLOGY_MATRIX

Module	adjmat
Description	Usage
Adjacency matrix in which two atoms are adjacent if they are connected topologically
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 they are connected topologically

The functionality in this action was developed as part of a project that attempted to study the nucleation of bubbles. The idea was to develop a criterion that would allow one to determine if two gas atoms $i$ and $j$ are part of the same bubble or not. This criterion would then be used to construct a adjacency matrix, which could be used in the same way that CONTACT_MATRIX is used in other methods.

The criterion that was developed to determine whether atom $i$ and $j$ are connected in this way works by determining if the density within a cylinder that is centered on the vector connecting atoms $i$ and $j$ is less than a certain threshold value. To make this determination we first determine whether any given atom $k$ is within the cylinder centered on the vector connecting atoms $i$ and $j$ by using the following expression

$f(\mathbf{r}_{ik}, \mathbf{r}_{ij}) = s_1\left( \sqrt{ |\mathbf{r}_{ij}|^2 - \left( \frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{ik}}{|\mathbf{r}_{ij} |} \right)^2} \right)s_2\left( -\frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{ik}}{|\mathbf{r}_{ij} |} \right) s_2\left( \frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{ik}}{|\mathbf{r}_{ij} |} - |\mathbf{r}_{ij}| \right)$

In this expression $s_1$ and $s_2$ are switching functions, while $\mathbf{r}_{lm}$ is used to indicate the vector connecting atoms $l$ and $m$ .

We then calculate the density for a grid of $M$ points along the vector connecting atom $i$ and atom $j$ using and find the maximum density on this grid using:

$\rho_{ij} = \max_{m \in M} \left[ \frac{M}{d_\textrm{max}} \right] \sum_k f(r_{ik}, r_{ij}) \int_{(m-1)d_{\textrm{max}}/M}^{ md_{\textrm{max}} /M } \textrm{d}x \quad K\left( \frac{x - r_{ks} \cdot r_{ij} }{ | r_{ks} | }\right)$

where $d_\textrm{max}$ is the D_MAX parameter of the switching function $s_3$ that appears in the next equation, $K$ is a kernal function and $s$ is used to represent a point in space that is $d_\textrm{max}$ from atom $j$ along the vector connecting atom $j$ to atom $i$ .

The final value that is stored in the $i, j$ element of the output matrix is:

$T_{ij} = s_3(|\mathbf{r}_{ij}|)s_4(\rho_{ij})$

where $s_3$ and $s_4$ are switching functions.

We ended up abandoning this method and the project (if you want drive bubble formation you are probably better off adding a bias on the volume of the simulation cell). However, if you would like to try this method an example input that uses this action would look like this:

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

matThe TOPOLOGY_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: TOPOLOGY_MATRIXAdjacency matrix in which two atoms are adjacent if they are connected topologically This action has hidden defaults. More details ...
    GROUPthe atoms for which you would like to calculate the adjacency matrix=1-85 BACKGROUND_ATOMSthe list of atoms that should be considered as part of the background density=86-210
    BIN_SIZEthe size to use for the bins=1.02 SIGMAthe width of the function to be used for kernel density estimation=0.17 KERNEL the type of kernel function to be used=triangular
    CYLINDER_SWITCHa switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.5 D_MAX=1.0}
    SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
    RADIUSswtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
    DENSITY_THRESHOLDa switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
...
mat: TOPOLOGY_MATRIXAdjacency matrix in which two atoms are adjacent if they are connected topologically This action uses the defaults shown here. More details ...
    GROUPthe atoms for which you would like to calculate the adjacency matrix=1-85 BACKGROUND_ATOMSthe list of atoms that should be considered as part of the background density=86-210
    BIN_SIZEthe size to use for the bins=1.02 SIGMAthe width of the function to be used for kernel density estimation=0.17 KERNEL the type of kernel function to be used=triangular
    CYLINDER_SWITCHa switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.5 D_MAX=1.0}
    SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
    RADIUSswtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
    DENSITY_THRESHOLDa switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
 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
...

The switching functions $s_1$ , $s_2$ , $s_3$ and $s_4$ are specified using the RADIUS, CYLINDER_SWITCH, SWITCH and DENSITY_THRESHOLD keywords respectively. We loop over the atoms in the group specified using the BACKGROUND_ATOMS keyword when looping over $k$ in the formulas above. An $85 \times 85$ matrix is output from the method as we are determining the connectivity between the atoms specified via the GROUP keyword.

The above example assumes that you want to calculate the connectivity within a single group of atoms. If you would calculate connectivity between two different groups of atoms you use the GROUPA and GROUPB keywords as shown below:

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

matThe TOPOLOGY_MATRIX action with label mat calculates the following quantities: Quantity    Type    Description  
mat.w matrix a matrix containing the weights for the bonds between each pair of atoms
mat.x matrix the projection of the bond on the x axis
mat.y matrix the projection of the bond on the y axis
mat.z matrix the projection of the bond on the z axis: TOPOLOGY_MATRIXAdjacency matrix in which two atoms are adjacent if they are connected topologically This action has hidden defaults. More details ...
    GROUPAwhen you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB=1-20 GROUPBwhen you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA=21-85 BACKGROUND_ATOMSthe list of atoms that should be considered as part of the background density=86-210
    BIN_SIZEthe size to use for the bins=1.02 SIGMAthe width of the function to be used for kernel density estimation=0.17 KERNEL the type of kernel function to be used=triangular
    CYLINDER_SWITCHa switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.5 D_MAX=1.0}
    SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
    RADIUSswtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
    DENSITY_THRESHOLDa switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
    COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix NOPBC don't use pbc
...
mat: TOPOLOGY_MATRIXAdjacency matrix in which two atoms are adjacent if they are connected topologically This action uses the defaults shown here. More details ...
    GROUPAwhen you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB=1-20 GROUPBwhen you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA=21-85 BACKGROUND_ATOMSthe list of atoms that should be considered as part of the background density=86-210
    BIN_SIZEthe size to use for the bins=1.02 SIGMAthe width of the function to be used for kernel density estimation=0.17 KERNEL the type of kernel function to be used=triangular
    CYLINDER_SWITCHa switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.5 D_MAX=1.0}
    SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
    RADIUSswtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
    DENSITY_THRESHOLDa switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
    COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix 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
...

Notice that we have also added the NOPBC and COMPONENTS keywords in this input. The action above thus outputs four matrices with the labels mat.w, mat.x, mat.y and mat.z. The matrix with the label mat.w is the adjacency matrix that would be output if you had not added the COMPONENTS flag. The $i,j$ component of the matrices mat.x, mat.y and mat.z contain the $x$ , $y$ and $z$ components of the vector connecting atoms $i$ and $k$ . Importantly, however, the components of these vectors are only stored in mat.x, mat.y and mat.z if the elements of mat.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.

The NOPBC flag, meanwhile, ensures that all distances are calculated in a way that does not take the periodic boundary conditions into account. By default, distances are calculated in a way that takes periodic boundary conditions into account.

The MASK keyword

You use the MASK keyword with TOPOLOGY_MATRIX in the same way that is used in CONTACT_MATRIX. This keyword thus expects a vector in input, which tells TOPOLOGY_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-85 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 TOPOLOGY_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: TOPOLOGY_MATRIXAdjacency matrix in which two atoms are adjacent if they are connected topologically This action has hidden defaults. More details ...
  GROUPthe atoms for which you would like to calculate the adjacency matrix=1-85 BACKGROUND_ATOMSthe list of atoms that should be considered as part of the background density=86-210
  BIN_SIZEthe size to use for the bins=1.02 SIGMAthe width of the function to be used for kernel density estimation=0.17 KERNEL the type of kernel function to be used=triangular
  CYLINDER_SWITCHa switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.5 D_MAX=1.0}
  SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
  RADIUSswtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
  DENSITY_THRESHOLDa switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
  MASKa vector that is used to used to determine which rows of the adjancency matrix to compute=sphere
...
cmap: TOPOLOGY_MATRIXAdjacency matrix in which two atoms are adjacent if they are connected topologically This action uses the defaults shown here. More details ...
  GROUPthe atoms for which you would like to calculate the adjacency matrix=1-85 BACKGROUND_ATOMSthe list of atoms that should be considered as part of the background density=86-210
  BIN_SIZEthe size to use for the bins=1.02 SIGMAthe width of the function to be used for kernel density estimation=0.17 KERNEL the type of kernel function to be used=triangular
  CYLINDER_SWITCHa switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.5 D_MAX=1.0}
  SWITCHThis keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
  RADIUSswtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
  DENSITY_THRESHOLDa switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.={RATIONAL R_0=0.1 D_MAX=0.5}
  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=85
# PLUMED interprets the command:
# ones: ONES SIZE=85
# 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,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 topological connections there are for each of the atoms that are within a spherical region that is centered on the point $(2.5,2.5,2.5)$ .

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
GROUP	atoms	the atoms for which you would like to calculate the adjacency matrix
GROUPA	atoms	when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB
GROUPB	atoms	when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA
BACKGROUND_ATOMS	atoms	the list of atoms that should be considered as part of the background density

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
GROUP	input	none	the atoms for which you would like to calculate the adjacency matrix
GROUPA	input	none	when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPB
GROUPB	input	none	when you are calculating the adjacency matrix between two sets of atoms this keyword is used to specify the atoms along with the keyword GROUPA
BACKGROUND_ATOMS	input	none	the list of atoms that should be considered as part of the background density
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
SWITCH	compulsory	none	This keyword is used if you want to employ an alternative to the continuous swiching function defined above. Options for this keyword are explained in the documentation for LESS_THAN.
RADIUS	compulsory	none	swtiching function that acts upon the radial distance of the atom from the center of the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.
CYLINDER_SWITCH	compulsory	none	a switching function on ( r_ij . Options for this keyword are explained in the documentation for LESS_THAN.
BIN_SIZE	compulsory	none	the size to use for the bins
DENSITY_THRESHOLD	compulsory	none	a switching function that acts upon the maximum density in the cylinder. Options for this keyword are explained in the documentation for LESS_THAN.
SIGMA	compulsory	none	the width of the function to be used for kernel density estimation
KERNEL	compulsory	gaussian	the type of kernel function to be used
COMPONENTS	optional	false	also calculate the components of the vector connecting the atoms in the contact matrix
NOPBC	optional	false	don't use pbc
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
ATOMS	You should use GROUP instead of this keyword which was used in older versions of PLUMED and is provided for back compatibility only

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