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 : ## 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 11 : void FuncSumHills::registerKeywords(Keywords& keys) {
214 11 : Function::registerKeywords(keys);
215 11 : keys.add("optional","HILLSFILES"," source file for hills creation(may be the same as HILLS)"); // this can be a vector!
216 11 : keys.add("optional","HISTOFILES"," source file for histogram creation(may be the same as HILLS)"); // also this can be a vector!
217 11 : keys.add("optional","HISTOSIGMA"," sigmas for binning when the histogram correction is needed ");
218 11 : keys.add("optional","PROJ"," only with sumhills: the projection on the CVs");
219 11 : keys.add("optional","KT"," only with sumhills: the kt factor when projection on CVs");
220 11 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
221 11 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
222 11 : keys.add("optional","GRID_BIN","the number of bins for the grid");
223 11 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
224 11 : keys.add("optional","INTERVAL","set one dimensional INTERVAL");
225 11 : keys.add("optional","OUTHILLS"," output file for hills ");
226 11 : keys.add("optional","OUTHISTO"," output file for histogram ");
227 11 : keys.add("optional","INITSTRIDE"," stride if you want an initial dump ");
228 11 : keys.add("optional","STRIDE"," stride when you do it on the fly ");
229 11 : keys.addFlag("ISCLTOOL",false,"use via plumed command line: calculate at read phase and then go");
230 11 : keys.addFlag("PARALLELREAD",false,"read parallel HILLS file");
231 11 : keys.addFlag("NEGBIAS",false,"dump negative bias ( -bias ) instead of the free energy: needed in well tempered with flexible hills ");
232 11 : keys.addFlag("NOHISTORY",false,"to be used with INITSTRIDE: it splits the bias/histogram in pieces without previous history ");
233 11 : keys.addFlag("MINTOZERO",false,"translate the resulting bias/histogram to have the minimum to zero ");
234 11 : keys.add("optional","FMT","the format that should be used to output real numbers");
235 22 : keys.setValueDescription("scalar","a scalar");
236 11 : }
237 :
238 9 : FuncSumHills::FuncSumHills(const ActionOptions&ao):
239 : Action(ao),
240 : Function(ao),
241 9 : initstride(-1),
242 9 : iscltool(false),
243 9 : integratehills(false),
244 9 : integratehisto(false),
245 9 : parallelread(false),
246 9 : negativebias(false),
247 9 : nohistory(false),
248 9 : minTOzero(false),
249 9 : doInt(false),
250 9 : lowI_(-1.),
251 9 : uppI_(-1.),
252 9 : beta(-1.),
253 9 : fmt("%14.9f") {
254 :
255 : // format
256 9 : parse("FMT",fmt);
257 9 : log<<" Output format is "<<fmt<<"\n";
258 : // here read
259 : // Grid Stuff
260 : std::vector<std::string> gmin;
261 18 : parseVector("GRID_MIN",gmin);
262 9 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) {
263 0 : error("not enough values for GRID_MIN");
264 : }
265 9 : plumed_massert(gmin.size()==getNumberOfArguments() || gmin.size()==0,"need GRID_MIN argument for this") ;
266 : std::vector<std::string> gmax;
267 18 : parseVector("GRID_MAX",gmax);
268 9 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) {
269 0 : error("not enough values for GRID_MAX");
270 : }
271 9 : plumed_massert(gmax.size()==getNumberOfArguments() || gmax.size()==0,"need GRID_MAX argument for this") ;
272 : std::vector<unsigned> gbin;
273 : std::vector<double> gspacing;
274 18 : parseVector("GRID_BIN",gbin);
275 9 : plumed_massert(gbin.size()==getNumberOfArguments() || gbin.size()==0,"need GRID_BIN argument for this") ;
276 9 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) {
277 0 : error("not enough values for GRID_BIN");
278 : }
279 18 : parseVector("GRID_SPACING",gspacing);
280 9 : plumed_massert(gspacing.size()==getNumberOfArguments() || gspacing.size()==0,"need GRID_SPACING argument for this") ;
281 9 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) {
282 0 : error("not enough values for GRID_SPACING");
283 : }
284 9 : if(gspacing.size()!=0 && gbin.size()==0) {
285 1 : log<<" The number of bins will be estimated from GRID_SPACING\n";
286 8 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
287 0 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
288 0 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
289 : }
290 9 : if(gspacing.size()!=0)
291 2 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
292 1 : if(gbin.size()==0) {
293 1 : gbin.assign(getNumberOfArguments(),1);
294 : }
295 : double a,b;
296 1 : Tools::convert(gmin[i],a);
297 1 : Tools::convert(gmax[i],b);
298 1 : unsigned n=std::ceil((b-a)/gspacing[i]);
299 1 : if(gbin[i]<n) {
300 1 : gbin[i]=n;
301 : }
302 : }
303 :
304 : // Inteval keyword
305 9 : std::vector<double> tmpI(2);
306 18 : parseVector("INTERVAL",tmpI);
307 9 : if(tmpI.size()!=2&&tmpI.size()!=0) {
308 0 : error("both a lower and an upper limits must be provided with INTERVAL");
309 9 : } else if(tmpI.size()==2) {
310 0 : lowI_=tmpI.at(0);
311 0 : uppI_=tmpI.at(1);
312 0 : if(getNumberOfArguments()!=1) {
313 0 : error("INTERVAL limits correction works only for monodimensional metadynamics!");
314 : }
315 0 : if(uppI_<lowI_) {
316 0 : error("The Upper limit must be greater than the Lower limit!");
317 : }
318 0 : doInt=true;
319 : }
320 9 : if(doInt) {
321 0 : log << " Upper and Lower limits boundaries for the bias are activated at " << lowI_ << " - " << uppI_<<"\n";
322 0 : log << " Using the same values as boundaries for the grid if not other value was defined (default: 200 bins)\n";
323 0 : std::ostringstream strsmin, strsmax;
324 0 : strsmin << lowI_;
325 0 : strsmax << uppI_;
326 0 : if(gmin.size()==0) {
327 0 : gmin.push_back(strsmin.str());
328 : }
329 0 : if(gmax.size()==0) {
330 0 : gmax.push_back(strsmax.str());
331 : }
332 0 : if(gbin.size()==0) {
333 0 : gbin.push_back(200);
334 : }
335 0 : }
336 :
337 :
338 : // hills file:
339 18 : parseVector("HILLSFILES",hillsFiles);
340 9 : if(hillsFiles.size()==0) {
341 1 : integratehills=false; // default behaviour
342 : } else {
343 8 : integratehills=true;
344 17 : for(unsigned i=0; i<hillsFiles.size(); i++) {
345 9 : log<<" hillsfile : "<<hillsFiles[i]<<"\n";
346 : }
347 : }
348 : // histo file:
349 18 : parseVector("HISTOFILES",histoFiles);
350 9 : if(histoFiles.size()==0) {
351 8 : integratehisto=false;
352 : } else {
353 1 : integratehisto=true;
354 2 : for(unsigned i=0; i<histoFiles.size(); i++) {
355 1 : log<<" histofile : "<<histoFiles[i]<<"\n";
356 : }
357 : }
358 : std::vector<double> histoSigma;
359 9 : if(integratehisto) {
360 1 : parseVector("HISTOSIGMA",histoSigma);
361 3 : for(unsigned i=0; i<histoSigma.size(); i++) {
362 2 : log<<" histosigma : "<<histoSigma[i]<<"\n";
363 : }
364 : }
365 :
366 : // needs a projection?
367 : proj.clear();
368 9 : parseVector("PROJ",proj);
369 9 : if(integratehills) {
370 8 : plumed_massert(proj.size()<getNumberOfArguments()," The number of projection must be less than the full list of arguments ");
371 : }
372 9 : if(integratehisto) {
373 1 : plumed_massert(proj.size()<=getNumberOfArguments()," The number of projection must be less or equal to the full list of arguments ");
374 : }
375 9 : if(integratehisto&&proj.size()==0) {
376 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
377 2 : proj.push_back(getPntrToArgument(i)->getName());
378 : }
379 : }
380 :
381 : // add some automatic hills width: not in case stride is defined
382 : // since when you start from zero the automatic size will be zero!
383 9 : if(gmin.size()==0 || gmax.size()==0) {
384 5 : log<<" \n";
385 5 : log<<" No boundaries defined: need to do a prescreening of hills \n";
386 : std::vector<Value*> tmphillsvalues, tmphistovalues;
387 5 : if(integratehills) {
388 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
389 8 : tmphillsvalues.push_back( getPntrToArgument(i) );
390 : }
391 : }
392 5 : if(integratehisto) {
393 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
394 2 : std::string ss = getPntrToArgument(i)->getName();
395 6 : for(unsigned j=0; j<proj.size(); j++) {
396 4 : if(proj[j]==ss) {
397 2 : tmphistovalues.push_back( getPntrToArgument(i) );
398 : }
399 : }
400 : }
401 : }
402 :
403 5 : if(integratehills) {
404 4 : FilesHandler hillsHandler(hillsFiles,parallelread,*this, log);
405 : std::vector<double> vmin,vmax;
406 : std::vector<unsigned> vbin;
407 4 : hillsHandler.getMinMaxBin(tmphillsvalues,comm,vmin,vmax,vbin);
408 4 : log<<" found boundaries from hillsfile: \n";
409 4 : gmin.resize(vmin.size());
410 4 : gmax.resize(vmax.size());
411 4 : if(gbin.size()==0) {
412 3 : gbin=vbin;
413 : } else {
414 1 : log<<" found nbins in input, this overrides the automatic choice \n";
415 : }
416 12 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
417 8 : Tools::convert(vmin[i],gmin[i]);
418 8 : Tools::convert(vmax[i],gmax[i]);
419 8 : log<<" variable "<< getPntrToArgument(i)->getName()<<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
420 : }
421 : }
422 : // if at this stage bins are not there then do it with histo
423 5 : if(gmin.size()==0) {
424 1 : FilesHandler histoHandler(histoFiles,parallelread,*this, log);
425 : std::vector<double> vmin,vmax;
426 : std::vector<unsigned> vbin;
427 1 : histoHandler.getMinMaxBin(tmphistovalues,comm,vmin,vmax,vbin,histoSigma);
428 1 : log<<" found boundaries from histofile: \n";
429 1 : gmin.resize(vmin.size());
430 1 : gmax.resize(vmax.size());
431 1 : if(gbin.size()==0) {
432 0 : gbin=vbin;
433 : } else {
434 1 : log<<" found nbins in input, this overrides the automatic choice \n";
435 : }
436 3 : for(unsigned i=0; i<proj.size(); i++) {
437 2 : Tools::convert(vmin[i],gmin[i]);
438 2 : Tools::convert(vmax[i],gmax[i]);
439 2 : log<<" variable "<< proj[i] <<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
440 : }
441 : }
442 5 : log<<" done!\n";
443 5 : log<<" \n";
444 : }
445 :
446 :
447 9 : if( proj.size() != 0 || integratehisto==true ) {
448 3 : parse("KT",beta);
449 7 : for(unsigned i=0; i<proj.size(); i++) {
450 4 : log<<" projection "<<i<<" : "<<proj[i]<<"\n";
451 : }
452 : // this should be only for projection or free energy from histograms
453 3 : plumed_massert(beta>0.,"if you make a projection or a histogram correction then you need KT flag!");
454 3 : beta=1./beta;
455 3 : log<<" beta is "<<beta<<"\n";
456 : }
457 : // is a cltool: then you start and then die
458 9 : parseFlag("ISCLTOOL",iscltool);
459 : //
460 9 : parseFlag("NEGBIAS",negativebias);
461 : //
462 9 : parseFlag("PARALLELREAD",parallelread);
463 : // stride
464 9 : parse("INITSTRIDE",initstride);
465 : // output suffix or names
466 9 : if(initstride<0) {
467 8 : log<<" Doing only one integration: no stride \n";
468 : outhills="fes.dat";
469 : outhisto="histo.dat";
470 : } else {
471 : outhills="fes_";
472 : outhisto="histo_";
473 1 : log<<" Doing integration slices every "<<initstride<<" kernels\n";
474 1 : parseFlag("NOHISTORY",nohistory);
475 1 : if(nohistory) {
476 0 : log<<" nohistory: each stride block has no memory of the previous block\n";
477 : }
478 : }
479 9 : parseFlag("MINTOZERO",minTOzero);
480 9 : if(minTOzero) {
481 0 : log<<" mintozero: bias/histogram will be translated to have the minimum value equal to zero\n";
482 : }
483 : //what might it be this?
484 : // here start
485 : // want something right now?? do it and return
486 : // your argument is a set of cvs
487 : // then you need: a hills / a colvar-like file (to do a histogram)
488 : // create a bias representation for this
489 9 : if(iscltool) {
490 :
491 : std::vector<Value*> tmphillsvalues, tmphistovalues;
492 9 : if(integratehills) {
493 23 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
494 : // allocate a new value from the old one: no deriv here
495 : // if we are summing hills then all the arguments are needed
496 15 : tmphillsvalues.push_back( getPntrToArgument(i) );
497 : }
498 : }
499 9 : if(integratehisto) {
500 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
501 2 : std::string ss = getPntrToArgument(i)->getName();
502 6 : for(unsigned j=0; j<proj.size(); j++) {
503 4 : if(proj[j]==ss) {
504 2 : tmphistovalues.push_back( getPntrToArgument(i) );
505 : }
506 : }
507 : }
508 : }
509 :
510 : // check if the files exists
511 9 : if(integratehills) {
512 8 : checkFilesAreExisting(hillsFiles);
513 16 : biasrep=Tools::make_unique<BiasRepresentation>(tmphillsvalues,comm, gmin, gmax, gbin, doInt, lowI_, uppI_);
514 8 : if(negativebias) {
515 1 : biasrep->setRescaledToBias(true);
516 1 : log<<" required the -bias instead of the free energy \n";
517 1 : if(initstride<0) {
518 : outhills="negativebias.dat";
519 : } else {
520 : outhills="negativebias_";
521 : }
522 : }
523 : }
524 :
525 9 : parse("OUTHILLS",outhills);
526 9 : parse("OUTHISTO",outhisto);
527 9 : if(integratehills) {
528 8 : log<<" output file for fes/bias is : "<<outhills<<"\n";
529 : }
530 9 : if(integratehisto) {
531 1 : log<<" output file for histogram is : "<<outhisto<<"\n";
532 : }
533 9 : checkRead();
534 :
535 9 : log<<"\n";
536 9 : log<<" Now calculating...\n";
537 9 : log<<"\n";
538 :
539 : // here it defines the column to be histogrammed, tmpvalues should be only
540 : // the list of the collective variable one want to consider
541 9 : if(integratehisto) {
542 1 : checkFilesAreExisting(histoFiles);
543 2 : historep=Tools::make_unique<BiasRepresentation>(tmphistovalues,comm,gmin,gmax,gbin,histoSigma);
544 : }
545 :
546 : // decide how to source hills ( serial/parallel )
547 : // here below the input control
548 : // say how many hills and it will read them from the
549 : // bunch of files provided, will update the representation
550 : // of hills (i.e. a list of hills and the associated grid)
551 :
552 : // decide how to source colvars ( serial parallel )
553 9 : std::unique_ptr<FilesHandler> hillsHandler;
554 9 : std::unique_ptr<FilesHandler> histoHandler;
555 :
556 9 : if(integratehills) {
557 16 : hillsHandler=Tools::make_unique<FilesHandler>(hillsFiles,parallelread,*this, log);
558 : }
559 9 : if(integratehisto) {
560 2 : histoHandler=Tools::make_unique<FilesHandler>(histoFiles,parallelread,*this, log);
561 : }
562 :
563 : // Stopwatch is logged when it goes out of scope
564 9 : Stopwatch sw(log);
565 :
566 : // Stopwatch is stopped when swh goes out of scope
567 9 : auto swh=sw.startStop("0 Summing hills");
568 :
569 : // read a number of hills and put in the bias representation
570 : int nfiles=0;
571 9 : bool ibias=integratehills;
572 9 : bool ihisto=integratehisto;
573 : while(true) {
574 10 : if( integratehills && ibias ) {
575 9 : if(nohistory) {
576 0 : biasrep->clear();
577 0 : log<<" clearing history before reading a new block\n";
578 : };
579 9 : log<<" reading hills: \n";
580 9 : ibias=hillsHandler->readBunch(biasrep.get(),initstride) ;
581 9 : log<<"\n";
582 : }
583 :
584 10 : if( integratehisto && ihisto ) {
585 1 : if(nohistory) {
586 0 : historep->clear();
587 0 : log<<" clearing history before reading a new block\n";
588 : };
589 1 : log<<" reading histogram: \n";
590 1 : ihisto=histoHandler->readBunch(historep.get(),initstride) ;
591 1 : log<<"\n";
592 : }
593 :
594 : // dump: need to project?
595 10 : if(proj.size()!=0) {
596 :
597 4 : if(integratehills) {
598 :
599 3 : log<<" Bias: Projecting on subgrid... \n";
600 3 : BiasWeight Bw(beta);
601 3 : Grid biasGrid=*(biasrep->getGridPtr());
602 3 : Grid smallGrid=biasGrid.project(proj,&Bw);
603 3 : OFile gridfile;
604 3 : gridfile.link(*this);
605 3 : std::ostringstream ostr;
606 3 : ostr<<nfiles;
607 : std::string myout;
608 3 : if(initstride>0) {
609 4 : myout=outhills+ostr.str()+".dat" ;
610 : } else {
611 : myout=outhills;
612 : }
613 3 : log<<" Bias: Writing subgrid on file "<<myout<<" \n";
614 3 : gridfile.open(myout);
615 3 : if(minTOzero) {
616 0 : smallGrid.setMinToZero();
617 : }
618 : smallGrid.setOutputFmt(fmt);
619 3 : smallGrid.writeToFile(gridfile);
620 3 : gridfile.close();
621 3 : if(!ibias) {
622 2 : integratehills=false; // once you get to the final bunch just give up
623 : }
624 3 : }
625 : // this should be removed
626 4 : if(integratehisto) {
627 :
628 1 : log<<" Histo: Projecting on subgrid... \n";
629 1 : Grid histoGrid=*(historep->getGridPtr());
630 :
631 1 : OFile gridfile;
632 1 : gridfile.link(*this);
633 1 : std::ostringstream ostr;
634 1 : ostr<<nfiles;
635 : std::string myout;
636 1 : if(initstride>0) {
637 0 : myout=outhisto+ostr.str()+".dat" ;
638 : } else {
639 : myout=outhisto;
640 : }
641 1 : log<<" Histo: Writing subgrid on file "<<myout<<" \n";
642 1 : gridfile.open(myout);
643 :
644 1 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
645 1 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
646 1 : if(minTOzero) {
647 0 : histoGrid.setMinToZero();
648 : }
649 : histoGrid.setOutputFmt(fmt);
650 1 : histoGrid.writeToFile(gridfile);
651 :
652 1 : if(!ihisto) {
653 1 : integratehisto=false; // once you get to the final bunch just give up
654 : }
655 1 : }
656 :
657 : } else {
658 :
659 6 : if(integratehills) {
660 :
661 6 : Grid biasGrid=*(biasrep->getGridPtr());
662 6 : biasGrid.scaleAllValuesAndDerivatives(-1.);
663 :
664 6 : OFile gridfile;
665 6 : gridfile.link(*this);
666 6 : std::ostringstream ostr;
667 6 : ostr<<nfiles;
668 : std::string myout;
669 6 : if(initstride>0) {
670 0 : myout=outhills+ostr.str()+".dat" ;
671 : } else {
672 : myout=outhills;
673 : }
674 6 : log<<" Writing full grid on file "<<myout<<" \n";
675 6 : gridfile.open(myout);
676 :
677 6 : if(minTOzero) {
678 0 : biasGrid.setMinToZero();
679 : }
680 : biasGrid.setOutputFmt(fmt);
681 6 : biasGrid.writeToFile(gridfile);
682 : // rescale back prior to accumulate
683 6 : if(!ibias) {
684 6 : integratehills=false; // once you get to the final bunch just give up
685 : }
686 6 : }
687 6 : if(integratehisto) {
688 :
689 0 : Grid histoGrid=*(historep->getGridPtr());
690 : // do this if you want a free energy from a grid, otherwise do not
691 0 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
692 0 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
693 :
694 0 : OFile gridfile;
695 0 : gridfile.link(*this);
696 0 : std::ostringstream ostr;
697 0 : ostr<<nfiles;
698 : std::string myout;
699 0 : if(initstride>0) {
700 0 : myout=outhisto+ostr.str()+".dat" ;
701 : } else {
702 : myout=outhisto;
703 : }
704 0 : log<<" Writing full grid on file "<<myout<<" \n";
705 0 : gridfile.open(myout);
706 :
707 : // also this is useful only for free energy
708 0 : if(minTOzero) {
709 0 : histoGrid.setMinToZero();
710 : }
711 : histoGrid.setOutputFmt(fmt);
712 0 : histoGrid.writeToFile(gridfile);
713 :
714 0 : if(!ihisto) {
715 0 : integratehisto=false; // once you get to the final bunch just give up
716 : }
717 0 : }
718 : }
719 10 : if ( !ibias && !ihisto) {
720 : break; //when both are over then just quit
721 : }
722 :
723 1 : nfiles++;
724 1 : }
725 :
726 : return;
727 9 : }
728 : // just an initialization but you need to do something on the fly?: need to connect with a metad run and its grid representation
729 : // your argument is a metad run
730 : // if the grid does not exist crash and say that you need some data
731 : // otherwise just link with it
732 :
733 9 : }
734 :
735 0 : void FuncSumHills::calculate() {
736 : // this should be connected only with a grid representation to metadynamics
737 : // at regular time just dump it
738 0 : plumed_merror("You should have never got here: this stuff is not yet implemented!");
739 : }
740 :
741 9 : bool FuncSumHills::checkFilesAreExisting(const std::vector<std::string> & hills ) {
742 9 : plumed_massert(hills.size()!=0,"the number of files provided should be at least one" );
743 : auto ifile=Tools::make_unique<IFile>();
744 9 : ifile->link(*this);
745 19 : for(unsigned i=0; i< hills.size(); i++) {
746 10 : plumed_massert(ifile->FileExist(hills[i]),"missing file "+hills[i]);
747 : }
748 9 : return true;
749 :
750 9 : }
751 :
752 : }
753 :
754 : }
755 :
756 :
|
