Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 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 calculated 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 : where 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 contours 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 analyze 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 : 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 override {
99 12 : return 2;
100 : }
101 0 : bool checkAllActive() const override {
102 0 : return gbuffer==0;
103 : }
104 : void clearAverage() override;
105 : void prepareForAveraging() override;
106 : void compute( const unsigned& current, MultiValue& myvals ) const override;
107 : void finishAveraging() override;
108 : };
109 :
110 13787 : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
111 :
112 5 : void FindContourSurface::registerKeywords( Keywords& keys ) {
113 5 : ContourFindingBase::registerKeywords( keys );
114 10 : keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
115 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");
116 5 : }
117 :
118 1 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
119 : Action(ao),
120 : ContourFindingBase(ao),
121 1 : firsttime(true),
122 1 : ones(ingrid->getDimension(),1) {
123 1 : if( ingrid->getDimension()<2 ) {
124 0 : error("cannot find dividing surface if input grid is one dimensional");
125 : }
126 :
127 : std::string dir;
128 1 : parse("SEARCHDIR",dir);
129 1 : parse("BUFFER",gbuffer);
130 1 : log.printf(" calculating location of contour on %d dimensional grid \n", ingrid->getDimension()-1 );
131 1 : if( gbuffer>0 ) {
132 1 : log.printf(" after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
133 : }
134 1 : checkRead();
135 :
136 : unsigned n=0;
137 1 : gdirs.resize( ingrid->getDimension()-1 );
138 4 : for(unsigned i=0; i<ingrid->getDimension(); ++i) {
139 3 : if( ingrid->getComponentName(i)==dir ) {
140 1 : dir_n=i;
141 : } else {
142 2 : if( n==gdirs.size() ) {
143 0 : error("could not find " + dir + " direction in input grid");
144 : }
145 2 : gdirs[n]=i;
146 2 : n++;
147 : }
148 : }
149 1 : if( n!=(ingrid->getDimension()-1) ) {
150 0 : error("output of grid is not understood");
151 : }
152 :
153 : // Create the input from the old string
154 2 : std::string vstring = "COMPONENTS=" + getLabel() + " COORDINATES=" + ingrid->getComponentName( gdirs[0] );
155 2 : for(unsigned i=1; i<gdirs.size(); ++i) {
156 2 : vstring += "," + ingrid->getComponentName( gdirs[i] );
157 : }
158 : vstring += " PBC=";
159 1 : if( ingrid->isPeriodic(gdirs[0]) ) {
160 : vstring+="T";
161 : } else {
162 : vstring+="F";
163 : }
164 2 : for(unsigned i=1; i<gdirs.size(); ++i) {
165 1 : if( ingrid->isPeriodic(gdirs[i]) ) {
166 : vstring+=",T";
167 : } else {
168 : vstring+=",F";
169 : }
170 : }
171 2 : auto grid=createGrid( "grid", vstring );
172 1 : grid->setNoDerivatives();
173 1 : setAveragingAction( std::move(grid), true );
174 2 : }
175 :
176 3 : void FindContourSurface::clearAverage() {
177 : // Set the boundaries of the output grid
178 : std::vector<double> fspacing;
179 3 : std::vector<unsigned> snbins( ingrid->getDimension()-1 );
180 3 : std::vector<std::string> smin( ingrid->getDimension()-1 ), smax( ingrid->getDimension()-1 );
181 9 : for(unsigned i=0; i<gdirs.size(); ++i) {
182 12 : smin[i]=ingrid->getMin()[gdirs[i]];
183 12 : smax[i]=ingrid->getMax()[gdirs[i]];
184 6 : snbins[i]=ingrid->getNbin()[gdirs[i]];
185 : }
186 3 : mygrid->setBounds( smin, smax, snbins, fspacing);
187 3 : resizeFunctions();
188 3 : ActionWithAveraging::clearAverage();
189 6 : }
190 :
191 3 : void FindContourSurface::prepareForAveraging() {
192 : // Create a task list if first time
193 3 : if( firsttime ) {
194 1 : std::vector<unsigned> find( ingrid->getDimension() );
195 1 : std::vector<unsigned> ind( mygrid->getDimension() );
196 197 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
197 196 : find.assign( find.size(), 0 );
198 196 : mygrid->getIndices( i, ind );
199 588 : for(unsigned j=0; j<gdirs.size(); ++j) {
200 392 : find[gdirs[j]]=ind[j];
201 : }
202 : // Current will be set equal to the start point for this grid index
203 196 : addTaskToList( ingrid->getIndex(find) );
204 : }
205 : // And prepare the task list
206 1 : deactivateAllTasks();
207 197 : for(unsigned i=0; i<getFullNumberOfTasks(); ++i) {
208 196 : taskFlags[i]=1;
209 : }
210 1 : lockContributors();
211 : // Set the direction in which to look for the contour
212 1 : direction.resize( ingrid->getDimension(), 0 );
213 1 : direction[dir_n] = 0.999999999*ingrid->getGridSpacing()[dir_n];
214 : }
215 3 : }
216 :
217 3 : void FindContourSurface::finishAveraging() {
218 : ContourFindingBase::finishAveraging();
219 : // And update the list of active grid points
220 3 : if( gbuffer>0 ) {
221 3 : std::vector<double> dx( ingrid->getGridSpacing() );
222 3 : std::vector<double> point( ingrid->getDimension() );
223 3 : std::vector<double> lpoint( mygrid->getDimension() );
224 : std::vector<unsigned> neighbours;
225 : unsigned num_neighbours;
226 3 : std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
227 3 : std::vector<bool> active( ingrid->getNumberOfPoints(), false );
228 3 : std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
229 591 : for(unsigned i=0; i<mygrid->getNumberOfPoints(); ++i) {
230 : // Retrieve the coordinates of this grid point
231 588 : mygrid->getGridPointCoordinates( i, lpoint );
232 588 : point[dir_n] = mygrid->getGridElement( i, 0 );
233 : // 0.5*dx added here to prevent problems with flooring of grid points
234 1764 : for(unsigned j=0; j<gdirs.size(); ++j) {
235 1176 : point[gdirs[j]]=lpoint[j] + 0.5*dx[gdirs[j]];
236 : }
237 588 : ingrid->getIndices( point, ugrid_indices );
238 : // Now activate buffer region
239 588 : ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
240 16464 : for(unsigned n=0; n<num_neighbours; ++n) {
241 15876 : active[ neighbours[n] ]=true;
242 : }
243 : }
244 3 : ingrid->activateThesePoints( active );
245 : }
246 3 : firsttime=false;
247 3 : }
248 :
249 588 : void FindContourSurface::compute( const unsigned& current, MultiValue& myvals ) const {
250 : std::vector<unsigned> neighbours;
251 : unsigned num_neighbours;
252 : unsigned nfound=0;
253 : double minp=0;
254 588 : std::vector<unsigned> bins_n( ingrid->getNbin() );
255 588 : unsigned shiftn=current;
256 588 : std::vector<unsigned> ind( ingrid->getDimension() );
257 588 : std::vector<double> point( ingrid->getDimension() );
258 : #ifndef NDEBUG
259 : std::vector<unsigned> oind( mygrid->getDimension() );
260 : mygrid->getIndices( current, oind );
261 : #endif
262 16749 : for(unsigned i=0; i<bins_n[dir_n]; ++i) {
263 : #ifndef NDEBUG
264 : std::vector<unsigned> base_ind( ingrid->getDimension() );
265 : ingrid->getIndices( shiftn, base_ind );
266 : for(unsigned j=0; j<gdirs.size(); ++j) {
267 : plumed_dbg_assert( base_ind[gdirs[j]]==oind[j] );
268 : }
269 : #endif
270 : // Ensure inactive grid points are ignored
271 16749 : if( ingrid->inactive( shiftn ) ) {
272 8892 : shiftn += ingrid->getStride()[dir_n];
273 8892 : continue;
274 : }
275 : // Get the index of the current grid point
276 7857 : ingrid->getIndices( shiftn, ind );
277 : // Exit if we are at the edge of the grid
278 7857 : if( !ingrid->isPeriodic(dir_n) && (ind[dir_n]+1)==bins_n[dir_n] ) {
279 0 : shiftn += ingrid->getStride()[dir_n];
280 0 : continue;
281 : }
282 :
283 : // Ensure points with inactive neighbours are ignored
284 7857 : ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
285 : bool cycle=false;
286 179020 : for(unsigned j=0; j<num_neighbours; ++j) {
287 172753 : if( ingrid->inactive( neighbours[j]) ) {
288 : cycle=true;
289 : break;
290 : }
291 : }
292 7857 : if( cycle ) {
293 1590 : shiftn += ingrid->getStride()[dir_n];
294 1590 : continue;
295 : }
296 :
297 : // Now get the function value at two points
298 6267 : double val1=getFunctionValue( shiftn ) - contour;
299 : double val2;
300 6267 : if( (ind[dir_n]+1)==bins_n[dir_n] ) {
301 0 : val2 = getFunctionValue( current ) - contour;
302 : } else {
303 6267 : val2=getFunctionValue( shiftn + ingrid->getStride()[dir_n] ) - contour;
304 : }
305 :
306 : // Check if the minimum is bracketed
307 6267 : if( val1*val2<0 ) {
308 588 : ingrid->getGridPointCoordinates( shiftn, point );
309 588 : findContour( direction, point );
310 588 : minp=point[dir_n];
311 : nfound++;
312 : break;
313 : }
314 :
315 :
316 : // This moves us on to the next point
317 5679 : shiftn += ingrid->getStride()[dir_n];
318 : }
319 : if( nfound==0 ) {
320 : std::string num;
321 0 : Tools::convert( getStep(), num );
322 0 : error("On step " + num + " failed to find required grid point");
323 : }
324 : myvals.setValue( 1, minp );
325 588 : }
326 :
327 : }
328 : }
|
