	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 in a three dimensional grid by searching over a Fibonacci sphere.

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 a set of equally spaced radii of a sphere that centered at on particular, user-specified atom or virtual atom. To ensure that these search radii are equally spaced on the surface of the sphere the search directions are generated by using a Fibonacci spiral projected on a sphere. In other words, the search directions are given by:

\[ \mathbf{r}_i = \left( \begin{matrix} \sqrt{1 - y^2} \cos(\phi) \\ \frac{2i}{n} - 1 + \frac{1}{n} \\ \sqrt{1 - y^2} \sin(\phi) \end{matrix} \right) \]

where \(y\) is the quantity second component of the vector defined above, \(n\) is the number of directions to look in and \(\phi\) is

\[ \phi = \mod(i + R, n) \pi ( 3 - \sqrt{5} ) \]

where \(R\) is a random variable between 0 and \(n-1\) that is generated during the read in of the input file and that is fixed during the whole calculation.

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 a sphere. If you are uncertain that the isocontours in your function have a spherical topology you should use FIND_CONTOUR in place of FIND_SPHERICAL_CONTOUR.

Examples

The following input demonstrates how this action can be used. The input here is used to study the shape of a droplet that has been formed during the condensation of Lennard Jones from the vapor. The input below achieves this by calculating the coordination numbers of all the atoms within the gas. Obviously, those atoms within the droplet will have a large value for the coordination number while the isolated atoms in the gas will have a low value. As such we can detect the sizes of the droplets by constructing a CONTACT_MATRIX whose \(ij\) element tells us whether atom \(i\) and atom \(j\) have coordination number that is greater that two. The atoms within the various droplets within the system can then be found by performing a DFSCLUSTERING on this matrix to detect the connected components. We can take the largest of these connected components and find the center of the droplet by exploiting the functionality within CENTER_OF_MULTICOLVAR. We can then construct a phase field based on the positions of the atoms in the largest cluster and the values of the coordination numbers of these atoms. The final line in the input then finds the a set of points on the dividing surface that separates the droplet from the surrounding gas. The value of the phase field on this isocontour is equal to 0.75.

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

# Calculate coordination numbers
c1: COORDINATIONNUMBER SPECIESthis keyword is used for colvars such as coordination number. 
=1-512 SWITCHthe switching function that it used in the construction of the contact matrix 
={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}   You cannot view the components that are calculated by each action for this input file. Sorry 
# Select coordination numbers that are more than 2.0
cf: MFILTER_MORE DATAcompulsory keyword 
the vector you wish to transform 
=c1 SWITCHcompulsory keyword 
the switching function that transform 
={RATIONAL D_0=2.0 R_0=0.1}  LOWMEM( default=off ) this flag does nothing and is present only to ensure back-compatibility
  You cannot view the components that are calculated by each action for this input file. Sorry 
# Build a contact matrix
mat: CONTACT_MATRIX ATOMSthe atoms for which you would like to calculate the adjacency matrix (basically equivalent
to GROUP) 
=cf SWITCHspecify the switching function to use between two sets of indistinguishable atoms.
={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}   You cannot view the components that are calculated by each action for this input file. Sorry 
# Find largest cluster
dfs: DFSCLUSTERING MATRIXthe input matrix (can use ARG instead) 
=mat LOWMEM( default=off ) this flag does nothing and is present only to ensure back-compatibility
  You cannot view the components that are calculated by each action for this input file. Sorry 
clust1: CLUSTER_PROPERTIES CLUSTERScompulsory keyword 
the label of the action that does the clustering 
=dfs CLUSTERcompulsory keyword ( default=1 )
which cluster would you like to look at 1 is the largest cluster, 2 is the second
largest, 3 is the the third largest and so on. 
=1  You cannot view the components that are calculated by each action for this input file. Sorry 
# Find center of largest cluster
trans1: MTRANSFORM_MORE DATA could not find this keyword 
=clust1 SWITCH could not find this keyword 
={RATIONAL D_0=2.0 R_0=0.1}  LOWMEM could not find this keyword 
  You cannot view the components that are calculated by each action for this input file. Sorry 
cent: CENTER_OF_MULTICOLVAR DATA could not find this keyword 
=trans1  You cannot view the components that are calculated by each action for this input file. Sorry 
# Calculate the phase field of the coordination
dens: MULTICOLVARDENS DATAthe multicolvar which you would like to calculate the density profile for 
=trans1 ORIGINcompulsory keyword 
we will use the position of this atom as the origin 
=cent DIRcompulsory keyword 
the direction in which to calculate the density profile 
=xyz NBINSthe number of bins to use in each direction (alternative to GRID_NBIN) 
=30,30,30 BANDWIDTHthe bandwidths for kernel density esimtation 
=2.0,2.0,2.0  You cannot view the components that are calculated by each action for this input file. Sorry 
# Find the isocontour around the nucleus
sc: FIND_SPHERICAL_CONTOUR GRID could not find this keyword 
=dens CONTOURcompulsory keyword 
the value we would like to draw the contour at in the space 
=0.85 INNER_RADIUScompulsory keyword 
the minimum radius on which to look for the contour 
=10.0 OUTER_RADIUScompulsory keyword 
the outer radius on which to look for the contour 
=40.0 NPOINTScompulsory keyword 
the number of points for which we are looking for the contour 
=100  You cannot view the components that are calculated by each action for this input file. Sorry 
# And print the grid to a file
GRID_TO_XYZ GRID could not find this keyword 
=sc FILE could not find this keyword 
=mysurface.xyz UNITS could not find this keyword 
=A  You cannot view the components that are calculated by each action for this input file. Sorry

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.
NPOINTS	the number of points for which we are looking for the contour
INNER_RADIUS	the minimum radius on which to look for the contour
OUTER_RADIUS	the outer radius on which to look for the contour
NBINS	( default=1 ) the number of discrete sections in which to divide the distance between the inner and outer radius when searching for a 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...