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 3 : static void registerKeywords( Keywords& keys ) {
121 3 : ContourFindingObject<gridtools::EvaluateGridFunction>::registerKeywords( keys );
122 3 : keys.add("compulsory","SEARCHDIR","In which directions do you wish to search for the contour.");
123 3 : }
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 : private:
159 : bool firststep;
160 : /// The parallel task manager
161 : PTM taskmanager;
162 : std::vector<std::string> gnames;
163 : gridtools::GridCoordinatesObject gridcoords;
164 : /// Get the input grid object for the grid that we are finding the contour in
165 : const gridtools::GridCoordinatesObject& getInputGridObject() const ;
166 : public:
167 : static void registerKeywords( Keywords& keys );
168 : explicit FindContourSurface(const ActionOptions&ao);
169 : unsigned getNumberOfDerivatives() override ;
170 : std::vector<std::string> getGridCoordinateNames() const override ;
171 : const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
172 : void calculate() override ;
173 : void getInputData( std::vector<double>& inputdata ) const override;
174 : static void performTask( std::size_t task_index,
175 : const FindContourSurfaceObject& actiondata,
176 : ParallelActionsInput& input,
177 : ParallelActionsOutput& output );
178 : };
179 :
180 : PLUMED_REGISTER_ACTION(FindContourSurface,"FIND_CONTOUR_SURFACE")
181 :
182 3 : void FindContourSurface::registerKeywords( Keywords& keys ) {
183 3 : ActionWithGrid::registerKeywords( keys );
184 6 : keys.addInputKeyword("compulsory","ARG","grid","the labels of the grid in which the contour will be found");
185 3 : FindContourSurfaceObject::registerKeywords( keys );
186 6 : keys.setValueDescription("grid","a grid containing the location of the points in the Willard-Chandler surface along the chosen direction");
187 3 : keys.addDOI("10.1088/1361-648X/aa893d");
188 3 : PTM::registerKeywords( keys );
189 3 : }
190 :
191 1 : FindContourSurface::FindContourSurface(const ActionOptions&ao):
192 : Action(ao),
193 : ActionWithGrid(ao),
194 1 : firststep(true),
195 1 : taskmanager(this) {
196 1 : if( getPntrToArgument(0)->getRank()<2 ) {
197 0 : error("cannot find dividing surface if input grid is one dimensional");
198 : }
199 :
200 : Value* gval=getPntrToArgument(0);
201 1 : gnames.resize( getPntrToArgument(0)->getRank()-1 );
202 :
203 1 : gridtools::ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( gval->getPntrToAction() );
204 1 : if( !ag ) {
205 0 : error("input argument must be a grid");
206 : }
207 1 : const gridtools::GridCoordinatesObject& mygridobj = ag->getGridCoordinatesObject();
208 2 : if( mygridobj.getGridType()=="fibonacci") {
209 0 : error("cannot search for contours in fibonacci grids");
210 : }
211 : // Now add a value
212 1 : std::vector<std::size_t> shape( mygridobj.getDimension()-1 );
213 1 : addValueWithDerivatives( shape );
214 1 : setNotPeriodic();
215 : // Read in the action input
216 : function::FunctionOptions options;
217 1 : FindContourSurfaceObject::read( taskmanager.getActionInput(), this, options );
218 : // Prepare the grid stuff for this action
219 1 : std::vector<bool> ipbc( mygridobj.getDimension()-1 );
220 3 : for(unsigned i=0; i<gnames.size(); ++i) {
221 2 : ipbc[i] = mygridobj.isPeriodic(taskmanager.getActionInput().gdirs[i]);
222 4 : gnames[i] = ag->getGridCoordinateNames()[taskmanager.getActionInput().gdirs[i]];
223 : }
224 1 : gridcoords.setup( "flat", ipbc, 0, 0.0 );
225 1 : log.printf(" calculating dividing surface along which function equals %f \n", taskmanager.getActionInput().cf.contour);
226 1 : taskmanager.setupParallelTaskManager( 0, 0 );
227 1 : }
228 :
229 8 : const gridtools::GridCoordinatesObject& FindContourSurface::getInputGridObject() const {
230 8 : return taskmanager.getActionInput().cf.function.getGridObject();
231 : }
232 :
233 3 : void FindContourSurface::calculate() {
234 3 : if( firststep ) {
235 : std::vector<double> fspacing;
236 1 : std::vector<std::size_t> snbins( gridcoords.getDimension() );
237 1 : std::vector<std::string> smin( gridcoords.getDimension() ), smax( gridcoords.getDimension() );
238 3 : for(unsigned i=0; i<taskmanager.getActionInput().gdirs.size(); ++i) {
239 4 : smin[i]=getInputGridObject().getMin()[taskmanager.getActionInput().gdirs[i]];
240 4 : smax[i]=getInputGridObject().getMax()[taskmanager.getActionInput().gdirs[i]];
241 2 : snbins[i]=getInputGridObject().getNbin(false)[taskmanager.getActionInput().gdirs[i]];
242 : }
243 1 : gridcoords.setBounds( smin, smax, snbins, fspacing );
244 2 : getPntrToComponent(0)->setShape( gridcoords.getNbin(true) );
245 :
246 1 : std::vector<unsigned> find( gridcoords.getDimension() );
247 1 : std::vector<unsigned> ind( gridcoords.getDimension() );
248 197 : for(unsigned i=0; i<gridcoords.getNumberOfPoints(); ++i) {
249 196 : find.assign( find.size(), 0 );
250 196 : gridcoords.getIndices( i, ind );
251 588 : for(unsigned j=0; j<taskmanager.getActionInput().gdirs.size(); ++j) {
252 392 : find[taskmanager.getActionInput().gdirs[j]]=ind[j];
253 : }
254 : }
255 :
256 : // Set the direction in which to look for the contour
257 1 : taskmanager.getActionInput().direction.resize( getInputGridObject().getDimension(), 0 );
258 1 : taskmanager.getActionInput().direction[taskmanager.getActionInput().dir_n] = 0.999999999*getInputGridObject().getGridSpacing()[taskmanager.getActionInput().dir_n];
259 1 : firststep=false;
260 1 : }
261 3 : taskmanager.runAllTasks();
262 3 : }
263 :
264 7 : unsigned FindContourSurface::getNumberOfDerivatives() {
265 7 : return gridcoords.getDimension();
266 : }
267 :
268 3 : std::vector<std::string> FindContourSurface::getGridCoordinateNames() const {
269 3 : return gnames;
270 : }
271 :
272 3 : const gridtools::GridCoordinatesObject& FindContourSurface::getGridCoordinatesObject() const {
273 3 : return gridcoords;
274 : }
275 :
276 3 : void FindContourSurface::getInputData( std::vector<double>& inputdata ) const {
277 : #ifndef DNDEBUG
278 3 : std::size_t rank = gridcoords.getDimension();
279 3 : std::size_t ndata = gridcoords.getNumberOfPoints();
280 3 : if( inputdata.size()!=ndata*rank ) {
281 1 : inputdata.resize( ndata*rank );
282 : }
283 3 : std::vector<double> point( rank );
284 591 : for(unsigned i=0; i<ndata; ++i) {
285 588 : gridcoords.getGridPointCoordinates( i, point );
286 1764 : for(unsigned j=0; j<rank; ++j) {
287 1176 : inputdata[i*rank + j] = point[j];
288 : }
289 : }
290 : #endif
291 3 : }
292 :
293 588 : void FindContourSurface::performTask( std::size_t task_index,
294 : const FindContourSurfaceObject& actiondata,
295 : ParallelActionsInput& input,
296 : ParallelActionsOutput& output ) {
297 : unsigned nfound=0;
298 : const gridtools::GridCoordinatesObject& gridobj = actiondata.cf.function.getGridObject();
299 588 : std::size_t nbins = gridobj.getNbin(false)[actiondata.dir_n];
300 : std::size_t shiftn=task_index;
301 588 : std::vector<unsigned> ind( gridobj.getDimension() );
302 588 : std::vector<double> point( gridobj.getDimension() );
303 : #ifndef DNDEBUG
304 : std::size_t rank = gridobj.getDimension()-1;
305 : View<const double> oind( input.inputdata+task_index*rank, rank );
306 : #endif
307 16842 : for(unsigned i=0; i<nbins; ++i) {
308 : #ifndef DNDEBUG
309 16842 : std::vector<double> base_ind( gridobj.getDimension() );
310 16842 : gridobj.getGridPointCoordinates( shiftn, base_ind );
311 50526 : for(unsigned j=0; j<actiondata.gdirs.size(); ++j) {
312 : plumed_dbg_assert( base_ind[actiondata.gdirs[j]]==oind[j] );
313 : }
314 : #endif
315 : // Get the index of the current grid point
316 16842 : gridobj.getIndices( shiftn, ind );
317 : // Exit if we are at the edge of the grid
318 16842 : if( !gridobj.isPeriodic(actiondata.dir_n) && (ind[actiondata.dir_n]+1)==nbins ) {
319 0 : shiftn += gridobj.getStride()[actiondata.dir_n];
320 : continue;
321 : }
322 :
323 : // Now get the function value at two points
324 16842 : double val1=actiondata.cf.function.function->get( shiftn ) - actiondata.cf.contour;
325 : double val2;
326 16842 : if( (ind[actiondata.dir_n]+1)==nbins ) {
327 0 : val2 = actiondata.cf.function.function->get( task_index ) - actiondata.cf.contour;
328 : } else {
329 16842 : val2 = actiondata.cf.function.function->get( shiftn + gridobj.getStride()[actiondata.dir_n] ) - actiondata.cf.contour;
330 : }
331 :
332 : // Check if the minimum is bracketed
333 16842 : if( val1*val2<0 ) {
334 588 : gridobj.getGridPointCoordinates( shiftn, point );
335 588 : ContourFindingObject<gridtools::EvaluateGridFunction>::findContour( actiondata.cf, actiondata.direction, point );
336 588 : output.values[0] = point[actiondata.dir_n];
337 : nfound++;
338 : break;
339 : }
340 :
341 : // This moves us on to the next point
342 16254 : shiftn += gridobj.getStride()[actiondata.dir_n];
343 : }
344 : if( nfound==0 ) {
345 0 : plumed_merror("failed to find required grid point");
346 : }
347 588 : }
348 :
349 : }
350 : }
|
