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 <algorithm>
31 : #include <cstdio>
32 : #include <string>
33 : #include <vector>
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 sums the Gaussians in a hills file and writes the resulting free energy on a grid.
44 :
45 : This tool is most frequently used after a [METAD](METAD.md) simulation 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 is `fes.dat`. It contains, for each point of the grid, the free energy
53 : and its derivatives with respect to the variables. The free energy is minus the sum of the Gaussians: with
54 : well-tempered metadynamics the heights written in the hills file are already scaled by $\gamma/(\gamma-1)$
55 : (where $\gamma$ is the bias factor), so that their sum is directly the estimate of minus the free energy.
56 : As always, the free energy is defined up to an additive constant; use `--mintozero` to shift its minimum to zero.
57 :
58 : sum_hills automatically detects the variables in the hills file and their periodicity. A column `X` of the file
59 : is a variable if the file also contains its width (`sigma_X`, or `sigma_X_X` for multivariate Gaussians),
60 : whatever the name of the variable. If you used flexible hills (multivariate Gaussian kernels), sum_hills
61 : understands it from the hills file.
62 :
63 : ## Choosing the grid
64 :
65 : By default the grid is chosen so as to contain all the Gaussians, with a spacing of about a third of the smallest width.
66 : You can instead choose the number of bins:
67 :
68 : ```plumed
69 : plumed sum_hills --bin 99,99 --hills PATHTOMYHILLSFILE
70 : ```
71 :
72 : In the command above the minimum and maximum of the grid are still detected for you. If you want to set them yourself
73 : you can do so as shown below:
74 :
75 : ```plumed
76 : plumed sum_hills --bin 99,99 --min -pi,-pi --max pi,pi --hills PATHTOMYHILLSFILE
77 : ```
78 :
79 : You can of course use numbers instead of -pi/pi. `--min` and `--max` must be given together. Instead of the number of bins
80 : you can give the grid spacing with `--spacing`, with or without `--min` and `--max`. If you give both `--bin` and `--spacing`,
81 : the larger number of bins is used for each variable.
82 :
83 : ## Multiple files and multiple walkers
84 :
85 : sum_hills can also accept multiple files:
86 :
87 : ```plumed
88 : plumed sum_hills --hills PATHTOMYHILLSFILE1,PATHTOMYHILLSFILE2,PATHTOMYHILLSFILE3
89 : ```
90 :
91 : The files are integrated one after the other. This does not change the final free energy, but when
92 : you use `--stride` (see below) with the hills files of multiple walkers the intermediate free energies would
93 : not follow the simulation. With `--sort time` the hills of all the files are merged in order of their time
94 : (hills with the same time are taken in the order of the files), so that each intermediate free energy
95 : contains the hills deposited by all the walkers up to a given time:
96 :
97 : ```plumed
98 : plumed sum_hills --hills HILLS.0,HILLS.1,HILLS.2 --sort time --stride 300
99 : ```
100 :
101 : With `--sort clock` the hills are merged in order of the `clock` field (the wall-clock time at which they were
102 : deposited, written with multiple walkers), which is more appropriate when the walkers run at different speeds.
103 : Each file should be ordered by the chosen field.
104 :
105 : ## Integrating out variables
106 :
107 : If you want to integrate out some variables from your free energy surface you can use `--idw`:
108 :
109 : ```plumed
110 : plumed sum_hills --hills PATHTOMYHILLSFILE --idw t1 --kt 0.6
111 : ```
112 :
113 : With `--idw` you list the variables that you want in the output free energy surface.
114 : All the other variables in the hills file are integrated out, i.e.
115 : $F(s_1) = -k_BT \log \sum_{s_2} e^{-F(s_1,s_2)/k_BT}$, where the sum runs over the points of the grid.
116 : `--kt` is the temperature of the system in energy units: be consistent with the units of your hills file,
117 : since plumed will not check this for you. If you need more variables you can list them separated by commas:
118 :
119 : ```plumed
120 : plumed sum_hills --hills PATHTOMYHILLSFILE --idw t1,t2 --kt 0.6
121 : ```
122 :
123 : ## Following the convergence
124 :
125 : With `--stride` you can write the free energy after each bunch of hills (the stride is a number of hills):
126 :
127 : ```plumed
128 : plumed sum_hills --stride 300 --hills PATHTOMYHILLSFILE
129 : ```
130 :
131 : The output files are then called `fes_0.dat`, `fes_1.dat`, etc. (see below), and the last one contains all the hills.
132 :
133 : With `--nohistory` each output contains only the hills of its own block of hills, instead of all the hills up to that point:
134 :
135 : ```plumed
136 : plumed sum_hills --stride 300 --hills PATHTOMYHILLSFILE --nohistory
137 : ```
138 :
139 : In this way you can monitor the bias deposited in blocks of time as a metadynamics simulation goes on:
140 : when the bias deposited in a block is constant over the whole domain, the simulation can be considered to be at convergence.
141 :
142 : To check the convergence you can also follow the free-energy difference between basins as the
143 : simulation goes on. With `--basins` you can define two or more basins as boxes, each given by the minimum and
144 : the maximum of each variable of the output, with the basins separated by colons. Along a periodic variable a
145 : basin can cross the boundary if its minimum is larger than its maximum. For example
146 :
147 : ```plumed
148 : plumed sum_hills --hills PATHTOMYHILLSFILE --stride 500 --kt 2.5 --basins -pi,0,-pi,pi:0,pi,-pi,pi
149 : ```
150 :
151 : writes, at each output, the time of the last hill read, the number of hills, and the free energy of each basin
152 : relative to the first one, computed as $-k_BT \log \sum_{s\in\mathrm{basin}} e^{-F(s)/k_BT}$ over the points of
153 : the grid, in the file `basins.dat` (the name can be set with `--outbasins`). With `--idw` the basins are defined on
154 : the projected variables. With `--histo` (see below) the same is done for the free energy from the histogram, in the file `histo_basins.dat`.
155 : A plateau of these differences is a necessary (not sufficient) indication of convergence.
156 :
157 : ## Time average of the bias
158 :
159 : With standard (not well-tempered) metadynamics the bias does not converge but fluctuates around the free energy.
160 : A better estimate is the time average of the bias after an initial transient, which you can obtain with `--average`:
161 :
162 : ```plumed
163 : plumed sum_hills --hills PATHTOMYHILLSFILE --average 1000
164 : ```
165 :
166 : computes the average of the bias between time 1000 (in the units of the time column of the hills file) and the time of
167 : the last hill. This is done by scaling the height of each hill deposited at time $t$ after 1000 by $(t_1-t)/(t_1-1000)$, where
168 : $t_1$ is the final time. You can also give the final time explicitly, e.g. `--average 1000,5000`: hills deposited after
169 : 5000 are ignored. The same average is useful for well-tempered metadynamics, where it reduces the fluctuations of the
170 : estimate of the free energy. Combined with `--stride`, each output is the average from the initial time up to the last hill read.
171 : This is the same average computed by the function `hills_time_average` of the plumed python module.
172 :
173 : ## Bias instead of free energy
174 :
175 : If you have run well-tempered metadynamics you can also ask sum_hills to output the negative bias, $-V(s)=(\gamma-1)/\gamma\, F(s)$,
176 : instead of the free energy by using `--negbias`:
177 :
178 : ```plumed
179 : plumed sum_hills --negbias --hills PATHTOMYHILLSFILE
180 : ```
181 :
182 : Here the default output file name is `negativebias.dat`. If the simulation was not well tempered the bias is minus the free energy,
183 : and the output is the same as without `--negbias`.
184 :
185 : ## Histograms
186 :
187 : From time to time you might need to use a hills or a colvar file as a set of points from which you want to build
188 : a free energy as $-k_BT \log P(s)$. You can do this with `--histo`:
189 :
190 : ```plumed
191 : plumed sum_hills --histo PATHTOMYCOLVARORHILLSFILE --sigma 0.2,0.2 --kt 0.6
192 : ```
193 :
194 : The histogram is computed with Gaussian kernels whose widths are given with `--sigma`, one for each variable,
195 : so it is a continuous histogram, and `--kt` is needed to convert it into a free energy. The default output file is `histo.dat`.
196 : With a hills file the variables are those of the hills; with a colvar file all the columns but time are variables, so you
197 : will usually select some of them with `--idw` (and give one width for each of them).
198 : Multiple input files can be given, separated by commas.
199 :
200 : All the points have the same weight: the result is the distribution sampled by the simulation, which is the
201 : free energy only if the simulation was not biased. **sum_hills does not reweight the points.** To compute a free
202 : energy from a biased simulation (e.g. from the COLVAR file of a metadynamics simulation) use [driver](driver.md) with
203 : [READ](READ.md), [REWEIGHT_BIAS](REWEIGHT_BIAS.md) (or [REWEIGHT_METAD](REWEIGHT_METAD.md)), [HISTOGRAM](HISTOGRAM.md)
204 : and [CONVERT_TO_FES](CONVERT_TO_FES.md), for example
205 :
206 : ```plumed
207 : #SETTINGS INPUTFILES=regtest/basic/rt-fametad-2-restart/COLVAR
208 : s: READ FILE=regtest/basic/rt-fametad-2-restart/COLVAR VALUES=sum_abs IGNORE_TIME IGNORE_FORCES
209 : metad: READ FILE=regtest/basic/rt-fametad-2-restart/COLVAR VALUES=MetaD.bias IGNORE_TIME IGNORE_FORCES
210 : w: REWEIGHT_BIAS TEMP=300 ARG=metad.bias
211 : hh: HISTOGRAM ARG=s GRID_MIN=0 GRID_MAX=1 GRID_BIN=100 BANDWIDTH=0.02 LOGWEIGHTS=w NORMALIZATION=true
212 : ff: CONVERT_TO_FES ARG=hh TEMP=300
213 : DUMPGRID ARG=ff FILE=fes.dat
214 : ```
215 :
216 : run with `plumed driver --plumed plumed.dat --noatoms`. With well-tempered metadynamics the bias changes with time:
217 : use the `rbias` component of [METAD](METAD.md) (with `CALC_RCT`) instead of the bias, or [REWEIGHT_METAD](REWEIGHT_METAD.md).
218 : Block averages, which give the statistical error of the free energy, are discussed in the documentation of [HISTOGRAM](HISTOGRAM.md)
219 : and in the tutorial linked there.
220 :
221 : You can compute the free energy from the hills and the histogram at the same time, e.g. from the same file:
222 :
223 : ```plumed
224 : plumed sum_hills --hills PATHTOMYHILLSFILE --histo PATHTOMYHILLSFILE --sigma 0.2,0.2 --kt 0.6
225 : ```
226 :
227 : `--nohistory` and `--stride` can also be used with `--histo`, to compute the histogram in blocks of points:
228 :
229 : ```plumed
230 : plumed sum_hills --histo PATHTOMYCOLVARORHILLSFILE --sigma 0.2,0.2 --kt 0.6 --stride 1000 --nohistory
231 : ```
232 :
233 : ## Output
234 :
235 : The output format can be controlled by using the `--fmt` option as shown below
236 :
237 : ```plumed
238 : plumed sum_hills --hills PATHTOMYHILLSFILE --fmt %8.3f
239 : ```
240 :
241 : where here we chose a float with length of 8 and 3 digits.
242 :
243 : The output file can be named in an arbitrary way by using the `--outfile` option as shown below:
244 :
245 : ```plumed
246 : plumed sum_hills --hills PATHTOMYHILLSFILE --outfile myfes.dat
247 : ```
248 :
249 : This command produces a file `myfes.dat` which contains the free energy surface.
250 :
251 : If you use `--stride`, the name specified with `--outfile` is the prefix of the files, so this command:
252 :
253 : ```plumed
254 : plumed sum_hills --hills PATHTOMYHILLSFILE --outfile myfes_ --stride 100
255 : ```
256 :
257 : produces output files called `myfes_0.dat`, `myfes_1.dat`, `myfes_2.dat` etc.
258 :
259 : The same is true for the output of the histogram, which is named with `--outhisto`. So this example
260 :
261 : ```plumed
262 : plumed sum_hills --histo HILLS --kt 2.5 --sigma 0.01 --outhisto myhisto.dat
263 : ```
264 :
265 : produces a file `myhisto.dat`, while this input
266 :
267 : ```plumed
268 : plumed sum_hills --histo HILLS --kt 2.5 --sigma 0.01 --outhisto myhisto_ --stride 100
269 : ```
270 :
271 : produces output files called `myhisto_0.dat`, `myhisto_1.dat`, `myhisto_2.dat` etc.
272 :
273 : */
274 : //+ENDPLUMEDOC
275 :
276 : class CLToolSumHills : public CLTool {
277 : public:
278 : static void registerKeywords( Keywords& keys );
279 : explicit CLToolSumHills(const CLToolOptions& co );
280 : int main(FILE* in,FILE*out,Communicator& pc) override;
281 : std::string description()const override;
282 : /// find a list of variables present, if they are periodic and which is the period
283 : /// return false if the file does not exist
284 : 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_);
285 : };
286 :
287 7321 : void CLToolSumHills::registerKeywords( Keywords& keys ) {
288 7321 : CLTool::registerKeywords( keys );
289 7321 : keys.addFlag("--help-debug",false,"print special options that can be used to create regtests");
290 7321 : keys.add("optional","--hills","the hills files to sum, separated by commas");
291 7321 : keys.add("optional","--histo","the files (hills or colvar) from which a histogram is computed, separated by commas; the points are not reweighted");
292 7321 : keys.add("optional","--stride","write the free energy (or histogram) every this number of hills (or points); by default only at the end");
293 7321 : keys.add("optional","--min","the lower bounds of the grid, one for each variable; by default they are detected from the hills");
294 7321 : keys.add("optional","--max","the upper bounds of the grid, one for each variable; by default they are detected from the hills");
295 7321 : keys.add("optional","--bin","the number of bins of the grid, one for each variable");
296 7321 : keys.add("optional","--spacing","the grid spacing, one for each variable, alternative to the number of bins");
297 7321 : keys.add("optional","--idw","the variables of the output (default all). With --hills the other variables are integrated out, with --histo the other variables are not considered");
298 7321 : keys.add("optional","--outfile","the output file for the free energy (fes.dat, or negativebias.dat with --negbias); with --stride the prefix of the files");
299 7321 : keys.add("optional","--outhisto","the output file for the histogram (histo.dat); with --stride the prefix of the files");
300 7321 : keys.add("optional","--kt","the temperature in energy units, needed with --idw, --histo and --basins");
301 7321 : keys.add("optional","--sigma","the widths of the Gaussian kernels of the histogram, one for each variable (needed with --histo)");
302 7321 : keys.addFlag("--negbias",false,"write minus the bias instead of the free energy (they differ only with well-tempered metadynamics)");
303 7321 : keys.addFlag("--nohistory",false,"with --stride, each output contains only the hills (or points) of its own block");
304 7321 : keys.addFlag("--mintozero",false,"shift the free energy (or histogram) so that its minimum is zero");
305 7321 : keys.add("optional","--fmt","the format of the numbers in the output, e.g. %14.9f");
306 7321 : keys.add("optional","--sort","time or clock: with multiple files (e.g. one per walker), merge the hills in order of time or clock instead of reading the files one after the other");
307 7321 : keys.add("optional","--average","time average of the bias from time t0 (t0) or between times t0 and t1 (t0,t1); without t1 the average ends at the last hill read");
308 7321 : keys.add("optional","--basins","basins separated by colons, each given by the minimum and maximum of each output variable (e.g. -3,-1:1,3); the free energy of each basin relative to the first one is written at each output");
309 7321 : keys.add("optional","--outbasins","the output file for the free energies of the basins (basins.dat)");
310 7321 : }
311 :
312 22 : CLToolSumHills::CLToolSumHills(const CLToolOptions& co ):
313 22 : CLTool(co) {
314 22 : inputdata=inputType::commandline;
315 22 : }
316 :
317 7 : std::string CLToolSumHills::description()const {
318 7 : return "sum the hills of metadynamics to compute the free energy";
319 : }
320 :
321 15 : int CLToolSumHills::main(FILE* in,FILE*out,Communicator& pc) {
322 :
323 : // Read the hills input file name
324 15 : std::vector<std::string> hillsFiles;
325 : bool dohills;
326 15 : dohills=parseVector("--hills",hillsFiles);
327 : // Read the histogram file
328 15 : std::vector<std::string> histoFiles;
329 : bool dohisto;
330 15 : dohisto=parseVector("--histo",histoFiles);
331 :
332 15 : plumed_massert(dohisto || dohills,"you should use --histo or/and --hills command");
333 :
334 15 : std::vector< std::vector<std::string> > vcvs;
335 15 : std::vector<std::string> vpmin;
336 15 : std::vector<std::string> vpmax;
337 : std::string lowI_, uppI_;
338 15 : if(dohills) {
339 : // parse it as it was a restart
340 : bool vmultivariate;
341 14 : if(!findCvsAndPeriodic(hillsFiles[0], vcvs, vpmin, vpmax, vmultivariate, lowI_, uppI_)) {
342 0 : plumed_merror("cannot find file "+hillsFiles[0]);
343 : }
344 14 : if(vcvs.empty()) {
345 0 : plumed_merror("no collective variable found in file "+hillsFiles[0]);
346 : }
347 : }
348 :
349 15 : std::vector< std::vector<std::string> > hcvs;
350 15 : std::vector<std::string> hpmin;
351 15 : std::vector<std::string> hpmax;
352 :
353 15 : std::vector<std::string> sigma;
354 15 : if(dohisto) {
355 : bool hmultivariate;
356 1 : if(!findCvsAndPeriodic(histoFiles[0], hcvs, hpmin, hpmax, hmultivariate, lowI_, uppI_)) {
357 0 : plumed_merror("cannot find file "+histoFiles[0]);
358 : }
359 1 : if(hcvs.empty()) {
360 0 : plumed_merror("no collective variable found in file "+histoFiles[0]);
361 : }
362 : // here need also the vector of sigmas
363 2 : parseVector("--sigma",sigma);
364 1 : if(sigma.size()==0) {
365 0 : plumed_merror("you should define --sigma vector when using histogram");
366 : }
367 : lowI_=uppI_="-1."; // Interval is not use for histograms
368 : }
369 :
370 1 : if(dohisto && dohills) {
371 0 : plumed_massert(vcvs==hcvs,"variables for histogram and bias should have the same labels");
372 0 : plumed_massert(hpmin==vpmin,"variables for histogram and bias should have the same min for periodicity");
373 0 : plumed_massert(hpmax==vpmax,"variables for histogram and bias should have the same max for periodicity");
374 : }
375 :
376 : // now put into a neutral vector
377 :
378 15 : std::vector< std::vector<std::string> > cvs;
379 15 : std::vector<std::string> pmin;
380 15 : std::vector<std::string> pmax;
381 :
382 15 : if(dohills) {
383 14 : cvs=vcvs;
384 14 : pmin=vpmin;
385 14 : pmax=vpmax;
386 : }
387 15 : if(dohisto) {
388 1 : cvs=hcvs;
389 1 : pmin=hpmin;
390 1 : pmax=hpmax;
391 : }
392 :
393 :
394 : // setup grids
395 : unsigned grid_check=0;
396 15 : std::vector<std::string> gmin(cvs.size());
397 30 : if(parseVector("--min",gmin)) {
398 7 : if(gmin.size()!=cvs.size() && gmin.size()!=0) {
399 0 : plumed_merror("not enough values for --min");
400 : }
401 : grid_check++;
402 : }
403 15 : std::vector<std::string> gmax(cvs.size() );
404 30 : if(parseVector("--max",gmax)) {
405 7 : if(gmax.size()!=cvs.size() && gmax.size()!=0) {
406 0 : plumed_merror("not enough values for --max");
407 : }
408 7 : grid_check++;
409 : }
410 15 : std::vector<std::string> gbin(cvs.size());
411 : bool grid_has_bin;
412 : grid_has_bin=false;
413 30 : if(parseVector("--bin",gbin)) {
414 9 : if(gbin.size()!=cvs.size() && gbin.size()!=0) {
415 0 : plumed_merror("not enough values for --bin");
416 : }
417 : grid_has_bin=true;
418 : }
419 15 : std::vector<std::string> gspacing(cvs.size());
420 : bool grid_has_spacing;
421 : grid_has_spacing=false;
422 30 : if(parseVector("--spacing",gspacing)) {
423 2 : if(gspacing.size()!=cvs.size() && gspacing.size()!=0) {
424 0 : plumed_merror("not enough values for --spacing");
425 : }
426 : grid_has_spacing=true;
427 : }
428 : // allowed: no grids only bin
429 : // not allowed: partial grid definition
430 15 : 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");
431 :
432 :
433 :
434 15 : PlumedMain plumed;
435 : std::string ss;
436 15 : unsigned nn=1;
437 : ss="setNatoms";
438 15 : plumed.cmd(ss,&nn);
439 15 : if(Communicator::initialized()) {
440 0 : plumed.cmd("setMPIComm",&pc.Get_comm());
441 : }
442 15 : plumed.cmd("init",&nn);
443 15 : std::vector <bool> isdone(cvs.size(),false);
444 41 : for(unsigned i=0; i<cvs.size(); i++) {
445 26 : if(!isdone[i]) {
446 : isdone[i]=true;
447 26 : std::vector<std::string> actioninput;
448 26 : std::vector <unsigned> inds;
449 26 : actioninput.push_back("FAKE");
450 52 : actioninput.push_back("ATOMS=1");
451 26 : actioninput.push_back("LABEL="+cvs[i][0]);
452 26 : std::vector<std::string> comps, periods;
453 26 : if(cvs[i].size()>1) {
454 1 : comps.push_back(cvs[i][1]);
455 1 : inds.push_back(i);
456 : }
457 26 : periods.push_back(pmin[i]);
458 26 : periods.push_back(pmax[i]);
459 37 : for(unsigned j=i+1; j<cvs.size(); j++) {
460 11 : if(cvs[i][0]==cvs[j][0] && !isdone[j]) {
461 0 : if(cvs[i].size()==1 || cvs[j].size()==1 ) {
462 0 : plumed_merror("you cannot have twice the same label and no components ");
463 : }
464 0 : if(cvs[j].size()>1) {
465 0 : comps.push_back(cvs[j][1]);
466 0 : periods.push_back(pmin[j]);
467 0 : periods.push_back(pmax[j]);
468 : isdone[j]=true;
469 0 : inds.push_back(j);
470 : }
471 : }
472 :
473 : }
474 : // drain all the components
475 : std::string addme;
476 26 : if(comps.size()>0) {
477 : addme="COMPONENTS=";
478 1 : for(unsigned ii=0; ii<comps.size()-1; ii++) {
479 0 : addme+=comps[ii]+",";
480 : }
481 : addme+=comps.back();
482 1 : actioninput.push_back(addme);
483 : }
484 : // periodicity (always explicit here)
485 : addme="PERIODIC=";
486 52 : for(unsigned j=0; j<periods.size()-1; j++) {
487 52 : addme+=periods[j]+",";
488 : }
489 : addme+=periods.back();
490 26 : actioninput.push_back(addme);
491 27 : for(unsigned j=0; j<inds.size(); j++) {
492 : unsigned jj;
493 1 : jj=inds[j];
494 1 : if(grid_check==2) {
495 : double gm;
496 : double pm;
497 1 : if(pmin[jj]!="none") {
498 1 : Tools::convert(gmin[jj],gm);
499 1 : Tools::convert(pmin[jj],pm);
500 1 : if( gm<pm ) {
501 0 : plumed_merror("Periodicity issue : GRID_MIN value ( "+gmin[jj]+" ) is less than periodicity in HILLS file in "+cvs[jj][0]+ " ( "+pmin[jj]+" ) ");
502 : }
503 : }
504 1 : if(pmax[jj]!="none") {
505 1 : Tools::convert(gmax[jj],gm);
506 1 : Tools::convert(pmax[jj],pm);
507 1 : if( gm>pm ) {
508 0 : plumed_merror("Periodicity issue : GRID_MAX value ( "+gmax[jj]+" ) is more than periodicity in HILLS file in "+cvs[jj][0]+ " ( "+pmax[jj]+" ) ");
509 : }
510 : }
511 : }
512 : }
513 :
514 26 : plumed.readInputWords(actioninput,false);
515 26 : }
516 :
517 : }
518 15 : unsigned ncv=cvs.size();
519 15 : std::vector<std::string> actioninput;
520 15 : std::vector<std::string> idw;
521 : // check if the variables to be used are correct
522 30 : if(parseVector("--idw",idw)) {
523 4 : for(unsigned i=0; i<idw.size(); i++) {
524 : bool found=false;
525 6 : for(unsigned j=0; j<cvs.size(); j++) {
526 4 : if(cvs[j].size()>1) {
527 0 : if(idw[i]==cvs[j][0]+"."+cvs[j][1]) {
528 : found=true;
529 : }
530 : } else {
531 4 : if(idw[i]==cvs[j][0]) {
532 : found=true;
533 : }
534 : }
535 : }
536 2 : if(!found) {
537 0 : plumed_merror("variable "+idw[i]+" is not found in the bunch of cvs: revise your --idw option" );
538 : }
539 : }
540 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 ");
541 : // in this case you need a beta factor!
542 : }
543 :
544 : std::string basins;
545 30 : parse("--basins",basins);
546 :
547 : std::string kt;
548 15 : kt=std::string("1.");// assign an arbitrary value just in case that idw.size()==cvs.size()
549 14 : if ( dohisto || idw.size()!=0 || basins.length()>0 ) {
550 8 : plumed_massert(parse("--kt",kt),"if you make a dimensionality reduction (--idw), a histogram (--histo) or use --basins then you need to define --kt ");
551 : }
552 :
553 : std::string addme;
554 :
555 15 : actioninput.push_back("FUNCSUMHILLS");
556 30 : actioninput.push_back("ISCLTOOL");
557 :
558 : // set names
559 : std::string outfile;
560 30 : if(parse("--outfile",outfile)) {
561 0 : actioninput.push_back("OUTHILLS="+outfile);
562 : }
563 : std::string outhisto;
564 30 : if(parse("--outhisto",outhisto)) {
565 2 : actioninput.push_back("OUTHISTO="+outhisto);
566 : }
567 :
568 :
569 : addme="ARG=";
570 26 : for(unsigned i=0; i<(ncv-1); i++) {
571 11 : if(cvs[i].size()==1) {
572 22 : addme+=std::string(cvs[i][0])+",";
573 : } else {
574 0 : addme+=std::string(cvs[i][0])+"."+std::string(cvs[i][1])+",";
575 : }
576 : }
577 15 : if(cvs[ncv-1].size()==1) {
578 28 : addme+=std::string(cvs[ncv-1][0]);
579 : } else {
580 3 : addme+=std::string(cvs[ncv-1][0])+"."+std::string(cvs[ncv-1][1]);
581 : }
582 15 : actioninput.push_back(addme);
583 15 : if(dohills) {
584 : addme="HILLSFILES=";
585 17 : for(unsigned i=0; i<hillsFiles.size()-1; i++) {
586 6 : addme+=hillsFiles[i]+",";
587 : }
588 : addme+=hillsFiles[hillsFiles.size()-1];
589 14 : actioninput.push_back(addme);
590 : // set the grid
591 : }
592 15 : if(grid_check==2) {
593 : addme="GRID_MAX=";
594 10 : for(unsigned i=0; i<(ncv-1); i++) {
595 6 : addme+=gmax[i]+",";
596 : }
597 : addme+=gmax[ncv-1];
598 7 : actioninput.push_back(addme);
599 : addme="GRID_MIN=";
600 10 : for(unsigned i=0; i<(ncv-1); i++) {
601 6 : addme+=gmin[i]+",";
602 : }
603 : addme+=gmin[ncv-1];
604 7 : actioninput.push_back(addme);
605 : }
606 15 : if(grid_has_bin) {
607 : addme="GRID_BIN=";
608 15 : for(unsigned i=0; i<(ncv-1); i++) {
609 12 : addme+=gbin[i]+",";
610 : }
611 : addme+=gbin[ncv-1];
612 9 : actioninput.push_back(addme);
613 : }
614 15 : if(grid_has_spacing) {
615 : addme="GRID_SPACING=";
616 3 : for(unsigned i=0; i<(ncv-1); i++) {
617 2 : addme+=gspacing[i]+",";
618 : }
619 : addme+=gspacing[ncv-1];
620 2 : actioninput.push_back(addme);
621 : }
622 : std::string stride;
623 : stride="";
624 30 : if(parse("--stride",stride)) {
625 4 : actioninput.push_back("INITSTRIDE="+stride);
626 : bool nohistory;
627 4 : parseFlag("--nohistory",nohistory);
628 4 : if(nohistory) {
629 0 : actioninput.push_back("NOHISTORY");
630 : }
631 : }
632 : bool mintozero;
633 15 : parseFlag("--mintozero",mintozero);
634 15 : if(mintozero) {
635 0 : actioninput.push_back("MINTOZERO");
636 : }
637 15 : if(idw.size()!=0) {
638 : addme="PROJ=";
639 2 : for(unsigned i=0; i<idw.size()-1; i++) {
640 0 : addme+=idw[i]+",";
641 : }
642 : addme+=idw.back();
643 2 : actioninput.push_back(addme);
644 : }
645 :
646 15 : if(dohisto) {
647 1 : if(idw.size()==0) {
648 1 : if(sigma.size()!=hcvs.size()) {
649 0 : plumed_merror("you should define as many --sigma vector as the number of collective variable used for the histogram ");
650 : }
651 : } else {
652 0 : if(idw.size()!=sigma.size()) {
653 0 : plumed_merror("you should define as many --sigma vector as the number of collective variable used for the histogram ");
654 : }
655 : }
656 : }
657 :
658 15 : if(idw.size()!=0 || dohisto || basins.length()>0) {
659 8 : actioninput.push_back("KT="+kt);
660 : }
661 15 : if(basins.length()>0) {
662 2 : actioninput.push_back("BASINS="+basins);
663 : std::string outbasins;
664 2 : if(parse("--outbasins",outbasins)) {
665 0 : actioninput.push_back("OUTBASINS="+outbasins);
666 : }
667 : }
668 15 : std::vector<std::string> average;
669 30 : if(parseVector("--average",average)) {
670 1 : if(average.size()!=1 && average.size()!=2) {
671 0 : plumed_merror("--average needs the initial time and optionally the final time");
672 : }
673 2 : actioninput.push_back("AVERAGE="+average[0]+(average.size()==2 ? ","+average[1] : ""));
674 : }
675 15 : if(dohisto) {
676 : addme="HISTOFILES=";
677 1 : for(unsigned i=0; i<histoFiles.size()-1; i++) {
678 0 : addme+=histoFiles[i]+",";
679 : }
680 : addme+=histoFiles[histoFiles.size()-1];
681 1 : actioninput.push_back(addme);
682 :
683 : addme="HISTOSIGMA=";
684 2 : for(unsigned i=0; i<sigma.size()-1; i++) {
685 2 : addme+=sigma[i]+",";
686 : }
687 : addme+=sigma.back();
688 1 : actioninput.push_back(addme);
689 : }
690 :
691 : bool negbias;
692 15 : parseFlag("--negbias",negbias);
693 15 : if(negbias) {
694 4 : actioninput.push_back("NEGBIAS");
695 : }
696 :
697 15 : if(lowI_!=uppI_) {
698 : addme="INTERVAL=";
699 0 : addme+=lowI_+",";
700 : addme+=uppI_;
701 0 : actioninput.push_back(addme);
702 : }
703 :
704 : std::string sortby;
705 30 : if(parse("--sort",sortby)) {
706 1 : if(sortby!="time" && sortby!="clock") {
707 0 : plumed_merror("--sort can only be time or clock");
708 : }
709 2 : actioninput.push_back("SORT="+sortby);
710 : }
711 :
712 : std::string fmt;
713 : fmt="";
714 15 : parse("--fmt",fmt);
715 15 : if(fmt!="") {
716 30 : actioninput.push_back("FMT="+fmt);
717 : }
718 :
719 :
720 15 : plumed.readInputWords(actioninput,false);
721 : // if not a grid, then set it up automatically
722 15 : return 0;
723 90 : }
724 :
725 15 : bool CLToolSumHills::findCvsAndPeriodic(const std::string & filename,
726 : std::vector< std::vector<std::string> > &cvs,
727 : std::vector<std::string> &pmin,
728 : std::vector<std::string> &pmax,
729 : bool &multivariate,
730 : std::string &lowI_,
731 : std::string &uppI_) {
732 15 : IFile ifile;
733 15 : ifile.allowIgnoredFields();
734 15 : if(!ifile.FileExist(filename)) {
735 : return false;
736 : }
737 15 : cvs.clear();
738 15 : pmin.clear();
739 15 : pmax.clear();
740 15 : ifile.open(filename);
741 15 : std::vector<std::string> fields;
742 15 : ifile.scanFieldList(fields);
743 307 : auto hasField=[&](const std::string & name) {
744 307 : return std::find(fields.begin(),fields.end(),name)!=fields.end();
745 15 : };
746 : // In a hills file a column X is a collective variable if its width is also
747 : // there (sigma_X, or sigma_X_X for multivariate kernels). Names are compared
748 : // exactly, so that variables such as sigma_1 or lifetime are accepted (see #1226).
749 : // In a file without widths (e.g. a COLVAR file used with --histo) all the
750 : // columns but time are variables.
751 : bool hasSigma=false;
752 30 : for(const auto & f : fields) {
753 69 : if(f!="time" && (hasField("sigma_"+f) || hasField("sigma_"+f+"_"+f))) {
754 : hasSigma=true;
755 : break;
756 : }
757 : }
758 180 : for(const auto & f : fields) {
759 165 : if(f=="time") {
760 15 : continue;
761 : }
762 856 : if(hasSigma && !hasField("sigma_"+f) && !hasField("sigma_"+f+"_"+f)) {
763 124 : continue;
764 : }
765 : // separate the components
766 : auto dot=f.find_first_of('.');
767 26 : if(dot!=std::string::npos) {
768 3 : cvs.emplace_back(std::vector<std::string> {f.substr(0,dot),f.substr(dot+1)});
769 : } else {
770 50 : cvs.emplace_back(std::vector<std::string> {f});
771 : }
772 : // get periodicity
773 26 : pmin.push_back("none");
774 26 : pmax.push_back("none");
775 52 : if(ifile.FieldExist("min_"+f)) {
776 32 : ifile.scanField("min_"+f,pmin.back());
777 : }
778 52 : if(ifile.FieldExist("max_"+f)) {
779 32 : ifile.scanField("max_"+f,pmax.back());
780 : }
781 : }
782 : // is multivariate ???
783 15 : multivariate=false;
784 30 : if(ifile.FieldExist("multivariate")) {
785 : std::string sss;
786 15 : ifile.scanField("multivariate",sss);
787 15 : multivariate=(sss=="true");
788 : }
789 : // do interval?
790 30 : if(ifile.FieldExist("lower_int")) {
791 0 : ifile.scanField("lower_int",lowI_);
792 0 : ifile.scanField("upper_int",uppI_);
793 : } else {
794 : lowI_="-1.";
795 : uppI_="-1.";
796 : }
797 15 : ifile.scanField();
798 : return true;
799 41 : }
800 :
801 :
802 21985 : PLUMED_REGISTER_CLTOOL(CLToolSumHills,"sum_hills")
803 :
804 :
805 :
806 : }
807 : }
|
