Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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_SURFACE
26 : /*
27 : Find an isocontour by searching along either the x, y or z direction.
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 dimensions it can be
33 : 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 : 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
37 : of points \f$\{x_c,y_c,z_c\}\f$ that have:
38 :
39 : \f[
40 : f(x_c,y_c,z_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 find by searching along lines
44 : that run parallel to the \f$x\f$, \f$y\f$ or \f$z\f$ axis of the simulation cell. The result is, therefore, a two dimensional
45 : function evaluated on a grid that gives us the height of the interface as a function of two coordinates.
46 :
47 : 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
48 : 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
49 : function have the appropriate topology you should use \ref FIND_CONTOUR in place of \ref FIND_CONTOUR_SURFACE.
50 :
51 :
52 : \par Examples
53 :
54 : The input shown below was used to analyse the results from a simulation of an interface between solid and molten Lennard Jones. The interface between
55 : the solid and the liquid was set up in the plane perpendicular to the \f$z\f$ direction of the simulation cell. The input below calculates something
56 : akin to a Willard-Chandler dividing surface \cite wcsurface between the solid phase and the liquid phase. There are two of these interfaces within the
57 : 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
58 : 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
59 : \ref FCCUBIC symmetry function using the \ref MULTICOLVARDENS action. Notice that we use the fact that we know roughly where the interface is when specifying
60 : how this phase field is to be calculated and specify the region over the \f$z\f$-axis in which we are going to search for the phase field in the line defining
61 : the \ref MULTICOLVARDENS. Once we have calculated the phase field we search for contour points on the lines that run parallel to the \f$z\f$-direction of the cell
62 : box using the FIND_CONTOUR_SURFACE command. The final result is a \f$14 \times 14\f$ grid of values for the height of the interface as a function of the \f$(x,y)\f$
63 : position. This grid is then output to a file called contour2.dat.
64 :
65 : 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
66 : on every step.
67 :
68 : \verbatim
69 : UNITS NATURAL
70 : FCCUBIC ...
71 : SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
72 : ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619 LABEL=fcc
73 : ... FCCUBIC
74 :
75 : dens2: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,50 ZREDUCED ZLOWER=6.0 ZUPPER=11.0 BANDWIDTH=1.0,1.0,1.0 CLEAR=1
76 :
77 : ss2: FIND_CONTOUR_SURFACE GRID=dens2 CONTOUR=0.42 SEARCHDIR=z STRIDE=1 CLEAR=1
78 : DUMPGRID GRID=ss2 FILE=contour2.dat FMT=%8.4f STRIDE=1
79 : \endverbatim
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 : namespace PLMD {
85 : namespace gridtools {
86 :
87 2 : class FindContourSurface : public ContourFindingBase {
88 : private:
89 : bool firsttime;
90 : unsigned dir_n;
91 : unsigned gbuffer;
92 : std::vector<unsigned> ones;
93 : std::vector<unsigned> gdirs;
94 : std::vector<double> direction;
95 : public:
96 : static void registerKeywords( Keywords& keys );
97 : explicit FindContourSurface(const ActionOptions&ao);
98 6 : unsigned getNumberOfQuantities() const { return 2; }
99 0 : bool checkAllActive() const { return gbuffer==0; }
100 : void clearAverage();
101 : void prepareForAveraging();
102 : void compute( const unsigned& current, MultiValue& myvals ) const ;
103 : void finishAveraging();
104 : };
105 :
106 2524 : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
107 :
108 2 : void FindContourSurface::registerKeywords( Keywords& keys ) {
109 2 : ContourFindingBase::registerKeywords( keys );
110 2 : keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
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 : keys.remove("FILE"); keys.remove("UNITS"); keys.remove("PRECISION");
113 2 : }
114 :
115 1 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
116 : Action(ao),
117 : ContourFindingBase(ao),
118 : firsttime(true),
119 1 : ones(ingrid->getDimension(),1)
120 : {
121 1 : if( ingrid->getDimension()<2 ) error("cannot find dividing surface if input grid is one dimensional");
122 :
123 1 : std::string dir; parse("SEARCHDIR",dir); parse("BUFFER",gbuffer);
124 1 : log.printf(" calculating location of contour on %d dimensional grid \n", ingrid->getDimension()-1 );
125 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);
126 1 : checkRead();
127 :
128 1 : unsigned n=0; gdirs.resize( ingrid->getDimension()-1 );
129 4 : for(unsigned i=0; i<ingrid->getDimension(); ++i) {
130 3 : if( ingrid->getComponentName(i)==dir ) {
131 1 : dir_n=i;
132 : } else {
133 2 : if( n==gdirs.size() ) error("could not find " + dir + " direction in input grid");
134 2 : gdirs[n]=i; n++;
135 : }
136 : }
137 1 : if( n!=(ingrid->getDimension()-1) ) error("output of grid is not understood");
138 :
139 : // Create the input from the old string
140 2 : std::string vstring = "COMPONENTS=" + getLabel() + " COORDINATES=" + ingrid->getComponentName( gdirs[0] );
141 1 : for(unsigned i=1; i<gdirs.size(); ++i) vstring += "," + ingrid->getComponentName( gdirs[i] );
142 1 : vstring += " PBC=";
143 1 : if( ingrid->isPeriodic(gdirs[0]) ) vstring+="T";
144 0 : else vstring+="F";
145 2 : for(unsigned i=1; i<gdirs.size(); ++i) {
146 1 : if( ingrid->isPeriodic(gdirs[i]) ) vstring+=",T"; else vstring+=",F";
147 : }
148 1 : createGrid( "grid", vstring ); mygrid->setNoDerivatives();
149 2 : setAveragingAction( mygrid, true );
150 1 : }
151 :
152 2 : void FindContourSurface::clearAverage() {
153 : // Set the boundaries of the output grid
154 4 : std::vector<double> fspacing; std::vector<unsigned> snbins( ingrid->getDimension()-1 );
155 4 : std::vector<std::string> smin( ingrid->getDimension()-1 ), smax( ingrid->getDimension()-1 );
156 6 : for(unsigned i=0; i<gdirs.size(); ++i) {
157 4 : smin[i]=ingrid->getMin()[gdirs[i]]; smax[i]=ingrid->getMax()[gdirs[i]];
158 4 : snbins[i]=ingrid->getNbin()[gdirs[i]];
159 : }
160 2 : mygrid->setBounds( smin, smax, snbins, fspacing); resizeFunctions();
161 4 : ActionWithAveraging::clearAverage();
162 2 : }
163 :
164 2 : void FindContourSurface::prepareForAveraging() {
165 : // Create a task list if first time
166 2 : if( firsttime ) {
167 1 : std::vector<unsigned> find( ingrid->getDimension() );
168 2 : std::vector<unsigned> ind( mygrid->getDimension() );
169 197 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
170 196 : find.assign( find.size(), 0 ); mygrid->getIndices( i, ind );
171 196 : for(unsigned j=0; j<gdirs.size(); ++j) find[gdirs[j]]=ind[j];
172 : // Current will be set equal to the start point for this grid index
173 196 : addTaskToList( ingrid->getIndex(find) );
174 : }
175 : // And prepare the task list
176 1 : deactivateAllTasks();
177 1 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
178 1 : lockContributors();
179 : // Set the direction in which to look for the contour
180 1 : direction.resize( ingrid->getDimension(), 0 );
181 2 : direction[dir_n] = 0.999999999*ingrid->getGridSpacing()[dir_n];
182 : }
183 2 : }
184 :
185 2 : void FindContourSurface::finishAveraging() {
186 2 : ContourFindingBase::finishAveraging();
187 : // And update the list of active grid points
188 2 : if( gbuffer>0 ) {
189 2 : std::vector<double> dx( ingrid->getGridSpacing() );
190 4 : std::vector<double> point( ingrid->getDimension() );
191 4 : std::vector<double> lpoint( mygrid->getDimension() );
192 4 : std::vector<unsigned> neighbours; unsigned num_neighbours;
193 4 : std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
194 4 : std::vector<bool> active( ingrid->getNumberOfPoints(), false );
195 4 : std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
196 394 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
197 : // Retrieve the coordinates of this grid point
198 392 : mygrid->getGridPointCoordinates( i, lpoint );
199 392 : point[dir_n] = mygrid->getGridElement( i, 0 );
200 : // 0.5*dx added here to prevent problems with flooring of grid points
201 392 : for(unsigned j=0; j<gdirs.size(); ++j) point[gdirs[j]]=lpoint[j] + 0.5*dx[gdirs[j]];
202 392 : ingrid->getIndices( point, ugrid_indices );
203 : // Now activate buffer region
204 392 : ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
205 392 : for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
206 : }
207 4 : ingrid->activateThesePoints( active );
208 : }
209 2 : firsttime=false;
210 2 : }
211 :
212 390 : void FindContourSurface::compute( const unsigned& current, MultiValue& myvals ) const {
213 390 : std::vector<unsigned> neighbours; unsigned num_neighbours; unsigned nfound=0; double minv=0, minp=0;
214 782 : std::vector<unsigned> bins_n( ingrid->getNbin() ); unsigned shiftn=current;
215 783 : std::vector<unsigned> ind( ingrid->getDimension() ); std::vector<double> point( ingrid->getDimension() );
216 : #ifndef DNDEBUG
217 783 : std::vector<unsigned> oind( mygrid->getDimension() ); mygrid->getIndices( current, oind );
218 : #endif
219 11166 : for(unsigned i=0; i<bins_n[dir_n]; ++i) {
220 : #ifndef DNDEBUG
221 11164 : std::vector<unsigned> base_ind( ingrid->getDimension() ); ingrid->getIndices( shiftn, base_ind );
222 11164 : for(unsigned j=0; j<gdirs.size(); ++j) plumed_dbg_assert( base_ind[gdirs[j]]==oind[j] );
223 : #endif
224 : // Ensure inactive grid points are ignored
225 11165 : if( ingrid->inactive( shiftn ) ) { shiftn += ingrid->getStride()[dir_n]; continue; }
226 : // Get the index of the current grid point
227 6720 : ingrid->getIndices( shiftn, ind );
228 : // Exit if we are at the edge of the grid
229 6720 : if( !ingrid->isPeriodic(dir_n) && (ind[dir_n]+1)==bins_n[dir_n] ) {
230 0 : shiftn += ingrid->getStride()[dir_n]; continue;
231 : }
232 :
233 : // Ensure points with inactive neighbours are ignored
234 6719 : ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
235 6707 : bool cycle=false;
236 166674 : for(unsigned j=0; j<num_neighbours; ++j) {
237 160749 : if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
238 : }
239 6720 : if( cycle ) { shiftn += ingrid->getStride()[dir_n]; continue; }
240 :
241 : // Now get the function value at two points
242 5925 : double val1=getFunctionValue( shiftn ) - contour; double val2;
243 5925 : if( (ind[dir_n]+1)==bins_n[dir_n] ) val2 = getFunctionValue( current ) - contour;
244 5924 : else val2=getFunctionValue( shiftn + ingrid->getStride()[dir_n] ) - contour;
245 :
246 : // Check if the minimum is bracketed
247 5925 : if( val1*val2<0 ) {
248 392 : ingrid->getGridPointCoordinates( shiftn, point ); findContour( direction, point );
249 392 : minp=point[dir_n]; nfound++; break;
250 : }
251 :
252 :
253 : // This moves us on to the next point
254 5533 : shiftn += ingrid->getStride()[dir_n];
255 5533 : }
256 392 : if( nfound==0 ) {
257 0 : std::string num; Tools::convert( getStep(), num );
258 0 : error("On step " + num + " failed to find required grid point");
259 : }
260 784 : myvals.setValue( 1, minp );
261 392 : }
262 :
263 : }
264 2523 : }
|
