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 :
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 : \plumedfile
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 : \endplumedfile
80 :
81 : */
82 : //+ENDPLUMEDOC
83 :
84 : namespace PLMD {
85 : namespace gridtools {
86 :
87 3 : 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 12 : 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 6453 : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
107 :
108 2 : void FindContourSurface::registerKeywords( Keywords& keys ) {
109 2 : ContourFindingBase::registerKeywords( keys );
110 8 : keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
111 10 : 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 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
115 : Action(ao),
116 : ContourFindingBase(ao),
117 : firsttime(true),
118 2 : ones(ingrid->getDimension(),1)
119 : {
120 2 : if( ingrid->getDimension()<2 ) error("cannot find dividing surface if input grid is one dimensional");
121 :
122 3 : std::string dir; parse("SEARCHDIR",dir); parse("BUFFER",gbuffer);
123 2 : log.printf(" calculating location of contour on %d dimensional grid \n", ingrid->getDimension()-1 );
124 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);
125 1 : checkRead();
126 :
127 2 : unsigned n=0; gdirs.resize( ingrid->getDimension()-1 );
128 11 : for(unsigned i=0; i<ingrid->getDimension(); ++i) {
129 6 : if( ingrid->getComponentName(i)==dir ) {
130 1 : dir_n=i;
131 : } else {
132 4 : if( n==gdirs.size() ) error("could not find " + dir + " direction in input grid");
133 2 : gdirs[n]=i; n++;
134 : }
135 : }
136 1 : if( n!=(ingrid->getDimension()-1) ) error("output of grid is not understood");
137 :
138 : // Create the input from the old string
139 5 : std::string vstring = "COMPONENTS=" + getLabel() + " COORDINATES=" + ingrid->getComponentName( gdirs[0] );
140 7 : for(unsigned i=1; i<gdirs.size(); ++i) vstring += "," + ingrid->getComponentName( gdirs[i] );
141 : vstring += " PBC=";
142 2 : if( ingrid->isPeriodic(gdirs[0]) ) vstring+="T";
143 : else vstring+="F";
144 5 : for(unsigned i=1; i<gdirs.size(); ++i) {
145 2 : if( ingrid->isPeriodic(gdirs[i]) ) vstring+=",T"; else vstring+=",F";
146 : }
147 2 : createGrid( "grid", vstring ); mygrid->setNoDerivatives();
148 1 : setAveragingAction( mygrid, true );
149 1 : }
150 :
151 3 : void FindContourSurface::clearAverage() {
152 : // Set the boundaries of the output grid
153 6 : std::vector<double> fspacing; std::vector<unsigned> snbins( ingrid->getDimension()-1 );
154 12 : std::vector<std::string> smin( ingrid->getDimension()-1 ), smax( ingrid->getDimension()-1 );
155 24 : for(unsigned i=0; i<gdirs.size(); ++i) {
156 24 : smin[i]=ingrid->getMin()[gdirs[i]]; smax[i]=ingrid->getMax()[gdirs[i]];
157 18 : snbins[i]=ingrid->getNbin()[gdirs[i]];
158 : }
159 3 : mygrid->setBounds( smin, smax, snbins, fspacing); resizeFunctions();
160 3 : ActionWithAveraging::clearAverage();
161 3 : }
162 :
163 3 : void FindContourSurface::prepareForAveraging() {
164 : // Create a task list if first time
165 3 : if( firsttime ) {
166 2 : std::vector<unsigned> find( ingrid->getDimension() );
167 2 : std::vector<unsigned> ind( mygrid->getDimension() );
168 394 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
169 392 : find.assign( find.size(), 0 ); mygrid->getIndices( i, ind );
170 1960 : for(unsigned j=0; j<gdirs.size(); ++j) find[gdirs[j]]=ind[j];
171 : // Current will be set equal to the start point for this grid index
172 196 : addTaskToList( ingrid->getIndex(find) );
173 : }
174 : // And prepare the task list
175 1 : deactivateAllTasks();
176 393 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) taskFlags[i]=1;
177 1 : lockContributors();
178 : // Set the direction in which to look for the contour
179 2 : direction.resize( ingrid->getDimension(), 0 );
180 3 : direction[dir_n] = 0.999999999*ingrid->getGridSpacing()[dir_n];
181 : }
182 3 : }
183 :
184 3 : void FindContourSurface::finishAveraging() {
185 : ContourFindingBase::finishAveraging();
186 : // And update the list of active grid points
187 3 : if( gbuffer>0 ) {
188 3 : std::vector<double> dx( ingrid->getGridSpacing() );
189 6 : std::vector<double> point( ingrid->getDimension() );
190 6 : std::vector<double> lpoint( mygrid->getDimension() );
191 : std::vector<unsigned> neighbours; unsigned num_neighbours;
192 6 : std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
193 6 : std::vector<bool> active( ingrid->getNumberOfPoints(), false );
194 6 : std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
195 1182 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
196 : // Retrieve the coordinates of this grid point
197 588 : mygrid->getGridPointCoordinates( i, lpoint );
198 1176 : point[dir_n] = mygrid->getGridElement( i, 0 );
199 : // 0.5*dx added here to prevent problems with flooring of grid points
200 8232 : for(unsigned j=0; j<gdirs.size(); ++j) point[gdirs[j]]=lpoint[j] + 0.5*dx[gdirs[j]];
201 588 : ingrid->getIndices( point, ugrid_indices );
202 : // Now activate buffer region
203 588 : ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
204 48216 : for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
205 : }
206 3 : ingrid->activateThesePoints( active );
207 : }
208 3 : firsttime=false;
209 3 : }
210 :
211 588 : void FindContourSurface::compute( const unsigned& current, MultiValue& myvals ) const {
212 : std::vector<unsigned> neighbours; unsigned num_neighbours; unsigned nfound=0; double minv=0, minp=0;
213 588 : std::vector<unsigned> bins_n( ingrid->getNbin() ); unsigned shiftn=current;
214 1764 : std::vector<unsigned> ind( ingrid->getDimension() ); std::vector<double> point( ingrid->getDimension() );
215 : #ifndef DNDEBUG
216 1176 : std::vector<unsigned> oind( mygrid->getDimension() ); mygrid->getIndices( current, oind );
217 : #endif
218 49659 : for(unsigned i=0; i<bins_n[dir_n]; ++i) {
219 : #ifndef DNDEBUG
220 33498 : std::vector<unsigned> base_ind( ingrid->getDimension() ); ingrid->getIndices( shiftn, base_ind );
221 100494 : for(unsigned j=0; j<gdirs.size(); ++j) plumed_dbg_assert( base_ind[gdirs[j]]==oind[j] );
222 : #endif
223 : // Ensure inactive grid points are ignored
224 42390 : if( ingrid->inactive( shiftn ) ) { shiftn += ingrid->getStride()[dir_n]; continue; }
225 : // Get the index of the current grid point
226 7857 : ingrid->getIndices( shiftn, ind );
227 : // Exit if we are at the edge of the grid
228 31428 : if( !ingrid->isPeriodic(dir_n) && (ind[dir_n]+1)==bins_n[dir_n] ) {
229 0 : shiftn += ingrid->getStride()[dir_n]; continue;
230 : }
231 :
232 : // Ensure points with inactive neighbours are ignored
233 7857 : ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
234 : bool cycle=false;
235 350183 : for(unsigned j=0; j<num_neighbours; ++j) {
236 345506 : if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
237 : }
238 9447 : if( cycle ) { shiftn += ingrid->getStride()[dir_n]; continue; }
239 :
240 : // Now get the function value at two points
241 6267 : double val1=getFunctionValue( shiftn ) - contour; double val2;
242 18801 : if( (ind[dir_n]+1)==bins_n[dir_n] ) val2 = getFunctionValue( current ) - contour;
243 18801 : else val2=getFunctionValue( shiftn + ingrid->getStride()[dir_n] ) - contour;
244 :
245 : // Check if the minimum is bracketed
246 6267 : if( val1*val2<0 ) {
247 588 : ingrid->getGridPointCoordinates( shiftn, point ); findContour( direction, point );
248 1176 : minp=point[dir_n]; nfound++; break;
249 : }
250 :
251 :
252 : // This moves us on to the next point
253 11358 : shiftn += ingrid->getStride()[dir_n];
254 : }
255 : if( nfound==0 ) {
256 0 : std::string num; Tools::convert( getStep(), num );
257 0 : error("On step " + num + " failed to find required grid point");
258 : }
259 : myvals.setValue( 1, minp );
260 588 : }
261 :
262 : }
263 4839 : }
|
