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 <limits>
30 : #include "tools/Tools.h"
31 : #include "tools/Stopwatch.h"
32 : #include "tools/Grid.h"
33 :
34 : namespace PLMD {
35 : namespace function {
36 :
37 :
38 : //+PLUMEDOC FUNCTION FUNCSUMHILLS
39 : /*
40 : This function is used by the command line tool [sum_hills](sum_hills.md) to sum the hills of a metadynamics simulation
41 : or to compute a histogram. It does all its work when it is created and it cannot be used in a simulation (it stops with an error).
42 :
43 : All its options are set by [sum_hills](sum_hills.md): see the documentation of [sum_hills](sum_hills.md) for their meaning.
44 :
45 : ## Examples
46 :
47 : There are currently no examples for this action: use [sum_hills](sum_hills.md) instead.
48 :
49 : */
50 : //+ENDPLUMEDOC
51 :
52 : class FilesHandler {
53 : std::vector <std::string> filenames;
54 : std::vector <std::unique_ptr<IFile>> ifiles;
55 : Action *action;
56 : Log *log;
57 : bool parallelread;
58 : unsigned beingread=0;
59 : bool isopen=false;
60 : /// key used to merge the files ("time" or "clock"); empty to read the files one after the other
61 : std::string sortby;
62 : /// a hill read from one of the files and not yet added (used when merging)
63 6 : struct PendingHill {
64 : std::unique_ptr<KernelFunctions> kernel;
65 : double biasf=1.0;
66 : double clock=0.0;
67 : double time=0.0;
68 : bool valid=false;
69 : };
70 : std::vector<PendingHill> pending;
71 : bool mergestarted=false;
72 : bool unsortedwarning=false;
73 : /// read the next hill of file i into pending[i] (merge mode)
74 : void readPending(BiasRepresentation *br, unsigned i);
75 : bool readBunchMerged(BiasRepresentation *br, int stride);
76 : public:
77 : FilesHandler(const std::vector<std::string> &filenames,
78 : bool parallelread,
79 : Action &myaction,
80 : Log &mylog,
81 : const std::string & sortby="");
82 : bool readBunch(BiasRepresentation *br, int stride=-1);
83 : bool scanOneHill(BiasRepresentation *br, IFile *ifile );
84 : void getMinMaxBin(const std::vector<Value*> & vals,
85 : Communicator &cc,
86 : std::vector<double> &vmin,
87 : std::vector<double> &vmax,
88 : std::vector<unsigned> &vbin);
89 : void getMinMaxBin(const std::vector<Value*> & vals,
90 : Communicator &cc,
91 : std::vector<double> &vmin,
92 : std::vector<double> &vmax,
93 : std::vector<unsigned> &vbin,
94 : const std::vector<double> &histosigma);
95 : };
96 23 : FilesHandler::FilesHandler(const std::vector<std::string> &fnames,
97 : const bool doparallelread,
98 : Action &act,
99 : Log &mylog,
100 23 : const std::string & sort)
101 23 : :filenames(fnames),
102 23 : log(&mylog),
103 23 : parallelread(doparallelread),
104 23 : sortby(sort) {
105 23 : plumed_massert(sortby=="" || sortby=="time" || sortby=="clock","files can only be sorted by time or clock");
106 23 : action=&act;
107 49 : for(unsigned i=0; i<filenames.size(); i++) {
108 : auto ifile=Tools::make_unique<IFile>();
109 26 : ifile->link(act);
110 26 : plumed_massert((ifile->FileExist(filenames[i])), "the file "+filenames[i]+" does not exist " );
111 26 : ifiles.emplace_back(std::move(ifile));
112 26 : }
113 :
114 23 : }
115 :
116 : // note that the FileHandler is completely transparent respect to the biasrepresentation
117 : // no check are made at this level
118 123 : void FilesHandler::readPending(BiasRepresentation *br, unsigned i) {
119 123 : IFile *ff=ifiles[i].get();
120 : PendingHill & p=pending[i];
121 : double t;
122 246 : if(!ff->scanField("time",t)) {
123 3 : (*log)<<" closing file "<<filenames[i]<<"\n";
124 3 : ff->close();
125 3 : p.valid=false;
126 : p.kernel.reset();
127 3 : return;
128 : }
129 120 : double c=0.0;
130 240 : if(ff->FieldExist("clock")) {
131 240 : ff->scanField("clock",c);
132 0 : } else if(sortby=="clock") {
133 0 : plumed_merror("cannot sort by clock: there is no clock field in file "+filenames[i]);
134 : }
135 120 : if(p.valid && !unsortedwarning) {
136 117 : const bool back = (sortby=="clock") ? (c<p.clock || (c==p.clock && t<p.time)) : (t<p.time);
137 117 : if(back) {
138 0 : (*log)<<" WARNING: file "<<filenames[i]<<" is not ordered by "<<sortby<<": hills are merged in the order in which they appear in each file\n";
139 0 : unsortedwarning=true;
140 : }
141 : }
142 120 : p.time=t;
143 120 : p.clock=c;
144 240 : p.kernel=br->readKernel(ff,p.biasf);
145 120 : if(br->hasSigmaInInput()) {
146 0 : ff->allowIgnoredFields();
147 : }
148 120 : ff->scanField();
149 120 : p.valid=true;
150 : }
151 :
152 5 : bool FilesHandler::readBunchMerged(BiasRepresentation *br, int stride) {
153 5 : if(!mergestarted) {
154 1 : (*log)<<" merging the hills of "<<filenames.size()<<" files in order of "<<sortby<<"\n";
155 1 : pending.resize(ifiles.size());
156 4 : for(unsigned i=0; i<ifiles.size(); i++) {
157 3 : (*log)<<" opening file "<<filenames[i]<<"\n";
158 3 : ifiles[i]->open(filenames[i]);
159 3 : readPending(br,i);
160 : }
161 1 : mergestarted=true;
162 : }
163 : while(true) {
164 : // the next hill is the earliest; ties are resolved with time (for clock) and then with the order of the files
165 : int next=-1;
166 484 : for(unsigned i=0; i<pending.size(); i++) {
167 363 : if(!pending[i].valid) {
168 6 : continue;
169 : }
170 357 : if(next<0) {
171 120 : next=i;
172 120 : continue;
173 : }
174 : const auto & a=pending[i];
175 237 : const auto & b=pending[next];
176 237 : const bool earlier = (sortby=="clock") ? (a.clock<b.clock || (a.clock==b.clock && a.time<b.time)) : (a.time<b.time);
177 237 : if(earlier) {
178 78 : next=i;
179 : }
180 : }
181 121 : if(next<0) {
182 1 : (*log)<<" final chunk: now with "<<br->getNumberOfKernels()<<" kernels \n";
183 5 : return false;
184 : }
185 120 : br->addKernel(std::move(pending[next].kernel),pending[next].biasf,pending[next].time);
186 120 : readPending(br,next);
187 120 : int n=br->getNumberOfKernels();
188 120 : if(stride>0 && n%stride==0) {
189 : // stop here unless all files are over
190 4 : for(const auto & p : pending)
191 4 : if(p.valid) {
192 4 : (*log)<<" done with this chunk: now with "<<n<<" kernels \n";
193 : return true;
194 : }
195 : }
196 116 : }
197 : }
198 :
199 32 : bool FilesHandler::readBunch(BiasRepresentation *br, int stride) {
200 32 : if(sortby!="" && ifiles.size()>1) {
201 5 : return readBunchMerged(br,stride);
202 : }
203 : bool morefiles=true;
204 27 : if(parallelread) {
205 0 : (*log)<<" doing parallelread \n";
206 0 : plumed_merror("parallelread is not yet implemented !!!");
207 : } else {
208 27 : (*log)<<" doing serialread \n";
209 : // read one by one hills
210 : // is the type defined? if not, assume it is a gaussian
211 : IFile *ff;
212 27 : ff=ifiles[beingread].get();
213 27 : if(!isopen) {
214 22 : (*log)<<" opening file "<<filenames[beingread]<<"\n";
215 22 : ff->open(filenames[beingread]);
216 22 : isopen=true;
217 : }
218 27 : int n=0;
219 : while(true) {
220 : bool fileisover=true;
221 6046 : while(scanOneHill(br,ff)) {
222 : // here do the dump if needed
223 6023 : n=br->getNumberOfKernels();
224 6023 : if(stride>0 && n%stride==0 && n!=0 ) {
225 5 : (*log)<<" done with this chunk: now with "<<n<<" kernels \n";
226 : fileisover=false;
227 : break;
228 : }
229 : }
230 : if(fileisover) {
231 23 : (*log)<<" closing file "<<filenames[beingread]<<"\n";
232 23 : ff->close();
233 23 : isopen=false;
234 23 : (*log)<<" now total "<<br->getNumberOfKernels()<<" kernels \n";
235 23 : beingread++;
236 23 : if(beingread<ifiles.size()) {
237 : ff=ifiles[beingread].get();
238 1 : ff->open(filenames[beingread]);
239 1 : (*log)<<" opening file "<<filenames[beingread]<<"\n";
240 1 : isopen=true;
241 : } else {
242 : morefiles=false;
243 22 : (*log)<<" final chunk: now with "<<n<<" kernels \n";
244 : break;
245 : }
246 : }
247 : // if there are no more files to read and this file is over then quit
248 : if(fileisover && !morefiles) {
249 : break;
250 : }
251 : // if you are in the middle of a file and you are here
252 : // then means that you read what you need to read
253 6 : if(!fileisover ) {
254 : break;
255 : }
256 : }
257 : }
258 27 : return morefiles;
259 : }
260 7 : void FilesHandler::getMinMaxBin(const std::vector<Value*> & vals,
261 : Communicator &cc,
262 : std::vector<double> &vmin,
263 : std::vector<double> &vmax,
264 : std::vector<unsigned> &vbin) {
265 : // create the representation (no grid)
266 7 : BiasRepresentation br(vals,cc);
267 : // read all the kernels
268 7 : readBunch(&br);
269 : // loop over the kernels and get the support
270 7 : br.getMinMaxBin(vmin,vmax,vbin);
271 7 : }
272 1 : void FilesHandler::getMinMaxBin(const std::vector<Value*> & vals,
273 : Communicator &cc,
274 : std::vector<double> &vmin,
275 : std::vector<double> &vmax,
276 : std::vector<unsigned> &vbin,
277 : const std::vector<double> &histosigma) {
278 1 : BiasRepresentation br(vals,cc,histosigma);
279 : // read all the kernels
280 1 : readBunch(&br);
281 : // loop over the kernels and get the support
282 1 : br.getMinMaxBin(vmin,vmax,vbin);
283 1 : }
284 6046 : bool FilesHandler::scanOneHill(BiasRepresentation *br, IFile *ifile ) {
285 : double time;
286 12092 : if(ifile->scanField("time",time)) {
287 : // the other fields (e.g. clock) are ignored
288 6023 : br->pushKernel(ifile,time);
289 6023 : if(br->hasSigmaInInput()) {
290 1092 : ifile->allowIgnoredFields();
291 : }
292 6023 : ifile->scanField();
293 6023 : return true;
294 : } else {
295 : return false;
296 : }
297 : }
298 :
299 :
300 900 : double mylog( double v1 ) {
301 900 : return std::log(v1);
302 : }
303 :
304 1800 : double mylogder( double v1 ) {
305 1800 : return 1./v1;
306 : }
307 :
308 :
309 :
310 : class FuncSumHills :
311 : public Function {
312 : std::vector<std::string> hillsFiles,histoFiles;
313 : std::vector<std::string> proj;
314 : int initstride=-1;
315 : bool iscltool=false;
316 : bool integratehills=false;
317 : bool integratehisto=false;
318 : bool parallelread=false;
319 : bool negativebias=false;
320 : bool nohistory=false;
321 : bool minTOzero=false;
322 : bool doInt=false;
323 : double lowI_=-1.0;
324 : double uppI_=-1.0;
325 : double beta=-1.0;
326 : std::string outhills;
327 : std::string outhisto;
328 : std::string fmt="%14.9f";
329 : std::string sortby;
330 : /// time average of the bias
331 : std::vector<double> average;
332 : /// boxes (min and max for each variable) defining the basins whose free energies are written
333 : std::vector<std::vector<double>> basinmin,basinmax;
334 : std::string outbasins;
335 : /// write the free energy of each basin relative to the first one
336 : void writeBasins(OFile & ofile, const Grid & fes, double time, int nkernels);
337 : std::unique_ptr<BiasRepresentation> biasrep;
338 : std::unique_ptr<BiasRepresentation> historep;
339 : public:
340 : explicit FuncSumHills(const ActionOptions&);
341 : void calculate() override; // this probably is not needed
342 : bool checkFilesAreExisting(const std::vector<std::string> & hills );
343 : static void registerKeywords(Keywords& keys);
344 : };
345 :
346 : PLUMED_REGISTER_ACTION(FuncSumHills,"FUNCSUMHILLS")
347 :
348 19 : void FuncSumHills::registerKeywords(Keywords& keys) {
349 19 : Function::registerKeywords(keys);
350 19 : keys.add("optional","HILLSFILES"," source file for hills creation(may be the same as HILLS)"); // this can be a vector!
351 19 : keys.add("optional","HISTOFILES"," source file for histogram creation(may be the same as HILLS)"); // also this can be a vector!
352 19 : keys.add("optional","HISTOSIGMA"," sigmas for binning when the histogram correction is needed ");
353 19 : keys.add("optional","PROJ"," only with sumhills: the projection on the CVs");
354 19 : keys.add("optional","KT"," only with sumhills: the kt factor when projection on CVs");
355 19 : keys.add("optional","GRID_MIN","the lower bounds for the grid");
356 19 : keys.add("optional","GRID_MAX","the upper bounds for the grid");
357 19 : keys.add("optional","GRID_BIN","the number of bins for the grid");
358 19 : keys.add("optional","GRID_SPACING","the approximate grid spacing (to be used as an alternative or together with GRID_BIN)");
359 19 : keys.add("optional","INTERVAL","set one dimensional INTERVAL");
360 19 : keys.add("optional","OUTHILLS"," output file for hills ");
361 19 : keys.add("optional","OUTHISTO"," output file for histogram ");
362 19 : keys.add("optional","INITSTRIDE"," stride if you want an initial dump ");
363 19 : keys.add("optional","STRIDE"," stride when you do it on the fly ");
364 19 : keys.addFlag("ISCLTOOL",false,"use via plumed command line: calculate at read phase and then go");
365 19 : keys.addFlag("PARALLELREAD",false,"read parallel HILLS file");
366 19 : keys.addFlag("NEGBIAS",false,"dump negative bias ( -bias ) instead of the free energy: needed in well tempered with flexible hills ");
367 19 : keys.addFlag("NOHISTORY",false,"to be used with INITSTRIDE: it splits the bias/histogram in pieces without previous history ");
368 19 : keys.addFlag("MINTOZERO",false,"translate the resulting bias/histogram to have the minimum to zero ");
369 19 : keys.add("optional","FMT","the format that should be used to output real numbers");
370 19 : keys.add("optional","SORT","merge the hills of multiple files in order of time or clock instead of reading the files one after the other");
371 19 : keys.add("optional","AVERAGE","time average of the bias: from time t0 (and until t1 if given) the hills are scaled by (t1-t)/(t1-t0), where t1 is by default the time of the last hill");
372 19 : keys.add("optional","BASINS","basins (separated by colons) given by the minimum and maximum of each variable of the output, e.g. -3,-1:1,3; their free energies relative to the first basin are written at each output");
373 19 : keys.add("optional","OUTBASINS","output file for the free energies of the basins");
374 38 : keys.setValueDescription("scalar","a scalar");
375 19 : }
376 :
377 15 : FuncSumHills::FuncSumHills(const ActionOptions&ao):
378 : Action(ao),
379 15 : Function(ao) {
380 : // format
381 15 : parse("FMT",fmt);
382 15 : log<<" Output format is "<<fmt<<"\n";
383 : // here read
384 : // Grid Stuff
385 15 : std::vector<std::string> gmin;
386 30 : parseVector("GRID_MIN",gmin);
387 15 : if(gmin.size()!=getNumberOfArguments() && gmin.size()!=0) {
388 0 : error("not enough values for GRID_MIN");
389 : }
390 15 : plumed_massert(gmin.size()==getNumberOfArguments() || gmin.size()==0,"need GRID_MIN argument for this") ;
391 15 : std::vector<std::string> gmax;
392 30 : parseVector("GRID_MAX",gmax);
393 15 : if(gmax.size()!=getNumberOfArguments() && gmax.size()!=0) {
394 0 : error("not enough values for GRID_MAX");
395 : }
396 15 : plumed_massert(gmax.size()==getNumberOfArguments() || gmax.size()==0,"need GRID_MAX argument for this") ;
397 15 : std::vector<unsigned> gbin;
398 15 : std::vector<double> gspacing;
399 30 : parseVector("GRID_BIN",gbin);
400 15 : plumed_massert(gbin.size()==getNumberOfArguments() || gbin.size()==0,"need GRID_BIN argument for this") ;
401 15 : if(gbin.size()!=getNumberOfArguments() && gbin.size()!=0) {
402 0 : error("not enough values for GRID_BIN");
403 : }
404 30 : parseVector("GRID_SPACING",gspacing);
405 15 : plumed_massert(gspacing.size()==getNumberOfArguments() || gspacing.size()==0,"need GRID_SPACING argument for this") ;
406 15 : if(gspacing.size()!=getNumberOfArguments() && gspacing.size()!=0) {
407 0 : error("not enough values for GRID_SPACING");
408 : }
409 15 : if(gspacing.size()!=0 && gbin.size()==0) {
410 2 : log<<" The number of bins will be estimated from GRID_SPACING\n";
411 13 : } else if(gspacing.size()!=0 && gbin.size()!=0) {
412 0 : log<<" You specified both GRID_BIN and GRID_SPACING\n";
413 0 : log<<" The more conservative (highest) number of bins will be used for each variable\n";
414 : }
415 : // the number of bins from GRID_SPACING is computed once the boundaries are known,
416 : // possibly after they have been detected from the hills
417 : const bool userbin=(gbin.size()!=0);
418 11 : auto applySpacing=[&]() {
419 11 : if(gspacing.size()==0) {
420 : return;
421 : }
422 2 : if(gbin.size()==0) {
423 2 : gbin.assign(gmin.size(),1);
424 : }
425 5 : for(unsigned i=0; i<gmin.size(); i++) {
426 : double a,b;
427 3 : Tools::convert(gmin[i],a);
428 3 : Tools::convert(gmax[i],b);
429 3 : unsigned n=std::ceil((b-a)/gspacing[i]);
430 3 : if(gbin[i]<n) {
431 3 : gbin[i]=n;
432 : }
433 : }
434 15 : };
435 15 : if(gmin.size()!=0 && gmax.size()!=0) {
436 7 : applySpacing();
437 : }
438 :
439 : // Inteval keyword
440 15 : std::vector<double> tmpI(2);
441 30 : parseVector("INTERVAL",tmpI);
442 15 : if(tmpI.size()!=2&&tmpI.size()!=0) {
443 0 : error("both a lower and an upper limits must be provided with INTERVAL");
444 15 : } else if(tmpI.size()==2) {
445 0 : lowI_=tmpI.at(0);
446 0 : uppI_=tmpI.at(1);
447 0 : if(getNumberOfArguments()!=1) {
448 0 : error("INTERVAL limits correction works only for monodimensional metadynamics!");
449 : }
450 0 : if(uppI_<lowI_) {
451 0 : error("The Upper limit must be greater than the Lower limit!");
452 : }
453 0 : doInt=true;
454 : }
455 15 : if(doInt) {
456 0 : log << " Upper and Lower limits boundaries for the bias are activated at " << lowI_ << " - " << uppI_<<"\n";
457 0 : log << " Using the same values as boundaries for the grid if not other value was defined (default: 200 bins)\n";
458 0 : std::ostringstream strsmin, strsmax;
459 0 : strsmin << lowI_;
460 0 : strsmax << uppI_;
461 0 : if(gmin.size()==0) {
462 0 : gmin.push_back(strsmin.str());
463 : }
464 0 : if(gmax.size()==0) {
465 0 : gmax.push_back(strsmax.str());
466 : }
467 0 : if(gbin.size()==0 && gspacing.size()==0) {
468 0 : gbin.push_back(200);
469 0 : } else if(gbin.size()==0) {
470 0 : applySpacing();
471 : }
472 0 : }
473 :
474 :
475 : // hills file:
476 30 : parseVector("HILLSFILES",hillsFiles);
477 15 : if(hillsFiles.size()==0) {
478 1 : integratehills=false; // default behaviour
479 : } else {
480 14 : integratehills=true;
481 31 : for(unsigned i=0; i<hillsFiles.size(); i++) {
482 17 : log<<" hillsfile : "<<hillsFiles[i]<<"\n";
483 : }
484 : }
485 : // histo file:
486 30 : parseVector("HISTOFILES",histoFiles);
487 15 : if(histoFiles.size()==0) {
488 14 : integratehisto=false;
489 : } else {
490 1 : integratehisto=true;
491 2 : for(unsigned i=0; i<histoFiles.size(); i++) {
492 1 : log<<" histofile : "<<histoFiles[i]<<"\n";
493 : }
494 : }
495 15 : std::vector<double> histoSigma;
496 15 : if(integratehisto) {
497 1 : parseVector("HISTOSIGMA",histoSigma);
498 3 : for(unsigned i=0; i<histoSigma.size(); i++) {
499 2 : log<<" histosigma : "<<histoSigma[i]<<"\n";
500 : }
501 : }
502 :
503 : // needs a projection?
504 15 : proj.clear();
505 15 : parseVector("PROJ",proj);
506 15 : if(integratehills) {
507 14 : plumed_massert(proj.size()<getNumberOfArguments()," The number of projection must be less than the full list of arguments ");
508 : }
509 15 : if(integratehisto) {
510 1 : plumed_massert(proj.size()<=getNumberOfArguments()," The number of projection must be less or equal to the full list of arguments ");
511 : }
512 15 : if(integratehisto&&proj.size()==0) {
513 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
514 2 : proj.push_back(getPntrToArgument(i)->getName());
515 : }
516 : }
517 :
518 : // add some automatic hills width: not in case stride is defined
519 : // since when you start from zero the automatic size will be zero!
520 15 : if(gmin.size()==0 || gmax.size()==0) {
521 8 : log<<" \n";
522 8 : log<<" No boundaries defined: need to do a prescreening of hills \n";
523 8 : std::vector<Value*> tmphillsvalues, tmphistovalues;
524 8 : if(integratehills) {
525 21 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
526 14 : tmphillsvalues.push_back( getPntrToArgument(i) );
527 : }
528 : }
529 8 : if(integratehisto) {
530 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
531 2 : std::string ss = getPntrToArgument(i)->getName();
532 6 : for(unsigned j=0; j<proj.size(); j++) {
533 4 : if(proj[j]==ss) {
534 2 : tmphistovalues.push_back( getPntrToArgument(i) );
535 : }
536 : }
537 : }
538 : }
539 :
540 8 : if(integratehills) {
541 7 : FilesHandler hillsHandler(hillsFiles,parallelread,*this, log);
542 7 : std::vector<double> vmin,vmax;
543 7 : std::vector<unsigned> vbin;
544 7 : hillsHandler.getMinMaxBin(tmphillsvalues,comm,vmin,vmax,vbin);
545 7 : log<<" found boundaries from hillsfile: \n";
546 7 : gmin.resize(vmin.size());
547 7 : gmax.resize(vmax.size());
548 21 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
549 14 : Tools::convert(vmin[i],gmin[i]);
550 14 : Tools::convert(vmax[i],gmax[i]);
551 : }
552 7 : if(userbin || gspacing.size()!=0) {
553 3 : log<<" found nbins or spacing in input, this overrides the automatic choice \n";
554 3 : applySpacing();
555 : } else {
556 4 : gbin=vbin;
557 : }
558 21 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
559 14 : log<<" variable "<< getPntrToArgument(i)->getName()<<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
560 : }
561 14 : }
562 : // if at this stage bins are not there then do it with histo
563 8 : if(gmin.size()==0) {
564 1 : FilesHandler histoHandler(histoFiles,parallelread,*this, log);
565 1 : std::vector<double> vmin,vmax;
566 1 : std::vector<unsigned> vbin;
567 1 : histoHandler.getMinMaxBin(tmphistovalues,comm,vmin,vmax,vbin,histoSigma);
568 1 : log<<" found boundaries from histofile: \n";
569 1 : gmin.resize(vmin.size());
570 1 : gmax.resize(vmax.size());
571 3 : for(unsigned i=0; i<proj.size(); i++) {
572 2 : Tools::convert(vmin[i],gmin[i]);
573 2 : Tools::convert(vmax[i],gmax[i]);
574 : }
575 1 : if(userbin || gspacing.size()!=0) {
576 1 : log<<" found nbins or spacing in input, this overrides the automatic choice \n";
577 1 : applySpacing();
578 : } else {
579 0 : gbin=vbin;
580 : }
581 3 : for(unsigned i=0; i<proj.size(); i++) {
582 2 : log<<" variable "<< proj[i] <<" min: "<<gmin[i]<<" max: "<<gmax[i]<<" nbin: "<<gbin[i]<<"\n";
583 : }
584 2 : }
585 8 : log<<" done!\n";
586 8 : log<<" \n";
587 : }
588 :
589 :
590 30 : parseVector("AVERAGE",average);
591 15 : if(average.size()!=0) {
592 1 : if(average.size()>2) {
593 0 : error("AVERAGE needs the initial time and optionally the final time");
594 : }
595 1 : if(!integratehills) {
596 0 : error("AVERAGE can only be used with hills files");
597 : }
598 1 : if(average.size()==2 && average[1]<=average[0]) {
599 0 : error("the final time of AVERAGE should be larger than the initial time");
600 : }
601 1 : log<<" time average of the bias from time "<<average[0];
602 1 : if(average.size()==2) {
603 1 : log<<" to time "<<average[1];
604 : }
605 1 : log<<"\n";
606 : }
607 : std::string basins;
608 30 : parse("BASINS",basins);
609 15 : if(basins.length()>0) {
610 : // the basins are defined on the variables of the output
611 1 : const unsigned nout=(proj.size()>0 ? proj.size() : getNumberOfArguments());
612 4 : for(const auto & b : Tools::getWords(basins,":")) {
613 3 : std::vector<std::string> words=Tools::getWords(b,",");
614 3 : if(words.size()!=2*nout) {
615 0 : error("each basin needs a minimum and a maximum for each of the "+std::to_string(nout)+" variables of the output: "+b);
616 : }
617 3 : std::vector<double> bmin(nout),bmax(nout);
618 9 : for(unsigned i=0; i<nout; i++) {
619 6 : if(!Tools::convertNoexcept(words[2*i],bmin[i]) || !Tools::convertNoexcept(words[2*i+1],bmax[i])) {
620 0 : error("cannot read basin "+b);
621 : }
622 : }
623 3 : basinmin.push_back(bmin);
624 3 : basinmax.push_back(bmax);
625 4 : }
626 1 : if(basinmin.size()<2) {
627 0 : error("BASINS needs at least two basins");
628 : }
629 1 : log<<" free energies of "<<basinmin.size()<<" basins relative to the first one are written at each output\n";
630 : }
631 15 : if( proj.size() != 0 || integratehisto==true || basinmin.size()>0 ) {
632 4 : parse("KT",beta);
633 8 : for(unsigned i=0; i<proj.size(); i++) {
634 4 : log<<" projection "<<i<<" : "<<proj[i]<<"\n";
635 : }
636 : // this should be only for projection or free energy from histograms
637 4 : plumed_massert(beta>0.,"if you make a projection, a histogram correction or use basins then you need KT flag!");
638 4 : beta=1./beta;
639 4 : log<<" beta is "<<beta<<"\n";
640 : }
641 : // is a cltool: then you start and then die
642 15 : parseFlag("ISCLTOOL",iscltool);
643 : //
644 15 : parseFlag("NEGBIAS",negativebias);
645 : //
646 15 : parseFlag("PARALLELREAD",parallelread);
647 : // stride
648 15 : parse("INITSTRIDE",initstride);
649 : // output suffix or names
650 15 : if(initstride<0) {
651 11 : log<<" Doing only one integration: no stride \n";
652 11 : outhills="fes.dat";
653 11 : outhisto="histo.dat";
654 : } else {
655 4 : outhills="fes_";
656 4 : outhisto="histo_";
657 4 : log<<" Doing integration slices every "<<initstride<<" kernels\n";
658 4 : parseFlag("NOHISTORY",nohistory);
659 4 : if(nohistory) {
660 0 : log<<" nohistory: each stride block has no memory of the previous block\n";
661 : }
662 : }
663 15 : parse("SORT",sortby);
664 15 : if(sortby!="" && sortby!="time" && sortby!="clock") {
665 0 : error("SORT can only be time or clock");
666 : }
667 15 : if(sortby!="") {
668 1 : log<<" hills from multiple files are merged in order of "<<sortby<<"\n";
669 : }
670 15 : parseFlag("MINTOZERO",minTOzero);
671 15 : if(minTOzero) {
672 0 : log<<" mintozero: bias/histogram will be translated to have the minimum value equal to zero\n";
673 : }
674 : //what might it be this?
675 : // here start
676 : // want something right now?? do it and return
677 : // your argument is a set of cvs
678 : // then you need: a hills / a colvar-like file (to do a histogram)
679 : // create a bias representation for this
680 15 : if(iscltool) {
681 :
682 15 : std::vector<Value*> tmphillsvalues, tmphistovalues;
683 15 : if(integratehills) {
684 38 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
685 : // allocate a new value from the old one: no deriv here
686 : // if we are summing hills then all the arguments are needed
687 24 : tmphillsvalues.push_back( getPntrToArgument(i) );
688 : }
689 : }
690 15 : if(integratehisto) {
691 3 : for(unsigned i=0; i<getNumberOfArguments(); i++) {
692 2 : std::string ss = getPntrToArgument(i)->getName();
693 6 : for(unsigned j=0; j<proj.size(); j++) {
694 4 : if(proj[j]==ss) {
695 2 : tmphistovalues.push_back( getPntrToArgument(i) );
696 : }
697 : }
698 : }
699 : }
700 :
701 : // check if the files exists
702 15 : if(integratehills) {
703 14 : checkFilesAreExisting(hillsFiles);
704 28 : biasrep=Tools::make_unique<BiasRepresentation>(tmphillsvalues,comm, gmin, gmax, gbin, doInt, lowI_, uppI_);
705 14 : if(average.size()>0) {
706 1 : biasrep->setTimeAverage(average[0],average.size()==2,(average.size()==2?average[1]:0.0));
707 : }
708 14 : if(negativebias && basinmin.size()>0) {
709 0 : error("BASINS requires a free energy and cannot be used with NEGBIAS");
710 : }
711 14 : if(negativebias) {
712 2 : biasrep->setRescaledToBias(true);
713 2 : log<<" required the -bias instead of the free energy \n";
714 2 : if(initstride<0) {
715 2 : outhills="negativebias.dat";
716 : } else {
717 0 : outhills="negativebias_";
718 : }
719 : }
720 : }
721 :
722 15 : parse("OUTHILLS",outhills);
723 15 : parse("OUTHISTO",outhisto);
724 15 : outbasins="basins.dat";
725 15 : parse("OUTBASINS",outbasins);
726 15 : if(integratehills) {
727 14 : log<<" output file for fes/bias is : "<<outhills<<"\n";
728 : }
729 15 : if(integratehisto) {
730 1 : log<<" output file for histogram is : "<<outhisto<<"\n";
731 : }
732 15 : checkRead();
733 :
734 15 : log<<"\n";
735 15 : log<<" Now calculating...\n";
736 15 : log<<"\n";
737 :
738 : // here it defines the column to be histogrammed, tmpvalues should be only
739 : // the list of the collective variable one want to consider
740 15 : if(integratehisto) {
741 1 : checkFilesAreExisting(histoFiles);
742 2 : historep=Tools::make_unique<BiasRepresentation>(tmphistovalues,comm,gmin,gmax,gbin,histoSigma);
743 : }
744 :
745 : // decide how to source hills ( serial/parallel )
746 : // here below the input control
747 : // say how many hills and it will read them from the
748 : // bunch of files provided, will update the representation
749 : // of hills (i.e. a list of hills and the associated grid)
750 :
751 : // decide how to source colvars ( serial parallel )
752 15 : std::unique_ptr<FilesHandler> hillsHandler;
753 15 : std::unique_ptr<FilesHandler> histoHandler;
754 :
755 15 : if(integratehills) {
756 28 : hillsHandler=Tools::make_unique<FilesHandler>(hillsFiles,parallelread,*this, log, sortby);
757 : }
758 15 : if(integratehisto) {
759 2 : histoHandler=Tools::make_unique<FilesHandler>(histoFiles,parallelread,*this, log, sortby);
760 : }
761 :
762 : // Stopwatch is logged when it goes out of scope
763 15 : OFile basinsfile, histobasinsfile;
764 15 : if(basinmin.size()>0) {
765 1 : if(integratehills) {
766 1 : basinsfile.link(*this);
767 1 : basinsfile.open(outbasins);
768 1 : basinsfile.fmtField(fmt);
769 1 : log<<" output file for the free energies of the basins is : "<<outbasins<<"\n";
770 : }
771 1 : if(integratehisto) {
772 0 : histobasinsfile.link(*this);
773 0 : histobasinsfile.open("histo_"+outbasins);
774 0 : histobasinsfile.fmtField(fmt);
775 0 : log<<" output file for the free energies of the basins from the histogram is : histo_"<<outbasins<<"\n";
776 : }
777 : }
778 :
779 15 : Stopwatch sw(log);
780 :
781 : // Stopwatch is stopped when swh goes out of scope
782 15 : auto swh=sw.startStop("0 Summing hills");
783 :
784 : // read a number of hills and put in the bias representation
785 : int nfiles=0;
786 15 : bool ibias=integratehills;
787 15 : bool ihisto=integratehisto;
788 : while(true) {
789 24 : if( integratehills && ibias ) {
790 23 : if(nohistory) {
791 0 : biasrep->clear();
792 0 : log<<" clearing history before reading a new block\n";
793 : };
794 23 : log<<" reading hills: \n";
795 23 : const int before=biasrep->getNumberOfKernels();
796 23 : ibias=hillsHandler->readBunch(biasrep.get(),initstride) ;
797 23 : log<<"\n";
798 : // when the number of hills is a multiple of the stride the last block is empty: do not write it again
799 23 : if(!ibias && nfiles>0 && biasrep->getNumberOfKernels()==before) {
800 1 : integratehills=false;
801 : }
802 : }
803 :
804 24 : if( integratehisto && ihisto ) {
805 1 : if(nohistory) {
806 0 : historep->clear();
807 0 : log<<" clearing history before reading a new block\n";
808 : };
809 1 : log<<" reading histogram: \n";
810 1 : const int before=historep->getNumberOfKernels();
811 1 : ihisto=histoHandler->readBunch(historep.get(),initstride) ;
812 1 : log<<"\n";
813 1 : if(!ihisto && nfiles>0 && historep->getNumberOfKernels()==before) {
814 0 : integratehisto=false;
815 : }
816 : }
817 :
818 : // dump: need to project?
819 24 : if(proj.size()!=0) {
820 :
821 4 : if(integratehills) {
822 :
823 3 : log<<" Bias: Projecting on subgrid... \n";
824 3 : BiasWeight Bw(beta);
825 3 : Grid biasGrid=*(biasrep->getGridPtr());
826 3 : Grid smallGrid=biasGrid.project(proj,&Bw);
827 3 : OFile gridfile;
828 3 : gridfile.link(*this);
829 3 : std::ostringstream ostr;
830 3 : ostr<<nfiles;
831 : std::string myout;
832 3 : if(initstride>0) {
833 6 : myout=outhills+ostr.str()+".dat" ;
834 : } else {
835 : myout=outhills;
836 : }
837 3 : log<<" Bias: Writing subgrid on file "<<myout<<" \n";
838 3 : gridfile.open(myout);
839 3 : if(minTOzero) {
840 0 : smallGrid.setMinToZero();
841 : }
842 : smallGrid.setOutputFmt(fmt);
843 3 : smallGrid.writeToFile(gridfile);
844 3 : if(basinmin.size()>0) {
845 0 : writeBasins(basinsfile,smallGrid,biasrep->getLastTime(),biasrep->getNumberOfKernels());
846 : }
847 3 : gridfile.close();
848 3 : if(!ibias) {
849 2 : integratehills=false; // once you get to the final bunch just give up
850 : }
851 3 : }
852 : // this should be removed
853 4 : if(integratehisto) {
854 :
855 1 : log<<" Histo: Projecting on subgrid... \n";
856 1 : Grid histoGrid=*(historep->getGridPtr());
857 :
858 1 : OFile gridfile;
859 1 : gridfile.link(*this);
860 1 : std::ostringstream ostr;
861 1 : ostr<<nfiles;
862 : std::string myout;
863 1 : if(initstride>0) {
864 0 : myout=outhisto+ostr.str()+".dat" ;
865 : } else {
866 : myout=outhisto;
867 : }
868 1 : log<<" Histo: Writing subgrid on file "<<myout<<" \n";
869 1 : gridfile.open(myout);
870 :
871 1 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
872 1 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
873 1 : if(minTOzero) {
874 0 : histoGrid.setMinToZero();
875 : }
876 : histoGrid.setOutputFmt(fmt);
877 1 : histoGrid.writeToFile(gridfile);
878 1 : if(basinmin.size()>0) {
879 0 : writeBasins(histobasinsfile,histoGrid,historep->getLastTime(),historep->getNumberOfKernels());
880 : }
881 :
882 1 : if(!ihisto) {
883 1 : integratehisto=false; // once you get to the final bunch just give up
884 : }
885 1 : }
886 :
887 : } else {
888 :
889 20 : if(integratehills) {
890 :
891 19 : Grid biasGrid=*(biasrep->getGridPtr());
892 19 : biasGrid.scaleAllValuesAndDerivatives(-1.);
893 :
894 19 : OFile gridfile;
895 19 : gridfile.link(*this);
896 19 : std::ostringstream ostr;
897 19 : ostr<<nfiles;
898 : std::string myout;
899 19 : if(initstride>0) {
900 30 : myout=outhills+ostr.str()+".dat" ;
901 : } else {
902 : myout=outhills;
903 : }
904 19 : log<<" Writing full grid on file "<<myout<<" \n";
905 19 : gridfile.open(myout);
906 :
907 19 : if(minTOzero) {
908 0 : biasGrid.setMinToZero();
909 : }
910 : biasGrid.setOutputFmt(fmt);
911 19 : biasGrid.writeToFile(gridfile);
912 19 : if(basinmin.size()>0) {
913 3 : writeBasins(basinsfile,biasGrid,biasrep->getLastTime(),biasrep->getNumberOfKernels());
914 : }
915 : // rescale back prior to accumulate
916 19 : if(!ibias) {
917 11 : integratehills=false; // once you get to the final bunch just give up
918 : }
919 19 : }
920 20 : if(integratehisto) {
921 :
922 0 : Grid histoGrid=*(historep->getGridPtr());
923 : // do this if you want a free energy from a grid, otherwise do not
924 0 : histoGrid.applyFunctionAllValuesAndDerivatives(&mylog,&mylogder);
925 0 : histoGrid.scaleAllValuesAndDerivatives(-1./beta);
926 :
927 0 : OFile gridfile;
928 0 : gridfile.link(*this);
929 0 : std::ostringstream ostr;
930 0 : ostr<<nfiles;
931 : std::string myout;
932 0 : if(initstride>0) {
933 0 : myout=outhisto+ostr.str()+".dat" ;
934 : } else {
935 : myout=outhisto;
936 : }
937 0 : log<<" Writing full grid on file "<<myout<<" \n";
938 0 : gridfile.open(myout);
939 :
940 : // also this is useful only for free energy
941 0 : if(minTOzero) {
942 0 : histoGrid.setMinToZero();
943 : }
944 : histoGrid.setOutputFmt(fmt);
945 0 : histoGrid.writeToFile(gridfile);
946 0 : if(basinmin.size()>0) {
947 0 : writeBasins(histobasinsfile,histoGrid,historep->getLastTime(),historep->getNumberOfKernels());
948 : }
949 :
950 0 : if(!ihisto) {
951 0 : integratehisto=false; // once you get to the final bunch just give up
952 : }
953 0 : }
954 : }
955 24 : if ( !ibias && !ihisto) {
956 : break; //when both are over then just quit
957 : }
958 :
959 9 : nfiles++;
960 9 : }
961 :
962 : return;
963 15 : }
964 : // just an initialization but you need to do something on the fly?: need to connect with a metad run and its grid representation
965 : // your argument is a metad run
966 : // if the grid does not exist crash and say that you need some data
967 : // otherwise just link with it
968 :
969 15 : }
970 :
971 3 : void FuncSumHills::writeBasins(OFile & ofile, const Grid & fes, double time, int nkernels) {
972 3 : const unsigned nb=basinmin.size();
973 3 : const unsigned nd=fes.getDimension();
974 3 : const auto pbc=fes.getIsPeriodic();
975 : // -kT log sum_s exp(-F(s)/kT) over the points of each basin, computed with a running maximum
976 3 : std::vector<double> fmax(nb,-std::numeric_limits<double>::max()), sum(nb,0.0);
977 3 : std::vector<double> x(nd);
978 10251 : for(Grid::index_t i=0; i<fes.getSize(); i++) {
979 10248 : const double mbf=-beta*fes.getValue(i);
980 10248 : if(!std::isfinite(mbf)) {
981 0 : continue;
982 : }
983 10248 : fes.getPoint(i,x);
984 40992 : for(unsigned b=0; b<nb; b++) {
985 : bool inside=true;
986 74592 : for(unsigned d=0; d<nd && inside; d++) {
987 43848 : const double lo=basinmin[b][d];
988 43848 : const double hi=basinmax[b][d];
989 43848 : if(pbc[d] && lo>hi) {
990 : // along a periodic variable a basin can cross the boundary
991 0 : inside=(x[d]>=lo || x[d]<=hi);
992 : } else {
993 43848 : inside=(x[d]>=lo && x[d]<=hi);
994 : }
995 : }
996 30744 : if(!inside) {
997 20208 : continue;
998 : }
999 10536 : if(mbf>fmax[b]) {
1000 349 : sum[b]=sum[b]*std::exp(fmax[b]-mbf)+1.0;
1001 349 : fmax[b]=mbf;
1002 : } else {
1003 10187 : sum[b]+=std::exp(mbf-fmax[b]);
1004 : }
1005 : }
1006 : }
1007 3 : std::vector<double> f(nb);
1008 12 : for(unsigned b=0; b<nb; b++) {
1009 9 : plumed_massert(sum[b]>0.0,"basin "+std::to_string(b)+" does not contain any point of the grid");
1010 9 : f[b]=-(std::log(sum[b])+fmax[b])/beta;
1011 : }
1012 3 : ofile.printField("time",time);
1013 3 : ofile.printField("nkernels",nkernels);
1014 9 : for(unsigned b=1; b<nb; b++) {
1015 12 : ofile.printField("deltaf_"+std::to_string(b),f[b]-f[0]);
1016 : }
1017 3 : ofile.printField();
1018 3 : }
1019 :
1020 0 : void FuncSumHills::calculate() {
1021 : // this should be connected only with a grid representation to metadynamics
1022 : // at regular time just dump it
1023 0 : plumed_merror("You should have never got here: this stuff is not yet implemented!");
1024 : }
1025 :
1026 15 : bool FuncSumHills::checkFilesAreExisting(const std::vector<std::string> & hills ) {
1027 15 : plumed_massert(hills.size()!=0,"the number of files provided should be at least one" );
1028 : auto ifile=Tools::make_unique<IFile>();
1029 15 : ifile->link(*this);
1030 33 : for(unsigned i=0; i< hills.size(); i++) {
1031 18 : plumed_massert(ifile->FileExist(hills[i]),"missing file "+hills[i]);
1032 : }
1033 15 : return true;
1034 :
1035 15 : }
1036 :
1037 : }
1038 :
1039 : }
1040 :
1041 :
|
