FIND_CONTOUR_SURFACE
This is part of the contour module
It is only available if you configure PLUMED with ./configure –enable-modules=contour . Furthermore, this feature is still being developed so take care when using it and report any problems on the mailing list.

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

As discussed in the part of the manual on Analysis 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 as a function of a few collective variables or it might be a phase field that has been calculated using MULTICOLVARDENS. 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 FIND_CONTOUR_SURFACE.

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 [116] 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 of the FCCUBIC symmetry function using the MULTICOLVARDENS action. 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 we are going to search for the phase field in the line defining the MULTICOLVARDENS. 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
tested on master
UNITS 
NATURAL
( default=off ) use natural units
fcc: FCCUBIC ...
SPECIES
this keyword is used for colvars such as coordination number.
=1-96000
SWITCH
the switching function that it used in the construction of the contact matrix
={CUBIC D_0=1.2 D_MAX=1.5}
ALPHA
compulsory keyword ( default=3.0 ) The alpha parameter of the angular function that is used for FCCUBIC
=27
PHI
compulsory keyword ( default=0.0 ) The Euler rotational angle phi
=0.0
THETA
compulsory keyword ( default=0.0 ) The Euler rotational angle theta
=-1.5708
PSI
compulsory keyword ( default=0.0 ) The Euler rotational angle psi
=-2.35619 ... dens2: MULTICOLVARDENS
DATA
the multicolvar which you would like to calculate the density profile for
=fcc
ORIGIN
compulsory keyword we will use the position of this atom as the origin
=1
DIR
compulsory keyword the direction in which to calculate the density profile
=xyz
NBINS
the number of bins to use in each direction (alternative to GRID_NBIN)
=14,14,50
ZREDUCED
could not find this keyword
ZLOWER
could not find this keyword
=6.0
ZUPPER
could not find this keyword
=11.0
BANDWIDTH
the bandwidths for kernel density esimtation
=1.0,1.0,1.0
CLEAR
compulsory keyword ( default=0 ) the frequency with which to clear the density
=1 ss2: FIND_CONTOUR_SURFACE
GRID
could not find this keyword
=dens2
CONTOUR
compulsory keyword the value we would like to draw the contour at in the space
=0.42
SEARCHDIR
compulsory keyword In which directions do you wish to search for the contour.
=z
STRIDE
could not find this keyword
=1
CLEAR
could not find this keyword
=1 DUMPGRID
GRID
the grid you would like to print (can also use ARG for specifying what is being printed)
=ss2
FILE
compulsory keyword ( default=density ) the file on which to write the grid.
=contour2.dat
FMT
the format that should be used to output real numbers
=%8.4f
STRIDE
compulsory keyword ( default=0 ) the frequency with which the grid should be output to the file.
=1
Glossary of keywords and components
Compulsory keywords
CONTOUR the value we would like to draw the contour at in the space
INTERPOLATION_TYPE ( default=spline ) the method to use for interpolation. Can be spline, linear, ceiling or floor.
SEARCHDIR In which directions do you wish to search for the contour.
Options
SERIAL ( default=off ) do the calculation in serial. Do not parallelize
ZERO_OUTSIDE_GRID_RANGE

( default=off ) if we are asked to evaluate the function for a number that is outside the range of the grid set it to zero

ARG the input for this action is the scalar output from one or more other actions. The particular scalars that you will use are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the scalars calculated by all the proceeding actions in the input file are taken. Some actions have multi-component outputs and each component of the output has a specific label. For example a DISTANCE action labelled dist may have three components x, y and z. To take just the x component you should use dist.x, if you wish to take all three components then use dist.*.More information on the referencing of Actions can be found in the section of the manual on the PLUMED Getting Started. Scalar values can also be referenced using POSIX regular expressions as detailed in the section on Regular Expressions. To use this feature you you must compile PLUMED with the appropriate flag.. You can use multiple instances of this keyword i.e. ARG1, ARG2, ARG3...