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