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 "CLTool.h"
23 : #include "core/CLToolRegister.h"
24 : #include "tools/Tools.h"
25 : #include "core/Action.h"
26 : #include "core/ActionRegister.h"
27 : #include "core/PlumedMain.h"
28 : #include "tools/Communicator.h"
29 : #include "tools/Random.h"
30 : #include <cstdio>
31 : #include <string>
32 : #include <vector>
33 : #include <iostream>
34 : #include "tools/File.h"
35 : #include "core/Value.h"
36 : #include "tools/Matrix.h"
37 :
38 : namespace PLMD {
39 : namespace cltools {
40 :
41 : //+PLUMEDOC TOOLS sum_hills
42 : /*
43 : sum_hills is a tool that allows one to to use plumed to post-process an existing hills/colvar file
44 :
45 : This tool is most frequently used as shown below after a [METAD](METAD.md) simulation has been performed to sum all the Gaussians in the hills file
46 : and output the free energy surface.
47 :
48 : ```plumed
49 : plumed sum_hills --hills PATHTOMYHILLSFILE
50 : ```
51 :
52 : The default name for the output file will be `fes.dat`
53 : Note that plumed automatically detects the
54 : number of the variables you have and their periodicity.
55 : Additionally, if you use flexible hills (multivariate Gaussian kernels), plumed will understand it from the HILLS file.
56 :
57 : The sum_hills tool can also accept multiple files that will be integrated one after the other as shown below
58 :
59 : ```plumed
60 : plumed sum_hills --hills PATHTOMYHILLSFILE1,PATHTOMYHILLSFILE2,PATHTOMYHILLSFILE3
61 : ```
62 :
63 : if you want to integrate out some variable from your output free energy surface you do
64 :
65 : ```plumed
66 : plumed sum_hills --hills PATHTOMYHILLSFILE --idw t1 --kt 0.6
67 : ```
68 :
69 : where with `--idw` you define the variables that you want in the output free energy surface.
70 : All other variables in the hills file will be integrated out. `--kt` defines the temperature of the system in energy units.
71 : (be consistent with the units you have in your hills: plumed will not check this for you)
72 : If you need more variables then you may use a comma separated syntax shown below:
73 :
74 : ```plumed
75 : plumed sum_hills --hills PATHTOMYHILLSFILE --idw t1,t2 --kt 0.6
76 : ```
77 :
78 : You can define the output grid with the number of bins you want by using a command like this one
79 :
80 : ```plumed
81 : plumed sum_hills --bin 99,99 --hills PATHTOMYHILLSFILE
82 : ```
83 :
84 : In the command above the min/max to use are detected for you. If you want to specify the min/max yourself you can do so by using
85 : a command like the one shown below:
86 :
87 : ```plumed
88 : plumed sum_hills --bin 99,99 --min -pi,-pi --max pi,pi --hills PATHTOMYHILLSFILE
89 : ```
90 :
91 : You can of course use numbers instead of -pi/pi.
92 :
93 : You can use a `--stride` keyword to have a dump of the free energy estimate after each bunch of hills you read as shown in the following command:
94 :
95 : ```plumed
96 : plumed sum_hills --stride 300 --hills PATHTOMYHILLSFILE
97 : ```
98 :
99 : If you have run well tempered metadynamics you can also ask `sum_hills` to output the negative bias instead of the free energy by using the keyword `--negbias`
100 : as shown below
101 :
102 : ```plumed
103 : plumed sum_hills --negbias --hills PATHTOMYHILLSFILE
104 : ```
105 :
106 : Here the default output file name will be negativebias.dat
107 :
108 : From time to time you might need to use HILLS or a COLVAR file
109 : as it was just a simple set of points from which you want to build
110 : a free energy by using -(1/beta)log(P). If you want to do this operation then
111 : you use the `--histo` option as shown below:
112 :
113 : ```plumed
114 : plumed sum_hills --histo PATHTOMYCOLVARORHILLSFILE --sigma 0.2,0.2 --kt 0.6
115 : ```
116 :
117 : in this case you need a `--kt` to do the reweighting and you
118 : need to set bandwidth parameters using the `--sigma` option for the histogram calculation as this histogram is computed using
119 : Gaussian kernels, so it will be a continuous histogram.
120 :
121 : For the command above the default output is called histo.dat.
122 : Note that also here you can have multiple input files separated by a comma.
123 :
124 : Additionally, if you want to compute the histogram and sum of the hills from the same file you can do this by using a command like this one:
125 :
126 : ```plumed
127 : plumed sum_hills --hills --histo PATHTOMYCOLVARORHILLSFILE --sigma 0.2,0.2 --kt 0.6
128 : ```
129 :
130 : The two files can be eventually the same
131 :
132 : Another interesting thing one can do is monitor the difference in blocks as a metadynamics goes on.
133 : When the bias deposited is constant over the whole domain one can consider to be at convergence.
134 : This can be done with the `--nohistory` keyword
135 :
136 : ```plumed
137 : plumed sum_hills --stride 300 --hills PATHTOMYHILLSFILE --nohistory
138 : ```
139 :
140 : and similarly one can do the same for an histogram file
141 :
142 : ```plumed
143 : plumed sum_hills --histo PATHTOMYCOLVARORHILLSFILE --sigma 0.2,0.2 --kt 0.6 --nohistory
144 : ```
145 :
146 : just to check the hypothetical free energy calculated in single blocks of time during a simulation
147 : and not in a cumulative way
148 :
149 : The output format can be controlled by using the `--fmt` option as shown below
150 :
151 : ```plumed
152 : plumed sum_hills --hills PATHTOMYHILLSFILE --fmt %8.3f
153 : ```
154 :
155 : where here we chose a float with length of 8 and 3 digits
156 :
157 : The output can be named in a arbitrary way by using the `--output` option as shown below:
158 :
159 : ```plumed
160 : plumed sum_hills --hills PATHTOMYHILLSFILE --outfile myfes.dat
161 : ```
162 :
163 : This command produces a file myfes.dat which contains the free energy surface.
164 :
165 : If you use stride, the name specied with `--output` is the suffix so this command:
166 :
167 : ```plumed
168 : plumed sum_hills --hills PATHTOMYHILLSFILE --outfile myfes_ --stride 100
169 : ```
170 :
171 : produces output files called myfes_0.dat, myfes_1.dat, myfes_2.dat etc.
172 :
173 : The same is true for the output coming from histogram that is used. So this example
174 :
175 : ```plumed
176 : plumed sum_hills --histo HILLS --kt 2.5 --sigma 0.01 --outhisto myhisto.dat
177 : ```
178 :
179 : produces a file myhisto.dat, while this input
180 :
181 : ```plumed
182 : plumed sum_hills --histo HILLS --kt 2.5 --sigma 0.01 --outhisto myhisto_ --stride 100
183 : ```
184 :
185 : produces output files called `myhisto_0.dat`, `myhisto_1.dat`, `myhisto_3.dat` etc..
186 :
187 : */
188 : //+ENDPLUMEDOC
189 :
190 : class CLToolSumHills : public CLTool {
191 : public:
192 : static void registerKeywords( Keywords& keys );
193 : explicit CLToolSumHills(const CLToolOptions& co );
194 : int main(FILE* in,FILE*out,Communicator& pc) override;
195 : std::string description()const override;
196 : /// find a list of variables present, if they are periodic and which is the period
197 : /// return false if the file does not exist
198 : static bool findCvsAndPeriodic(const std::string & filename, std::vector< std::vector <std::string> > &cvs,std::vector<std::string> &pmin,std::vector<std::string> &pmax, bool &multivariate, std::string &lowI_, std::string &uppI_);
199 : };
200 :
201 5418 : void CLToolSumHills::registerKeywords( Keywords& keys ) {
202 5418 : CLTool::registerKeywords( keys );
203 5418 : keys.addFlag("--help-debug",false,"print special options that can be used to create regtests");
204 5418 : keys.add("optional","--hills","specify the name of the hills file");
205 5418 : keys.add("optional","--histo","specify the name of the file for histogram a colvar/hills file is good");
206 5418 : keys.add("optional","--stride","specify the stride for integrating hills file (default 0=never)");
207 5418 : keys.add("optional","--min","the lower bounds for the grid");
208 5418 : keys.add("optional","--max","the upper bounds for the grid");
209 5418 : keys.add("optional","--bin","the number of bins for the grid");
210 5418 : keys.add("optional","--spacing","grid spacing, alternative to the number of bins");
211 5418 : keys.add("optional","--idw","specify the variables to be used for the free-energy/histogram (default is all). With --hills the other variables will be integrated out, with --histo the other variables won't be considered");
212 5418 : keys.add("optional","--outfile","specify the output file for sumhills");
213 5418 : keys.add("optional","--outhisto","specify the output file for the histogram");
214 5418 : keys.add("optional","--kt","specify temperature in energy units for integrating out variables");
215 5418 : keys.add("optional","--sigma"," a vector that specify the sigma for binning (only needed when doing histogram ");
216 5418 : keys.addFlag("--negbias",false," print the negative bias instead of the free energy (only needed with well tempered runs and flexible hills) ");
217 5418 : keys.addFlag("--nohistory",false," to be used with --stride: it splits the bias/histogram in pieces without previous history ");
218 5418 : keys.addFlag("--mintozero",false," it translate all the minimum value in bias/histogram to zero (useful to compare results) ");
219 5418 : keys.add("optional","--fmt","specify the output format");
220 5418 : }
221 :
222 14 : CLToolSumHills::CLToolSumHills(const CLToolOptions& co ):
223 14 : CLTool(co) {
224 14 : inputdata=commandline;
225 14 : }
226 :
227 5 : std::string CLToolSumHills::description()const {
228 5 : return "sum the hills with plumed";
229 : }
230 :
231 9 : int CLToolSumHills::main(FILE* in,FILE*out,Communicator& pc) {
232 :
233 : // Read the hills input file name
234 : std::vector<std::string> hillsFiles;
235 : bool dohills;
236 18 : dohills=parseVector("--hills",hillsFiles);
237 : // Read the histogram file
238 : std::vector<std::string> histoFiles;
239 : bool dohisto;
240 9 : dohisto=parseVector("--histo",histoFiles);
241 :
242 9 : plumed_massert(dohisto || dohills,"you should use --histo or/and --hills command");
243 :
244 : std::vector< std::vector<std::string> > vcvs;
245 : std::vector<std::string> vpmin;
246 : std::vector<std::string> vpmax;
247 : std::string lowI_, uppI_;
248 9 : if(dohills) {
249 : // parse it as it was a restart
250 : bool vmultivariate;
251 8 : findCvsAndPeriodic(hillsFiles[0], vcvs, vpmin, vpmax, vmultivariate, lowI_, uppI_);
252 : }
253 :
254 : std::vector< std::vector<std::string> > hcvs;
255 : std::vector<std::string> hpmin;
256 : std::vector<std::string> hpmax;
257 :
258 : std::vector<std::string> sigma;
259 9 : if(dohisto) {
260 : bool hmultivariate;
261 1 : findCvsAndPeriodic(histoFiles[0], hcvs, hpmin, hpmax, hmultivariate, lowI_, uppI_);
262 : // here need also the vector of sigmas
263 2 : parseVector("--sigma",sigma);
264 1 : if(sigma.size()==0) {
265 0 : plumed_merror("you should define --sigma vector when using histogram");
266 : }
267 : lowI_=uppI_="-1."; // Interval is not use for histograms
268 : }
269 :
270 9 : if(dohisto && dohills) {
271 0 : plumed_massert(vcvs==hcvs,"variables for histogram and bias should have the same labels");
272 0 : plumed_massert(hpmin==vpmin,"variables for histogram and bias should have the same min for periodicity");
273 0 : plumed_massert(hpmax==vpmax,"variables for histogram and bias should have the same max for periodicity");
274 : }
275 :
276 : // now put into a neutral vector
277 :
278 : std::vector< std::vector<std::string> > cvs;
279 : std::vector<std::string> pmin;
280 : std::vector<std::string> pmax;
281 :
282 9 : if(dohills) {
283 8 : cvs=vcvs;
284 8 : pmin=vpmin;
285 8 : pmax=vpmax;
286 : }
287 9 : if(dohisto) {
288 1 : cvs=hcvs;
289 1 : pmin=hpmin;
290 1 : pmax=hpmax;
291 : }
292 :
293 :
294 : // setup grids
295 : unsigned grid_check=0;
296 9 : std::vector<std::string> gmin(cvs.size());
297 18 : if(parseVector("--min",gmin)) {
298 4 : if(gmin.size()!=cvs.size() && gmin.size()!=0) {
299 0 : plumed_merror("not enough values for --min");
300 : }
301 : grid_check++;
302 : }
303 9 : std::vector<std::string> gmax(cvs.size() );
304 18 : if(parseVector("--max",gmax)) {
305 4 : if(gmax.size()!=cvs.size() && gmax.size()!=0) {
306 0 : plumed_merror("not enough values for --max");
307 : }
308 4 : grid_check++;
309 : }
310 9 : std::vector<std::string> gbin(cvs.size());
311 : bool grid_has_bin;
312 : grid_has_bin=false;
313 18 : if(parseVector("--bin",gbin)) {
314 5 : if(gbin.size()!=cvs.size() && gbin.size()!=0) {
315 0 : plumed_merror("not enough values for --bin");
316 : }
317 : grid_has_bin=true;
318 : }
319 9 : std::vector<std::string> gspacing(cvs.size());
320 : bool grid_has_spacing;
321 : grid_has_spacing=false;
322 18 : if(parseVector("--spacing",gspacing)) {
323 1 : if(gspacing.size()!=cvs.size() && gspacing.size()!=0) {
324 0 : plumed_merror("not enough values for --spacing");
325 : }
326 : grid_has_spacing=true;
327 : }
328 : // allowed: no grids only bin
329 : // not allowed: partial grid definition
330 9 : plumed_massert( gmin.size()==gmax.size() && (gmin.size()==0 || gmin.size()==cvs.size() ),"you should specify --min and --max together with same number of components");
331 :
332 :
333 :
334 9 : PlumedMain plumed;
335 : std::string ss;
336 9 : unsigned nn=1;
337 : ss="setNatoms";
338 9 : plumed.cmd(ss,&nn);
339 9 : if(Communicator::initialized()) {
340 0 : plumed.cmd("setMPIComm",&pc.Get_comm());
341 : }
342 9 : plumed.cmd("init",&nn);
343 9 : std::vector <bool> isdone(cvs.size(),false);
344 26 : for(unsigned i=0; i<cvs.size(); i++) {
345 17 : if(!isdone[i]) {
346 : isdone[i]=true;
347 : std::vector<std::string> actioninput;
348 : std::vector <unsigned> inds;
349 17 : actioninput.push_back("FAKE");
350 34 : actioninput.push_back("ATOMS=1");
351 34 : actioninput.push_back("LABEL="+cvs[i][0]);
352 : std::vector<std::string> comps, periods;
353 17 : if(cvs[i].size()>1) {
354 1 : comps.push_back(cvs[i][1]);
355 1 : inds.push_back(i);
356 : }
357 17 : periods.push_back(pmin[i]);
358 17 : periods.push_back(pmax[i]);
359 25 : for(unsigned j=i+1; j<cvs.size(); j++) {
360 8 : if(cvs[i][0]==cvs[j][0] && !isdone[j]) {
361 0 : if(cvs[i].size()==1 || cvs[j].size()==1 ) {
362 0 : plumed_merror("you cannot have twice the same label and no components ");
363 : }
364 0 : if(cvs[j].size()>1) {
365 0 : comps.push_back(cvs[j][1]);
366 0 : periods.push_back(pmin[j]);
367 0 : periods.push_back(pmax[j]);
368 : isdone[j]=true;
369 0 : inds.push_back(j);
370 : }
371 : }
372 :
373 : }
374 : // drain all the components
375 : std::string addme;
376 17 : if(comps.size()>0) {
377 : addme="COMPONENTS=";
378 1 : for(unsigned i=0; i<comps.size()-1; i++) {
379 0 : addme+=comps[i]+",";
380 : }
381 : addme+=comps.back();
382 1 : actioninput.push_back(addme);
383 : }
384 : // periodicity (always explicit here)
385 : addme="PERIODIC=";
386 34 : for(unsigned j=0; j<periods.size()-1; j++) {
387 34 : addme+=periods[j]+",";
388 : }
389 : addme+=periods.back();
390 17 : actioninput.push_back(addme);
391 18 : for(unsigned j=0; j<inds.size(); j++) {
392 : unsigned jj;
393 1 : jj=inds[j];
394 1 : if(grid_check==2) {
395 : double gm;
396 : double pm;
397 1 : if(pmin[jj]!="none") {
398 1 : Tools::convert(gmin[jj],gm);
399 1 : Tools::convert(pmin[jj],pm);
400 1 : if( gm<pm ) {
401 0 : plumed_merror("Periodicity issue : GRID_MIN value ( "+gmin[jj]+" ) is less than periodicity in HILLS file in "+cvs[jj][0]+ " ( "+pmin[jj]+" ) ");
402 : }
403 : }
404 1 : if(pmax[jj]!="none") {
405 1 : Tools::convert(gmax[jj],gm);
406 1 : Tools::convert(pmax[jj],pm);
407 1 : if( gm>pm ) {
408 0 : plumed_merror("Periodicity issue : GRID_MAX value ( "+gmax[jj]+" ) is more than periodicity in HILLS file in "+cvs[jj][0]+ " ( "+pmax[jj]+" ) ");
409 : }
410 : }
411 : }
412 : }
413 :
414 : // for(unsigned i=0;i< actioninput.size();i++){
415 : // cerr<<"AA "<<actioninput[i]<<endl;
416 : // }
417 17 : plumed.readInputWords(actioninput,false);
418 34 : }
419 :
420 : }
421 9 : unsigned ncv=cvs.size();
422 : std::vector<std::string> actioninput;
423 : std::vector<std::string> idw;
424 : // check if the variables to be used are correct
425 18 : if(parseVector("--idw",idw)) {
426 4 : for(unsigned i=0; i<idw.size(); i++) {
427 : bool found=false;
428 6 : for(unsigned j=0; j<cvs.size(); j++) {
429 4 : if(cvs[j].size()>1) {
430 0 : if(idw[i]==cvs[j][0]+"."+cvs[j][1]) {
431 : found=true;
432 : }
433 : } else {
434 4 : if(idw[i]==cvs[j][0]) {
435 : found=true;
436 : }
437 : }
438 : }
439 2 : if(!found) {
440 0 : plumed_merror("variable "+idw[i]+" is not found in the bunch of cvs: revise your --idw option" );
441 : }
442 : }
443 2 : plumed_massert( idw.size()<=cvs.size(),"the number of variables to be integrated should be at most equal to the total number of cvs ");
444 : // in this case you need a beta factor!
445 : }
446 :
447 : std::string kt;
448 9 : kt=std::string("1.");// assign an arbitrary value just in case that idw.size()==cvs.size()
449 9 : if ( dohisto || idw.size()!=0 ) {
450 6 : plumed_massert(parse("--kt",kt),"if you make a dimensionality reduction (--idw) or a histogram (--histo) then you need to define --kt ");
451 : }
452 :
453 : std::string addme;
454 :
455 9 : actioninput.push_back("FUNCSUMHILLS");
456 18 : actioninput.push_back("ISCLTOOL");
457 :
458 : // set names
459 : std::string outfile;
460 18 : if(parse("--outfile",outfile)) {
461 0 : actioninput.push_back("OUTHILLS="+outfile);
462 : }
463 : std::string outhisto;
464 18 : if(parse("--outhisto",outhisto)) {
465 2 : actioninput.push_back("OUTHISTO="+outhisto);
466 : }
467 :
468 :
469 : addme="ARG=";
470 17 : for(unsigned i=0; i<(ncv-1); i++) {
471 8 : if(cvs[i].size()==1) {
472 16 : addme+=std::string(cvs[i][0])+",";
473 : } else {
474 0 : addme+=std::string(cvs[i][0])+"."+std::string(cvs[i][1])+",";
475 : }
476 : }
477 9 : if(cvs[ncv-1].size()==1) {
478 16 : addme+=std::string(cvs[ncv-1][0]);
479 : } else {
480 2 : addme+=std::string(cvs[ncv-1][0])+"."+std::string(cvs[ncv-1][1]);
481 : }
482 9 : actioninput.push_back(addme);
483 : //for(unsigned i=0;i< actioninput.size();i++){
484 : // cerr<<"AA "<<actioninput[i]<<endl;
485 : //}
486 9 : if(dohills) {
487 : addme="HILLSFILES=";
488 9 : for(unsigned i=0; i<hillsFiles.size()-1; i++) {
489 2 : addme+=hillsFiles[i]+",";
490 : }
491 : addme+=hillsFiles[hillsFiles.size()-1];
492 8 : actioninput.push_back(addme);
493 : // set the grid
494 : }
495 9 : if(grid_check==2) {
496 : addme="GRID_MAX=";
497 7 : for(unsigned i=0; i<(ncv-1); i++) {
498 6 : addme+=gmax[i]+",";
499 : }
500 : addme+=gmax[ncv-1];
501 4 : actioninput.push_back(addme);
502 : addme="GRID_MIN=";
503 7 : for(unsigned i=0; i<(ncv-1); i++) {
504 6 : addme+=gmin[i]+",";
505 : }
506 : addme+=gmin[ncv-1];
507 4 : actioninput.push_back(addme);
508 : }
509 9 : if(grid_has_bin) {
510 : addme="GRID_BIN=";
511 10 : for(unsigned i=0; i<(ncv-1); i++) {
512 10 : addme+=gbin[i]+",";
513 : }
514 : addme+=gbin[ncv-1];
515 5 : actioninput.push_back(addme);
516 : }
517 9 : if(grid_has_spacing) {
518 : addme="GRID_SPACING=";
519 1 : for(unsigned i=0; i<(ncv-1); i++) {
520 0 : addme+=gspacing[i]+",";
521 : }
522 : addme+=gspacing[ncv-1];
523 1 : actioninput.push_back(addme);
524 : }
525 : std::string stride;
526 : stride="";
527 18 : if(parse("--stride",stride)) {
528 1 : actioninput.push_back("INITSTRIDE="+stride);
529 : bool nohistory;
530 1 : parseFlag("--nohistory",nohistory);
531 1 : if(nohistory) {
532 0 : actioninput.push_back("NOHISTORY");
533 : }
534 : }
535 : bool mintozero;
536 9 : parseFlag("--mintozero",mintozero);
537 9 : if(mintozero) {
538 0 : actioninput.push_back("MINTOZERO");
539 : }
540 9 : if(idw.size()!=0) {
541 : addme="PROJ=";
542 2 : for(unsigned i=0; i<idw.size()-1; i++) {
543 0 : addme+=idw[i]+",";
544 : }
545 : addme+=idw.back();
546 2 : actioninput.push_back(addme);
547 : }
548 :
549 9 : if(dohisto) {
550 1 : if(idw.size()==0) {
551 1 : if(sigma.size()!=hcvs.size()) {
552 0 : plumed_merror("you should define as many --sigma vector as the number of collective variable used for the histogram ");
553 : }
554 : } else {
555 0 : if(idw.size()!=sigma.size()) {
556 0 : plumed_merror("you should define as many --sigma vector as the number of collective variable used for the histogram ");
557 : }
558 : }
559 : }
560 :
561 9 : if(idw.size()!=0 || dohisto) {
562 6 : actioninput.push_back("KT="+kt);
563 : }
564 9 : if(dohisto) {
565 : addme="HISTOFILES=";
566 1 : for(unsigned i=0; i<histoFiles.size()-1; i++) {
567 0 : addme+=histoFiles[i]+",";
568 : }
569 : addme+=histoFiles[histoFiles.size()-1];
570 1 : actioninput.push_back(addme);
571 :
572 : addme="HISTOSIGMA=";
573 2 : for(unsigned i=0; i<sigma.size()-1; i++) {
574 2 : addme+=sigma[i]+",";
575 : }
576 : addme+=sigma.back();
577 1 : actioninput.push_back(addme);
578 : }
579 :
580 : bool negbias;
581 9 : parseFlag("--negbias",negbias);
582 9 : if(negbias) {
583 2 : actioninput.push_back("NEGBIAS");
584 : }
585 :
586 9 : if(lowI_!=uppI_) {
587 : addme="INTERVAL=";
588 0 : addme+=lowI_+",";
589 : addme+=uppI_;
590 0 : actioninput.push_back(addme);
591 : }
592 :
593 : std::string fmt;
594 : fmt="";
595 18 : parse("--fmt",fmt);
596 9 : if(fmt!="") {
597 18 : actioninput.push_back("FMT="+fmt);
598 : }
599 :
600 :
601 : // for(unsigned i=0;i< actioninput.size();i++){
602 : // cerr<<"AA "<<actioninput[i]<<endl;
603 : // }
604 9 : plumed.readInputWords(actioninput,false);
605 : // if not a grid, then set it up automatically
606 9 : return 0;
607 27 : }
608 :
609 9 : bool CLToolSumHills::findCvsAndPeriodic(const std::string & filename, std::vector< std::vector<std::string> > &cvs, std::vector<std::string> &pmin,std::vector<std::string> &pmax, bool &multivariate, std::string &lowI_, std::string &uppI_) {
610 9 : IFile ifile;
611 9 : ifile.allowIgnoredFields();
612 : std::vector<std::string> fields;
613 9 : if(ifile.FileExist(filename)) {
614 : cvs.clear();
615 : pmin.clear();
616 : pmax.clear();
617 9 : ifile.open(filename);
618 9 : ifile.scanFieldList(fields);
619 : bool before_sigma=true;
620 121 : for(unsigned i=0; i<fields.size(); i++) {
621 : size_t pos = 0;
622 : size_t founds,foundm,foundp;
623 : //found=(fields[i].find("sigma_", pos) || fields[i].find("min_", pos) || fields[i].find("max_", pos) ) ;
624 112 : founds=fields[i].find("sigma_", pos) ;
625 112 : foundm=fields[i].find("min_", pos) ;
626 112 : foundp=fields[i].find("max_", pos) ;
627 112 : if (founds!=std::string::npos || foundm!=std::string::npos || foundp!=std::string::npos ) {
628 : before_sigma=false;
629 : }
630 : // cvs are after time and before sigmas
631 : size_t found;
632 112 : found=fields[i].find("time", pos);
633 112 : if( found==std::string::npos && before_sigma) {
634 : // separate the components
635 : size_t dot=fields[i].find_first_of('.');
636 : std::vector<std::string> ss;
637 : // this loop does not take into account repetitions
638 17 : if(dot!=std::string::npos) {
639 1 : std::string a=fields[i].substr(0,dot);
640 1 : std::string name=fields[i].substr(dot+1);
641 1 : ss.push_back(a);
642 1 : ss.push_back(name);
643 1 : cvs.push_back(ss);
644 : } else {
645 : std::vector<std::string> ss;
646 16 : ss.push_back(fields[i]);
647 16 : cvs.push_back(ss);
648 16 : }
649 : //std::cerr<<"found variable number "<<cvs.size()<<" : "<<cvs.back()[0]<<std::endl;
650 : //if((cvs.back()).size()!=1){
651 : // std::cerr<<"component "<<(cvs.back()).back()<<std::endl;
652 : //}
653 : // get periodicity
654 17 : pmin.push_back("none");
655 34 : pmax.push_back("none");
656 : std::string mm;
657 17 : if((cvs.back()).size()>1) {
658 2 : mm=cvs.back()[0]+"."+cvs.back()[1];
659 : } else {
660 : mm=cvs.back()[0];
661 : }
662 34 : if(ifile.FieldExist("min_"+mm)) {
663 : std::string val;
664 28 : ifile.scanField("min_"+mm,val);
665 : pmin[pmin.size()-1]=val;
666 : // std::cerr<<"found min : "<<pmin.back()<<std::endl;
667 : }
668 : //std::cerr<<"found min : "<<pmin.back()<<std::endl;
669 34 : if(ifile.FieldExist("max_"+mm)) {
670 : std::string val;
671 28 : ifile.scanField("max_"+mm,val);
672 : pmax[pmax.size()-1]=val;
673 : // std::cerr<<"found max : "<<pmax.back()<<std::endl;
674 : }
675 : //std::cerr<<"found max : "<<pmax.back()<<std::endl;
676 17 : }
677 : }
678 : // is multivariate ???
679 : std::string sss;
680 9 : multivariate=false;
681 18 : if(ifile.FieldExist("multivariate")) {
682 : ;
683 18 : ifile.scanField("multivariate",sss);
684 9 : if(sss=="true") {
685 6 : multivariate=true;
686 3 : } else if(sss=="false") {
687 3 : multivariate=false;
688 : }
689 : }
690 : // do interval?
691 18 : if(ifile.FieldExist("lower_int")) {
692 0 : ifile.scanField("lower_int",lowI_);
693 0 : ifile.scanField("upper_int",uppI_);
694 : } else {
695 : lowI_="-1.";
696 : uppI_="-1.";
697 : }
698 9 : ifile.scanField();
699 : return true;
700 : } else {
701 : return false;
702 : }
703 9 : }
704 :
705 :
706 16268 : PLUMED_REGISTER_CLTOOL(CLToolSumHills,"sum_hills")
707 :
708 :
709 :
710 : }
711 : }
|
