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/ActionWithArguments.h"
23 : #include "core/ActionPilot.h"
24 : #include "core/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "ActionWithGrid.h"
27 : #include "tools/OFile.h"
28 :
29 : namespace PLMD {
30 : namespace gridtools {
31 :
32 : //+PLUMEDOC GRIDANALYSIS DUMPCUBE
33 : /*
34 : Output a three dimensional grid using the Gaussian cube file format.
35 :
36 : Suppose you have calculated the value of a function on a three dimensional grid.
37 : This function might be a \ref HISTOGRAM or it might be a free energy energy surface
38 : that was calculated from this histogram by using \ref CONVERT_TO_FES. Alternatively,
39 : your function might be a phase-field that was calculated using \ref MULTICOLVARDENS.
40 : Whatever the function is, however, you obviously cannot show it using a typical contour
41 : plotting program such as gnuplot as you have three input variables.
42 :
43 : Tools like VMD have nice features for plotting these types of three dimensional functions
44 : but typically you are required to use a Gaussian cube file format to input the data. This
45 : action thus allows you to output a function evaluated on a grid to a Gaussian cube file format.
46 :
47 : \par Examples
48 :
49 : The input below can be used to post process a trajectory. A histogram as a function of the distance
50 : between atoms 1 and 2, the distance between atom 1 and 3 and the angle between the vector connecting atoms 1 and
51 : 2 and 1 and 3 is computed using kernel density estimation. Once all the data contained in the trajectory has been read in and
52 : all the kernels have been added the resulting histogram is output to a file called histoA1.cube. This file has the
53 : Gaussian cube file format. The histogram can thus be visualized using tools such as VMD.
54 :
55 : \plumedfile
56 : x1: DISTANCE ATOMS=1,2
57 : x2: DISTANCE ATOMS=1,3
58 : x3: ANGLE ATOMS=1,2,3
59 :
60 : hA1: HISTOGRAM ARG=x1,x2,x3 GRID_MIN=0.0,0.0,0.0 GRID_MAX=3.0,3.0,3.0 GRID_BIN=10,10,10 BANDWIDTH=1.0,1.0,1.0
61 : DUMPCUBE GRID=hA1 FILE=histoA1.cube
62 : \endplumedfile
63 :
64 : */
65 : //+ENDPLUMEDOC
66 :
67 : //+PLUMEDOC GRIDANALYSIS DUMPGRID
68 : /*
69 : Output the function on the grid to a file with the PLUMED grid format.
70 :
71 : PLUMED provides a number of actions that calculate the values of functions on grids.
72 : For instance, whenever you calculate a free energy as a function of a collective variable using
73 : \ref HISTOGRAM and \ref CONVERT_TO_FES you will generally want to output the value of the free energy at a number of points on a
74 : discrete grid that covers the CV space uniformly. Alternatively you may want to calculate
75 : what value some symmetry function takes at different points inside your simulation cell using \ref MULTICOLVARDENS.
76 :
77 : This action allows you to output these functions calculated on a grid using a format that can be read in using gnuplot
78 : and other such plotting programs. The file output using this action will have a header that contains some essential
79 : information about the function plotted and that looks something like this:
80 :
81 : \verbatim
82 : #! FIELDS x y hA1 dhA1_x dhA1_x
83 : #! SET normalisation 2.0000
84 : #! SET min_x 0.0
85 : #! SET max_x 3.0
86 : #! SET nbins_x 100
87 : #! SET periodic_x false
88 : #! SET min_y 0.0
89 : #! SET max_y 3.0
90 : #! SET nbins_y 100
91 : #! SET periodic_y false
92 : \endverbatim
93 :
94 : The header shown here tells us that we have grid showing the values that a function with two arguments x and y
95 : takes at various points in our cell. The lines beneath the first line then tell us a little bit about these two
96 : input arguments.
97 :
98 : The remaining lines of the file give us information on the positions of our grid points and the value the function and
99 : its partial derivatives with respect to x and y. If the header is as above a list of values of the function that have
100 : x=0 and 100 values of y between 0.0 and 3.0 will be provided. This block of data will be followed with a blank line.
101 : There will then be a second block of values which will all have been evaluated the same value of x and all possible values
102 : for y. This block is then followed by a blank line again and this pattern continues until all points of the grid have been covered.
103 :
104 : \par Examples
105 :
106 : The following input monitors two torsional angles during a simulation
107 : and outputs a continuous histogram as a function of them at the end of the simulation.
108 : \plumedfile
109 : TORSION ATOMS=1,2,3,4 LABEL=r1
110 : TORSION ATOMS=2,3,4,5 LABEL=r2
111 : HISTOGRAM ...
112 : ARG=r1,r2
113 : GRID_MIN=-3.14,-3.14
114 : GRID_MAX=3.14,3.14
115 : GRID_BIN=200,200
116 : BANDWIDTH=0.05,0.05
117 : LABEL=hh
118 : ... HISTOGRAM
119 :
120 : DUMPGRID GRID=hh FILE=histo
121 : \endplumedfile
122 :
123 : The following input monitors two torsional angles during a simulation
124 : and outputs a discrete histogram as a function of them at the end of the simulation.
125 : \plumedfile
126 : TORSION ATOMS=1,2,3,4 LABEL=r1
127 : TORSION ATOMS=2,3,4,5 LABEL=r2
128 : HISTOGRAM ...
129 : ARG=r1,r2
130 : KERNEL=DISCRETE
131 : GRID_MIN=-3.14,-3.14
132 : GRID_MAX=3.14,3.14
133 : GRID_BIN=200,200
134 : LABEL=hh
135 : ... HISTOGRAM
136 :
137 : DUMPGRID GRID=hh FILE=histo
138 : \endplumedfile
139 :
140 : The following input monitors two torsional angles during a simulation
141 : and outputs the histogram accumulated thus far every 100000 steps.
142 : \plumedfile
143 : TORSION ATOMS=1,2,3,4 LABEL=r1
144 : TORSION ATOMS=2,3,4,5 LABEL=r2
145 : HISTOGRAM ...
146 : ARG=r1,r2
147 : GRID_MIN=-3.14,-3.14
148 : GRID_MAX=3.14,3.14
149 : GRID_BIN=200,200
150 : BANDWIDTH=0.05,0.05
151 : LABEL=hh
152 : ... HISTOGRAM
153 :
154 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
155 : \endplumedfile
156 :
157 : The following input monitors two torsional angles during a simulation
158 : and outputs a separate histogram for each 100000 steps worth of trajectory.
159 : Notice how the CLEAR keyword is used here and how it is not used in the
160 : previous example.
161 :
162 : \plumedfile
163 : TORSION ATOMS=1,2,3,4 LABEL=r1
164 : TORSION ATOMS=2,3,4,5 LABEL=r2
165 : HISTOGRAM ...
166 : ARG=r1,r2 CLEAR=100000
167 : GRID_MIN=-3.14,-3.14
168 : GRID_MAX=3.14,3.14
169 : GRID_BIN=200,200
170 : BANDWIDTH=0.05,0.05
171 : LABEL=hh
172 : ... HISTOGRAM
173 :
174 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
175 : \endplumedfile
176 :
177 : */
178 : //+ENDPLUMEDOC
179 :
180 : class DumpGrid :
181 : public ActionWithArguments,
182 : public ActionPilot {
183 : private:
184 : std::string fmt, filename;
185 : bool onefile, xyzfile;
186 : public:
187 : static void registerKeywords( Keywords& keys );
188 : explicit DumpGrid(const ActionOptions&ao);
189 132 : ~DumpGrid() {}
190 105 : void calculate() override {}
191 105 : void apply() override {}
192 : void update() override ;
193 : };
194 :
195 : PLUMED_REGISTER_ACTION(DumpGrid,"DUMPCUBE")
196 : PLUMED_REGISTER_ACTION(DumpGrid,"DUMPGRID")
197 :
198 70 : void DumpGrid::registerKeywords( Keywords& keys ) {
199 70 : Action::registerKeywords( keys );
200 70 : ActionPilot::registerKeywords( keys );
201 70 : ActionWithArguments::registerKeywords( keys );
202 140 : keys.addInputKeyword("compulsory","ARG","grid","the label for the grid that you would like to output");
203 70 : keys.add("optional","GRID","the grid you would like to print (can also use ARG for specifying what is being printed)");
204 70 : keys.add("compulsory","STRIDE","0","the frequency with which the grid should be output to the file. Default of zero means dump at end of calculation");
205 70 : keys.add("compulsory","FILE","density","the file on which to write the grid.");
206 70 : keys.add("optional","FMT","the format that should be used to output real numbers");
207 70 : keys.addFlag("PRINT_XYZ",false,"output coordinates on fibonacci grid to xyz file");
208 70 : keys.addFlag("PRINT_ONE_FILE",false,"output grids one after the other in a single file");
209 70 : }
210 :
211 66 : DumpGrid::DumpGrid(const ActionOptions&ao):
212 : Action(ao),
213 : ActionWithArguments(ao),
214 : ActionPilot(ao),
215 66 : fmt("%f") {
216 66 : if( getNumberOfArguments()==0 ) {
217 : std::vector<Value*> grids;
218 0 : parseArgumentList("GRID",grids);
219 0 : requestArguments(grids);
220 : }
221 66 : if( getNumberOfArguments()!=1 ) {
222 0 : error("should only be one argument");
223 : }
224 66 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) {
225 0 : error("input should be a grid");
226 : }
227 66 : if( getName()=="DUMPCUBE" && getPntrToArgument(0)->getRank()!=3 ) {
228 0 : error("input should be a three dimensional grid");
229 : }
230 132 : parse("FILE",filename);
231 66 : if(filename.length()==0) {
232 0 : error("name out output file was not specified");
233 : }
234 66 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
235 66 : if( !ag ) {
236 0 : error( getPntrToArgument(0)->getName() + " is not grid");
237 : }
238 :
239 66 : log.printf(" outputting grid with label %s to file named %s",getPntrToArgument(0)->getName().c_str(), filename.c_str() );
240 66 : parse("FMT",fmt);
241 66 : log.printf(" with format %s \n", fmt.c_str() );
242 66 : if( getName()=="DUMPGRID" ) {
243 122 : fmt = " " + fmt;
244 5 : } else if( getName()=="DUMPCUBE" ) {
245 10 : fmt = fmt + " ";
246 : }
247 66 : parseFlag("PRINT_ONE_FILE", onefile);
248 66 : parseFlag("PRINT_XYZ",xyzfile);
249 66 : if( xyzfile ) {
250 4 : if( getName()=="DUMPCUBE" ) {
251 0 : error("PRINT_XYZ flag not compatible with DUMPCUBE");
252 : }
253 8 : if( ag->getGridCoordinatesObject().getGridType()=="flat" ) {
254 0 : error("can only use PRINT_XYZ option for fibonacci grids");
255 : }
256 4 : log.printf(" outputting grid to xyzfile\n");
257 : }
258 66 : if( onefile ) {
259 2 : log.printf(" printing all grids on a single file \n");
260 : } else {
261 64 : log.printf(" printing all grids on separate files \n");
262 : }
263 66 : }
264 :
265 137 : void DumpGrid::update() {
266 137 : OFile ofile;
267 137 : ofile.link(*this);
268 137 : if( onefile ) {
269 10 : ofile.enforceRestart();
270 10 : ofile.open( filename );
271 : } else {
272 127 : ofile.setBackupString("analysis");
273 127 : ofile.open( filename );
274 : }
275 :
276 137 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
277 137 : plumed_assert( ag );
278 137 : const GridCoordinatesObject & mygrid = ag->getGridCoordinatesObject();
279 :
280 137 : if( getName()=="DUMPCUBE" ) {
281 : double lunit=1.0;
282 11 : std::vector<std::string> argnames( ag->getGridCoordinateNames() );
283 :
284 11 : ofile.printf("PLUMED CUBE FILE\n");
285 11 : ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
286 : // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
287 : bool isdists=true;
288 11 : std::vector<double> extent(3);
289 11 : std::vector<std::string> min( mygrid.getMin() ), max( mygrid.getMax() );
290 44 : for(unsigned j=0; j<3; ++j) {
291 33 : if( argnames[j].find(".")!=std::string::npos ) {
292 3 : std::size_t dot = argnames[j].find(".");
293 3 : std::string name = argnames[j].substr(dot+1);
294 6 : if( name!="x" && name!="y" && name!="z" ) {
295 : isdists=false;
296 : }
297 : } else {
298 : isdists=false;
299 : }
300 :
301 : double mind, maxd;
302 33 : Tools::convert( min[j], mind );
303 33 : Tools::convert( max[j], maxd );
304 33 : if( mygrid.isPeriodic(j) ) {
305 3 : extent[j]=maxd-mind;
306 : } else {
307 30 : extent[j]=maxd-mind+mygrid.getGridSpacing()[j];
308 : }
309 : }
310 11 : if( isdists ) {
311 1 : if( plumed.usingNaturalUnits() ) {
312 : lunit = 1.0/0.5292;
313 : } else {
314 0 : lunit = plumed.getUnits().getLength()/.05929;
315 : }
316 : }
317 22 : std::string ostr = "%d " + fmt + fmt + fmt + "\n";
318 : Value* gval=getPntrToArgument(0);
319 11 : ofile.printf(ostr.c_str(),1,-0.5*lunit*extent[0],-0.5*lunit*extent[1],-0.5*lunit*extent[2]);
320 11 : ofile.printf(ostr.c_str(),mygrid.getNbin(true)[0],lunit*mygrid.getGridSpacing()[0],0.0,0.0); // Number of bins in each direction followed by
321 11 : ofile.printf(ostr.c_str(),mygrid.getNbin(true)[1],0.0,lunit*mygrid.getGridSpacing()[1],0.0); // shape of voxel
322 22 : ofile.printf(ostr.c_str(),mygrid.getNbin(true)[2],0.0,0.0,lunit*mygrid.getGridSpacing()[2]);
323 11 : ofile.printf(ostr.c_str(),1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
324 11 : std::vector<unsigned> pp(3);
325 11 : std::vector<unsigned> nbin( mygrid.getNbin(true) );
326 135 : for(pp[0]=0; pp[0]<nbin[0]; ++pp[0]) {
327 1530 : for(pp[1]=0; pp[1]<nbin[1]; ++pp[1]) {
328 20204 : for(pp[2]=0; pp[2]<nbin[2]; ++pp[2]) {
329 18798 : unsigned ival=pp[pp.size()-1];
330 56394 : for(unsigned i=pp.size()-1; i>0; --i) {
331 37596 : ival=ival*nbin[i-1]+pp[i-1];
332 : }
333 18798 : ofile.printf(fmt.c_str(), gval->get(ival) );
334 18798 : if(pp[2]%6==5) {
335 1994 : ofile.printf("\n");
336 : }
337 : }
338 1406 : ofile.printf("\n");
339 : }
340 : }
341 148 : } else if( xyzfile ) {
342 13 : std::vector<double> coords(3);
343 : Value* myarg = getPntrToArgument(0);
344 13 : unsigned nvals = myarg->getNumberOfValues();
345 13 : ofile.printf("%d\n\n", nvals);
346 1797 : for(unsigned i=0; i<nvals; ++i) {
347 1784 : double val = myarg->get(i);
348 1784 : mygrid.getGridPointCoordinates( i, coords );
349 1784 : ofile.printf("X");
350 7136 : for(unsigned j=0; j<3; ++j) {
351 10704 : ofile.printf((" " + fmt).c_str(), val*coords[j] );
352 : }
353 1784 : ofile.printf("\n");
354 : }
355 : } else {
356 : Value* gval=getPntrToArgument(0);
357 113 : std::vector<double> xx( gval->getRank() );
358 113 : std::vector<unsigned> ind( gval->getRank() );
359 113 : std::vector<std::string> argn( ag->getGridCoordinateNames() );
360 226 : if( mygrid.getGridType()=="fibonacci" ) {
361 12 : ofile.addConstantField("nbins");
362 : } else {
363 214 : plumed_assert( mygrid.getGridType()=="flat" );
364 229 : for(unsigned i=0; i<gval->getRank(); ++i) {
365 244 : ofile.addConstantField("min_" + argn[i] );
366 244 : ofile.addConstantField("max_" + argn[i] );
367 244 : ofile.addConstantField("nbins_" + argn[i] );
368 244 : ofile.addConstantField("periodic_" + argn[i] );
369 : }
370 : }
371 :
372 213260 : for(unsigned i=0; i<gval->getNumberOfValues(); ++i) {
373 : // Retrieve and print the grid coordinates
374 213147 : mygrid.getGridPointCoordinates( i, ind, xx );
375 213147 : if(i>0 && gval->getRank()==2 && ind[gval->getRank()-2]==0) {
376 986 : ofile.printf("\n");
377 : }
378 213147 : ofile.fmtField(fmt);
379 426294 : if( mygrid.getGridType()=="fibonacci" ) {
380 1200 : ofile.printField("nbins", static_cast<int>(gval->getNumberOfValues()) );
381 : } else {
382 682787 : for(unsigned j=0; j<gval->getRank(); ++j) {
383 940480 : ofile.printField("min_" + argn[j], mygrid.getMin()[j] );
384 940480 : ofile.printField("max_" + argn[j], mygrid.getMax()[j] );
385 940480 : ofile.printField("nbins_" + argn[j], static_cast<int>(mygrid.getNbin(false)[j]) );
386 470240 : if( mygrid.isPeriodic(j) ) {
387 360452 : ofile.printField("periodic_" + argn[j], "true" );
388 : } else {
389 580028 : ofile.printField("periodic_" + argn[j], "false" );
390 : }
391 : }
392 : }
393 : // Print the grid coordinates
394 685187 : for(unsigned j=0; j<gval->getRank(); ++j) {
395 472040 : ofile.fmtField(fmt);
396 472040 : ofile.printField(argn[j],xx[j]);
397 : }
398 : // Print value
399 213147 : ofile.fmtField(fmt);
400 213147 : ofile.printField( gval->getName(), gval->get(i) );
401 : // Print the derivatives
402 685187 : for(unsigned j=0; j<gval->getRank(); ++j) {
403 472040 : ofile.fmtField(fmt);
404 944080 : ofile.printField( "d" + gval->getName() + "_" + argn[j], gval->getGridDerivative(i,j) );
405 : }
406 213147 : ofile.printField();
407 : }
408 113 : }
409 137 : }
410 :
411 :
412 : }
413 : }
|
