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