Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-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 "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 11 : 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 23 : 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 12 : plumed_massert((ifile->FileExist(filenames[i])), "the file "+filenames[i]+" does not exist " );
78 : }
79 :
80 11 : }
81 :
82 22 : FilesHandler::~FilesHandler() {
83 11 : 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 12 : 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 12 : ff=ifiles[beingread];
99 12 : if(!isopen) {
100 11 : (*log)<<" opening file "<<filenames[beingread]<<"\n";
101 11 : ff->open(filenames[beingread]); isopen=true;
102 : }
103 : int n;
104 : while(true) {
105 13 : bool fileisover=true;
106 5579 : 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 1 : fileisover=false;
112 1 : break;
113 : }
114 : }
115 13 : if(fileisover) {
116 12 : (*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 12 : if(beingread<ifiles.size()) {
122 1 : ff=ifiles[beingread]; ff->open(filenames[beingread]);
123 1 : (*log)<<" opening file "<<filenames[beingread]<<"\n";
124 1 : isopen=true;
125 : } else {
126 11 : morefiles=false;
127 11 : (*log)<<" final chunk: now with "<<n<<" kernels \n";
128 11 : break;
129 : }
130 : }
131 : // if there are no more files to read and this file is over then quit
132 2 : 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 1 : }
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 3 : 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 1 : 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 5566 : if(ifile->scanField("time",dummy)) {
159 : //(*log)<<" scanning one hill: "<<dummy<<" \n";
160 5554 : if(ifile->FieldExist("biasf")) ifile->scanField("biasf",dummy);
161 5554 : 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 5554 : return true;
168 : } else {
169 12 : 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 2530 : PLUMED_REGISTER_ACTION(FuncSumHills,"FUNCSUMHILLS")
210 :
211 8 : void FuncSumHills::registerKeywords(Keywords& keys) {
212 8 : Function::registerKeywords(keys);
213 8 : keys.use("ARG");
214 8 : keys.add("optional","HILLSFILES"," source file for hills creation(may be the same as HILLS)"); // this can be a vector!
215 8 : keys.add("optional","HISTOFILES"," source file for histogram creation(may be the same as HILLS)"); // also this can be a vector!
216 8 : keys.add("optional","HISTOSIGMA"," sigmas for binning when the histogram correction is needed ");
217 8 : keys.add("optional","PROJ"," only with sumhills: the projection on the cvs");
218 8 : keys.add("optional","KT"," only with sumhills: the kt factor when projection on cvs");
219 8 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
220 8 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
221 8 : keys.add("optional","GRID_BIN","the number of bins for the grid");
222 8 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
223 8 : keys.add("optional","INTERVAL","set monodimensional INTERVAL");
224 8 : keys.add("optional","OUTHILLS"," output file for hills ");
225 8 : keys.add("optional","OUTHISTO"," output file for histogram ");
226 8 : keys.add("optional","INITSTRIDE"," stride if you want an initial dump ");
227 8 : keys.add("optional","STRIDE"," stride when you do it on the fly ");
228 8 : keys.addFlag("ISCLTOOL",true,"use via plumed commandline: calculate at read phase and then go");
229 8 : keys.addFlag("PARALLELREAD",false,"read parallel HILLS file");
230 8 : keys.addFlag("NEGBIAS",false,"dump negative bias ( -bias ) instead of the free energy: needed in welltempered with flexible hills ");
231 8 : keys.addFlag("NOHISTORY",false,"to be used with INITSTRIDE: it splits the bias/histogram in pieces without previous history ");
232 8 : keys.addFlag("MINTOZERO",false,"translate the resulting bias/histogram to have the minimum to zero ");
233 8 : 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 7 : historep(NULL)
254 : {
255 :
256 : // format
257 7 : parse("FMT",fmt);
258 7 : log<<" Output format is "<<fmt<<"\n";
259 : // here read
260 : // Grid Stuff
261 7 : vector<std::string> gmin;
262 7 : 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 7 : vector<std::string> gmax;
266 7 : 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 7 : vector<unsigned> gbin;
270 7 : vector<double> gspacing;
271 7 : parseVector("GRID_BIN",gbin);
272 7 : plumed_massert(gbin.size()==getNumberOfArguments() || gbin.size()==0,"need GRID_BIN argument for this") ;
273 7 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) error("not enough values for GRID_BIN");
274 7 : parseVector("GRID_SPACING",gspacing);
275 7 : plumed_massert(gspacing.size()==getNumberOfArguments() || gspacing.size()==0,"need GRID_SPACING argument for this") ;
276 7 : 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 7 : 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 7 : parseVector("HILLSFILES",hillsFiles);
317 7 : if(hillsFiles.size()==0) {
318 1 : integratehills=false; // default behaviour
319 : } else {
320 6 : integratehills=true;
321 6 : for(unsigned i=0; i<hillsFiles.size(); i++) log<<" hillsfile : "<<hillsFiles[i]<<"\n";
322 : }
323 : // histo file:
324 7 : parseVector("HISTOFILES",histoFiles);
325 7 : if(histoFiles.size()==0) {
326 6 : integratehisto=false;
327 : } else {
328 1 : integratehisto=true;
329 1 : for(unsigned i=0; i<histoFiles.size(); i++) log<<" histofile : "<<histoFiles[i]<<"\n";
330 : }
331 7 : vector<double> histoSigma;
332 7 : if(integratehisto) {
333 1 : parseVector("HISTOSIGMA",histoSigma);
334 1 : for(unsigned i=0; i<histoSigma.size(); i++) log<<" histosigma : "<<histoSigma[i]<<"\n";
335 : }
336 :
337 : // needs a projection?
338 7 : proj.clear();
339 7 : 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 7 : if(integratehisto&&proj.size()==0) {
347 1 : 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 7 : 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 8 : std::vector<Value*> tmphillsvalues, tmphistovalues;
356 4 : if(integratehills) {
357 3 : for(unsigned i=0; i<getNumberOfArguments(); i++)tmphillsvalues.push_back( getPntrToArgument(i) );
358 : }
359 4 : if(integratehisto) {
360 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
361 2 : std::string ss = getPntrToArgument(i)->getName();
362 6 : for(unsigned j=0; j<proj.size(); j++) {
363 4 : if(proj[j]==ss) tmphistovalues.push_back( getPntrToArgument(i) );
364 : }
365 2 : }
366 : }
367 :
368 4 : if(integratehills) {
369 : FilesHandler *hillsHandler;
370 3 : hillsHandler=new FilesHandler(hillsFiles,parallelread,*this, log);
371 6 : vector<double> vmin,vmax;
372 6 : vector<unsigned> vbin;
373 3 : hillsHandler->getMinMaxBin(tmphillsvalues,comm,vmin,vmax,vbin);
374 3 : log<<" found boundaries from hillsfile: \n";
375 3 : gmin.resize(vmin.size());
376 3 : gmax.resize(vmax.size());
377 3 : if(gbin.size()==0) {
378 2 : gbin=vbin;
379 : } else {
380 1 : log<<" found nbins in input, this overrides the automatic choice \n";
381 : }
382 9 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
383 6 : Tools::convert(vmin[i],gmin[i]);
384 6 : Tools::convert(vmax[i],gmax[i]);
385 6 : log<<" variable "<< getPntrToArgument(i)->getName()<<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
386 : }
387 6 : delete hillsHandler;
388 : }
389 : // if at this stage bins are not there then do it with histo
390 4 : if(gmin.size()==0) {
391 : FilesHandler *histoHandler;
392 1 : histoHandler=new FilesHandler(histoFiles,parallelread,*this, log);
393 2 : vector<double> vmin,vmax;
394 2 : vector<unsigned> vbin;
395 1 : histoHandler->getMinMaxBin(tmphistovalues,comm,vmin,vmax,vbin,histoSigma);
396 1 : log<<" found boundaries from histofile: \n";
397 1 : gmin.resize(vmin.size());
398 1 : gmax.resize(vmax.size());
399 1 : if(gbin.size()==0) {
400 0 : gbin=vbin;
401 : } else {
402 1 : log<<" found nbins in input, this overrides the automatic choice \n";
403 : }
404 3 : for(unsigned i=0; i<proj.size(); i++) {
405 2 : Tools::convert(vmin[i],gmin[i]);
406 2 : Tools::convert(vmax[i],gmax[i]);
407 2 : log<<" variable "<< proj[i] <<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
408 : }
409 2 : delete histoHandler;
410 : }
411 4 : log<<" done!\n";
412 8 : log<<" \n";
413 : }
414 :
415 :
416 7 : if( proj.size() != 0 || integratehisto==true ) {
417 2 : parse("KT",beta);
418 2 : for(unsigned i=0; i<proj.size(); i++) log<<" projection "<<i<<" : "<<proj[i]<<"\n";
419 : // this should be only for projection or free energy from histograms
420 2 : plumed_massert(beta>0.,"if you make a projection or a histogram correction then you need KT flag!");
421 2 : beta=1./beta;
422 2 : log<<" beta is "<<beta<<"\n";
423 : }
424 : // is a cltool: then you start and then die
425 7 : parseFlag("ISCLTOOL",iscltool);
426 : //
427 7 : parseFlag("NEGBIAS",negativebias);
428 : //
429 7 : parseFlag("PARALLELREAD",parallelread);
430 : // stride
431 7 : parse("INITSTRIDE",initstride);
432 : // output suffix or names
433 7 : if(initstride<0) {
434 6 : log<<" Doing only one integration: no stride \n";
435 6 : outhills="fes.dat"; outhisto="histo.dat";
436 : }
437 : else {
438 1 : outhills="fes_"; outhisto="histo_";
439 1 : log<<" Doing integration slices every "<<initstride<<" kernels\n";
440 1 : parseFlag("NOHISTORY",nohistory);
441 1 : if(nohistory)log<<" nohistory: each stride block has no memory of the previous block\n";
442 : }
443 7 : parseFlag("MINTOZERO",minTOzero);
444 7 : if(minTOzero)log<<" mintozero: bias/histogram will be translated to have the minimum value equal to zero\n";
445 : //what might it be this?
446 : // here start
447 : // want something right now?? do it and return
448 : // your argument is a set of cvs
449 : // then you need: a hills / a colvar-like file (to do a histogram)
450 : // create a bias representation for this
451 7 : if(iscltool) {
452 :
453 14 : std::vector<Value*> tmphillsvalues, tmphistovalues;
454 7 : if(integratehills) {
455 18 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
456 : // allocate a new value from the old one: no deriv here
457 : // if we are summing hills then all the arguments are needed
458 12 : tmphillsvalues.push_back( getPntrToArgument(i) );
459 : }
460 : }
461 7 : if(integratehisto) {
462 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
463 2 : std::string ss = getPntrToArgument(i)->getName();
464 6 : for(unsigned j=0; j<proj.size(); j++) {
465 4 : if(proj[j]==ss) tmphistovalues.push_back( getPntrToArgument(i) );
466 : }
467 2 : }
468 : }
469 :
470 : // check if the files exists
471 7 : if(integratehills) {
472 6 : checkFilesAreExisting(hillsFiles);
473 6 : biasrep=new BiasRepresentation(tmphillsvalues,comm, gmin, gmax, gbin, doInt, lowI_, uppI_);
474 6 : if(negativebias) {
475 1 : biasrep->setRescaledToBias(true);
476 1 : log<<" required the -bias instead of the free energy \n";
477 1 : if(initstride<0) {outhills="negativebias.dat";}
478 0 : else {outhills="negativebias_";}
479 : }
480 : }
481 :
482 7 : parse("OUTHILLS",outhills);
483 7 : parse("OUTHISTO",outhisto);
484 7 : if(integratehills)log<<" output file for fes/bias is : "<<outhills<<"\n";
485 7 : if(integratehisto)log<<" output file for histogram is : "<<outhisto<<"\n";
486 7 : checkRead();
487 :
488 7 : log<<"\n";
489 7 : log<<" Now calculating...\n";
490 7 : log<<"\n";
491 :
492 : // here it defines the column to be histogrammed, tmpvalues should be only
493 : // the list of the collective variable one want to consider
494 7 : if(integratehisto) {
495 1 : checkFilesAreExisting(histoFiles);
496 1 : historep=new BiasRepresentation(tmphistovalues,comm,gmin,gmax,gbin,histoSigma);
497 : }
498 :
499 : // decide how to source hills ( serial/parallel )
500 : // here below the input control
501 : // say how many hills and it will read them from the
502 : // bunch of files provided, will update the representation
503 : // of hills (i.e. a list of hills and the associated grid)
504 :
505 : // decide how to source colvars ( serial parallel )
506 : FilesHandler *hillsHandler;
507 : FilesHandler *histoHandler;
508 :
509 7 : hillsHandler=NULL;
510 7 : histoHandler=NULL;
511 :
512 7 : if(integratehills) hillsHandler=new FilesHandler(hillsFiles,parallelread,*this, log);
513 7 : if(integratehisto) histoHandler=new FilesHandler(histoFiles,parallelread,*this, log);
514 :
515 14 : Stopwatch sw;
516 :
517 7 : sw.start("0 Summing hills");
518 :
519 : // read a number of hills and put in the bias representation
520 7 : int nfiles=0;
521 7 : bool ibias=integratehills; bool ihisto=integratehisto;
522 : while(true) {
523 8 : if( integratehills && ibias ) {
524 7 : if(nohistory) {biasrep->clear(); log<<" clearing history before reading a new block\n";};
525 7 : log<<" reading hills: \n";
526 7 : ibias=hillsHandler->readBunch(biasrep,initstride) ; log<<"\n";
527 : }
528 :
529 8 : if( integratehisto && ihisto ) {
530 1 : if(nohistory) {historep->clear(); log<<" clearing history before reading a new block\n";};
531 1 : log<<" reading histogram: \n";
532 1 : ihisto=histoHandler->readBunch(historep,initstride) ; log<<"\n";
533 : }
534 :
535 : // dump: need to project?
536 8 : if(proj.size()!=0) {
537 :
538 3 : if(integratehills) {
539 :
540 2 : log<<" Bias: Projecting on subgrid... \n";
541 2 : BiasWeight *Bw=new BiasWeight(beta);
542 2 : Grid biasGrid=*(biasrep->getGridPtr());
543 4 : Grid smallGrid=biasGrid.project(proj,Bw);
544 4 : OFile gridfile; gridfile.link(*this);
545 4 : std::ostringstream ostr; ostr<<nfiles;
546 4 : string myout;
547 2 : if(initstride>0) { myout=outhills+ostr.str()+".dat" ;} else {myout=outhills;}
548 2 : log<<" Bias: Writing subgrid on file "<<myout<<" \n";
549 2 : gridfile.open(myout);
550 2 : if(minTOzero) smallGrid.setMinToZero();
551 2 : smallGrid.setOutputFmt(fmt);
552 2 : smallGrid.writeToFile(gridfile);
553 2 : gridfile.close();
554 2 : if(!ibias)integratehills=false;// once you get to the final bunch just give up
555 4 : delete Bw;
556 : }
557 : // this should be removed
558 3 : if(integratehisto) {
559 :
560 1 : log<<" Histo: Projecting on subgrid... \n";
561 1 : Grid histoGrid=*(historep->getGridPtr());
562 :
563 2 : OFile gridfile; gridfile.link(*this);
564 2 : std::ostringstream ostr; ostr<<nfiles;
565 2 : string myout;
566 1 : if(initstride>0) { myout=outhisto+ostr.str()+".dat" ;} else {myout=outhisto;}
567 1 : log<<" Histo: Writing subgrid on file "<<myout<<" \n";
568 1 : gridfile.open(myout);
569 :
570 1 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
571 1 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
572 1 : if(minTOzero) histoGrid.setMinToZero();
573 1 : histoGrid.setOutputFmt(fmt);
574 1 : histoGrid.writeToFile(gridfile);
575 1 : gridfile.close();
576 :
577 2 : if(!ihisto)integratehisto=false;// once you get to the final bunch just give up
578 : }
579 :
580 : } else {
581 :
582 5 : if(integratehills) {
583 :
584 5 : Grid biasGrid=*(biasrep->getGridPtr());
585 5 : biasGrid.scaleAllValuesAndDerivatives(-1.);
586 :
587 10 : OFile gridfile; gridfile.link(*this);
588 10 : std::ostringstream ostr; ostr<<nfiles;
589 10 : string myout;
590 5 : if(initstride>0) { myout=outhills+ostr.str()+".dat" ;} else {myout=outhills;}
591 5 : log<<" Writing full grid on file "<<myout<<" \n";
592 5 : gridfile.open(myout);
593 :
594 5 : if(minTOzero) biasGrid.setMinToZero();
595 5 : biasGrid.setOutputFmt(fmt);
596 5 : biasGrid.writeToFile(gridfile);
597 5 : gridfile.close();
598 : // rescale back prior to accumulate
599 10 : if(!ibias)integratehills=false;// once you get to the final bunch just give up
600 : }
601 5 : if(integratehisto) {
602 :
603 0 : Grid histoGrid=*(historep->getGridPtr());
604 : // do this if you want a free energy from a grid, otherwise do not
605 0 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
606 0 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
607 :
608 0 : OFile gridfile; gridfile.link(*this);
609 0 : std::ostringstream ostr; ostr<<nfiles;
610 0 : string myout;
611 0 : if(initstride>0) { myout=outhisto+ostr.str()+".dat" ;} else {myout=outhisto;}
612 0 : log<<" Writing full grid on file "<<myout<<" \n";
613 0 : gridfile.open(myout);
614 :
615 : // also this is usefull only for free energy
616 0 : if(minTOzero) histoGrid.setMinToZero();
617 0 : histoGrid.setOutputFmt(fmt);
618 0 : histoGrid.writeToFile(gridfile);
619 0 : gridfile.close();
620 :
621 0 : if(!ihisto)integratehisto=false; // once you get to the final bunch just give up
622 : }
623 : }
624 8 : if ( !ibias && !ihisto) break; //when both are over then just quit
625 :
626 1 : nfiles++;
627 : }
628 7 : if(hillsHandler) delete hillsHandler;
629 7 : if(histoHandler) delete histoHandler;
630 :
631 7 : sw.stop("0 Summing hills");
632 :
633 7 : log<<sw;
634 :
635 21 : return;
636 0 : }
637 : // just an initialization but you need to do something on the fly?: need to connect with a metad run and its grid representation
638 : // your argument is a metad run
639 : // if the grid does not exist crash and say that you need some data
640 : // otherwise just link with it
641 :
642 : }
643 :
644 0 : void FuncSumHills::calculate() {
645 : // this should be connected only with a grid representation to metadynamics
646 : // at regular time just dump it
647 0 : plumed_merror("You should have never got here: this stuff is not yet implemented!");
648 : }
649 :
650 28 : FuncSumHills::~FuncSumHills() {
651 7 : if(historep) delete historep;
652 7 : if(biasrep) delete biasrep;
653 21 : }
654 :
655 7 : bool FuncSumHills::checkFilesAreExisting(const vector<string> & hills ) {
656 7 : plumed_massert(hills.size()!=0,"the number of files provided should be at least one" );
657 7 : IFile *ifile = new IFile();
658 7 : ifile->link(*this);
659 15 : for(unsigned i=0; i< hills.size(); i++) {
660 8 : plumed_massert(ifile->FileExist(hills[i]),"missing file "+hills[i]);
661 : }
662 7 : delete ifile;
663 7 : return true;
664 :
665 : }
666 :
667 : }
668 :
669 2523 : }
670 :
671 :
|
