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 "ContourFindingObject.h"
24 : #include "gridtools/ActionWithGrid.h"
25 : #include "gridtools/EvaluateGridFunction.h"
26 : #include "core/ParallelTaskManager.h"
27 :
28 : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR_SURFACE
29 : /*
30 : Find an isocontour by searching along either the x, y or z direction.
31 :
32 : As discussed in the documentation for the [gridtools](module_gridtools.md), PLUMED contains a number of tools that allow you to calculate
33 : a function on a grid. The function on this grid might be a [HISTOGRAM](HISTOGRAM.md) or it might be one of the phase fields that are
34 : discussed [here](module_contour.md). If this function has one or two input
35 : arguments it is relatively straightforward to plot the function. If by contrast the data has a three dimensions it can be
36 : difficult to visualize.
37 :
38 : This action provides one tool for visualizing these functions. It can be used to search for a set of points on a contour
39 : where the function takes a particular value. In other words, for the function $f(x,y,z)$ this action would find a set
40 : of points $\{x_c,y_c,z_c\}$ that have:
41 :
42 : $$
43 : f(x_c,y_c,z_c) - c = 0
44 : $$
45 :
46 : where $c$ is some constant value that is specified by the user. The points on this contour are find by searching along lines
47 : that run parallel to the $x$, $y$ or $z$ axis of the simulation cell. The result is, therefore, a two dimensional
48 : function evaluated on a grid that gives us the height of the interface as a function of two coordinates.
49 :
50 : 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
51 : 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
52 : function have the appropriate topology you should use [FIND_CONTOUR](FIND_CONTOUR.md) in place of this action.
53 :
54 : ## Examples
55 :
56 : 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
57 : the solid and the liquid was set up in the plane perpendicular to the $z$ direction of the simulation cell. The input below calculates something
58 : akin to a Willard-Chandler dividing surface (see [contour](module_contour.md)) between the solid phase and the liquid phase. There are two of these interfaces within the
59 : 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
60 : simulation box. The input below detects the height profile of one of these two interfaces. It does so by computing a phase field average from the values, $s_i$, of the
61 : [FCCUBIC](FCCUBIC.md) symmetry functions for each of the atoms using the following expression.
62 :
63 : $$
64 : \rho'(x,y,z) = \frac{ \sum_{i=1}^N s_i K\left(\frac{x-x_i}{\lambda}, \frac{y-y_i}{\lambda}, \frac{z-z_i}{\lambda}\right) }{ \sum_{i=1}^N K\left(\frac{x-x_i}{\lambda}, \frac{y-y_i}{\lambda}, \frac{z-z_i}{\lambda}\right) }
65 : $$
66 :
67 : where $(x_i, y_i, z_i)$ is the position of atom $i$ relative to the position of atom 1, $K$ is a Gaussian kernel function and $\lambda=1.0$.
68 :
69 : Notice that we use the fact that we know roughly where the interface is when specifying how this phase field is to be calculated and specify the region over the $z$-axis
70 : in which the [KDE](KDE.md) is computed. Once we have calculated the phase field we search for contour points on the lines that run parallel to the $z$-direction of the cell
71 : box using the FIND_CONTOUR_SURFACE command. The final result is a $14 \times 14$ grid of values for the height of the interface as a function of the $(x,y)$
72 : position. This grid is then output to a file called `contour2.dat`.
73 :
74 : 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
75 : on every step.
76 :
77 : ```plumed
78 : UNITS NATURAL
79 :
80 : # This calculates the value of a set of symmetry functions for the atoms of interest
81 : fcc: FCCUBIC ...
82 : SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
83 : ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619
84 : ...
85 :
86 : # This determines the positions of the atoms of interest relative to the position of atom 1
87 : dens2_dist: DISTANCES ORIGIN=1 ATOMS=fcc COMPONENTS
88 : # This computes the numerator in the expression above for the phase field
89 : dens2_numer: KDE ...
90 : VOLUMES=fcc_n ARG=dens2_dist.x,dens2_dist.y,dens2_dist.z
91 : GRID_BIN=14,14,50 GRID_MIN=auto,auto,6.0
92 : GRID_MAX=auto,auto,11.0 BANDWIDTH=1.0,1.0,1.0
93 : ...
94 : # This computes the denominator
95 : dens2_denom: KDE ...
96 : ARG=dens2_dist.x,dens2_dist.y,dens2_dist.z
97 : GRID_BIN=14,14,50 GRID_MIN=auto,auto,6.0
98 : GRID_MAX=auto,auto,11.0 BANDWIDTH=1.0,1.0,1.0
99 : ...
100 : # This computes the final phase field
101 : dens2: CUSTOM ARG=dens2_numer,dens2_denom FUNC=x/y PERIODIC=NO
102 :
103 : # We can now find and print the location of the two dimensional contour surface
104 : ss2: FIND_CONTOUR_SURFACE ARG=dens2 CONTOUR=0.42 SEARCHDIR=dens2_dist.z
105 : DUMPGRID ARG=ss2 FILE=contour2.dat STRIDE=1
106 : ```
107 :
108 : */
109 : //+ENDPLUMEDOC
110 :
111 : namespace PLMD {
112 : namespace contour {
113 :
114 1 : class FindContourSurfaceObject {
115 : public:
116 : unsigned dir_n;
117 : std::vector<unsigned> gdirs;
118 : std::vector<double> direction;
119 : ContourFindingObject<gridtools::EvaluateGridFunction> cf;
120 5 : static void registerKeywords( Keywords& keys ) {
121 5 : ContourFindingObject<gridtools::EvaluateGridFunction>::registerKeywords( keys );
122 5 : keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
123 5 : }
124 1 : static void read( FindContourSurfaceObject& func, ActionWithArguments* action, function::FunctionOptions& options ) {
125 1 : ContourFindingObject<gridtools::EvaluateGridFunction>::read( func.cf, action, options );
126 : std::string dir;
127 1 : action->parse("SEARCHDIR",dir);
128 1 : action->log.printf(" calculating location of contour on %d dimensional grid \n", (action->getPntrToArgument(0))->getRank()-1 );
129 1 : Value* gval=func.cf.function.function;
130 1 : gridtools::ActionWithGrid* ag = gridtools::ActionWithGrid::getInputActionWithGrid( gval->getPntrToAction() );
131 1 : if( !ag ) {
132 0 : action->error("input argument must be a grid");
133 : }
134 1 : std::vector<std::string> argn( ag->getGridCoordinateNames() );
135 1 : func.gdirs.resize( gval->getRank()-1 );
136 : unsigned n=0;
137 4 : for(unsigned i=0; i<gval->getRank(); ++i) {
138 3 : if( argn[i]==dir ) {
139 1 : func.dir_n=i;
140 : } else {
141 2 : if( n==func.gdirs.size() ) {
142 0 : action->error("could not find " + dir + " direction in input grid");
143 : }
144 2 : func.gdirs[n]=i;
145 2 : n++;
146 : }
147 : }
148 1 : if( n!=(gval->getRank()-1) ) {
149 0 : action->error("output of grid is not understood");
150 : }
151 2 : }
152 : };
153 :
154 : class FindContourSurface : public gridtools::ActionWithGrid {
155 : public:
156 : using input_type = FindContourSurfaceObject;
157 : using PTM = ParallelTaskManager<FindContourSurface>;
158 : typedef typename PTM::ParallelActionsInput ParallelActionsInput;
159 : typedef typename PTM::ParallelActionsOutput ParallelActionsOutput;
160 : private:
161 : bool firststep;
162 : /// The parallel task manager
163 : PTM taskmanager;
164 : std::vector<std::string> gnames;
165 : gridtools::GridCoordinatesObject gridcoords;
166 : /// Get the input grid object for the grid that we are finding the contour in
167 : const gridtools::GridCoordinatesObject& getInputGridObject() const ;
168 : public:
169 : static void registerKeywords( Keywords& keys );
170 : explicit FindContourSurface(const ActionOptions&ao);
171 : unsigned getNumberOfDerivatives() override ;
172 : std::vector<std::string> getGridCoordinateNames() const override ;
173 : const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
174 : void calculate() override ;
175 : void getArgumentInputData( std::vector<double>& inputdata ) const override;
176 : static void performTask( std::size_t task_index,
177 : const FindContourSurfaceObject& actiondata,
178 : ParallelActionsInput& input,
179 : ParallelActionsOutput& output );
180 : };
181 :
182 : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
183 :
184 5 : void FindContourSurface::registerKeywords( Keywords& keys ) {
185 5 : ActionWithGrid::registerKeywords( keys );
186 10 : keys.addInputKeyword("compulsory","ARG","grid","the labels of the grid in which the contour will be found");
187 5 : FindContourSurfaceObject::registerKeywords( keys );
188 10 : keys.setValueDescription("grid","a grid containing the location of the points in the Willard-Chandler surface along the chosen direction");
189 5 : keys.addDOI("10.1088/1361-648X/aa893d");
190 5 : PTM::registerKeywords( keys );
191 5 : }
192 :
193 1 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
194 : Action(ao),
195 : ActionWithGrid(ao),
196 1 : firststep(true),
197 1 : taskmanager(this) {
198 1 : if( getPntrToArgument(0)->getRank()<2 ) {
199 0 : error("cannot find dividing surface if input grid is one dimensional");
200 : }
201 :
202 : Value* gval=getPntrToArgument(0);
203 1 : gnames.resize( getPntrToArgument(0)->getRank()-1 );
204 :
205 1 : gridtools::ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( gval->getPntrToAction() );
206 1 : if( !ag ) {
207 0 : error("input argument must be a grid");
208 : }
209 1 : const gridtools::GridCoordinatesObject& mygridobj = ag->getGridCoordinatesObject();
210 1 : if( mygridobj.getGridType()=="fibonacci") {
211 0 : error("cannot search for contours in fibonacci grids");
212 : }
213 : // Now add a value
214 1 : std::vector<std::size_t> shape( mygridobj.getDimension()-1 );
215 1 : addValueWithDerivatives( shape );
216 1 : setNotPeriodic();
217 : // Read in the action input
218 1 : function::FunctionOptions options;
219 1 : FindContourSurfaceObject::read( taskmanager.getActionInput(), this, options );
220 : // Prepare the grid stuff for this action
221 1 : std::vector<bool> ipbc( mygridobj.getDimension()-1 );
222 3 : for(unsigned i=0; i<gnames.size(); ++i) {
223 2 : ipbc[i] = mygridobj.isPeriodic(taskmanager.getActionInput().gdirs[i]);
224 4 : gnames[i] = ag->getGridCoordinateNames()[taskmanager.getActionInput().gdirs[i]];
225 : }
226 1 : gridcoords.setup( "flat", ipbc, 0, 0.0 );
227 1 : log.printf(" calculating dividing surface along which function equals %f \n", taskmanager.getActionInput().cf.contour);
228 1 : taskmanager.setupParallelTaskManager( 0, 0 );
229 1 : }
230 :
231 8 : const gridtools::GridCoordinatesObject& FindContourSurface::getInputGridObject() const {
232 8 : return taskmanager.getActionInput().cf.function.getGridObject();
233 : }
234 :
235 3 : void FindContourSurface::calculate() {
236 3 : if( firststep ) {
237 1 : std::vector<double> fspacing;
238 1 : std::vector<std::size_t> snbins( gridcoords.getDimension() );
239 1 : std::vector<std::string> smin( gridcoords.getDimension() ), smax( gridcoords.getDimension() );
240 3 : for(unsigned i=0; i<taskmanager.getActionInput().gdirs.size(); ++i) {
241 4 : smin[i]=getInputGridObject().getMin()[taskmanager.getActionInput().gdirs[i]];
242 4 : smax[i]=getInputGridObject().getMax()[taskmanager.getActionInput().gdirs[i]];
243 2 : snbins[i]=getInputGridObject().getNbin(false)[taskmanager.getActionInput().gdirs[i]];
244 : }
245 1 : gridcoords.setBounds( smin, smax, snbins, fspacing );
246 1 : getPntrToComponent(0)->setShape( gridcoords.getNbin(true) );
247 :
248 1 : std::vector<unsigned> find( gridcoords.getDimension() );
249 1 : std::vector<unsigned> ind( gridcoords.getDimension() );
250 197 : for(unsigned i=0; i<gridcoords.getNumberOfPoints(); ++i) {
251 196 : find.assign( find.size(), 0 );
252 196 : gridcoords.getIndices( i, ind );
253 588 : for(unsigned j=0; j<taskmanager.getActionInput().gdirs.size(); ++j) {
254 392 : find[taskmanager.getActionInput().gdirs[j]]=ind[j];
255 : }
256 : }
257 :
258 : // Set the direction in which to look for the contour
259 1 : taskmanager.getActionInput().direction.resize( getInputGridObject().getDimension(), 0 );
260 1 : taskmanager.getActionInput().direction[taskmanager.getActionInput().dir_n] = 0.999999999*getInputGridObject().getGridSpacing()[taskmanager.getActionInput().dir_n];
261 1 : firststep=false;
262 1 : }
263 3 : taskmanager.runAllTasks();
264 3 : }
265 :
266 7 : unsigned FindContourSurface::getNumberOfDerivatives() {
267 7 : return gridcoords.getDimension();
268 : }
269 :
270 3 : std::vector<std::string> FindContourSurface::getGridCoordinateNames() const {
271 3 : return gnames;
272 : }
273 :
274 3 : const gridtools::GridCoordinatesObject& FindContourSurface::getGridCoordinatesObject() const {
275 3 : return gridcoords;
276 : }
277 :
278 3 : void FindContourSurface::getArgumentInputData( std::vector<double>& inputdata ) const {
279 : #ifndef DNDEBUG
280 3 : std::size_t rank = gridcoords.getDimension();
281 3 : std::size_t ndata = gridcoords.getNumberOfPoints();
282 3 : if( inputdata.size()!=ndata*rank ) {
283 1 : inputdata.resize( ndata*rank );
284 : }
285 3 : std::vector<double> point( rank );
286 591 : for(unsigned i=0; i<ndata; ++i) {
287 588 : gridcoords.getGridPointCoordinates( i, point );
288 1764 : for(unsigned j=0; j<rank; ++j) {
289 1176 : inputdata[i*rank + j] = point[j];
290 : }
291 : }
292 : #endif
293 3 : }
294 :
295 588 : void FindContourSurface::performTask( std::size_t task_index,
296 : const FindContourSurfaceObject& actiondata,
297 : ParallelActionsInput& input,
298 : ParallelActionsOutput& output ) {
299 : unsigned nfound=0;
300 : const gridtools::GridCoordinatesObject& gridobj = actiondata.cf.function.getGridObject();
301 588 : std::size_t nbins = gridobj.getNbin(false)[actiondata.dir_n];
302 : std::size_t shiftn=task_index;
303 588 : std::vector<unsigned> ind( gridobj.getDimension() );
304 588 : std::vector<double> point( gridobj.getDimension() );
305 : #ifndef DNDEBUG
306 : std::size_t rank = gridobj.getDimension()-1;
307 : View<const double> oind( input.inputdata+task_index*rank, rank );
308 : #endif
309 16842 : for(unsigned i=0; i<nbins; ++i) {
310 : #ifndef DNDEBUG
311 16842 : std::vector<double> base_ind( gridobj.getDimension() );
312 16842 : gridobj.getGridPointCoordinates( shiftn, base_ind );
313 50526 : for(unsigned j=0; j<actiondata.gdirs.size(); ++j) {
314 : plumed_dbg_assert( base_ind[actiondata.gdirs[j]]==oind[j] );
315 : }
316 : #endif
317 : // Get the index of the current grid point
318 16842 : gridobj.getIndices( shiftn, ind );
319 : // Exit if we are at the edge of the grid
320 16842 : if( !gridobj.isPeriodic(actiondata.dir_n) && (ind[actiondata.dir_n]+1)==nbins ) {
321 0 : shiftn += gridobj.getStride()[actiondata.dir_n];
322 : continue;
323 : }
324 :
325 : // Now get the function value at two points
326 16842 : double val1=actiondata.cf.function.function->get( shiftn ) - actiondata.cf.contour;
327 : double val2;
328 16842 : if( (ind[actiondata.dir_n]+1)==nbins ) {
329 0 : val2 = actiondata.cf.function.function->get( task_index ) - actiondata.cf.contour;
330 : } else {
331 16842 : val2 = actiondata.cf.function.function->get( shiftn + gridobj.getStride()[actiondata.dir_n] ) - actiondata.cf.contour;
332 : }
333 :
334 : // Check if the minimum is bracketed
335 16842 : if( val1*val2<0 ) {
336 588 : gridobj.getGridPointCoordinates( shiftn, point );
337 588 : ContourFindingObject<gridtools::EvaluateGridFunction>::findContour( actiondata.cf, actiondata.direction, point );
338 588 : output.values[0] = point[actiondata.dir_n];
339 : nfound++;
340 : break;
341 : }
342 :
343 : // This moves us on to the next point
344 16254 : shiftn += gridobj.getStride()[actiondata.dir_n];
345 : }
346 : if( nfound==0 ) {
347 0 : plumed_merror("failed to find required grid point");
348 : }
349 588 : }
350 :
351 : }
352 : }
|
