Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "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 calculated 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 LOWMEM SWITCH={SMAP R_0=0.5 A=8 B=8}
78 : center: CENTER_OF_MULTICOLVAR DATA=tfcc
79 :
80 : dens: MULTICOLVARDENS ...
81 : DATA=tfcc ORIGIN=center DIR=xyz
82 : NBINS=80,80,80 BANDWIDTH=1.0,1.0,1.0 STRIDE=1 CLEAR=1
83 : ...
84 :
85 : FIND_CONTOUR GRID=dens CONTOUR=0.5 FILE=mycontour.xyz
86 : \endplumedfile
87 :
88 : */
89 : //+ENDPLUMEDOC
90 :
91 : namespace PLMD {
92 : namespace gridtools {
93 :
94 : class FindContour : public ContourFindingBase {
95 : private:
96 : bool firsttime;
97 : unsigned gbuffer;
98 : /// Stuff for output
99 : OFile of;
100 : double lenunit;
101 : std::string fmt_xyz;
102 : public:
103 : static void registerKeywords( Keywords& keys );
104 : explicit FindContour(const ActionOptions&ao);
105 0 : bool checkAllActive() const override {
106 0 : return gbuffer==0;
107 : }
108 : void prepareForAveraging() override;
109 0 : bool isPeriodic() override {
110 0 : return false;
111 : }
112 7 : unsigned getNumberOfQuantities() const override {
113 7 : return 1 + ingrid->getDimension();
114 : }
115 : void compute( const unsigned& current, MultiValue& myvals ) const override;
116 : void finishAveraging() override;
117 : };
118 :
119 13787 : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
120 :
121 5 : void FindContour::registerKeywords( Keywords& keys ) {
122 5 : ContourFindingBase::registerKeywords( keys );
123 : // We want a better way of doing this bit
124 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");
125 10 : keys.add("compulsory","FILE","file on which to output coordinates");
126 10 : keys.add("compulsory","UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
127 10 : keys.add("optional", "PRECISION","The number of digits in trajectory file");
128 5 : }
129 :
130 1 : FindContour::FindContour(const ActionOptions&ao):
131 : Action(ao),
132 : ContourFindingBase(ao),
133 1 : firsttime(true) {
134 :
135 1 : parse("BUFFER",gbuffer);
136 1 : if( gbuffer>0 ) {
137 0 : log.printf(" after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
138 : }
139 :
140 : std::string file;
141 2 : parse("FILE",file);
142 1 : if( file.length()==0 ) {
143 0 : error("name out output file was not specified");
144 : }
145 1 : std::string type=Tools::extension(file);
146 1 : log<<" file name "<<file<<"\n";
147 1 : if(type!="xyz") {
148 0 : error("can only print xyz file type with contour finding");
149 : }
150 :
151 : fmt_xyz="%f";
152 : std::string precision;
153 2 : parse("PRECISION",precision);
154 1 : if(precision.length()>0) {
155 : int p;
156 1 : Tools::convert(precision,p);
157 1 : log<<" with precision "<<p<<"\n";
158 : std::string a,b;
159 1 : Tools::convert(p+5,a);
160 1 : Tools::convert(p,b);
161 2 : fmt_xyz="%"+a+"."+b+"f";
162 : }
163 : std::string unitname;
164 2 : parse("UNITS",unitname);
165 1 : if(unitname!="PLUMED") {
166 0 : Units myunit;
167 0 : myunit.setLength(unitname);
168 0 : lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
169 0 : } else {
170 1 : lenunit=1.0;
171 : }
172 1 : of.link(*this);
173 1 : of.open(file);
174 1 : checkRead();
175 1 : mydata=buildDataStashes( NULL );
176 1 : }
177 :
178 2 : void FindContour::prepareForAveraging() {
179 : // Create a task list if first time
180 2 : if( firsttime ) {
181 16465 : for(unsigned i=0; i<ingrid->getDimension()*ingrid->getNumberOfPoints(); ++i) {
182 16464 : addTaskToList( i );
183 : }
184 : }
185 2 : firsttime=false;
186 2 : deactivateAllTasks();
187 :
188 : // We now need to identify the grid points that we need to search through
189 2 : std::vector<unsigned> nbin( ingrid->getNbin() );
190 2 : std::vector<unsigned> ind( ingrid->getDimension() );
191 2 : std::vector<unsigned> ones( ingrid->getDimension(), 1 );
192 : unsigned num_neighbours;
193 : std::vector<unsigned> neighbours;
194 10978 : for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
195 : // Ensure inactive grid points are ignored
196 10976 : if( ingrid->inactive(i) ) {
197 0 : continue;
198 : }
199 :
200 : // Get the index of the current grid point
201 10976 : ingrid->getIndices( i, ind );
202 10976 : ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
203 : bool cycle=false;
204 307328 : for(unsigned j=0; j<num_neighbours; ++j) {
205 296352 : if( ingrid->inactive( neighbours[j]) ) {
206 : cycle=true;
207 : break;
208 : }
209 : }
210 10976 : if( cycle ) {
211 0 : continue;
212 : }
213 :
214 : // Get the value of a point on the grid
215 10976 : double val1=getFunctionValue( i ) - contour;
216 : bool edge=false;
217 43904 : for(unsigned j=0; j<ingrid->getDimension(); ++j) {
218 : // Make sure we don't search at the edge of the grid
219 32928 : if( !ingrid->isPeriodic(j) && (ind[j]+1)==nbin[j] ) {
220 0 : continue;
221 32928 : } else if( (ind[j]+1)==nbin[j] ) {
222 : edge=true;
223 1960 : ind[j]=0;
224 : } else {
225 30968 : ind[j]+=1;
226 : }
227 32928 : double val2=getFunctionValue( ind ) - contour;
228 32928 : if( val1*val2<0 ) {
229 1108 : taskFlags[ ingrid->getDimension()*i + j ] = 1;
230 : }
231 32928 : if( ingrid->isPeriodic(j) && edge ) {
232 : edge=false;
233 1960 : ind[j]=nbin[j]-1;
234 : } else {
235 30968 : ind[j]-=1;
236 : }
237 : }
238 : }
239 2 : lockContributors();
240 2 : }
241 :
242 554 : void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
243 : // Retrieve the initial grid point coordinates
244 554 : unsigned gpoint = std::floor( current / ingrid->getDimension() );
245 554 : std::vector<double> point( ingrid->getDimension() );
246 554 : ingrid->getGridPointCoordinates( gpoint, point );
247 :
248 : // Retrieve the direction we are searching for the contour
249 554 : unsigned gdir = current%(ingrid->getDimension() );
250 554 : std::vector<double> direction( ingrid->getDimension(), 0 );
251 554 : direction[gdir] = 0.999999999*ingrid->getGridSpacing()[gdir];
252 :
253 : // Now find the contour
254 : findContour( direction, point );
255 : // And transfer to the store data vessel
256 2216 : for(unsigned i=0; i<ingrid->getDimension(); ++i) {
257 1662 : myvals.setValue( 1+i, point[i] );
258 : }
259 554 : }
260 :
261 1 : void FindContour::finishAveraging() {
262 : // And update the list of active grid points
263 1 : if( gbuffer>0 ) {
264 : std::vector<unsigned> neighbours;
265 : unsigned num_neighbours;
266 0 : std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
267 0 : std::vector<bool> active( ingrid->getNumberOfPoints(), false );
268 0 : std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
269 0 : for(unsigned i=0; i<getCurrentNumberOfActiveTasks(); ++i) {
270 : // Get the point we are operating on
271 0 : unsigned ipoint = std::floor( getActiveTask(i) / ingrid->getDimension() );
272 : // Get the indices of this point
273 0 : ingrid->getIndices( ipoint, ugrid_indices );
274 : // Now activate buffer region
275 0 : ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
276 0 : for(unsigned n=0; n<num_neighbours; ++n) {
277 0 : active[ neighbours[n] ]=true;
278 : }
279 : }
280 0 : ingrid->activateThesePoints( active );
281 : }
282 1 : std::vector<double> point( 1 + ingrid->getDimension() );
283 1 : of.printf("%u\n",mydata->getNumberOfStoredValues());
284 1 : of.printf("Points found on isocontour\n");
285 555 : for(unsigned i=0; i<mydata->getNumberOfStoredValues(); ++i) {
286 554 : mydata->retrieveSequentialValue( i, false, point );
287 554 : of.printf("X");
288 2216 : for(unsigned j=0; j<ingrid->getDimension(); ++j) {
289 3324 : of.printf( (" " + fmt_xyz).c_str(), lenunit*point[1+j] );
290 : }
291 554 : of.printf("\n");
292 : }
293 1 : }
294 :
295 : }
296 : }
|
