Action: FIND_CONTOUR_SURFACE

Module	contour
Description	Usage
Find an isocontour by searching along either the x, y or z direction.
output value	type
a grid containing the location of the points in the Willard-Chandler surface along the chosen direction	grid

Details and examples

Find an isocontour by searching along either the x, y or z direction.

As discussed in the documentation for the gridtools, PLUMED contains a number of tools that allow you to calculate a function on a grid. The function on this grid might be a HISTOGRAM or it might be one of the phase fields that are discussed here. If this function has one or two input arguments it is relatively straightforward to plot the function. If by contrast the data has a three dimensions it can be difficult to visualize.

This action provides one tool for visualizing these functions. It can be used to search for a set of points on a contour where the function takes a particular value. In other words, for the function $f(x,y,z)$ this action would find a set of points $\{x_c,y_c,z_c\}$ that have:

$f(x_c,y_c,z_c) - c = 0$

where $c$ is some constant value that is specified by the user. The points on this contour are find by searching along lines that run parallel to the $x$ , $y$ or $z$ axis of the simulation cell. The result is, therefore, a two dimensional function evaluated on a grid that gives us the height of the interface as a function of two coordinates.

It is important to note that this action can only be used to detect contours in three dimensional functions. In addition, this action will fail to find the full set of contour points if the contour does not have the same topology as an infinite plane. If you are uncertain that the isocontours in your function have the appropriate topology you should use FIND_CONTOUR in place of this action.

Examples

The input shown below was used to analyze the results from a simulation of an interface between solid and molten Lennard Jones. The interface between the solid and the liquid was set up in the plane perpendicular to the $z$ direction of the simulation cell. The input below calculates something akin to a Willard-Chandler dividing surface (see contour) between the solid phase and the liquid phase. There are two of these interfaces within the simulation box because of the periodic boundary conditions but we were able to determine that one of these two surfaces lies in a particular part of the simulation box. The input below detects the height profile of one of these two interfaces. It does so by computing a phase field average from the values, $s_i$ , of the FCCUBIC symmetry functions for each of the atoms using the following expression.

$\rho'(x,y,z) = \frac{ \sum_{i=1}^N s_i K\left(\frac{x-x_i}{\lambda}, \frac{y-y_i}{\lambda}, \frac{z-z_i}{\lambda}\right) }{ \sum_{i=1}^N K\left(\frac{x-x_i}{\lambda}, \frac{y-y_i}{\lambda}, \frac{z-z_i}{\lambda}\right) }$

where $(x_i, y_i, z_i)$ is the position of atom $i$ relative to the position of atom 1, $K$ is a Gaussian kernel function and $\lambda=1.0$ .

Notice that we use the fact that we know roughly where the interface is when specifying how this phase field is to be calculated and specify the region over the $z$ -axis in which the KDE is computed. Once we have calculated the phase field we search for contour points on the lines that run parallel to the $z$ -direction of the cell box using the FIND_CONTOUR_SURFACE command. The final result is a $14 \times 14$ grid of values for the height of the interface as a function of the $(x,y)$ position. This grid is then output to a file called contour2.dat.

Notice that the commands below calculate the instantaneous position of the surface separating the solid and liquid and that as such the accumulated average is cleared on every step.

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

UNITSThis command sets the internal units for the code. More details NATURAL use natural units

# This calculates the value of a set of symmetry functions for the atoms of interest
The UNITS action with label  calculates somethingfccThe FCCUBIC action with label fcc calculates the following quantities: Quantity    Type    Description  
fcc vector the symmetry function for each of the specified atoms: FCCUBICMeasure how similar the environment around atoms is to that found in a FCC structure. This action is a shortcut. More details ...
  SPECIESthe list of atoms for which the symmetry function is being calculated and the atoms that can be in the environments=1-96000 SWITCHthe switching function that it used in the construction of the contact matrix. Options for this keyword are explained in the documentation for LESS_THAN.={CUBIC D_0=1.2 D_MAX=1.5}
  ALPHA The alpha parameter of the angular function that is used for FCCUBIC=27 PHI The Euler rotational angle phi=0.0 THETA The Euler rotational angle theta=-1.5708 PSI The Euler rotational angle psi=-2.35619
...
# PLUMED interprets the command:
# fcc: FCCUBIC ...
#   SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
#   ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619
# ...
# as follows (Click the red comment above to revert to the short version of the input):
fcc_grpThe GROUP action with label fcc_grp calculates the following quantities: Quantity    Type    Description  
fcc_grp atoms indices of atoms specified in GROUP: GROUPDefine a group of atoms so that a particular list of atoms can be referenced with a single label in definitions of CVs or virtual atoms. More details ATOMSthe numerical indexes for the set of atoms in the group=1-96000
fcc_matThe CONTACT_MATRIX action with label fcc_mat calculates the following quantities: Quantity    Type    Description  
fcc_mat.w matrix a matrix containing the weights for the bonds between each pair of atoms
fcc_mat.x matrix the projection of the bond on the x axis
fcc_mat.y matrix the projection of the bond on the y axis
fcc_mat.z matrix the projection of the bond on the z axis: CONTACT_MATRIXAdjacency matrix in which two atoms are adjacent if they are within a certain cutoff. More details GROUPspecifies the list of atoms that should be assumed indistinguishable=1-96000 SWITCHthe input for the switching function that acts upon the distance between each pair of atoms. Options for this keyword are explained in the documentation for LESS_THAN.={CUBIC D_0=1.2 D_MAX=1.5} COMPONENTS also calculate the components of the vector connecting the atoms in the contact matrix
fcc_xrotThe CUSTOM action with label fcc_xrot calculates the following quantities: Quantity    Type    Description  
fcc_xrot 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=fcc_mat.x,fcc_mat.y,fcc_mat.z FUNCthe function you wish to evaluate=-0.707104*x+2.59736e-06*y+0.70711*z PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
fcc_yrotThe CUSTOM action with label fcc_yrot calculates the following quantities: Quantity    Type    Description  
fcc_yrot 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=fcc_mat.x,fcc_mat.y,fcc_mat.z FUNCthe function you wish to evaluate=0.70711*x+2.59734e-06*y+0.707104*z PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
fcc_zrotThe CUSTOM action with label fcc_zrot calculates the following quantities: Quantity    Type    Description  
fcc_zrot 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=fcc_mat.x,fcc_mat.y,fcc_mat.z FUNCthe function you wish to evaluate=-0*x+1*y+-3.67321e-06*z PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
fcc_vfuncThe FCCUBIC_FUNC action with label fcc_vfunc calculates the following quantities: Quantity    Type    Description  
fcc_vfunc matrix the matrix obtained by doing an element-wise application of a function that measures the similarity with an fcc environment to the input matrix: FCCUBIC_FUNCMeasure how similar the environment around atoms is to that found in a FCC structure. More details ARGthe values input to this function=fcc_mat.x,fcc_mat.y,fcc_mat.z ALPHA The alpha parameter of the angular function=27
fcc_wvfuncThe CUSTOM action with label fcc_wvfunc calculates the following quantities: Quantity    Type    Description  
fcc_wvfunc 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=fcc_vfunc,fcc_mat.w 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
fcc_onesThe CONSTANT action with label fcc_ones calculates the following quantities: Quantity    Type    Description  
fcc_ones vector the constant value that was read from the plumed input: ONESCreate a constant vector with all elements equal to one More details SIZEthe number of ones that you would like to create=96000
fccThe MATRIX_VECTOR_PRODUCT action with label fcc calculates the following quantities: Quantity    Type    Description  
fcc 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=fcc_wvfunc,fcc_ones
fcc_denomThe MATRIX_VECTOR_PRODUCT action with label fcc_denom calculates the following quantities: Quantity    Type    Description  
fcc_denom 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=fcc_mat.w,fcc_ones
fcc_nThe CUSTOM action with label fcc_n calculates the following quantities: Quantity    Type    Description  
fcc_n 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=fcc,fcc_denom FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO
# --- End of included input --- 
# This determines the positions of the atoms of interest relative to the position of atom 1
dens2_distThe DISTANCE action with label dens2_dist calculates the following quantities: Quantity    Type    Description  
dens2_dist.x vector the x-component of the vector connecting the two atoms
dens2_dist.y vector the y-component of the vector connecting the two atoms
dens2_dist.z vector the z-component of the vector connecting the two atoms: DISTANCESCalculate the distances between multiple piars of atoms This action is a shortcut. More details ORIGINcalculate the distance of all the atoms specified using the ATOMS keyword from this point=1 ATOMSthe pairs of atoms that you would like to calculate the angles for=fcc COMPONENTS calculate the x, y and z components of the distance separately and store them as label
# PLUMED interprets the command:
# dens2_dist: DISTANCES ORIGIN=1 ATOMS=fcc COMPONENTS
# as follows (Click the red comment above to revert to the short version of the input):
dens2_dist: DISTANCECalculate the distance/s between pairs of atoms. More details COMPONENTS calculate the x, y and z components of the distance separately and store them as label ATOMS1the pair of atom that we are calculating the distance between=1,1 ATOMS2the pair of atom that we are calculating the distance between=1,2 ATOMS3the pair of atom that we are calculating the distance between=1,3 ATOMS4the pair of atom that we are calculating the distance between=1,4 ATOMS5the pair of atom that we are calculating the distance between=1,5     # Action input conctinues with 95995 further ATOMSn keywords, 
# --- End of included input --- # This computes the numerator in the expression above for the phase field
dens2_numerThe KDE action with label dens2_numer calculates the following quantities: Quantity    Type    Description  
dens2_numer grid a function on a grid that was obtained by doing a Kernel Density Estimation using the input arguments: KDECreate a histogram from the input scalar/vector/matrix using KDE This action has hidden defaults. More details ...
   VOLUMESthis keyword take the label of an action that calculates a vector of values=fcc_n ARGthe label for the value that should be used to construct the histogram=dens2_dist.x,dens2_dist.y,dens2_dist.z
   GRID_BINthe number of bins for the grid=14,14,50 GRID_MIN the lower bounds for the grid=auto,auto,6.0
   GRID_MAX the upper bounds for the grid=auto,auto,11.0 BANDWIDTHthe bandwidths for kernel density esimtation=1.0,1.0,1.0
...
dens2_numer: KDECreate a histogram from the input scalar/vector/matrix using KDE This action uses the defaults shown here. More details ...
   VOLUMESthis keyword take the label of an action that calculates a vector of values=fcc_n ARGthe label for the value that should be used to construct the histogram=dens2_dist.x,dens2_dist.y,dens2_dist.z
   GRID_BINthe number of bins for the grid=14,14,50 GRID_MIN the lower bounds for the grid=auto,auto,6.0
   GRID_MAX the upper bounds for the grid=auto,auto,11.0 BANDWIDTHthe bandwidths for kernel density esimtation=1.0,1.0,1.0
 KERNEL the kernel function you are using=GAUSSIAN
...
# This computes the denominator
dens2_denomThe KDE action with label dens2_denom calculates the following quantities: Quantity    Type    Description  
dens2_denom grid a function on a grid that was obtained by doing a Kernel Density Estimation using the input arguments: KDECreate a histogram from the input scalar/vector/matrix using KDE This action has hidden defaults. More details ...
   ARGthe label for the value that should be used to construct the histogram=dens2_dist.x,dens2_dist.y,dens2_dist.z
   GRID_BINthe number of bins for the grid=14,14,50 GRID_MIN the lower bounds for the grid=auto,auto,6.0
   GRID_MAX the upper bounds for the grid=auto,auto,11.0 BANDWIDTHthe bandwidths for kernel density esimtation=1.0,1.0,1.0
...
dens2_denom: KDECreate a histogram from the input scalar/vector/matrix using KDE This action uses the defaults shown here. More details ...
   ARGthe label for the value that should be used to construct the histogram=dens2_dist.x,dens2_dist.y,dens2_dist.z
   GRID_BINthe number of bins for the grid=14,14,50 GRID_MIN the lower bounds for the grid=auto,auto,6.0
   GRID_MAX the upper bounds for the grid=auto,auto,11.0 BANDWIDTHthe bandwidths for kernel density esimtation=1.0,1.0,1.0
 KERNEL the kernel function you are using=GAUSSIAN
...
# This computes the final phase field
dens2The CUSTOM action with label dens2 calculates the following quantities: Quantity    Type    Description  
dens2 grid the grid obtained by doing an element-wise application of an arbitrary function to the input grid: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=dens2_numer,dens2_denom FUNCthe function you wish to evaluate=x/y PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO

# We can now find and print the location of the two dimensional contour surface
ss2The FIND_CONTOUR_SURFACE action with label ss2 calculates the following quantities: Quantity    Type    Description  
ss2 grid a grid containing the location of the points in the Willard-Chandler surface along the chosen direction: FIND_CONTOUR_SURFACEFind an isocontour by searching along either the x, y or z direction. This action has hidden defaults. More details ARGthe labels of the grid in which the contour will be found=dens2 CONTOURthe value we would like to draw the contour at in the space=0.42 SEARCHDIRIn which directions do you wish to search for the contour=dens2_dist.z
ss2: FIND_CONTOUR_SURFACEFind an isocontour by searching along either the x, y or z direction. This action uses the defaults shown here. More details ARGthe labels of the grid in which the contour will be found=dens2 CONTOURthe value we would like to draw the contour at in the space=0.42 SEARCHDIRIn which directions do you wish to search for the contour=dens2_dist.z  INTERPOLATION_TYPE the method to use for interpolation=spline
DUMPGRIDOutput the function on the grid to a file with the PLUMED grid format. More details ARGthe label for the grid that you would like to output=ss2 FILE the file on which to write the grid=contour2.dat STRIDE the frequency with which the grid should be output to the file=1

Input

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

Keyword	Type	Description
ARG	grid	the labels of the grid in which the contour will be found

Full list of keywords

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

Keyword	Type	Default	Description
ARG	input	none	the labels of the grid in which the contour will be found
CONTOUR	compulsory	none	the value we would like to draw the contour at in the space
INTERPOLATION_TYPE	compulsory	spline	the method to use for interpolation
SEARCHDIR	compulsory	none	In which directions do you wish to search for the contour
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.

References

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

E. Baldi, M. Ceriotti, G. A. Tribello, Extracting the interfacial free energy and anisotropy from a smooth fluctuating dividing surface. Journal of Physics: Condensed Matter. 29, 445001 (2017)

Quantity	Type	Description
fcc	vector	the symmetry function for each of the specified atoms