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 140 : keys.add("optional","GRID","the grid you would like to print (can also use ARG for specifying what is being printed)");
204 140 : 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 140 : keys.add("compulsory","FILE","density","the file on which to write the grid.");
206 140 : keys.add("optional","FMT","the format that should be used to output real numbers");
207 140 : keys.addFlag("PRINT_XYZ",false,"output coordinates on fibonacci grid to xyz file");
208 140 : 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 : {
217 66 : if( getNumberOfArguments()==0 ) {
218 0 : std::vector<Value*> grids; parseArgumentList("GRID",grids); requestArguments(grids);
219 : }
220 66 : if( getNumberOfArguments()!=1 ) error("should only be one argument");
221 66 : if( getPntrToArgument(0)->getRank()==0 || !getPntrToArgument(0)->hasDerivatives() ) error("input should be a grid");
222 66 : if( getName()=="DUMPCUBE" && getPntrToArgument(0)->getRank()!=3 ) error("input should be a three dimensional grid");
223 132 : parse("FILE",filename);
224 66 : if(filename.length()==0) error("name out output file was not specified");
225 66 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
226 66 : if( !ag ) error( getPntrToArgument(0)->getName() + " is not grid");
227 :
228 66 : log.printf(" outputting grid with label %s to file named %s",getPntrToArgument(0)->getName().c_str(), filename.c_str() );
229 132 : parse("FMT",fmt); log.printf(" with format %s \n", fmt.c_str() );
230 137 : if( getName()=="DUMPGRID" ) fmt = " " + fmt; else if( getName()=="DUMPCUBE" ) fmt = fmt + " ";
231 132 : parseFlag("PRINT_ONE_FILE", onefile); parseFlag("PRINT_XYZ",xyzfile);
232 66 : if( xyzfile ) {
233 4 : if( getName()=="DUMPCUBE" ) error("PRINT_XYZ flag not compatible with DUMPCUBE");
234 8 : if( ag->getGridCoordinatesObject().getGridType()=="flat" ) error("can only use PRINT_XYZ option for fibonacci grids");
235 4 : log.printf(" outputting grid to xyzfile\n");
236 : }
237 66 : if( onefile ) log.printf(" printing all grids on a single file \n");
238 64 : else log.printf(" printing all grids on separate files \n");
239 66 : }
240 :
241 137 : void DumpGrid::update() {
242 137 : OFile ofile; ofile.link(*this);
243 137 : if( onefile ) {
244 10 : ofile.enforceRestart();
245 10 : ofile.open( filename );
246 : } else {
247 127 : ofile.setBackupString("analysis");
248 127 : ofile.open( filename );
249 : }
250 :
251 137 : ActionWithGrid* ag = ActionWithGrid::getInputActionWithGrid( getPntrToArgument(0)->getPntrToAction() );
252 137 : plumed_assert( ag ); const GridCoordinatesObject & mygrid = ag->getGridCoordinatesObject();
253 :
254 137 : if( getName()=="DUMPCUBE" ) {
255 11 : double lunit=1.0; std::vector<std::string> argnames( ag->getGridCoordinateNames() );
256 :
257 11 : ofile.printf("PLUMED CUBE FILE\n");
258 11 : ofile.printf("OUTER LOOP: X, MIDDLE LOOP: Y, INNER LOOP: Z\n");
259 : // Number of atoms followed by position of origin (origin set so that center of grid is in center of cell)
260 11 : bool isdists=true; std::vector<double> extent(3); std::vector<std::string> min( mygrid.getMin() ), max( mygrid.getMax() );
261 44 : for(unsigned j=0; j<3; ++j) {
262 33 : if( argnames[j].find(".")!=std::string::npos ) {
263 3 : std::size_t dot = argnames[j].find(".");
264 3 : std::string name = argnames[j].substr(dot+1);
265 6 : if( name!="x" && name!="y" && name!="z" ) isdists=false;
266 : } else isdists=false;
267 :
268 33 : double mind, maxd; Tools::convert( min[j], mind ); Tools::convert( max[j], maxd );
269 33 : if( mygrid.isPeriodic(j) ) extent[j]=maxd-mind;
270 30 : else { extent[j]=maxd-mind+mygrid.getGridSpacing()[j]; }
271 : }
272 11 : if( isdists ) {
273 1 : if( plumed.usingNaturalUnits() ) lunit = 1.0/0.5292;
274 0 : else lunit = plumed.getUnits().getLength()/.05929;
275 : }
276 22 : std::string ostr = "%d " + fmt + fmt + fmt + "\n"; Value* gval=getPntrToArgument(0);
277 11 : ofile.printf(ostr.c_str(),1,-0.5*lunit*extent[0],-0.5*lunit*extent[1],-0.5*lunit*extent[2]);
278 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
279 11 : ofile.printf(ostr.c_str(),mygrid.getNbin(true)[1],0.0,lunit*mygrid.getGridSpacing()[1],0.0); // shape of voxel
280 22 : ofile.printf(ostr.c_str(),mygrid.getNbin(true)[2],0.0,0.0,lunit*mygrid.getGridSpacing()[2]);
281 11 : ofile.printf(ostr.c_str(),1,0.0,0.0,0.0); // Fake atom otherwise VMD doesn't work
282 11 : std::vector<unsigned> pp(3); std::vector<unsigned> nbin( mygrid.getNbin(true) );
283 135 : for(pp[0]=0; pp[0]<nbin[0]; ++pp[0]) {
284 1530 : for(pp[1]=0; pp[1]<nbin[1]; ++pp[1]) {
285 20204 : for(pp[2]=0; pp[2]<nbin[2]; ++pp[2]) {
286 18798 : unsigned ival=pp[pp.size()-1];
287 56394 : for(unsigned i=pp.size()-1; i>0; --i) ival=ival*nbin[i-1]+pp[i-1];
288 18798 : ofile.printf(fmt.c_str(), gval->get(ival) );
289 18798 : if(pp[2]%6==5) ofile.printf("\n");
290 : }
291 1406 : ofile.printf("\n");
292 : }
293 : }
294 148 : } else if( xyzfile ) {
295 13 : std::vector<double> coords(3);
296 : Value* myarg = getPntrToArgument(0);
297 13 : unsigned nvals = myarg->getNumberOfValues(); ofile.printf("%d\n\n", nvals);
298 1797 : for(unsigned i=0; i<nvals; ++i) {
299 1784 : double val = myarg->get(i); mygrid.getGridPointCoordinates( i, coords );
300 1784 : ofile.printf("X");
301 7136 : for(unsigned j=0; j<3; ++j) ofile.printf((" " + fmt).c_str(), val*coords[j] );
302 1784 : ofile.printf("\n");
303 : }
304 : } else {
305 : Value* gval=getPntrToArgument(0);
306 113 : std::vector<double> xx( gval->getRank() );
307 113 : std::vector<unsigned> ind( gval->getRank() );
308 113 : std::vector<std::string> argn( ag->getGridCoordinateNames() );
309 226 : if( mygrid.getGridType()=="fibonacci" ) {
310 12 : ofile.addConstantField("nbins");
311 : } else {
312 214 : plumed_assert( mygrid.getGridType()=="flat" );
313 229 : for(unsigned i=0; i<gval->getRank(); ++i) {
314 244 : ofile.addConstantField("min_" + argn[i] );
315 244 : ofile.addConstantField("max_" + argn[i] );
316 244 : ofile.addConstantField("nbins_" + argn[i] );
317 244 : ofile.addConstantField("periodic_" + argn[i] );
318 : }
319 : }
320 :
321 213260 : for(unsigned i=0; i<gval->getNumberOfValues(); ++i) {
322 : // Retrieve and print the grid coordinates
323 213147 : mygrid.getGridPointCoordinates( i, ind, xx );
324 213147 : if(i>0 && gval->getRank()==2 && ind[gval->getRank()-2]==0) ofile.printf("\n");
325 213147 : ofile.fmtField(fmt);
326 426294 : if( mygrid.getGridType()=="fibonacci" ) {
327 1200 : ofile.printField("nbins", static_cast<int>(gval->getNumberOfValues()) );
328 : } else {
329 682787 : for(unsigned j=0; j<gval->getRank(); ++j) {
330 940480 : ofile.printField("min_" + argn[j], mygrid.getMin()[j] );
331 940480 : ofile.printField("max_" + argn[j], mygrid.getMax()[j] );
332 940480 : ofile.printField("nbins_" + argn[j], static_cast<int>(mygrid.getNbin(false)[j]) );
333 650466 : if( mygrid.isPeriodic(j) ) ofile.printField("periodic_" + argn[j], "true" );
334 580028 : else ofile.printField("periodic_" + argn[j], "false" );
335 : }
336 : }
337 : // Print the grid coordinates
338 685187 : for(unsigned j=0; j<gval->getRank(); ++j) { ofile.fmtField(fmt); ofile.printField(argn[j],xx[j]); }
339 : // Print value
340 213147 : ofile.fmtField(fmt); ofile.printField( gval->getName(), gval->get(i) );
341 : // Print the derivatives
342 685187 : for(unsigned j=0; j<gval->getRank(); ++j) { ofile.fmtField(fmt); ofile.printField( "d" + gval->getName() + "_" + argn[j], gval->getGridDerivative(i,j) ); }
343 213147 : ofile.printField();
344 : }
345 113 : }
346 137 : }
347 :
348 :
349 : }
350 : }
|
