FIND_SPHERICAL_CONTOUR

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 numbersc1:COORDINATIONNUMBERSPECIES=1-512this keyword is used for colvars such as coordination number.SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0} # Select coordination numbers that are more than 2.0This keyword is used if you want to employ an alternative to the continuous switching function defined above.cf:MFILTER_MOREDATA=compulsory keywordThe multicolvar that calculates the set of base quantities that we are interested inc1SWITCH={RATIONAL D_0=2.0 R_0=0.1}This keyword is used if you want to employ an alternative to the continuous switching function defined above.LOWMEM# Build a contact matrix( default=off ) lower the memory requirementsmat:CONTACT_MATRIXATOMS=The list of atoms for which you would like to calculate the contact matrix.cfSWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0} # Find largest clusterThis keyword is used if you want to employ an alternative to the continuous switching function defined above.dfs:DFSCLUSTERINGMATRIX=compulsory keywordthe action that calculates the adjacency matrix vessel we would like to analyzematLOWMEM( default=off ) lower the memory requirementsclust1:CLUSTER_PROPERTIESCLUSTERS=compulsory keywordthe label of the action that does the clusteringdfsCLUSTER=1 # Find center of largest 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.trans1:MTRANSFORM_MOREDATA=compulsory keywordThe multicolvar that calculates the set of base quantities that we are interested inclust1LOWMEM( default=off ) lower the memory requirementscent:CENTER_OF_MULTICOLVARDATA=compulsory keywordfind the average value for a multicolvartrans1# Calculate the phase field of the coordinationdens:MULTICOLVARDENSDATA=compulsory keywordthe multicolvar which you would like to calculate the density profile fortrans1ORIGIN=we will use the position of this atom as the origin.centDIR=xyzcompulsory keywordthe direction in which to calculate the density profileNBINS=30,30,30the number of bins to use to represent the density profileBANDWIDTH=2.0,2.0,2.0 # Find the isocontour around the nucleuscompulsory keywordthe bandwidths for kernel density estimationsc:FIND_SPHERICAL_CONTOURGRID=compulsory keywordthe action that creates the input grid you would like to usedensCONTOUR=0.85compulsory keywordthe value we would like to draw the contour at in the spaceINNER_RADIUS=10.0compulsory keywordthe minimum radius on which to look for the contourOUTER_RADIUS=40.0compulsory keywordthe outer radius on which to look for the contourNPOINTS=100 # And print the grid to a file GRID_TO_XYZcompulsory keywordthe number of points for which we are looking for the contourGRID=compulsory keywordthe action that creates the grid you would like to outputscFILE=mysurface.xyzcompulsory keyword ( default=density )the file on which to write the grid.UNITS=Acompulsory keyword ( default=PLUMED )the units in which to print out the coordinates.

- Glossary of keywords and components

- Compulsory keywords

STRIDE | ( default=1 ) the frequency with which the data should be collected and added to the quantity being averaged |

CLEAR | ( default=0 ) the frequency with which to clear all the accumulated data. The default value of 0 implies that all the data will be used and that the grid will never be cleared |

NORMALIZATION | ( default=true ) This controls how the data is normalized it can be set equal to true, false or ndata. The differences between these options are explained in the manual page for HISTOGRAM |

GRID | the action that creates the input grid you would like to use |

CONTOUR | the value we would like to draw the contour at in the space |

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 use MPI |

LOWMEM | ( default=off ) lower the memory requirements |

TIMINGS | ( default=off ) output information on the timings of the various parts of the calculation |

LOGWEIGHTS | list of actions that calculates log weights that should be used to weight configurations when calculating averages |

CONCENTRATION | the concentration parameter for Von Mises-Fisher distributions |

COMPONENT | if your input is a vector field use this to specify the component of the input vector field for which you wish to use |