Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 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 :
25 : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR
26 : /*
27 : Find an isocontour in a smooth function.
28 :
29 : As discussed in the part of the manual on \ref Analysis PLUMED contains a number of tools that allow you to calculate
30 : a function on a grid. The function on this grid might be a \ref HISTOGRAM as a function of a few collective variables
31 : or it might be a phase field that has been calcualted using \ref MULTICOLVARDENS. If this function has one or two input
32 : arguments it is relatively straightforward to plot the function. If by contrast the data has a three or more dimensions
33 : it can be difficult to visualize.
34 :
35 : This action provides one tool for visualizing these functions. It can be used to search for a set of points on a contour
36 : where the function takes a particular values. In other words, for the function \f$f(x,y)\f$ this action would find a set
37 : of points \f$\{x_c,y_c\}\f$ that have:
38 :
39 : \f[
40 : f(x_c,y_c) - c = 0
41 : \f]
42 :
43 : where \f$c\f$ is some constant value that is specified by the user. The points on this contour are detected using a variant
44 : on the marching squares or marching cubes algorithm, which you can find information on here:
45 :
46 : https://en.wikipedia.org/wiki/Marching_squares
47 : https://en.wikipedia.org/wiki/Marching_cubes
48 :
49 : As such, and unlike \ref FIND_CONTOUR_SURFACE or \ref FIND_SPHERICAL_CONTOUR, the function input to this action can have any dimension.
50 : Furthermore, the topology of the contour will be determined by the algorithm and does not need to be specified by the user.
51 :
52 : \par Examples
53 :
54 : The input below allows you to calculate something akin to a Willard-Chandler dividing surface \cite wcsurface.
55 : The simulation cell in this case contains a solid phase and a liquid phase. The Willard-Chandler surface is the
56 : surface that separates the parts of the box containing the solid from the parts containing the liquid. To compute the position
57 : of this surface the \ref FCCUBIC symmetry function is calculated for each of the atoms in the system from on the geometry of the
58 : atoms in the first coordination sphere of each of the atoms. These quantities are then transformed using a switching function.
59 : This procedure generates a single number for each atom in the system and this quantity has a value of one for atoms that are in
60 : parts of the box that resemble the solid structure and zero for atoms that are in parts of the box that resemble the liquid.
61 : The position of a virtual atom is then computed using \ref CENTER_OF_MULTICOLVAR and a phase field model is constructed using
62 : \ref MULTICOLVARDENS. These procedure ensures that we have a continuous function that gives a measure of the average degree of
63 : solidness at each point in the simulation cell. The Willard-Chandler dividing surface is calculated by finding a a set of points
64 : at which the value of this phase field is equal to 0.5. This set of points is output to file called mycontour.dat. A new contour
65 : is found on every single step for each frame that is read in.
66 :
67 : \verbatim
68 : UNITS NATURAL
69 : FCCUBIC ...
70 : SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
71 : ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619 LABEL=fcc
72 : ... FCCUBIC
73 :
74 : tfcc: MTRANSFORM_MORE DATA=fcc SWITCH={SMAP R_0=0.5 A=8 B=8}
75 : center: CENTER_OF_MULTICOLVAR DATA=tfcc
76 :
77 : MULTICOLVARDENS ...
78 : DATA=tfcc ORIGIN=center DIR=xyz LABEL=dens
79 : NBINS=80,80,80 BANDWIDTH=1.0,1.0,1.0 STRIDE=25
80 : LABEL=dens STRIDE=1 CLEAR=1
81 : ... MULTICOLVARDENS
82 :
83 : FIND_CONTOUR GRID=dens CONTOUR=0.5 FILE=mycontour.dat
84 : \endverbatim
85 :
86 : */
87 : //+ENDPLUMEDOC
88 :
89 : namespace PLMD {
90 : namespace gridtools {
91 :
92 2 : class FindContour : public ContourFindingBase {
93 : private:
94 : bool firsttime;
95 : unsigned gbuffer;
96 : public:
97 : static void registerKeywords( Keywords& keys );
98 : explicit FindContour(const ActionOptions&ao);
99 0 : bool checkAllActive() const { return gbuffer==0; }
100 : void prepareForAveraging();
101 0 : bool isPeriodic() { return false; }
102 : void compute( const unsigned& current, MultiValue& myvals ) const ;
103 : void finishAveraging();
104 : };
105 :
106 2524 : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
107 :
108 2 : void FindContour::registerKeywords( Keywords& keys ) {
109 2 : ContourFindingBase::registerKeywords( keys );
110 : // We want a better way of doing this bit
111 2 : keys.add("compulsory","BUFFER","0","number of buffer grid points around location where grid was found on last step. If this is zero the full grid is calculated on each step");
112 2 : }
113 :
114 1 : FindContour::FindContour(const ActionOptions&ao):
115 : Action(ao),
116 : ContourFindingBase(ao),
117 1 : firsttime(true)
118 : {
119 :
120 1 : parse("BUFFER",gbuffer);
121 1 : if( gbuffer>0 ) log.printf(" after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
122 1 : checkRead();
123 1 : }
124 :
125 2 : void FindContour::prepareForAveraging() {
126 : // Create a task list if first time
127 2 : if( firsttime ) {
128 1 : for(unsigned i=0; i<ingrid->getDimension()*ingrid->getNumberOfPoints(); ++i) addTaskToList( i );
129 : }
130 2 : firsttime=false; deactivateAllTasks();
131 :
132 : // We now need to identify the grid points that we need to search through
133 2 : std::vector<unsigned> nbin( ingrid->getNbin() );
134 4 : std::vector<unsigned> ind( ingrid->getDimension() );
135 4 : std::vector<unsigned> ones( ingrid->getDimension(), 1 );
136 4 : unsigned num_neighbours; std::vector<unsigned> neighbours;
137 10978 : for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
138 : // Ensure inactive grid points are ignored
139 10976 : if( ingrid->inactive(i) ) continue;
140 :
141 : // Get the index of the current grid point
142 10976 : ingrid->getIndices( i, ind );
143 10976 : ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
144 10976 : bool cycle=false;
145 307328 : for(unsigned j=0; j<num_neighbours; ++j) {
146 296352 : if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
147 : }
148 10976 : if( cycle ) continue;
149 :
150 : // Get the value of a point on the grid
151 10976 : double val1=getFunctionValue( i ) - contour;
152 10976 : bool edge=false;
153 43904 : for(unsigned j=0; j<ingrid->getDimension(); ++j) {
154 : // Make sure we don't search at the edge of the grid
155 32928 : if( !ingrid->isPeriodic(j) && (ind[j]+1)==nbin[j] ) continue;
156 32928 : else if( (ind[j]+1)==nbin[j] ) { edge=true; ind[j]=0; }
157 30968 : else ind[j]+=1;
158 32928 : double val2=getFunctionValue( ind ) - contour;
159 32928 : if( val1*val2<0 ) taskFlags[ ingrid->getDimension()*i + j ] = 1;
160 32928 : if( ingrid->isPeriodic(j) && edge ) { edge=false; ind[j]=nbin[j]-1; }
161 30968 : else ind[j]-=1;
162 : }
163 : }
164 4 : lockContributors();
165 2 : }
166 :
167 554 : void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
168 : // Retrieve the initial grid point coordinates
169 554 : unsigned gpoint = std::floor( current / ingrid->getDimension() );
170 554 : std::vector<double> point( ingrid->getDimension() );
171 554 : ingrid->getGridPointCoordinates( gpoint, point );
172 :
173 : // Retrieve the direction we are searching for the contour
174 554 : unsigned gdir = current%(ingrid->getDimension() );
175 1108 : std::vector<double> direction( ingrid->getDimension(), 0 );
176 554 : direction[gdir] = 0.999999999*ingrid->getGridSpacing()[gdir];
177 :
178 : // Now find the contour
179 554 : findContour( direction, point );
180 : // And transfer to the store data vessel
181 1108 : for(unsigned i=0; i<ingrid->getDimension(); ++i) myvals.setValue( 1+i, point[i] );
182 554 : }
183 :
184 1 : void FindContour::finishAveraging() {
185 1 : ContourFindingBase::finishAveraging();
186 : // And update the list of active grid points
187 1 : if( gbuffer>0 ) {
188 0 : std::vector<unsigned> neighbours; unsigned num_neighbours;
189 0 : std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
190 0 : std::vector<bool> active( ingrid->getNumberOfPoints(), false );
191 0 : std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
192 0 : for(unsigned i=0; i<getCurrentNumberOfActiveTasks(); ++i) {
193 : // Get the point we are operating on
194 0 : unsigned ipoint = std::floor( getActiveTask(i) / ingrid->getDimension() );
195 : // Get the indices of this point
196 0 : ingrid->getIndices( ipoint, ugrid_indices );
197 : // Now activate buffer region
198 0 : ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
199 0 : for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
200 : }
201 0 : ingrid->activateThesePoints( active );
202 : }
203 1 : }
204 :
205 : }
206 2523 : }
|
