Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "vesselbase/StoreDataVessel.h"
24 : #include "ContourFindingBase.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/Atoms.h"
27 :
28 : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR
29 : /*
30 : Find an isocontour in a smooth function.
31 :
32 : As discussed in the part of the manual on \ref Analysis 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 \ref HISTOGRAM as a function of a few collective variables
34 : or it might be a phase field that has been calcualted using \ref MULTICOLVARDENS. 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 or more dimensions
36 : it can be 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 values. In other words, for the function \f$f(x,y)\f$ this action would find a set
40 : of points \f$\{x_c,y_c\}\f$ that have:
41 :
42 : \f[
43 : f(x_c,y_c) - c = 0
44 : \f]
45 :
46 : where \f$c\f$ is some constant value that is specified by the user. The points on this contour are detected using a variant
47 : on the marching squares or marching cubes algorithm, which you can find information on here:
48 :
49 : https://en.wikipedia.org/wiki/Marching_squares
50 : https://en.wikipedia.org/wiki/Marching_cubes
51 :
52 : As such, and unlike \ref FIND_CONTOUR_SURFACE or \ref FIND_SPHERICAL_CONTOUR, the function input to this action can have any dimension.
53 : Furthermore, the topology of the contour will be determined by the algorithm and does not need to be specified by the user.
54 :
55 : \par Examples
56 :
57 : The input below allows you to calculate something akin to a Willard-Chandler dividing surface \cite wcsurface.
58 : The simulation cell in this case contains a solid phase and a liquid phase. The Willard-Chandler surface is the
59 : surface that separates the parts of the box containing the solid from the parts containing the liquid. To compute the position
60 : of this surface the \ref FCCUBIC symmetry function is calculated for each of the atoms in the system from on the geometry of the
61 : atoms in the first coordination sphere of each of the atoms. These quantities are then transformed using a switching function.
62 : This procedure generates a single number for each atom in the system and this quantity has a value of one for atoms that are in
63 : parts of the box that resemble the solid structure and zero for atoms that are in parts of the box that resemble the liquid.
64 : The position of a virtual atom is then computed using \ref CENTER_OF_MULTICOLVAR and a phase field model is constructed using
65 : \ref MULTICOLVARDENS. These procedure ensures that we have a continuous function that gives a measure of the average degree of
66 : 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 : \plumedfile
71 : UNITS NATURAL
72 : FCCUBIC ...
73 : SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
74 : ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619 LABEL=fcc
75 : ... FCCUBIC
76 :
77 : tfcc: MTRANSFORM_MORE DATA=fcc SWITCH={SMAP R_0=0.5 A=8 B=8}
78 : center: CENTER_OF_MULTICOLVAR DATA=tfcc
79 :
80 : MULTICOLVARDENS ...
81 : DATA=tfcc ORIGIN=center DIR=xyz LABEL=dens
82 : NBINS=80,80,80 BANDWIDTH=1.0,1.0,1.0 STRIDE=25
83 : LABEL=dens STRIDE=1 CLEAR=1
84 : ... MULTICOLVARDENS
85 :
86 : FIND_CONTOUR GRID=dens CONTOUR=0.5 FILE=mycontour.dat
87 : \endplumedfile
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : namespace PLMD {
93 : namespace gridtools {
94 :
95 3 : class FindContour : public ContourFindingBase {
96 : private:
97 : bool firsttime;
98 : unsigned gbuffer;
99 : /// Stuff for output
100 : OFile of;
101 : double lenunit;
102 : std::string fmt_xyz;
103 : /// The data is stored in a grid
104 : vesselbase::StoreDataVessel* mydata;
105 : public:
106 : static void registerKeywords( Keywords& keys );
107 : explicit FindContour(const ActionOptions&ao);
108 0 : bool checkAllActive() const { return gbuffer==0; }
109 : void prepareForAveraging();
110 0 : bool isPeriodic() { return false; }
111 14 : unsigned getNumberOfQuantities() const { return 1 + ingrid->getDimension(); }
112 : void compute( const unsigned& current, MultiValue& myvals ) const ;
113 : void finishAveraging();
114 : };
115 :
116 6453 : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
117 :
118 2 : void FindContour::registerKeywords( Keywords& keys ) {
119 2 : ContourFindingBase::registerKeywords( keys );
120 : // We want a better way of doing this bit
121 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");
122 8 : keys.add("compulsory","FILE","file on which to output coordinates");
123 10 : keys.add("compulsory","UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
124 8 : keys.add("optional", "PRECISION","The number of digits in trajectory file");
125 2 : }
126 :
127 1 : FindContour::FindContour(const ActionOptions&ao):
128 : Action(ao),
129 : ContourFindingBase(ao),
130 1 : firsttime(true)
131 : {
132 :
133 2 : parse("BUFFER",gbuffer);
134 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);
135 :
136 2 : std::string file; parse("FILE",file);
137 1 : if( file.length()==0 ) error("name out output file was not specified");
138 1 : std::string type=Tools::extension(file);
139 1 : log<<" file name "<<file<<"\n";
140 1 : if(type!="xyz") error("can only print xyz file type with contour finding");
141 :
142 : fmt_xyz="%f";
143 2 : std::string precision; parse("PRECISION",precision);
144 1 : if(precision.length()>0) {
145 1 : int p; Tools::convert(precision,p);
146 1 : log<<" with precision "<<p<<"\n";
147 : std::string a,b;
148 1 : Tools::convert(p+5,a);
149 1 : Tools::convert(p,b);
150 5 : fmt_xyz="%"+a+"."+b+"f";
151 : }
152 2 : std::string unitname; parse("UNITS",unitname);
153 1 : if(unitname!="PLUMED") {
154 0 : Units myunit; myunit.setLength(unitname);
155 0 : lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
156 : }
157 1 : else lenunit=1.0;
158 1 : of.link(*this); of.open(file);
159 1 : checkRead(); mydata=buildDataStashes( NULL );
160 1 : }
161 :
162 2 : void FindContour::prepareForAveraging() {
163 : // Create a task list if first time
164 2 : if( firsttime ) {
165 16466 : for(unsigned i=0; i<ingrid->getDimension()*ingrid->getNumberOfPoints(); ++i) addTaskToList( i );
166 : }
167 2 : firsttime=false; deactivateAllTasks();
168 :
169 : // We now need to identify the grid points that we need to search through
170 2 : std::vector<unsigned> nbin( ingrid->getNbin() );
171 4 : std::vector<unsigned> ind( ingrid->getDimension() );
172 4 : std::vector<unsigned> ones( ingrid->getDimension(), 1 );
173 : unsigned num_neighbours; std::vector<unsigned> neighbours;
174 21956 : for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
175 : // Ensure inactive grid points are ignored
176 10976 : if( ingrid->inactive(i) ) continue;
177 :
178 : // Get the index of the current grid point
179 10976 : ingrid->getIndices( i, ind );
180 10976 : ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
181 : bool cycle=false;
182 603680 : for(unsigned j=0; j<num_neighbours; ++j) {
183 592704 : if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
184 : }
185 10976 : if( cycle ) continue;
186 :
187 : // Get the value of a point on the grid
188 10976 : double val1=getFunctionValue( i ) - contour;
189 : bool edge=false;
190 120736 : for(unsigned j=0; j<ingrid->getDimension(); ++j) {
191 : // Make sure we don't search at the edge of the grid
192 32928 : if( !ingrid->isPeriodic(j) && (ind[j]+1)==nbin[j] ) continue;
193 65856 : else if( (ind[j]+1)==nbin[j] ) { edge=true; ind[j]=0; }
194 30968 : else ind[j]+=1;
195 32928 : double val2=getFunctionValue( ind ) - contour;
196 35144 : if( val1*val2<0 ) taskFlags[ ingrid->getDimension()*i + j ] = 1;
197 69776 : if( ingrid->isPeriodic(j) && edge ) { edge=false; ind[j]=nbin[j]-1; }
198 30968 : else ind[j]-=1;
199 : }
200 : }
201 2 : lockContributors();
202 2 : }
203 :
204 554 : void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
205 : // Retrieve the initial grid point coordinates
206 1108 : unsigned gpoint = std::floor( current / ingrid->getDimension() );
207 554 : std::vector<double> point( ingrid->getDimension() );
208 554 : ingrid->getGridPointCoordinates( gpoint, point );
209 :
210 : // Retrieve the direction we are searching for the contour
211 1108 : unsigned gdir = current%(ingrid->getDimension() );
212 554 : std::vector<double> direction( ingrid->getDimension(), 0 );
213 1662 : direction[gdir] = 0.999999999*ingrid->getGridSpacing()[gdir];
214 :
215 : // Now find the contour
216 : findContour( direction, point );
217 : // And transfer to the store data vessel
218 7756 : for(unsigned i=0; i<ingrid->getDimension(); ++i) myvals.setValue( 1+i, point[i] );
219 554 : }
220 :
221 1 : void FindContour::finishAveraging() {
222 : // And update the list of active grid points
223 1 : if( gbuffer>0 ) {
224 : std::vector<unsigned> neighbours; unsigned num_neighbours;
225 0 : std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
226 0 : std::vector<bool> active( ingrid->getNumberOfPoints(), false );
227 0 : std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
228 0 : for(unsigned i=0; i<getCurrentNumberOfActiveTasks(); ++i) {
229 : // Get the point we are operating on
230 0 : unsigned ipoint = std::floor( getActiveTask(i) / ingrid->getDimension() );
231 : // Get the indices of this point
232 0 : ingrid->getIndices( ipoint, ugrid_indices );
233 : // Now activate buffer region
234 0 : ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
235 0 : for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
236 : }
237 0 : ingrid->activateThesePoints( active );
238 : }
239 2 : std::vector<double> point( 1 + ingrid->getDimension() );
240 1 : of.printf("%u\n",mydata->getNumberOfStoredValues());
241 1 : of.printf("Points found on isocontour\n");
242 555 : for(unsigned i=0; i<mydata->getNumberOfStoredValues(); ++i) {
243 554 : mydata->retrieveSequentialValue( i, false, point ); of.printf("X");
244 9418 : for(unsigned j=0; j<ingrid->getDimension(); ++j) of.printf( (" " + fmt_xyz).c_str(), lenunit*point[1+j] );
245 554 : of.printf("\n");
246 : }
247 1 : }
248 :
249 : }
250 4839 : }
|
