Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2019 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "core/ActionRegister.h"
23 : #include "ContourFindingBase.h"
24 : #include "tools/Random.h"
25 :
26 : //+PLUMEDOC GRIDANALYSIS FIND_SPHERICAL_CONTOUR
27 : /*
28 : Find an isocontour in a three dimensional grid by searching over a Fibonacci sphere.
29 :
30 : As discussed in the part of the manual on \ref Analysis PLUMED contains a number of tools that allow you to calculate
31 : a function on a grid. The function on this grid might be a \ref HISTOGRAM as a function of a few collective variables
32 : or it might be a phase field that has been calcualted using \ref MULTICOLVARDENS. If this function has one or two input
33 : arguments it is relatively straightforward to plot the function. If by contrast the data has a three dimensions it can be
34 : difficult to visualize.
35 :
36 : This action provides one tool for visualizing these functions. It can be used to search for a set of points on a contour
37 : wher the function takes a particular value. In other words, for the function \f$f(x,y,z)\f$ this action would find a set
38 : of points \f$\{x_c,y_c,z_c\}\f$ that have:
39 :
40 : \f[
41 : f(x_c,y_c,z_c) - c = 0
42 : \f]
43 :
44 : where \f$c\f$ is some constant value that is specified by the user. The points on this contour are find by searching along a
45 : set of equally spaced radii of a sphere that centered at on particular, user-speciified atom or virtual atom. To ensure that
46 : these search radii are equally spaced on the surface of the sphere the search directions are generated by using a fibonacci
47 : spiral projected on a sphere. In other words, the search directions are given by:
48 :
49 : \f[
50 : \mathbf{r}_i = \left(
51 : \begin{matrix}
52 : \sqrt{1 - y^2} \cos(\phi) \\
53 : \frac{2i}{n} - 1 + \frac{1}{n} \\
54 : \sqrt{1 - y^2} \sin(\phi)
55 : \end{matrix}
56 : \right)
57 : \f]
58 :
59 : where \f$y\f$ is the quantity second component of the vector defined above, \f$n\f$ is the number of directions to look in and \f$\phi\f$ is
60 :
61 : \f[
62 : \phi = \mod(i + R, n) \pi ( 3 - \sqrt{5} )
63 : \f]
64 :
65 : where \f$R\f$ is a random variable between 0 and \f$n-1\f$ that is generated during the read in of the input file and that is fixed during
66 : the whole calculation.
67 :
68 : It is important to note that this action can only be used to detect countours in three dimensional functions. In addition, this action will fail to
69 : 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
70 : function have a spherical topology you should use \ref FIND_CONTOUR in place of \ref FIND_SPHERICAL_CONTOUR.
71 :
72 : \par Examples
73 :
74 : 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
75 : condensation of Lennard Jones from the vapour. The input below achieves this by calculating the coordination numbers of all the atoms within the gas.
76 : 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.
77 : As such we can detect the sizes of the droplets by constructing a \ref CONTACT_MATRIX whose \f$ij\f$ element tells us whether atom \f$i\f$ and atom \f$j\f$
78 : have coordination number that is greater that two. The atoms within the various droplets within the system can then be found by performing a
79 : \ref 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
80 : by exploiting the functionality within \ref CENTER_OF_MULTICOLVAR. We can then construct a phase field based on the positions of the atoms in the largest
81 : 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
82 : teh droplet from the surrounding gas. The value of the phase field on this isocontour is equal to 0.75.
83 :
84 : \plumedfile
85 : # Calculate coordination numbers
86 : c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
87 : # Select coordination numbers that are more than 2.0
88 : cf: MFILTER_MORE DATA=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
89 : # Build a contact matrix
90 : mat: CONTACT_MATRIX ATOMS=cf SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
91 : # Find largest cluster
92 : dfs: DFSCLUSTERING MATRIX=mat
93 : clust1: CLUSTER_PROPERTIES CLUSTERS=dfs CLUSTER=1
94 : # Find center of largest cluster
95 : trans1: MTRANSFORM_MORE DATA=clust1 SWITCH={RATIONAL D_0=2.0 R_0=0.1} LOWMEM
96 : cent: CENTER_OF_MULTICOLVAR DATA=trans1
97 : # Calculate the phase field of the coordination
98 : dens: MULTICOLVARDENS DATA=trans1 ORIGIN=cent DIR=xyz NBINS=30,30,30 BANDWIDTH=2.0,2.0,2.0
99 : # Find the isocontour around the nucleus
100 : FIND_SPHERICAL_CONTOUR GRID=dens CONTOUR=0.85 INNER_RADIUS=10.0 OUTER_RADIUS=40.0 FILE=mysurface.xyz UNITS=A PRECISION=4 NPOINTS=100
101 : \endplumedfile
102 :
103 : */
104 : //+ENDPLUMEDOC
105 :
106 : namespace PLMD {
107 : namespace gridtools {
108 :
109 2 : class FindSphericalContour : public ContourFindingBase {
110 : private:
111 : unsigned nbins;
112 : double min, max;
113 : public:
114 : static void registerKeywords( Keywords& keys );
115 : explicit FindSphericalContour(const ActionOptions&ao);
116 4 : unsigned getNumberOfQuantities() const { return 2; }
117 : void compute( const unsigned& current, MultiValue& myvals ) const ;
118 : };
119 :
120 6453 : PLUMED_REGISTER_ACTION(FindSphericalContour,"FIND_SPHERICAL_CONTOUR")
121 :
122 2 : void FindSphericalContour::registerKeywords( Keywords& keys ) {
123 2 : ContourFindingBase::registerKeywords( keys );
124 8 : keys.add("compulsory","NPOINTS","the number of points for which we are looking for the contour");
125 8 : keys.add("compulsory","INNER_RADIUS","the minimum radius on which to look for the contour");
126 8 : keys.add("compulsory","OUTER_RADIUS","the outer radius on which to look for the contour");
127 10 : keys.add("compulsory","NBINS","1","the number of discrete sections in which to divide the distance between the inner and outer radius when searching for a contour");
128 2 : }
129 :
130 1 : FindSphericalContour::FindSphericalContour(const ActionOptions&ao):
131 : Action(ao),
132 1 : ContourFindingBase(ao)
133 : {
134 2 : if( ingrid->getDimension()!=3 ) error("input grid must be three dimensional");
135 :
136 2 : unsigned npoints; parse("NPOINTS",npoints);
137 1 : log.printf(" searching for %u points on dividing surface \n",npoints);
138 4 : parse("INNER_RADIUS",min); parse("OUTER_RADIUS",max); parse("NBINS",nbins);
139 1 : log.printf(" expecting to find dividing surface at radii between %f and %f \n",min,max);
140 1 : log.printf(" looking for contour in windows of length %f \n", (max-min)/nbins);
141 : // Set this here so the same set of grid points are used on every turn
142 2 : std::string vstring = "TYPE=fibonacci COMPONENTS=" + getLabel() + " COORDINATES=x,y,z PBC=F,F,F";
143 2 : createGrid( "grid", vstring ); mygrid->setNoDerivatives();
144 1 : setAveragingAction( mygrid, true ); mygrid->setupFibonacciGrid( npoints );
145 :
146 1 : checkRead();
147 : // Create a task list
148 1 : for(unsigned i=0; i<npoints; ++i) addTaskToList( i );
149 1 : deactivateAllTasks();
150 201 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
151 1 : lockContributors();
152 1 : }
153 :
154 100 : void FindSphericalContour::compute( const unsigned& current, MultiValue& myvals ) const {
155 : // Generate contour point on inner sphere
156 100 : std::vector<double> contour_point(3), direction(3), der(3), tmp(3);
157 : // Retrieve this contour point from grid
158 100 : mygrid->getGridPointCoordinates( current, direction );
159 : // Now setup contour point on inner sphere
160 700 : for(unsigned j=0; j<3; ++j) {
161 900 : contour_point[j] = min*direction[j];
162 600 : direction[j] = (max-min)*direction[j] / static_cast<double>(nbins);
163 : }
164 : bool found=false;
165 100 : for(unsigned k=0; k<nbins; ++k) {
166 1300 : for(unsigned j=0; j<3; ++j) tmp[j] = contour_point[j] + direction[j];
167 : double val1 = getDifferenceFromContour( contour_point, der );
168 : double val2 = getDifferenceFromContour( tmp, der );
169 100 : if( val1*val2<0 ) {
170 : findContour( direction, contour_point );
171 700 : double norm=0; for(unsigned j=0; j<3; ++j) norm += contour_point[j]*contour_point[j];
172 100 : myvals.setValue( 1, sqrt(norm) ); found=true; break;
173 : }
174 0 : for(unsigned j=0; j<3; ++j) contour_point[j] = tmp[j];
175 : }
176 0 : if( !found ) error("range does not bracket the dividing surface");
177 100 : }
178 :
179 : }
180 4839 : }
|
