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 "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, whenever 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 beneath 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 continuous 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 : LABEL=hh
134 : ... HISTOGRAM
135 :
136 : DUMPGRID GRID=hh FILE=histo STRIDE=100000
137 : \endplumedfile
138 :
139 : */
140 : //+ENDPLUMEDOC
141 :
142 : class DumpGrid : public GridPrintingBase {
143 : public:
144 : static void registerKeywords( Keywords& keys );
145 : explicit DumpGrid(const ActionOptions&ao);
146 : void printGrid( OFile& ofile ) const override;
147 : };
148 :
149 10517 : PLUMED_REGISTER_ACTION(DumpGrid,"DUMPGRID")
150 :
151 50 : void DumpGrid::registerKeywords( Keywords& keys ) {
152 50 : GridPrintingBase::registerKeywords( keys );
153 50 : }
154 :
155 49 : DumpGrid::DumpGrid(const ActionOptions&ao):
156 : Action(ao),
157 49 : GridPrintingBase(ao)
158 : {
159 98 : if( ingrid->getType()!="flat" ) error("cannot dump grid of type " + ingrid->getType() + " using DUMPGRID");
160 98 : fmt = " " + fmt; checkRead();
161 49 : }
162 :
163 60 : void DumpGrid::printGrid( OFile& ofile ) const {
164 60 : ofile.addConstantField("normalisation");
165 128 : for(unsigned i=0; i<ingrid->getDimension(); ++i) {
166 136 : ofile.addConstantField("min_" + ingrid->getComponentName(i) );
167 136 : ofile.addConstantField("max_" + ingrid->getComponentName(i) );
168 136 : ofile.addConstantField("nbins_" + ingrid->getComponentName(i) );
169 136 : ofile.addConstantField("periodic_" + ingrid->getComponentName(i) );
170 : }
171 :
172 60 : std::vector<double> xx( ingrid->getDimension() );
173 60 : std::vector<unsigned> ind( ingrid->getDimension() );
174 103636 : for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
175 103576 : ingrid->getIndices( i, ind );
176 103576 : if(i>0 && ingrid->getDimension()==2 && ind[ingrid->getDimension()-2]==0) ofile.printf("\n");
177 103576 : ofile.fmtField(fmt); ofile.printField("normalisation", ingrid->getNorm() );
178 306047 : for(unsigned j=0; j<ingrid->getDimension(); ++j) {
179 404942 : ofile.printField("min_" + ingrid->getComponentName(j), ingrid->getMin()[j] );
180 404942 : ofile.printField("max_" + ingrid->getComponentName(j), ingrid->getMax()[j] );
181 404942 : ofile.printField("nbins_" + ingrid->getComponentName(j), static_cast<int>(ingrid->getNbin()[j]) );
182 219655 : if( ingrid->isPeriodic(j) ) ofile.printField("periodic_" + ingrid->getComponentName(j), "true" );
183 370574 : else ofile.printField("periodic_" + ingrid->getComponentName(j), "false" );
184 : }
185 : // Retrieve and print the grid coordinates
186 103576 : ingrid->getGridPointCoordinates(i, xx );
187 306047 : for(unsigned j=0; j<ingrid->getDimension(); ++j) { ofile.fmtField(fmt); ofile.printField(ingrid->getComponentName(j),xx[j]); }
188 440087 : for(unsigned j=0; j<ingrid->getNumberOfQuantities(); ++j) {
189 336511 : ofile.fmtField(fmt); ofile.printField(ingrid->arg_names[ingrid->dimension+j], ingrid->getGridElement( i, j ) );
190 : }
191 103576 : ofile.printField();
192 : }
193 60 : }
194 :
195 : }
196 : }
|
