Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2011-2018 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 "ActionWithGrid.h"
23 : #include "core/ActionRegister.h"
24 : #include "lepton/Lepton.h"
25 : #include "tools/IFile.h"
26 :
27 : //+PLUMEDOC GRIDCALC REFERENCE_GRID
28 : /*
29 : Setup a constant grid by either reading values from a file or definining a function in input
30 :
31 : This action allows you to create a constant function and store the values that this function takes
32 : at a grid of points. There are two ways that you can create the constant function. You can either
33 : read in the value of the function at a set of grid points from a file as is done in this example input:
34 :
35 : ```plumed
36 : #SETTINGS INPUTFILES=regtest/gridtools/rt-test-fib-read/kde.grid
37 : ffg: REFERENCE_GRID FILE=regtest/gridtools/rt-test-fib-read/kde.grid VALUE=h
38 : DUMPGRID ARG=ffg FILE=output.grid
39 : ```
40 :
41 : The input file in this case has the format for grids that is discussed in the documentation for [DUMPGRID](DUMPGRID.md).
42 :
43 : The second way you can create these constant functions is illustrated below.
44 :
45 : ```plumed
46 : d2: REFERENCE_GRID GRID_MIN=0 GRID_MAX=10 GRID_BIN=20 FUNC=d1*d1 VAR=d1 PERIODIC=NO
47 : d1: DISTANCE ATOMS=1,2
48 : ff: EVALUATE_FUNCTION_FROM_GRID GRID=d2
49 : PRINT ARG=d1,ff FMT=%8.4f FILE=colvar
50 : ```
51 :
52 : This input illustrates a rather elaborate (and approximate) method for evaluating the square of the distance between atoms 1 and 2.
53 : As you can see, the option involves using the GRID_MIN, GRID_MAX and GRID_BIN options to create a grid of points. You can then use the FUNC option
54 : to specify the function that should be evaluated at each of those grid points.
55 : If you use this option the lepton library that is discussed in the documentation for [CUSTOM](CUSTOM.md) is used to evaluate the
56 : function at the various grid points.
57 :
58 : N.B. This method with lepton was implemented to facilitate the implementation of the normalisation in the implementation of the [RDF](RDF.md) shortcut.
59 :
60 : Lastly note that you can specify the grid spacing in the input to this action rather than the number of bins as shown below:
61 :
62 : ```plumed
63 : d2: REFERENCE_GRID GRID_MIN=0 GRID_MAX=10 GRID_SPACING=0.5 FUNC=d1*d1 VAR=d1 PERIODIC=NO
64 : d1: DISTANCE ATOMS=1,2
65 : ff: EVALUATE_FUNCTION_FROM_GRID GRID=d2
66 : PRINT ARG=d1,ff FMT=%8.4f FILE=colvar
67 : ```
68 :
69 : */
70 : //+ENDPLUMEDOC
71 :
72 : namespace PLMD {
73 : namespace gridtools {
74 :
75 : static std::map<std::string, double> leptonConstants= {
76 : {"e", std::exp(1.0)},
77 : {"log2e", 1.0/std::log(2.0)},
78 : {"log10e", 1.0/std::log(10.0)},
79 : {"ln2", std::log(2.0)},
80 : {"ln10", std::log(10.0)},
81 : {"pi", pi},
82 : {"pi_2", pi*0.5},
83 : {"pi_4", pi*0.25},
84 : // {"1_pi", 1.0/pi},
85 : // {"2_pi", 2.0/pi},
86 : // {"2_sqrtpi", 2.0/std::sqrt(pi)},
87 : {"sqrt2", std::sqrt(2.0)},
88 : {"sqrt1_2", std::sqrt(0.5)}
89 : };
90 :
91 : class ReadGridInSetup : public ActionWithGrid {
92 : private:
93 : GridCoordinatesObject gridobject;
94 : std::vector<std::string> dernames;
95 : void createGridAndValue( const std::string& gtype, const std::vector<bool>& ipbc, const unsigned& nfermi,
96 : const std::vector<std::string>& gmin, const std::vector<std::string>& gmax,
97 : const std::vector<std::size_t>& gbin, std::vector<double>& gspacing );
98 : public:
99 : static void registerKeywords( Keywords& keys );
100 : explicit ReadGridInSetup(const ActionOptions&ao);
101 : unsigned getNumberOfDerivatives() override ;
102 : std::vector<std::string> getGridCoordinateNames() const override ;
103 : const GridCoordinatesObject& getGridCoordinatesObject() const override ;
104 0 : void calculate() override {}
105 : };
106 :
107 : PLUMED_REGISTER_ACTION(ReadGridInSetup,"REFERENCE_GRID")
108 :
109 41 : void ReadGridInSetup::registerKeywords( Keywords& keys ) {
110 41 : ActionWithGrid::registerKeywords(keys);
111 41 : keys.add("optional","FUNC","the function to compute on the grid");
112 41 : keys.add("compulsory","GRID_MIN","auto","the lower bounds for the grid");
113 41 : keys.add("compulsory","GRID_MAX","auto","the upper bounds for the grid");
114 41 : keys.add("compulsory","PERIODIC","are the grid directions periodic");
115 41 : keys.add("optional","GRID_BIN","the number of bins for the grid");
116 41 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
117 41 : keys.add("optional","VAR","the names to give each of the grid directions in the function. If you have up to three grid coordinates in your function you can use x, y and z to refer to them. Otherwise you must use this flag to give your variables names.");
118 41 : keys.add("compulsory","FILE","the name of the file that contains the reference data");
119 41 : keys.add("compulsory","VALUE","the name of the value that should be read from the grid");
120 82 : keys.setValueDescription("grid","the constant function on the grid that was specified in input");
121 41 : }
122 :
123 23 : ReadGridInSetup::ReadGridInSetup(const ActionOptions&ao):
124 : Action(ao),
125 23 : ActionWithGrid(ao) {
126 : std::string func;
127 48 : parse("FUNC",func);
128 23 : if( func.length()>0 ) {
129 : // Read in stuff for grid
130 16 : std::vector<std::string> gmin;
131 32 : parseVector("GRID_MIN",gmin);
132 16 : std::vector<std::string> gmax(gmin.size());
133 32 : parseVector("GRID_MAX",gmax);
134 16 : std::vector<std::size_t> gbin(gmin.size());
135 32 : parseVector("GRID_BIN",gbin);
136 16 : std::vector<double> spacing(gmin.size());
137 32 : parseVector("GRID_SPACING",spacing);
138 16 : if( gbin.size()!=gmin.size() && spacing.size()!=gmin.size() ) {
139 1 : error("GRID_BIN or GRID_SPACING must be set");
140 : }
141 15 : std::vector<std::string> pbc(gmin.size());
142 30 : parseVector("PERIODIC",pbc);
143 16 : std::vector<bool> ipbc( pbc.size() );
144 32 : for(unsigned i=0; i<ipbc.size(); ++i) {
145 17 : if( pbc[i]=="YES" ) {
146 : ipbc[i]=true;
147 17 : } else if( pbc[i]=="NO" ) {
148 : ipbc[i]=false;
149 : } else {
150 0 : error( pbc[i] + " is not a valid instruction to the PERIODIC keyword");
151 : }
152 : }
153 :
154 : // Read in the variables
155 30 : parseVector("VAR",dernames);
156 15 : if(dernames.size()==0) {
157 10 : dernames.resize(gmin.size());
158 10 : if(gmin.size()>3) {
159 0 : error("Using more than 3 arguments you should explicitly write their names with VAR");
160 : }
161 10 : if(dernames.size()>0) {
162 : dernames[0]="x";
163 : }
164 10 : if(dernames.size()>1) {
165 : dernames[1]="y";
166 : }
167 10 : if(dernames.size()>2) {
168 : dernames[2]="z";
169 : }
170 : }
171 15 : if(dernames.size()!=gmin.size()) {
172 0 : error("Size of VAR array should be the same as number of grid dimensions");
173 : }
174 :
175 : // Create the grid and the value of the grid
176 15 : createGridAndValue( "flat", ipbc, 0, gmin, gmax, gbin, spacing );
177 :
178 : // Read in stuff for function
179 14 : log.printf(" evaluating function : %s\n",func.c_str());
180 14 : log.printf(" with variables :");
181 30 : for(unsigned i=0; i<dernames.size(); i++) {
182 16 : log.printf(" %s",dernames[i].c_str());
183 : }
184 14 : log.printf("\n");
185 28 : gbin=gridobject.getNbin(false);
186 14 : log.printf(" on %zu", gbin[0]);
187 16 : for(unsigned i=1; i<gbin.size(); ++i) {
188 2 : log.printf(" by %zu \n", gbin[i]);
189 : }
190 14 : log.printf(" grid of points between (%s", gmin[0].c_str() );
191 16 : for(unsigned i=1; i<gmin.size(); ++i) {
192 2 : log.printf(", %s", gmin[i].c_str() );
193 : }
194 14 : log.printf(") and (%s", gmax[0].c_str() );
195 16 : for(unsigned i=1; i<gmax.size(); ++i) {
196 2 : log.printf(", %s", gmax[i].c_str() );
197 : }
198 14 : log.printf(")\n");
199 :
200 43 : lepton::CompiledExpression expression=[&](const lepton::ParsedExpression& pe) {
201 14 : log<<" function as parsed by lepton: "<<pe<<"\n";
202 14 : return pe.createCompiledExpression();
203 : }
204 15 : (lepton::Parser::parse(func).optimize(leptonConstants));
205 :
206 28 : for(auto &p: expression.getVariables()) {
207 14 : if(std::find(dernames.begin(),dernames.end(),p)==dernames.end()) {
208 0 : error("variable " + p + " is not defined");
209 : }
210 : }
211 14 : log<<" derivatives as computed by lepton:\n";
212 14 : std::vector<lepton::CompiledExpression> expression_deriv( dernames.size() );
213 30 : for(unsigned i=0; i<dernames.size(); i++) {
214 16 : lepton::ParsedExpression pe=lepton::Parser::parse(func).differentiate(dernames[i]).optimize(leptonConstants);
215 16 : log<<" "<<pe<<"\n";
216 16 : expression_deriv[i]=pe.createCompiledExpression();
217 : }
218 : // And finally calculate all the grid points
219 14 : std::vector<double> dder( dernames.size() ), xx( dernames.size() );
220 14 : Value* valout=getPntrToComponent(0);
221 1662 : for(unsigned index=0; index<valout->getNumberOfValues(); ++index) {
222 1648 : gridobject.getGridPointCoordinates( index, xx );
223 4673 : for(unsigned j=0; j<xx.size(); ++j) {
224 : try {
225 3025 : expression.getVariableReference(dernames[j])=xx[j];
226 1377 : } catch(PLMD::lepton::Exception& exc) {
227 : // this is necessary since in some cases lepton things a variable is not present even though it is present
228 : // e.g. func=0*x
229 1377 : }
230 : }
231 1648 : valout->set( index, expression.evaluate() );
232 4673 : for(unsigned k=0; k<xx.size(); ++k) {
233 8804 : for(unsigned j=0; j<xx.size(); ++j) {
234 : try {
235 5779 : expression_deriv[k].getVariableReference(dernames[j])=xx[j];
236 4131 : } catch(PLMD::lepton::Exception& exc) {
237 : // this is necessary since in some cases lepton things a variable is not present even though it is present
238 : // e.g. func=0*x
239 4131 : }
240 : }
241 3025 : valout->addGridDerivatives( index, k, expression_deriv[k].evaluate() );
242 : }
243 : }
244 37 : } else {
245 : std::string valuestr;
246 14 : parse("VALUE",valuestr);
247 : std::string tstyle, filen;
248 7 : parse("FILE",filen);
249 7 : log.printf(" reading function %s on grid from file %s \n", valuestr.c_str(), filen.c_str() );
250 7 : IFile ifile;
251 7 : ifile.open(filen);
252 7 : if( !ifile.FieldExist( valuestr ) ) {
253 0 : error("could not find grid value in input file");
254 : }
255 7 : std::vector<std::string> fieldnames;
256 7 : ifile.scanFieldList( fieldnames );
257 :
258 : // Retrieve the names of the variables the grid is computed over
259 : bool flatgrid=false;
260 68 : for(unsigned i=0; i<fieldnames.size(); ++i) {
261 61 : if( fieldnames[i].find("min_")!=std::string::npos ) {
262 : flatgrid=true;
263 : }
264 122 : std::size_t dot = fieldnames[i].find_first_of("d" + valuestr + "_" );
265 122 : if( fieldnames[i].find("d" + valuestr + "_")!=std::string::npos ) {
266 40 : dernames.push_back( fieldnames[i].substr(dot+2+valuestr.length()) );
267 : }
268 : }
269 7 : if( flatgrid && dernames.size()==0 ) {
270 0 : error("could not find any derivatives for value " + valuestr + " in input file. Header should contain at least columns with a name starting d" + valuestr + "_");
271 : }
272 : // Now get all the header data for the grid
273 7 : std::vector<std::string> gmin( dernames.size() ), gmax( dernames.size() );
274 : std::string pstring;
275 : int gbin1;
276 7 : std::vector<std::size_t> gbin( dernames.size() );
277 7 : std::vector<bool> ipbc( dernames.size() );
278 7 : if( !flatgrid ) {
279 6 : ifile.scanField( "nbins", gbin1);
280 6 : gbin[0]=gbin1;
281 6 : std::vector<double> spacing;
282 6 : createGridAndValue( "fibonacci", ipbc, gbin[0], gmin, gmax, gbin, spacing );
283 6 : } else {
284 3 : for(unsigned i=0; i<dernames.size(); ++i) {
285 4 : ifile.scanField( "min_" + dernames[i], gmin[i]);
286 4 : ifile.scanField( "max_" + dernames[i], gmax[i]);
287 4 : ifile.scanField( "periodic_" + dernames[i], pstring );
288 2 : ifile.scanField( "nbins_" + dernames[i], gbin1);
289 2 : gbin[i]=gbin1;
290 2 : if( pstring=="true" ) {
291 1 : log.printf(" for periodic coordinate %s minimum is %s maximum is %s and number of bins is %zu \n",dernames[i].c_str(),gmin[i].c_str(),gmax[i].c_str(),gbin[i]);
292 : ipbc[i]=true;
293 1 : } else if( pstring=="false" ) {
294 1 : log.printf(" for coordinate %s minimum is %s maximum is %s and number of bins is %zu \n",dernames[i].c_str(),gmin[i].c_str(),gmax[i].c_str(),gbin[i]);
295 : ipbc[i]=false;
296 : } else {
297 0 : error("do not understand periodidicy of " + dernames[i] );
298 : }
299 :
300 4 : bool hasder=ifile.FieldExist( "d" + valuestr + "_" + dernames[i] );
301 2 : if( !hasder ) {
302 0 : plumed_merror("missing derivatives from grid file");
303 : }
304 : }
305 1 : std::vector<double> spacing;
306 1 : createGridAndValue( "flat", ipbc, 0, gmin, gmax, gbin, spacing );
307 1 : }
308 : // And finally read all the grid points
309 7 : Value* valout=getPntrToComponent(0);
310 7 : std::vector<double> dder( dernames.size() ), xx( dernames.size() );
311 2527 : for(unsigned i=0; i<valout->getNumberOfValues(); ++i) {
312 : double x, val;
313 2520 : ifile.scanField( valuestr, val );
314 9660 : for(unsigned j=0; j<dernames.size(); ++j) {
315 7140 : ifile.scanField(dernames[j],x);
316 7140 : if( !flatgrid ) {
317 12600 : ifile.scanField("nbins", gbin1);
318 : } else {
319 840 : xx[j]=x+gridobject.getGridSpacing()[j]/2.0;
320 1680 : ifile.scanField( "min_" + dernames[j], gmin[j]);
321 1680 : ifile.scanField( "max_" + dernames[j], gmax[j]);
322 1680 : ifile.scanField( "nbins_" + dernames[j], gbin1);
323 1680 : ifile.scanField( "periodic_" + dernames[j], pstring );
324 : }
325 : }
326 9660 : for(unsigned j=0; j<dernames.size(); ++j) {
327 21420 : ifile.scanField( "d" + valuestr + "_" + dernames[j], dder[j] );
328 : }
329 :
330 2520 : unsigned index=gridobject.getIndex(xx);
331 2520 : if( !flatgrid ) {
332 2100 : index=i;
333 : }
334 2520 : valout->set( index, val );
335 9660 : for(unsigned j=0; j<dernames.size(); ++j) {
336 7140 : valout->addGridDerivatives( index, j, dder[j] );
337 : }
338 2520 : ifile.scanField();
339 : }
340 7 : ifile.close();
341 21 : }
342 29 : }
343 :
344 22 : void ReadGridInSetup::createGridAndValue( const std::string& gtype, const std::vector<bool>& ipbc, const unsigned& nfermi,
345 : const std::vector<std::string>& gmin, const std::vector<std::string>& gmax,
346 : const std::vector<std::size_t>& gbin, std::vector<double>& gspacing ) {
347 22 : gridobject.setup( gtype, ipbc, nfermi, 0.0 );
348 22 : if( gtype=="flat" ) {
349 16 : gridobject.setBounds( gmin, gmax, gbin, gspacing );
350 : // Now create the value
351 15 : std::vector<std::size_t> shape( gridobject.getNbin(true) );
352 15 : ActionWithValue::addValueWithDerivatives( shape );
353 15 : setNotPeriodic();
354 15 : } else {
355 6 : std::vector<std::size_t> shape( 3 );
356 6 : shape[0]=gbin[0];
357 6 : shape[1]=shape[2]=1;
358 6 : ActionWithValue::addValueWithDerivatives( shape );
359 6 : setNotPeriodic();
360 6 : }
361 42 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
362 21 : getPntrToComponent(i)->setConstant();
363 21 : getPntrToComponent(i)->setDerivativeIsZeroWhenValueIsZero();
364 : }
365 : // This ensures we set the flag to never active the action. We can say we have atoms here as we don't need them
366 : // to calculate the CV
367 21 : setupConstantValues( true );
368 21 : }
369 :
370 80 : unsigned ReadGridInSetup::getNumberOfDerivatives() {
371 80 : return dernames.size();
372 : }
373 :
374 6 : std::vector<std::string> ReadGridInSetup::getGridCoordinateNames() const {
375 6 : return dernames;
376 : }
377 :
378 18491 : const GridCoordinatesObject& ReadGridInSetup::getGridCoordinatesObject() const {
379 18491 : return gridobject;
380 : }
381 :
382 : }
383 : }
384 :
|
