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