Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2020 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 "FindContour.h"
23 : #include "core/ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 :
26 : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR
27 : /*
28 : Find an isocontour in a smooth function.
29 :
30 : As discussed in the documentation for the [gridtools](module_gridtools.md), PLUMED contains a number of tools that allow you to calculate
31 : 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
32 : discussed [here](module_contour.md). If this function has one or two input
33 : arguments it is relatively straightforward to plot the function. If by contrast the data has a three dimensions it can be
34 : difficult to visualize.
35 :
36 : This action provides one tool for visualizing these functions. It can be used to search for a set of points on a contour
37 : where the function takes a particular values. In other words, for the function $f(x,y)$ this action would find a set
38 : of points $\{x_c,y_c\}$ that have:
39 :
40 : $$
41 : f(x_c,y_c) - c = 0
42 : $$
43 :
44 : where $c$ is some constant value that is specified by the user. The points on this contour are detected using a variant
45 : on the [marching squares](https://en.wikipedia.org/wiki/Marching_squares) or [marching cubes](https://en.wikipedia.org/wiki/Marching_cubes) algorithm.
46 : As such, and unlike [FIND_CONTOUR_SURFACE](FIND_CONTOUR_SURFACE.md) or [FIND_SPHERICAL_CONTOUR](FIND_SPHERICAL_CONTOUR.md), the function input to this action is not required to have three arguments.
47 : You can find a contour in any function with 2 or more arguments. Furthermore, the topology of the contour will be determined by the algorithm and does not need to be specified by the user.
48 :
49 : ## Examples
50 :
51 : 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
52 : 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
53 : 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
54 : 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
55 : 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 transformed values, $f(s_i)$, of the
56 : [FCCUBIC](FCCUBIC.md) symmetry functions for each of the atoms using the following expression.
57 :
58 : $$
59 : \rho'(x,y,z) = \frac{ \sum_{i=1}^N f(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) }
60 : $$
61 :
62 : where $(x_i, y_i, z_i)$ is the position of atom $i$ relative to the position of atom 1, $f$ is a switching function, $K$ is a Gaussian kernel function and $\lambda=1.0$.
63 :
64 : The virtual atom that is computed using [CENTER](CENTER.md) is located in the center of the region where the atoms are solid like and the positions of the atoms
65 : $(x_i, y_i, z_i)$ are given relative to this position when computing the phase field in the expression above. This phase field provides a continuous function that
66 : gives a measure of the average degree of solidness at each point in the simulation cell. The Willard-Chandler dividing surface is calculated by finding a a set of points
67 : at which the value of this phase field is equal to 0.5. This set of points is output to file called `mycontour.dat`. A new contour
68 : is found on every single step for each frame that is read in.
69 :
70 : ```plumed
71 : UNITS NATURAL
72 :
73 : # This calculates the value of a set of symmetry functions for the atoms of interest
74 : fcc: FCCUBIC ...
75 : SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
76 : ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619
77 : ...
78 :
79 : # Transform the symmetry functions with a switching function
80 : tfcc: LESS_THAN ARG=fcc SWITCH={SMAP R_0=0.5 A=8 B=8}
81 :
82 : # Now compute the center of the solid like region
83 : center: CENTER ATOMS=1-96000 WEIGHTS=tfcc
84 :
85 : # This determines the positions of the atoms of interest relative to the center of the solid region
86 : dens_dist: DISTANCES ORIGIN=center ATOMS=1-96000 COMPONENTS
87 : # This computes the numerator in the expression above for the phase field
88 : dens_numer: KDE ...
89 : VOLUMES=tfcc ARG=dens_dist.x,dens_dist.y,dens_dist.z
90 : GRID_BIN=80,80,80 BANDWIDTH=1.0,1.0,1.0
91 : ...
92 : # This computes the denominator
93 : dens_denom: KDE ...
94 : ARG=dens_dist.x,dens_dist.y,dens_dist.z
95 : GRID_BIN=80,80,80 BANDWIDTH=1.0,1.0,1.0
96 : ...
97 : # This computes the final phase field
98 : dens: CUSTOM ARG=dens_numer,dens_denom FUNC=x/y PERIODIC=NO
99 :
100 : # Find the isocontour
101 : cont: FIND_CONTOUR ARG=dens CONTOUR=0.5
102 : # Use the special method for outputting the contour to a file
103 : DUMPCONTOUR ARG=cont FILE=surface.xyz
104 : ```
105 :
106 : Notice that with this method you have to use the [DUMPCONTOUR](DUMPCONTOUR.md) action to output any contour you find to a file.
107 :
108 : */
109 : //+ENDPLUMEDOC
110 :
111 : namespace PLMD {
112 : namespace contour {
113 :
114 : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
115 :
116 5 : void FindContour::registerKeywords( Keywords& keys ) {
117 5 : ActionWithVector::registerKeywords( keys );
118 5 : ActionWithValue::useCustomisableComponents(keys);
119 10 : keys.addInputKeyword("compulsory","ARG","grid","the labels of the grid in which the contour will be found");
120 5 : ContourFindingObject<gridtools::EvaluateGridFunction>::registerKeywords( keys );
121 5 : 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");
122 5 : keys.addDOI("10.1039/B805786A");
123 5 : keys.addDOI("10.1021/jp909219k");
124 5 : keys.addDOI("10.1088/1361-648X/aa893d");
125 5 : keys.addDOI("10.1063/1.5134461");
126 5 : PTM::registerKeywords( keys );
127 5 : }
128 :
129 2 : FindContour::FindContour(const ActionOptions&ao):
130 : Action(ao),
131 : ActionWithVector(ao),
132 2 : firststep(true),
133 2 : taskmanager(this) {
134 2 : parse("BUFFER",gbuffer);
135 2 : if( gbuffer>0 ) {
136 0 : log.printf(" after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
137 : }
138 :
139 2 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
140 2 : std::vector<std::string> argn( ag->getGridCoordinateNames() );
141 :
142 2 : std::vector<std::size_t> shape(1);
143 2 : shape[0]=0;
144 8 : for(unsigned i=0; i<argn.size(); ++i ) {
145 6 : addComponent( argn[i], shape );
146 6 : componentIsNotPeriodic( argn[i] );
147 : }
148 : function::FunctionOptions options;
149 2 : ContourFindingObject<gridtools::EvaluateGridFunction>::read( taskmanager.getActionInput(), this, options );
150 2 : log.printf(" calculating dividing surface along which function equals %f \n", taskmanager.getActionInput().contour);
151 2 : }
152 :
153 3 : std::string FindContour::getOutputComponentDescription( const std::string& cname, const Keywords& keys ) const {
154 6 : return "a vector of coordinates for the contour along the " + cname + " direction";
155 : }
156 :
157 2 : void FindContour::calculate() {
158 2 : if( firststep ) {
159 1 : std::vector<std::size_t> shape(1);
160 1 : shape[0] = getPntrToArgument(0)->getRank()*getPntrToArgument(0)->getNumberOfValues();
161 4 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
162 3 : getPntrToComponent(i)->setShape( shape );
163 : }
164 1 : active_cells.resize( shape[0] );
165 1 : taskmanager.setupParallelTaskManager( 0, 0 );
166 1 : firststep = false;
167 : }
168 2 : taskmanager.runAllTasks();
169 2 : }
170 :
171 0 : unsigned FindContour::getNumberOfDerivatives() {
172 0 : return 0;
173 : }
174 :
175 3 : void FindContour::getNumberOfTasks( unsigned& ntasks ) {
176 3 : ntasks = active_cells.size();
177 :
178 : Value* gval=getPntrToArgument(0);
179 3 : unsigned npoints = gval->getNumberOfValues();
180 3 : std::vector<unsigned> ind( gval->getRank() );
181 3 : std::vector<unsigned> ones( gval->getRank(), 1 );
182 3 : std::vector<std::size_t> nbin( getInputGridObject().getNbin( false ) );
183 : unsigned num_neighbours;
184 : std::vector<unsigned> neighbours;
185 :
186 : std::fill( active_cells.begin(), active_cells.end(), 0 );
187 16467 : for(unsigned i=0; i<npoints; ++i) {
188 : // Get the index of the current grid point
189 16464 : getInputGridObject().getIndices( i, ind );
190 16464 : getInputGridObject().getNeighbors( ind, ones, num_neighbours, neighbours );
191 : // Get the value of a point on the grid
192 16464 : double val1=gval->get( i ) - taskmanager.getActionInput().contour;
193 : bool edge=false;
194 65856 : for(unsigned j=0; j<gval->getRank(); ++j) {
195 : // Make sure we don't search at the edge of the grid
196 98784 : if( !getInputGridObject().isPeriodic(j) && (ind[j]+1)==nbin[j] ) {
197 0 : continue;
198 49392 : } else if( (ind[j]+1)==nbin[j] ) {
199 : edge=true;
200 2940 : ind[j]=0;
201 : } else {
202 46452 : ind[j]+=1;
203 : }
204 49392 : double val2=gval->get( getInputGridObject().getIndex(ind) ) - taskmanager.getActionInput().contour;
205 49392 : if( val1*val2<0 ) {
206 1662 : active_cells[gval->getRank()*i + j] = 1;
207 : }
208 98784 : if( getInputGridObject().isPeriodic(j) && edge ) {
209 : edge=false;
210 2940 : ind[j]=nbin[j]-1;
211 : } else {
212 46452 : ind[j]-=1;
213 : }
214 : }
215 : }
216 3 : }
217 :
218 32928 : int FindContour::checkTaskIsActive( const unsigned& taskno ) const {
219 32928 : if( active_cells[taskno]>0 ) {
220 1108 : return 1;
221 : }
222 : return -1;
223 : }
224 :
225 2 : void FindContour::getInputData( std::vector<double>& inputdata ) const {
226 2 : std::size_t rank = getPntrToArgument(0)->getRank();
227 2 : std::size_t ndata = getInputGridObject().getNumberOfPoints()*rank;
228 2 : if( inputdata.size()!=2*rank*ndata ) {
229 1 : inputdata.resize( 2*rank*ndata );
230 : }
231 2 : std::vector<double> point( getPntrToArgument(0)->getRank() );
232 10978 : for(unsigned i=0; i<getInputGridObject().getNumberOfPoints(); ++i) {
233 10976 : getInputGridObject().getGridPointCoordinates( i, point );
234 43904 : for(unsigned j=0; j<rank; ++j) {
235 131712 : for(unsigned k=0; k<rank; ++k) {
236 98784 : inputdata[i*rank*2*rank + j*2*rank + k] = point[k];
237 98784 : inputdata[i*rank*2*rank + j*2*rank + rank + k] = 0;
238 : }
239 32928 : inputdata[i*rank*2*rank + j*2*rank + rank + j] = 0.999999999*getInputGridObject().getGridSpacing()[j];
240 : }
241 : }
242 2 : }
243 :
244 1108 : void FindContour::performTask( std::size_t task_index,
245 : const ContourFindingObject<gridtools::EvaluateGridFunction>& actiondata,
246 : ParallelActionsInput& input,
247 : ParallelActionsOutput& output ) {
248 :
249 1108 : std::size_t rank = actiondata.function.getGridObject().getDimension();
250 1108 : std::vector<double> direction( rank ), point( rank );
251 4432 : for(unsigned i=0; i<rank; ++i) {
252 3324 : point[i] = input.inputdata[ 2*rank*task_index + i];
253 3324 : direction[i] = input.inputdata[ 2*rank*task_index + rank + i];
254 : }
255 : // Now find the contour
256 1108 : ContourFindingObject<gridtools::EvaluateGridFunction>::findContour( actiondata, direction, point );
257 4432 : for(unsigned i=0; i<rank; ++i) {
258 3324 : output.values[i] = point[i];
259 : }
260 1108 : }
261 :
262 : }
263 : }
|
