Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "GridPrintingBase.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/OFile.h"
25 :
26 : namespace PLMD {
27 : namespace gridtools {
28 :
29 : //+PLUMEDOC GRIDANALYSIS DUMPGRID
30 : /*
31 : Output the function on the grid to a file with the PLUMED grid format.
32 :
33 : PLUMED provides a number of actions that calculate the values of functions on grids.
34 : For instance, whenver you calculate a free energy as a function of a collective variable using
35 : \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
36 : discrete grid that covers the CV space uniformly. Alternatively you may want to calculate
37 : what value some symmetry function takes at different points inside your simulation cell using \ref MULTICOLVARDENS.
38 :
39 : This action allows you to output these functions calculated on a grid using a format that can be read in using gnuplot
40 : and other such plotting programs. The file output using this action will have a header that contains some essential
41 : information about the function plotted and that looks something like this:
42 :
43 : \verbatim
44 : #! FIELDS x y hA1 dhA1_x dhA1_x
45 : #! SET normalisation 2.0000
46 : #! SET min_x 0.0
47 : #! SET max_x 3.0
48 : #! SET nbins_x 100
49 : #! SET periodic_x false
50 : #! SET min_y 0.0
51 : #! SET max_y 3.0
52 : #! SET nbins_y 100
53 : #! SET periodic_y false
54 : \endverbatim
55 :
56 : The header shown here tells us that we have grid showing the values that a function with two arguments x and y
57 : takes at various points in our cell. The lines beheath the first line then tell us a little bit about these two
58 : input arguments.
59 :
60 : The remaining lines of the file give us information on the positions of our grid points and the value the function and
61 : its partial derivatives with respect to x and y. If the header is as above a list of values of the function that have
62 : 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.
63 : There will then be a second block of values which will all have been evaluated the same value of x and all possible values
64 : 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.
65 :
66 : \par Examples
67 :
68 : The following input monitors two torsional angles during a simulation
69 : and outputs a continuos histogram as a function of them at the end of the simulation.
70 : \plumedfile
71 : TORSION ATOMS=1,2,3,4 LABEL=r1
72 : TORSION ATOMS=2,3,4,5 LABEL=r2
73 : HISTOGRAM ...
74 : ARG=r1,r2
75 : GRID_MIN=-3.14,-3.14
76 : GRID_MAX=3.14,3.14
77 : GRID_BIN=200,200
78 : BANDWIDTH=0.05,0.05
79 : LABEL=hh
80 : ... HISTOGRAM
81 :
82 : DUMPGRID GRID=hh FILE=histo
83 : \endplumedfile
84 :
85 : The following input monitors two torsional angles during a simulation
86 : and outputs a discrete histogram as a function of them at the end of the simulation.
87 : \plumedfile
88 : TORSION ATOMS=1,2,3,4 LABEL=r1
89 : TORSION ATOMS=2,3,4,5 LABEL=r2
90 : HISTOGRAM ...
91 : ARG=r1,r2
92 : KERNEL=DISCRETE
93 : GRID_MIN=-3.14,-3.14
94 : GRID_MAX=3.14,3.14
95 : GRID_BIN=200,200
96 : LABEL=hh
97 : ... HISTOGRAM
98 :
99 : DUMPGRID GRID=hh FILE=histo
100 : \endplumedfile
101 :
102 : The following input monitors two torsional angles during a simulation
103 : and outputs the histogram accumulated thus far every 100000 steps.
104 : \plumedfile
105 : TORSION ATOMS=1,2,3,4 LABEL=r1
106 : TORSION ATOMS=2,3,4,5 LABEL=r2
107 : HISTOGRAM ...
108 : ARG=r1,r2
109 : GRID_MIN=-3.14,-3.14
110 : GRID_MAX=3.14,3.14
111 : GRID_BIN=200,200
112 : BANDWIDTH=0.05,0.05
113 : LABEL=hh
114 : ... HISTOGRAM
115 :
116 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
117 : \endplumedfile
118 :
119 : The following input monitors two torsional angles during a simulation
120 : and outputs a separate histogram for each 100000 steps worth of trajectory.
121 : Notice how the CLEAR keyword is used here and how it is not used in the
122 : previous example.
123 :
124 : \plumedfile
125 : TORSION ATOMS=1,2,3,4 LABEL=r1
126 : TORSION ATOMS=2,3,4,5 LABEL=r2
127 : HISTOGRAM ...
128 : ARG=r1,r2 CLEAR=100000
129 : GRID_MIN=-3.14,-3.14
130 : GRID_MAX=3.14,3.14
131 : GRID_BIN=200,200
132 : BANDWIDTH=0.05,0.05
133 : GRID_WFILE=histo
134 : LABEL=hh
135 : ... HISTOGRAM
136 :
137 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
138 : \endplumedfile
139 :
140 : */
141 : //+ENDPLUMEDOC
142 :
143 38 : class DumpGrid : public GridPrintingBase {
144 : public:
145 : static void registerKeywords( Keywords& keys );
146 : explicit DumpGrid(const ActionOptions&ao);
147 : void printGrid( OFile& ofile ) const ;
148 : };
149 :
150 6490 : PLUMED_REGISTER_ACTION(DumpGrid,"DUMPGRID")
151 :
152 39 : void DumpGrid::registerKeywords( Keywords& keys ) {
153 39 : GridPrintingBase::registerKeywords( keys );
154 39 : }
155 :
156 38 : DumpGrid::DumpGrid(const ActionOptions&ao):
157 : Action(ao),
158 38 : GridPrintingBase(ao)
159 : {
160 76 : if( ingrid->getType()!="flat" ) error("cannot dump grid of type " + ingrid->getType() + " using DUMPGRID");
161 76 : fmt = " " + fmt; checkRead();
162 38 : }
163 :
164 53 : void DumpGrid::printGrid( OFile& ofile ) const {
165 106 : ofile.addConstantField("normalisation");
166 283 : for(unsigned i=0; i<ingrid->getDimension(); ++i) {
167 177 : ofile.addConstantField("min_" + ingrid->getComponentName(i) );
168 236 : ofile.addConstantField("max_" + ingrid->getComponentName(i) );
169 236 : ofile.addConstantField("nbins_" + ingrid->getComponentName(i) );
170 236 : ofile.addConstantField("periodic_" + ingrid->getComponentName(i) );
171 : }
172 :
173 53 : std::vector<double> xx( ingrid->getDimension() );
174 106 : std::vector<unsigned> ind( ingrid->getDimension() );
175 180948 : for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
176 90421 : ingrid->getIndices( i, ind );
177 266977 : if(i>0 && ingrid->getDimension()==2 && ind[ingrid->getDimension()-2]==0) ofile.printf("\n");
178 271263 : ofile.fmtField(fmt); ofile.printField("normalisation", ingrid->getNorm() );
179 710687 : for(unsigned j=0; j<ingrid->getDimension(); ++j) {
180 883075 : ofile.printField("min_" + ingrid->getComponentName(j), ingrid->getMin()[j] );
181 883075 : ofile.printField("max_" + ingrid->getComponentName(j), ingrid->getMax()[j] );
182 1059690 : ofile.printField("nbins_" + ingrid->getComponentName(j), static_cast<int>(ingrid->getNbin()[j]) );
183 401566 : if( ingrid->isPeriodic(j) ) ofile.printField("periodic_" + ingrid->getComponentName(j), "true" );
184 822655 : else ofile.printField("periodic_" + ingrid->getComponentName(j), "false" );
185 : }
186 : // Retrieve and print the grid coordinates
187 90421 : ingrid->getGridPointCoordinates(i, xx );
188 1063917 : for(unsigned j=0; j<ingrid->getDimension(); ++j) { ofile.fmtField(fmt); ofile.printField(ingrid->getComponentName(j),xx[j]); }
189 714738 : for(unsigned j=0; j<ingrid->getNumberOfQuantities(); ++j) {
190 533896 : ofile.fmtField(fmt); ofile.printField(ingrid->arg_names[ingrid->dimension+j], ingrid->getGridElement( i, j ) );
191 : }
192 90421 : ofile.printField();
193 : }
194 53 : }
195 :
196 : }
197 4839 : }
|
