Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "ActionRegister.h"
23 : #include "Function.h"
24 : #include "tools/Exception.h"
25 : #include "tools/Communicator.h"
26 : #include "tools/BiasRepresentation.h"
27 : #include "tools/KernelFunctions.h"
28 : #include "tools/File.h"
29 : #include "tools/Tools.h"
30 : #include "tools/Stopwatch.h"
31 : #include "tools/Grid.h"
32 : #include <iostream>
33 : #include <memory>
34 :
35 : using namespace std;
36 :
37 : namespace PLMD {
38 : namespace function {
39 :
40 :
41 : //+PLUMEDOC FUNCTION FUNCSUMHILLS
42 : /*
43 : This function is intended to be called by the command line tool sum_hills. It is meant to integrate a HILLS file or an HILLS file interpreted as a histogram in a variety of ways. It is, therefore, not expected that you use this during your dynamics (it will crash!)
44 :
45 : In the future one could implement periodic integration during the metadynamics
46 : or straightforward MD as a tool to check convergence
47 :
48 : \par Examples
49 :
50 : There are currently no examples for this keyword.
51 :
52 : */
53 : //+ENDPLUMEDOC
54 :
55 11 : class FilesHandler {
56 : vector <string> filenames;
57 : vector <std::unique_ptr<IFile>> ifiles;
58 : Action *action;
59 : Log *log;
60 : bool parallelread;
61 : unsigned beingread;
62 : bool isopen;
63 : public:
64 : FilesHandler(const vector<string> &filenames, const bool ¶llelread, Action &myaction, Log &mylog);
65 : bool readBunch(BiasRepresentation *br, int stride);
66 : bool scanOneHill(BiasRepresentation *br, IFile *ifile );
67 : void getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin);
68 : void getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin, vector<double> &histosigma);
69 : };
70 22 : FilesHandler::FilesHandler(const vector<string> &filenames, const bool ¶llelread, Action &action, Log &mylog ):filenames(filenames),log(&mylog),parallelread(parallelread),beingread(0),isopen(false) {
71 11 : this->action=&action;
72 46 : for(unsigned i=0; i<filenames.size(); i++) {
73 12 : std::unique_ptr<IFile> ifile(new IFile());
74 12 : ifile->link(action);
75 12 : plumed_massert((ifile->FileExist(filenames[i])), "the file "+filenames[i]+" does not exist " );
76 12 : ifiles.emplace_back(std::move(ifile));
77 : }
78 :
79 11 : }
80 :
81 : // note that the FileHandler is completely transparent respect to the biasrepresentation
82 : // no check are made at this level
83 12 : bool FilesHandler::readBunch(BiasRepresentation *br, int stride = -1) {
84 : bool morefiles; morefiles=true;
85 12 : if(parallelread) {
86 0 : (*log)<<" doing parallelread \n";
87 0 : plumed_merror("parallelread is not yet implemented !!!");
88 : } else {
89 12 : (*log)<<" doing serialread \n";
90 : // read one by one hills
91 : // is the type defined? if not, assume it is a gaussian
92 : IFile *ff;
93 12 : ff=ifiles[beingread].get();
94 12 : if(!isopen) {
95 11 : (*log)<<" opening file "<<filenames[beingread]<<"\n";
96 22 : ff->open(filenames[beingread]); isopen=true;
97 : }
98 : int n;
99 : while(true) {
100 : bool fileisover=true;
101 5579 : while(scanOneHill(br,ff)) {
102 : // here do the dump if needed
103 5554 : n=br->getNumberOfKernels();
104 5554 : if(stride>0 && n%stride==0 && n!=0 ) {
105 1 : (*log)<<" done with this chunk: now with "<<n<<" kernels \n";
106 : fileisover=false;
107 : break;
108 : }
109 : }
110 13 : if(fileisover) {
111 24 : (*log)<<" closing file "<<filenames[beingread]<<"\n";
112 12 : ff->close();
113 12 : isopen=false;
114 12 : (*log)<<" now total "<<br->getNumberOfKernels()<<" kernels \n";
115 12 : beingread++;
116 24 : if(beingread<ifiles.size()) {
117 2 : ff=ifiles[beingread].get(); ff->open(filenames[beingread]);
118 2 : (*log)<<" opening file "<<filenames[beingread]<<"\n";
119 1 : isopen=true;
120 : } else {
121 : morefiles=false;
122 11 : (*log)<<" final chunk: now with "<<n<<" kernels \n";
123 : break;
124 : }
125 : }
126 : // if there are no more files to read and this file is over then quit
127 : if(fileisover && !morefiles) {break;}
128 : // if you are in the middle of a file and you are here
129 : // then means that you read what you need to read
130 2 : if(!fileisover ) {break;}
131 : }
132 : }
133 12 : return morefiles;
134 : }
135 3 : void FilesHandler::getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin) {
136 : // create the representation (no grid)
137 3 : BiasRepresentation br(vals,cc);
138 : // read all the kernels
139 3 : readBunch(&br);
140 : // loop over the kernels and get the support
141 3 : br.getMinMaxBin(vmin,vmax,vbin);
142 3 : }
143 1 : void FilesHandler::getMinMaxBin(vector<Value*> vals, Communicator &cc, vector<double> &vmin, vector<double> &vmax, vector<unsigned> &vbin, vector<double> &histosigma) {
144 1 : BiasRepresentation br(vals,cc,histosigma);
145 : // read all the kernels
146 1 : readBunch(&br);
147 : // loop over the kernels and get the support
148 1 : br.getMinMaxBin(vmin,vmax,vbin);
149 : //for(unsigned i=0;i<vals.size();i++){cerr<<"XXX "<<vmin[i]<<" "<<vmax[i]<<" "<<vbin[i]<<"\n";}
150 1 : }
151 5566 : bool FilesHandler::scanOneHill(BiasRepresentation *br, IFile *ifile ) {
152 : double dummy;
153 11132 : if(ifile->scanField("time",dummy)) {
154 : //(*log)<<" scanning one hill: "<<dummy<<" \n";
155 16662 : if(ifile->FieldExist("biasf")) ifile->scanField("biasf",dummy);
156 11108 : if(ifile->FieldExist("clock")) ifile->scanField("clock",dummy);
157 : // keep this intermediate function in case you need to parse more data in the future
158 5554 : br->pushKernel(ifile);
159 : //(*log)<<" read hill\n";
160 5554 : if(br->hasSigmaInInput())ifile->allowIgnoredFields();
161 5554 : ifile->scanField();
162 : return true;
163 : } else {
164 : return false;
165 : }
166 : }
167 :
168 :
169 900 : double mylog( double v1 ) {
170 900 : return log(v1);
171 : }
172 :
173 1800 : double mylogder( double v1 ) {
174 1800 : return 1./v1;
175 : }
176 :
177 :
178 :
179 21 : class FuncSumHills :
180 : public Function
181 : {
182 : vector<string> hillsFiles,histoFiles;
183 : vector<string> proj;
184 : int initstride;
185 : bool iscltool,integratehills,integratehisto,parallelread;
186 : bool negativebias;
187 : bool nohistory;
188 : bool minTOzero;
189 : bool doInt;
190 : double lowI_;
191 : double uppI_;
192 : double beta;
193 : string outhills,outhisto,fmt;
194 : std::unique_ptr<BiasRepresentation> biasrep;
195 : std::unique_ptr<BiasRepresentation> historep;
196 : public:
197 : explicit FuncSumHills(const ActionOptions&);
198 : void calculate() override; // this probably is not needed
199 : bool checkFilesAreExisting(const vector<string> & hills );
200 : static void registerKeywords(Keywords& keys);
201 : };
202 :
203 7846 : PLUMED_REGISTER_ACTION(FuncSumHills,"FUNCSUMHILLS")
204 :
205 8 : void FuncSumHills::registerKeywords(Keywords& keys) {
206 8 : Function::registerKeywords(keys);
207 16 : keys.use("ARG");
208 32 : keys.add("optional","HILLSFILES"," source file for hills creation(may be the same as HILLS)"); // this can be a vector!
209 32 : keys.add("optional","HISTOFILES"," source file for histogram creation(may be the same as HILLS)"); // also this can be a vector!
210 32 : keys.add("optional","HISTOSIGMA"," sigmas for binning when the histogram correction is needed ");
211 32 : keys.add("optional","PROJ"," only with sumhills: the projection on the CVs");
212 32 : keys.add("optional","KT"," only with sumhills: the kt factor when projection on CVs");
213 32 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
214 32 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
215 32 : keys.add("optional","GRID_BIN","the number of bins for the grid");
216 32 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
217 32 : keys.add("optional","INTERVAL","set one dimensional INTERVAL");
218 32 : keys.add("optional","OUTHILLS"," output file for hills ");
219 32 : keys.add("optional","OUTHISTO"," output file for histogram ");
220 32 : keys.add("optional","INITSTRIDE"," stride if you want an initial dump ");
221 32 : keys.add("optional","STRIDE"," stride when you do it on the fly ");
222 24 : keys.addFlag("ISCLTOOL",true,"use via plumed command line: calculate at read phase and then go");
223 24 : keys.addFlag("PARALLELREAD",false,"read parallel HILLS file");
224 24 : keys.addFlag("NEGBIAS",false,"dump negative bias ( -bias ) instead of the free energy: needed in well tempered with flexible hills ");
225 24 : keys.addFlag("NOHISTORY",false,"to be used with INITSTRIDE: it splits the bias/histogram in pieces without previous history ");
226 24 : keys.addFlag("MINTOZERO",false,"translate the resulting bias/histogram to have the minimum to zero ");
227 32 : keys.add("optional","FMT","the format that should be used to output real numbers");
228 8 : }
229 :
230 7 : FuncSumHills::FuncSumHills(const ActionOptions&ao):
231 : Action(ao),
232 : Function(ao),
233 : initstride(-1),
234 : iscltool(false),
235 : integratehills(false),
236 : integratehisto(false),
237 : parallelread(false),
238 : negativebias(false),
239 : nohistory(false),
240 : minTOzero(false),
241 : doInt(false),
242 : lowI_(-1.),
243 : uppI_(-1.),
244 : beta(-1.),
245 21 : fmt("%14.9f")
246 : {
247 :
248 : // format
249 14 : parse("FMT",fmt);
250 7 : log<<" Output format is "<<fmt<<"\n";
251 : // here read
252 : // Grid Stuff
253 : vector<std::string> gmin;
254 14 : parseVector("GRID_MIN",gmin);
255 7 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) error("not enough values for GRID_MIN");
256 7 : plumed_massert(gmin.size()==getNumberOfArguments() || gmin.size()==0,"need GRID_MIN argument for this") ;
257 0 : vector<std::string> gmax;
258 14 : parseVector("GRID_MAX",gmax);
259 7 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) error("not enough values for GRID_MAX");
260 7 : plumed_massert(gmax.size()==getNumberOfArguments() || gmax.size()==0,"need GRID_MAX argument for this") ;
261 : vector<unsigned> gbin;
262 : vector<double> gspacing;
263 14 : parseVector("GRID_BIN",gbin);
264 7 : plumed_massert(gbin.size()==getNumberOfArguments() || gbin.size()==0,"need GRID_BIN argument for this") ;
265 7 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
266 14 : parseVector("GRID_SPACING",gspacing);
267 7 : plumed_massert(gspacing.size()==getNumberOfArguments() || gspacing.size()==0,"need GRID_SPACING argument for this") ;
268 7 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) error("not enough values for GRID_SPACING");
269 7 : if(gspacing.size()!=0 && gbin.size()==0) {
270 0 : log<<" The number of bins will be estimated from GRID_SPACING\n";
271 7 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
272 0 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
273 0 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
274 : }
275 7 : if(gspacing.size()!=0) for(unsigned i=0; i<getNumberOfArguments(); i++) {
276 0 : if(gbin.size()==0) gbin.assign(getNumberOfArguments(),1);
277 : double a,b;
278 0 : Tools::convert(gmin[i],a);
279 0 : Tools::convert(gmax[i],b);
280 0 : unsigned n=((b-a)/gspacing[i])+1;
281 0 : if(gbin[i]<n) gbin[i]=n;
282 : }
283 :
284 : // Inteval keyword
285 7 : vector<double> tmpI(2);
286 14 : parseVector("INTERVAL",tmpI);
287 7 : if(tmpI.size()!=2&&tmpI.size()!=0) error("both a lower and an upper limits must be provided with INTERVAL");
288 7 : else if(tmpI.size()==2) {
289 0 : lowI_=tmpI.at(0);
290 0 : uppI_=tmpI.at(1);
291 0 : if(getNumberOfArguments()!=1) error("INTERVAL limits correction works only for monodimensional metadynamics!");
292 0 : if(uppI_<lowI_) error("The Upper limit must be greater than the Lower limit!");
293 0 : doInt=true;
294 : }
295 7 : if(doInt) {
296 0 : log << " Upper and Lower limits boundaries for the bias are activated at " << lowI_ << " - " << uppI_<<"\n";
297 0 : log << " Using the same values as boundaries for the grid if not other value was defined (default: 200 bins)\n";
298 0 : std::ostringstream strsmin, strsmax;
299 0 : strsmin << lowI_;
300 0 : strsmax << uppI_;
301 0 : if(gmin.size()==0) gmin.push_back(strsmin.str());
302 0 : if(gmax.size()==0) gmax.push_back(strsmax.str());
303 0 : if(gbin.size()==0) gbin.push_back(200);
304 : }
305 :
306 :
307 : // hills file:
308 14 : parseVector("HILLSFILES",hillsFiles);
309 7 : if(hillsFiles.size()==0) {
310 1 : integratehills=false; // default behaviour
311 : } else {
312 6 : integratehills=true;
313 26 : for(unsigned i=0; i<hillsFiles.size(); i++) log<<" hillsfile : "<<hillsFiles[i]<<"\n";
314 : }
315 : // histo file:
316 14 : parseVector("HISTOFILES",histoFiles);
317 7 : if(histoFiles.size()==0) {
318 6 : integratehisto=false;
319 : } else {
320 1 : integratehisto=true;
321 4 : for(unsigned i=0; i<histoFiles.size(); i++) log<<" histofile : "<<histoFiles[i]<<"\n";
322 : }
323 : vector<double> histoSigma;
324 7 : if(integratehisto) {
325 2 : parseVector("HISTOSIGMA",histoSigma);
326 6 : for(unsigned i=0; i<histoSigma.size(); i++) log<<" histosigma : "<<histoSigma[i]<<"\n";
327 : }
328 :
329 : // needs a projection?
330 7 : proj.clear();
331 14 : parseVector("PROJ",proj);
332 7 : if(integratehills) {
333 6 : plumed_massert(proj.size()<getNumberOfArguments()," The number of projection must be less than the full list of arguments ");
334 : }
335 7 : if(integratehisto) {
336 1 : plumed_massert(proj.size()<=getNumberOfArguments()," The number of projection must be less or equal to the full list of arguments ");
337 : }
338 8 : if(integratehisto&&proj.size()==0) {
339 5 : for(unsigned i=0; i<getNumberOfArguments(); i++) proj.push_back(getPntrToArgument(i)->getName());
340 : }
341 :
342 : // add some automatic hills width: not in case stride is defined
343 : // since when you start from zero the automatic size will be zero!
344 10 : if(gmin.size()==0 || gmax.size()==0) {
345 4 : log<<" \n";
346 4 : log<<" No boundaries defined: need to do a prescreening of hills \n";
347 : std::vector<Value*> tmphillsvalues, tmphistovalues;
348 4 : if(integratehills) {
349 21 : for(unsigned i=0; i<getNumberOfArguments(); i++)tmphillsvalues.push_back( getPntrToArgument(i) );
350 : }
351 4 : if(integratehisto) {
352 5 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
353 2 : std::string ss = getPntrToArgument(i)->getName();
354 10 : for(unsigned j=0; j<proj.size(); j++) {
355 8 : if(proj[j]==ss) tmphistovalues.push_back( getPntrToArgument(i) );
356 : }
357 : }
358 : }
359 :
360 4 : if(integratehills) {
361 3 : FilesHandler hillsHandler(hillsFiles,parallelread,*this, log);
362 : vector<double> vmin,vmax;
363 : vector<unsigned> vbin;
364 6 : hillsHandler.getMinMaxBin(tmphillsvalues,comm,vmin,vmax,vbin);
365 3 : log<<" found boundaries from hillsfile: \n";
366 3 : gmin.resize(vmin.size());
367 3 : gmax.resize(vmax.size());
368 3 : if(gbin.size()==0) {
369 2 : gbin=vbin;
370 : } else {
371 1 : log<<" found nbins in input, this overrides the automatic choice \n";
372 : }
373 15 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
374 12 : Tools::convert(vmin[i],gmin[i]);
375 6 : Tools::convert(vmax[i],gmax[i]);
376 6 : log<<" variable "<< getPntrToArgument(i)->getName()<<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
377 : }
378 : }
379 : // if at this stage bins are not there then do it with histo
380 4 : if(gmin.size()==0) {
381 1 : FilesHandler histoHandler(histoFiles,parallelread,*this, log);
382 : vector<double> vmin,vmax;
383 : vector<unsigned> vbin;
384 2 : histoHandler.getMinMaxBin(tmphistovalues,comm,vmin,vmax,vbin,histoSigma);
385 1 : log<<" found boundaries from histofile: \n";
386 1 : gmin.resize(vmin.size());
387 1 : gmax.resize(vmax.size());
388 1 : if(gbin.size()==0) {
389 0 : gbin=vbin;
390 : } else {
391 1 : log<<" found nbins in input, this overrides the automatic choice \n";
392 : }
393 5 : for(unsigned i=0; i<proj.size(); i++) {
394 2 : Tools::convert(vmin[i],gmin[i]);
395 2 : Tools::convert(vmax[i],gmax[i]);
396 2 : log<<" variable "<< proj[i] <<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
397 : }
398 : }
399 4 : log<<" done!\n";
400 4 : log<<" \n";
401 : }
402 :
403 :
404 7 : if( proj.size() != 0 || integratehisto==true ) {
405 4 : parse("KT",beta);
406 10 : for(unsigned i=0; i<proj.size(); i++) log<<" projection "<<i<<" : "<<proj[i]<<"\n";
407 : // this should be only for projection or free energy from histograms
408 2 : plumed_massert(beta>0.,"if you make a projection or a histogram correction then you need KT flag!");
409 2 : beta=1./beta;
410 2 : log<<" beta is "<<beta<<"\n";
411 : }
412 : // is a cltool: then you start and then die
413 14 : parseFlag("ISCLTOOL",iscltool);
414 : //
415 14 : parseFlag("NEGBIAS",negativebias);
416 : //
417 14 : parseFlag("PARALLELREAD",parallelread);
418 : // stride
419 14 : parse("INITSTRIDE",initstride);
420 : // output suffix or names
421 7 : if(initstride<0) {
422 6 : log<<" Doing only one integration: no stride \n";
423 12 : outhills="fes.dat"; outhisto="histo.dat";
424 : }
425 : else {
426 2 : outhills="fes_"; outhisto="histo_";
427 1 : log<<" Doing integration slices every "<<initstride<<" kernels\n";
428 2 : parseFlag("NOHISTORY",nohistory);
429 1 : if(nohistory)log<<" nohistory: each stride block has no memory of the previous block\n";
430 : }
431 14 : parseFlag("MINTOZERO",minTOzero);
432 7 : if(minTOzero)log<<" mintozero: bias/histogram will be translated to have the minimum value equal to zero\n";
433 : //what might it be this?
434 : // here start
435 : // want something right now?? do it and return
436 : // your argument is a set of cvs
437 : // then you need: a hills / a colvar-like file (to do a histogram)
438 : // create a bias representation for this
439 7 : if(iscltool) {
440 :
441 : std::vector<Value*> tmphillsvalues, tmphistovalues;
442 7 : if(integratehills) {
443 30 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
444 : // allocate a new value from the old one: no deriv here
445 : // if we are summing hills then all the arguments are needed
446 24 : tmphillsvalues.push_back( getPntrToArgument(i) );
447 : }
448 : }
449 7 : if(integratehisto) {
450 5 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
451 2 : std::string ss = getPntrToArgument(i)->getName();
452 10 : for(unsigned j=0; j<proj.size(); j++) {
453 8 : if(proj[j]==ss) tmphistovalues.push_back( getPntrToArgument(i) );
454 : }
455 : }
456 : }
457 :
458 : // check if the files exists
459 7 : if(integratehills) {
460 6 : checkFilesAreExisting(hillsFiles);
461 6 : biasrep.reset(new BiasRepresentation(tmphillsvalues,comm, gmin, gmax, gbin, doInt, lowI_, uppI_));
462 6 : if(negativebias) {
463 1 : biasrep->setRescaledToBias(true);
464 1 : log<<" required the -bias instead of the free energy \n";
465 1 : if(initstride<0) {outhills="negativebias.dat";}
466 0 : else {outhills="negativebias_";}
467 : }
468 : }
469 :
470 14 : parse("OUTHILLS",outhills);
471 14 : parse("OUTHISTO",outhisto);
472 7 : if(integratehills)log<<" output file for fes/bias is : "<<outhills<<"\n";
473 7 : if(integratehisto)log<<" output file for histogram is : "<<outhisto<<"\n";
474 7 : checkRead();
475 :
476 7 : log<<"\n";
477 7 : log<<" Now calculating...\n";
478 7 : log<<"\n";
479 :
480 : // here it defines the column to be histogrammed, tmpvalues should be only
481 : // the list of the collective variable one want to consider
482 7 : if(integratehisto) {
483 1 : checkFilesAreExisting(histoFiles);
484 1 : historep.reset(new BiasRepresentation(tmphistovalues,comm,gmin,gmax,gbin,histoSigma));
485 : }
486 :
487 : // decide how to source hills ( serial/parallel )
488 : // here below the input control
489 : // say how many hills and it will read them from the
490 : // bunch of files provided, will update the representation
491 : // of hills (i.e. a list of hills and the associated grid)
492 :
493 : // decide how to source colvars ( serial parallel )
494 7 : std::unique_ptr<FilesHandler> hillsHandler;
495 7 : std::unique_ptr<FilesHandler> histoHandler;
496 :
497 7 : if(integratehills) hillsHandler.reset(new FilesHandler(hillsFiles,parallelread,*this, log));
498 7 : if(integratehisto) histoHandler.reset(new FilesHandler(histoFiles,parallelread,*this, log));
499 :
500 : // Stopwatch is logged when it goes out of scope
501 14 : Stopwatch sw(log);
502 :
503 : // Stopwatch is stopped when swh goes out of scope
504 21 : auto swh=sw.startStop("0 Summing hills");
505 :
506 : // read a number of hills and put in the bias representation
507 : int nfiles=0;
508 7 : bool ibias=integratehills; bool ihisto=integratehisto;
509 : while(true) {
510 8 : if( integratehills && ibias ) {
511 7 : if(nohistory) {biasrep->clear(); log<<" clearing history before reading a new block\n";};
512 7 : log<<" reading hills: \n";
513 14 : ibias=hillsHandler->readBunch(biasrep.get(),initstride) ; log<<"\n";
514 : }
515 :
516 8 : if( integratehisto && ihisto ) {
517 1 : if(nohistory) {historep->clear(); log<<" clearing history before reading a new block\n";};
518 1 : log<<" reading histogram: \n";
519 2 : ihisto=histoHandler->readBunch(historep.get(),initstride) ; log<<"\n";
520 : }
521 :
522 : // dump: need to project?
523 8 : if(proj.size()!=0) {
524 :
525 3 : if(integratehills) {
526 :
527 2 : log<<" Bias: Projecting on subgrid... \n";
528 2 : BiasWeight Bw(beta);
529 4 : Grid biasGrid=*(biasrep->getGridPtr());
530 4 : Grid smallGrid=biasGrid.project(proj,&Bw);
531 4 : OFile gridfile; gridfile.link(*this);
532 4 : std::ostringstream ostr; ostr<<nfiles;
533 : string myout;
534 10 : if(initstride>0) { myout=outhills+ostr.str()+".dat" ;} else {myout=outhills;}
535 2 : log<<" Bias: Writing subgrid on file "<<myout<<" \n";
536 2 : gridfile.open(myout);
537 2 : if(minTOzero) smallGrid.setMinToZero();
538 : smallGrid.setOutputFmt(fmt);
539 2 : smallGrid.writeToFile(gridfile);
540 2 : gridfile.close();
541 2 : if(!ibias)integratehills=false;// once you get to the final bunch just give up
542 : }
543 : // this should be removed
544 3 : if(integratehisto) {
545 :
546 1 : log<<" Histo: Projecting on subgrid... \n";
547 1 : Grid histoGrid=*(historep->getGridPtr());
548 :
549 2 : OFile gridfile; gridfile.link(*this);
550 2 : std::ostringstream ostr; ostr<<nfiles;
551 : string myout;
552 1 : if(initstride>0) { myout=outhisto+ostr.str()+".dat" ;} else {myout=outhisto;}
553 1 : log<<" Histo: Writing subgrid on file "<<myout<<" \n";
554 1 : gridfile.open(myout);
555 :
556 1 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
557 1 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
558 1 : if(minTOzero) histoGrid.setMinToZero();
559 : histoGrid.setOutputFmt(fmt);
560 1 : histoGrid.writeToFile(gridfile);
561 :
562 2 : if(!ihisto)integratehisto=false;// once you get to the final bunch just give up
563 : }
564 :
565 : } else {
566 :
567 5 : if(integratehills) {
568 :
569 5 : Grid biasGrid=*(biasrep->getGridPtr());
570 5 : biasGrid.scaleAllValuesAndDerivatives(-1.);
571 :
572 10 : OFile gridfile; gridfile.link(*this);
573 10 : std::ostringstream ostr; ostr<<nfiles;
574 : string myout;
575 5 : if(initstride>0) { myout=outhills+ostr.str()+".dat" ;} else {myout=outhills;}
576 5 : log<<" Writing full grid on file "<<myout<<" \n";
577 5 : gridfile.open(myout);
578 :
579 5 : if(minTOzero) biasGrid.setMinToZero();
580 : biasGrid.setOutputFmt(fmt);
581 5 : biasGrid.writeToFile(gridfile);
582 : // rescale back prior to accumulate
583 10 : if(!ibias)integratehills=false;// once you get to the final bunch just give up
584 : }
585 5 : if(integratehisto) {
586 :
587 0 : Grid histoGrid=*(historep->getGridPtr());
588 : // do this if you want a free energy from a grid, otherwise do not
589 0 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
590 0 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
591 :
592 0 : OFile gridfile; gridfile.link(*this);
593 0 : std::ostringstream ostr; ostr<<nfiles;
594 : string myout;
595 0 : if(initstride>0) { myout=outhisto+ostr.str()+".dat" ;} else {myout=outhisto;}
596 0 : log<<" Writing full grid on file "<<myout<<" \n";
597 0 : gridfile.open(myout);
598 :
599 : // also this is usefull only for free energy
600 0 : if(minTOzero) histoGrid.setMinToZero();
601 : histoGrid.setOutputFmt(fmt);
602 0 : histoGrid.writeToFile(gridfile);
603 :
604 0 : if(!ihisto)integratehisto=false; // once you get to the final bunch just give up
605 : }
606 : }
607 8 : if ( !ibias && !ihisto) break; //when both are over then just quit
608 :
609 1 : nfiles++;
610 : }
611 :
612 7 : return;
613 0 : }
614 : // just an initialization but you need to do something on the fly?: need to connect with a metad run and its grid representation
615 : // your argument is a metad run
616 : // if the grid does not exist crash and say that you need some data
617 : // otherwise just link with it
618 :
619 : }
620 :
621 0 : void FuncSumHills::calculate() {
622 : // this should be connected only with a grid representation to metadynamics
623 : // at regular time just dump it
624 0 : plumed_merror("You should have never got here: this stuff is not yet implemented!");
625 : }
626 :
627 7 : bool FuncSumHills::checkFilesAreExisting(const vector<string> & hills ) {
628 7 : plumed_massert(hills.size()!=0,"the number of files provided should be at least one" );
629 7 : std::unique_ptr<IFile> ifile(new IFile());
630 7 : ifile->link(*this);
631 23 : for(unsigned i=0; i< hills.size(); i++) {
632 8 : plumed_massert(ifile->FileExist(hills[i]),"missing file "+hills[i]);
633 : }
634 7 : return true;
635 :
636 : }
637 :
638 : }
639 :
640 5874 : }
641 :
642 :
|
