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