1 #include <iostream>
2 #include <sstream>
3 #include <fstream>
4 #include <iomanip>
5 #include <cmath>
6 #include <algorithm>
7
8 // used to set the absolute path of a replica file
9 #if defined (WIN32) && !defined(__CYGWIN__)
10 #include <direct.h>
11 #define CHDIR ::_chdir
12 #define GETCWD ::_getcwd
13 #define PATHSEP "\\"
14 #else
15 #include <unistd.h>
16 #define CHDIR ::chdir
17 #define GETCWD ::getcwd
18 #define PATHSEP "/"
19 #endif
20
21
22 #include "colvar.h"
23 #include "colvarbias_meta.h"
24
25
colvarbias_meta()26 colvarbias_meta::colvarbias_meta()
27 : colvarbias(),
28 state_file_step (0),
29 new_hills_begin (hills.end())
30 {
31 }
32
33
colvarbias_meta(std::string const & conf,char const * key)34 colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
35 : colvarbias (conf, key),
36 state_file_step (0),
37 new_hills_begin (hills.end())
38 {
39 if (cvm::n_abf_biases > 0)
40 cvm::log ("Warning: ABF and metadynamics running at the "
41 "same time can lead to inconsistent results.\n");
42
43 get_keyval (conf, "hillWeight", hill_weight, 0.01);
44 if (hill_weight == 0.0)
45 cvm::log ("Warning: hillWeight has been set to zero, "
46 "this bias will have no effect.\n");
47
48 get_keyval (conf, "newHillFrequency", new_hill_freq, 1000);
49
50 get_keyval (conf, "hillWidth", hill_width, std::sqrt (2.0 * PI) / 2.0);
51
52 {
53 bool b_replicas = false;
54 get_keyval (conf, "multipleReplicas", b_replicas, false);
55 if (b_replicas)
56 comm = multiple_replicas;
57 else
58 comm = single_replica;
59 }
60
61 get_keyval (conf, "useGrids", use_grids, true);
62
63 if (use_grids) {
64 get_keyval (conf, "gridsUpdateFrequency", grids_freq, new_hill_freq);
65 get_keyval (conf, "rebinGrids", rebin_grids, false);
66
67 expand_grids = false;
68 for (size_t i = 0; i < colvars.size(); i++) {
69 if (colvars[i]->expand_boundaries) {
70 expand_grids = true;
71 cvm::log ("Metadynamics bias \""+this->name+"\""+
72 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
73 ": Will expand grids when the colvar \""+
74 colvars[i]->name+"\" approaches its boundaries.\n");
75 }
76 }
77
78 get_keyval (conf, "keepHills", keep_hills, false);
79 if (! get_keyval (conf, "writeFreeEnergyFile", dump_fes, true))
80 get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
81 get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false);
82
83 for (size_t i = 0; i < colvars.size(); i++) {
84 colvars[i]->enable (colvar::task_grid);
85 }
86
87 hills_energy = new colvar_grid_scalar (colvars);
88 hills_energy_gradients = new colvar_grid_gradient (colvars);
89 } else {
90 rebin_grids = false;
91 keep_hills = false;
92 dump_fes = false;
93 dump_fes_save = false;
94 dump_replica_fes = false;
95 }
96
97 if (comm != single_replica) {
98
99 if (expand_grids)
100 cvm::fatal_error ("Error: expandBoundaries is not supported when "
101 "using more than one replicas; please allocate "
102 "wide enough boundaries for each colvar"
103 "ahead of time.\n");
104
105 if (get_keyval (conf, "dumpPartialFreeEnergyFile", dump_replica_fes, false)) {
106 if (dump_replica_fes && (! dump_fes)) {
107 cvm::log ("Enabling \"dumpFreeEnergyFile\".\n");
108 }
109 }
110
111 get_keyval (conf, "replicaID", replica_id, std::string (""));
112 if (!replica_id.size())
113 cvm::fatal_error ("Error: replicaID must be defined "
114 "when using more than one replica.\n");
115
116 get_keyval (conf, "replicasRegistry",
117 replicas_registry_file,
118 (this->name+".replicas.registry.txt"));
119
120 get_keyval (conf, "replicaUpdateFrequency",
121 replica_update_freq, new_hill_freq);
122
123 if (keep_hills)
124 cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+
125 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
126 ": keepHills with more than one replica can lead to a very "
127 "large amount input/output and slow down your calculations. "
128 "Please consider disabling it.\n");
129
130
131 {
132 // TODO: one may want to specify the path manually for intricated filesystems?
133 char *pwd = new char[3001];
134 if (GETCWD (pwd, 3000) == NULL)
135 cvm::fatal_error ("Error: cannot get the path of the current working directory.\n");
136 replica_list_file =
137 (std::string (pwd)+std::string (PATHSEP)+
138 this->name+"."+replica_id+".files.txt");
139 // replica_hills_file and replica_state_file are those written
140 // by the current replica; within the mirror biases, they are
141 // those by another replica
142 replica_hills_file =
143 (std::string (pwd)+std::string (PATHSEP)+
144 cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills");
145 replica_state_file =
146 (std::string (pwd)+std::string (PATHSEP)+
147 cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state");
148 delete pwd;
149 }
150
151 // now register this replica
152
153 // first check that it isn't already there
154 bool registered_replica = false;
155 std::ifstream reg_is (replicas_registry_file.c_str());
156 if (reg_is.good()) { // the file may not be there yet
157 std::string existing_replica ("");
158 std::string existing_replica_file ("");
159 while ((reg_is >> existing_replica) && existing_replica.size() &&
160 (reg_is >> existing_replica_file) && existing_replica_file.size()) {
161 if (existing_replica == replica_id) {
162 // this replica was already registered
163 replica_list_file = existing_replica_file;
164 reg_is.close();
165 registered_replica = true;
166 break;
167 }
168 }
169 reg_is.close();
170 }
171
172 // if this replica was not included yet, we should generate a
173 // new record for it: but first, we write this replica's files,
174 // for the others to read
175
176 // open the "hills" buffer file
177 replica_hills_os.open (replica_hills_file.c_str());
178 if (!replica_hills_os.good())
179 cvm::fatal_error ("Error: in opening file \""+
180 replica_hills_file+"\" for writing.\n");
181 replica_hills_os.setf (std::ios::scientific, std::ios::floatfield);
182
183 // write the state file (so that there is always one available)
184 write_replica_state_file();
185 // schedule to read the state files of the other replicas
186 for (size_t ir = 0; ir < replicas.size(); ir++) {
187 (replicas[ir])->replica_state_file_in_sync = false;
188 }
189
190 // if we're running without grids, use a growing list of "hills" files
191 // otherwise, just one state file and one "hills" file as buffer
192 std::ofstream list_os (replica_list_file.c_str(),
193 (use_grids ? std::ios::trunc : std::ios::app));
194 if (! list_os.good())
195 cvm::fatal_error ("Error: in opening file \""+
196 replica_list_file+"\" for writing.\n");
197 list_os << "stateFile " << replica_state_file << "\n";
198 list_os << "hillsFile " << replica_hills_file << "\n";
199 list_os.close();
200
201 // finally, if add a new record for this replica to the registry
202 if (! registered_replica) {
203 std::ofstream reg_os (replicas_registry_file.c_str(), std::ios::app);
204 if (! reg_os.good())
205 cvm::fatal_error ("Error: in opening file \""+
206 replicas_registry_file+"\" for writing.\n");
207 reg_os << replica_id << " " << replica_list_file << "\n";
208 reg_os.close();
209 }
210 }
211
212 get_keyval (conf, "writeHillsTrajectory", b_hills_traj, false);
213 if (b_hills_traj) {
214 std::string const traj_file_name (cvm::output_prefix+
215 ".colvars."+this->name+
216 ( (comm != single_replica) ?
217 ("."+replica_id) :
218 ("") )+
219 ".hills.traj");
220 hills_traj_os.open (traj_file_name.c_str());
221 if (!hills_traj_os.good())
222 cvm::fatal_error ("Error: in opening hills output file \"" +
223 traj_file_name + "\".\n");
224 }
225
226 // for well-tempered metadynamics
227 get_keyval (conf, "wellTempered", well_tempered, false);
228 get_keyval (conf, "biasTemperature", bias_temperature, -1.0);
229 if ((bias_temperature == -1.0) && well_tempered) {
230 cvm::fatal_error ("Error: biasTemperature is not set.\n");
231 }
232 if (well_tempered) {
233 cvm::log("Well-tempered metadynamics is used.\n");
234 cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
235 }
236
237 if (cvm::debug())
238 cvm::log ("Done initializing the metadynamics bias \""+this->name+"\""+
239 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
240
241 save_delimiters = false;
242 }
243
244
~colvarbias_meta()245 colvarbias_meta::~colvarbias_meta()
246 {
247 if (hills_energy) {
248 delete hills_energy;
249 hills_energy = NULL;
250 }
251
252 if (hills_energy_gradients) {
253 delete hills_energy_gradients;
254 hills_energy_gradients = NULL;
255 }
256
257 if (replica_hills_os.good())
258 replica_hills_os.close();
259
260 if (hills_traj_os.good())
261 hills_traj_os.close();
262
263 if (cvm::n_meta_biases > 0)
264 cvm::n_meta_biases -= 1;
265 }
266
267
268
269 // **********************************************************************
270 // Hill management member functions
271 // **********************************************************************
272
273 std::list<colvarbias_meta::hill>::const_iterator
create_hill(colvarbias_meta::hill const & h)274 colvarbias_meta::create_hill (colvarbias_meta::hill const &h)
275 {
276 hill_iter const hills_end = hills.end();
277 hills.push_back (h);
278 if (new_hills_begin == hills_end) {
279 // if new_hills_begin is unset, set it for the first time
280 new_hills_begin = hills.end();
281 new_hills_begin--;
282 }
283
284 if (use_grids) {
285
286 // also add it to the list of hills that are off-grid, which may
287 // need to be computed analytically when the colvar returns
288 // off-grid
289 cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers, true);
290 if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
291 hills_off_grid.push_back (h);
292 }
293 }
294
295 // output to trajectory (if specified)
296 if (hills_traj_os.good()) {
297 hills_traj_os << (hills.back()).output_traj();
298 if (cvm::debug()) {
299 hills_traj_os.flush();
300 }
301 }
302
303 has_data = true;
304 return hills.end();
305 }
306
307
308 std::list<colvarbias_meta::hill>::const_iterator
delete_hill(hill_iter & h)309 colvarbias_meta::delete_hill (hill_iter &h)
310 {
311 if (cvm::debug()) {
312 cvm::log ("Deleting hill from the metadynamics bias \""+this->name+"\""+
313 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
314 ", with step number "+
315 cvm::to_str (h->it)+(h->replica.size() ?
316 ", replica id \""+h->replica :
317 "")+".\n");
318 }
319
320 if (use_grids && hills_off_grid.size()) {
321 for (hill_iter hoff = hills_off_grid.begin();
322 hoff != hills_off_grid.end(); hoff++) {
323 if (*h == *hoff) {
324 hills_off_grid.erase (hoff);
325 break;
326 }
327 }
328 }
329
330 if (hills_traj_os.good()) {
331 // output to the trajectory
332 hills_traj_os << "# DELETED this hill: "
333 << (hills.back()).output_traj()
334 << "\n";
335 if (cvm::debug())
336 hills_traj_os.flush();
337 }
338
339 return hills.erase (h);
340 }
341
342
update()343 cvm::real colvarbias_meta::update()
344 {
345 if (cvm::debug())
346 cvm::log ("Updating the metadynamics bias \""+this->name+"\""+
347 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
348
349 if (use_grids) {
350
351 std::vector<int> curr_bin = hills_energy->get_colvars_index();
352
353 if (expand_grids) {
354
355 // first of all, expand the grids, if specified
356 if (cvm::debug())
357 cvm::log ("Metadynamics bias \""+this->name+"\""+
358 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
359 ": current coordinates on the grid: "+
360 cvm::to_str (curr_bin)+".\n");
361
362 bool changed_grids = false;
363 int const min_buffer =
364 (3 * (size_t) std::floor (hill_width)) + 1;
365
366 std::vector<int> new_sizes (hills_energy->sizes());
367 std::vector<colvarvalue> new_lower_boundaries (hills_energy->lower_boundaries);
368 std::vector<colvarvalue> new_upper_boundaries (hills_energy->upper_boundaries);
369
370 for (size_t i = 0; i < colvars.size(); i++) {
371
372 if (! colvars[i]->expand_boundaries)
373 continue;
374
375 cvm::real &new_lb = new_lower_boundaries[i].real_value;
376 cvm::real &new_ub = new_upper_boundaries[i].real_value;
377 int &new_size = new_sizes[i];
378 bool changed_lb = false, changed_ub = false;
379
380 if (!colvars[i]->hard_lower_boundary)
381 if (curr_bin[i] < min_buffer) {
382 int const extra_points = (min_buffer - curr_bin[i]);
383 new_lb -= extra_points * colvars[i]->width;
384 new_size += extra_points;
385 // changed offset in this direction => the pointer needs to
386 // be changed, too
387 curr_bin[i] += extra_points;
388
389 changed_lb = true;
390 cvm::log ("Metadynamics bias \""+this->name+"\""+
391 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
392 ": new lower boundary for colvar \""+
393 colvars[i]->name+"\", at "+
394 cvm::to_str (new_lower_boundaries[i])+".\n");
395 }
396
397 if (!colvars[i]->hard_upper_boundary)
398 if (curr_bin[i] > new_size - min_buffer - 1) {
399 int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
400 new_ub += extra_points * colvars[i]->width;
401 new_size += extra_points;
402
403 changed_ub = true;
404 cvm::log ("Metadynamics bias \""+this->name+"\""+
405 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
406 ": new upper boundary for colvar \""+
407 colvars[i]->name+"\", at "+
408 cvm::to_str (new_upper_boundaries[i])+".\n");
409 }
410
411 if (changed_lb || changed_ub)
412 changed_grids = true;
413 }
414
415 if (changed_grids) {
416
417 // map everything into new grids
418
419 colvar_grid_scalar *new_hills_energy =
420 new colvar_grid_scalar (*hills_energy);
421 colvar_grid_gradient *new_hills_energy_gradients =
422 new colvar_grid_gradient (*hills_energy_gradients);
423
424 // supply new boundaries to the new grids
425
426 new_hills_energy->lower_boundaries = new_lower_boundaries;
427 new_hills_energy->upper_boundaries = new_upper_boundaries;
428 new_hills_energy->create (new_sizes, 0.0, 1);
429
430 new_hills_energy_gradients->lower_boundaries = new_lower_boundaries;
431 new_hills_energy_gradients->upper_boundaries = new_upper_boundaries;
432 new_hills_energy_gradients->create (new_sizes, 0.0, colvars.size());
433
434 new_hills_energy->map_grid (*hills_energy);
435 new_hills_energy_gradients->map_grid (*hills_energy_gradients);
436
437 delete hills_energy;
438 delete hills_energy_gradients;
439 hills_energy = new_hills_energy;
440 hills_energy_gradients = new_hills_energy_gradients;
441
442 curr_bin = hills_energy->get_colvars_index();
443 if (cvm::debug())
444 cvm::log ("Coordinates on the new grid: "+
445 cvm::to_str (curr_bin)+".\n");
446 }
447 }
448 }
449
450 // add a new hill if the required time interval has passed
451 if ((cvm::step_absolute() % new_hill_freq) == 0) {
452
453 if (cvm::debug())
454 cvm::log ("Metadynamics bias \""+this->name+"\""+
455 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
456 ": adding a new hill at step "+cvm::to_str (cvm::step_absolute())+".\n");
457
458 switch (comm) {
459
460 case single_replica:
461 if (well_tempered) {
462 std::vector<int> curr_bin = hills_energy->get_colvars_index();
463 cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin);
464 cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
465 create_hill (hill ((hill_weight*exp_weight), colvars, hill_width));
466 } else {
467 create_hill (hill (hill_weight, colvars, hill_width));
468 }
469 break;
470
471 case multiple_replicas:
472 if (well_tempered) {
473 std::vector<int> curr_bin = hills_energy->get_colvars_index();
474 cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin);
475 cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
476 create_hill (hill ((hill_weight*exp_weight), colvars, hill_width, replica_id));
477 } else {
478 create_hill (hill (hill_weight, colvars, hill_width, replica_id));
479 }
480 if (replica_hills_os.good()) {
481 replica_hills_os << hills.back();
482 } else {
483 cvm::fatal_error ("Error: in metadynamics bias \""+this->name+"\""+
484 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
485 " while writing hills for the other replicas.\n");
486 }
487 break;
488 }
489 }
490
491 // sync with the other replicas (if needed)
492 if (comm != single_replica) {
493
494 // reread the replicas registry
495 if ((cvm::step_absolute() % replica_update_freq) == 0) {
496 update_replicas_registry();
497 // empty the output buffer
498 replica_hills_os.flush();
499
500 read_replica_files();
501 }
502 }
503
504 // calculate the biasing energy and forces
505 bias_energy = 0.0;
506 for (size_t ir = 0; ir < colvars.size(); ir++) {
507 colvar_forces[ir].reset();
508 }
509 if (comm == multiple_replicas)
510 for (size_t ir = 0; ir < replicas.size(); ir++) {
511 replicas[ir]->bias_energy = 0.0;
512 for (size_t ic = 0; ic < colvars.size(); ic++) {
513 replicas[ir]->colvar_forces[ic].reset();
514 }
515 }
516
517 if (use_grids) {
518
519 // get the forces from the grid
520
521 if ((cvm::step_absolute() % grids_freq) == 0) {
522 // map the most recent gaussians to the grids
523 project_hills (new_hills_begin, hills.end(),
524 hills_energy, hills_energy_gradients);
525 new_hills_begin = hills.end();
526
527 // TODO: we may want to condense all into one replicas array,
528 // including "this" as the first element
529 if (comm == multiple_replicas) {
530 for (size_t ir = 0; ir < replicas.size(); ir++) {
531 replicas[ir]->project_hills (replicas[ir]->new_hills_begin,
532 replicas[ir]->hills.end(),
533 replicas[ir]->hills_energy,
534 replicas[ir]->hills_energy_gradients);
535 replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
536 }
537 }
538 }
539
540 std::vector<int> curr_bin = hills_energy->get_colvars_index();
541 if (cvm::debug())
542 cvm::log ("Metadynamics bias \""+this->name+"\""+
543 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
544 ": current coordinates on the grid: "+
545 cvm::to_str (curr_bin)+".\n");
546
547 if (hills_energy->index_ok (curr_bin)) {
548
549 // within the grid: add the energy and the forces from there
550
551 bias_energy += hills_energy->value (curr_bin);
552 for (size_t ic = 0; ic < colvars.size(); ic++) {
553 cvm::real const *f = &(hills_energy_gradients->value (curr_bin));
554 colvar_forces[ic].real_value += -1.0 * f[ic];
555 // the gradients are stored, not the forces
556 }
557 if (comm == multiple_replicas)
558 for (size_t ir = 0; ir < replicas.size(); ir++) {
559 bias_energy += replicas[ir]->hills_energy->value (curr_bin);
560 cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value (curr_bin));
561 for (size_t ic = 0; ic < colvars.size(); ic++) {
562 colvar_forces[ic].real_value += -1.0 * f[ic];
563 }
564 }
565
566 } else {
567
568 // off the grid: compute analytically only the hills at the grid's edges
569
570 calc_hills (hills_off_grid.begin(), hills_off_grid.end(), bias_energy);
571 for (size_t ic = 0; ic < colvars.size(); ic++) {
572 calc_hills_force (ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces);
573 }
574
575 if (comm == multiple_replicas)
576 for (size_t ir = 0; ir < replicas.size(); ir++) {
577 calc_hills (replicas[ir]->hills_off_grid.begin(),
578 replicas[ir]->hills_off_grid.end(),
579 bias_energy);
580 for (size_t ic = 0; ic < colvars.size(); ic++) {
581 calc_hills_force (ic,
582 replicas[ir]->hills_off_grid.begin(),
583 replicas[ir]->hills_off_grid.end(),
584 colvar_forces);
585 }
586 }
587 }
588 }
589
590 // now include the hills that have not been binned yet (starting
591 // from new_hills_begin)
592
593 calc_hills (new_hills_begin, hills.end(), bias_energy);
594 for (size_t ic = 0; ic < colvars.size(); ic++) {
595 calc_hills_force (ic, new_hills_begin, hills.end(), colvar_forces);
596 }
597
598 if (cvm::debug())
599 cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
600 ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
601
602 if (cvm::debug())
603 cvm::log ("Metadynamics bias \""+this->name+"\""+
604 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
605 ": adding the forces from the other replicas.\n");
606
607 if (comm == multiple_replicas)
608 for (size_t ir = 0; ir < replicas.size(); ir++) {
609 calc_hills (replicas[ir]->new_hills_begin,
610 replicas[ir]->hills.end(),
611 bias_energy);
612 for (size_t ic = 0; ic < colvars.size(); ic++) {
613 calc_hills_force (ic,
614 replicas[ir]->new_hills_begin,
615 replicas[ir]->hills.end(),
616 colvar_forces);
617 }
618 if (cvm::debug())
619 cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
620 ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
621 }
622
623 return bias_energy;
624 }
625
626
calc_hills(colvarbias_meta::hill_iter h_first,colvarbias_meta::hill_iter h_last,cvm::real & energy,std::vector<colvarvalue> const & colvar_values)627 void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter h_first,
628 colvarbias_meta::hill_iter h_last,
629 cvm::real &energy,
630 std::vector<colvarvalue> const &colvar_values)
631 {
632 std::vector<colvarvalue> curr_values (colvars.size());
633 for (size_t i = 0; i < colvars.size(); i++) {
634 curr_values[i].type (colvars[i]->type());
635 }
636
637 if (colvar_values.size()) {
638 for (size_t i = 0; i < colvars.size(); i++) {
639 curr_values[i] = colvar_values[i];
640 }
641 } else {
642 for (size_t i = 0; i < colvars.size(); i++) {
643 curr_values[i] = colvars[i]->value();
644 }
645 }
646
647 for (hill_iter h = h_first; h != h_last; h++) {
648
649 // compute the gaussian exponent
650 cvm::real cv_sqdev = 0.0;
651 for (size_t i = 0; i < colvars.size(); i++) {
652 colvarvalue const &x = curr_values[i];
653 colvarvalue const ¢er = h->centers[i];
654 cvm::real const half_width = 0.5 * h->widths[i];
655 cv_sqdev += (colvars[i]->dist2 (x, center)) / (half_width*half_width);
656 }
657
658 // compute the gaussian
659 if (cv_sqdev > 23.0) {
660 // set it to zero if the exponent is more negative than log(1.0E-05)
661 h->value (0.0);
662 } else {
663 h->value (std::exp (-0.5*cv_sqdev));
664 }
665 energy += h->energy();
666 }
667 }
668
669
calc_hills_force(size_t const & i,colvarbias_meta::hill_iter h_first,colvarbias_meta::hill_iter h_last,std::vector<colvarvalue> & forces,std::vector<colvarvalue> const & values)670 void colvarbias_meta::calc_hills_force (size_t const &i,
671 colvarbias_meta::hill_iter h_first,
672 colvarbias_meta::hill_iter h_last,
673 std::vector<colvarvalue> &forces,
674 std::vector<colvarvalue> const &values)
675 {
676 // Retrieve the value of the colvar
677 colvarvalue x (values.size() ? values[i].type() : colvars[i]->type());
678 x = (values.size() ? values[i] : colvars[i]->value());
679
680 // do the type check only once (all colvarvalues in the hills series
681 // were already saved with their types matching those in the
682 // colvars)
683
684 switch (colvars[i]->type()) {
685
686 case colvarvalue::type_scalar:
687 for (hill_iter h = h_first; h != h_last; h++) {
688 if (h->value() == 0.0) continue;
689 colvarvalue const ¢er = h->centers[i];
690 cvm::real const half_width = 0.5 * h->widths[i];
691 forces[i].real_value +=
692 ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
693 (colvars[i]->dist2_lgrad (x, center)).real_value );
694 }
695 break;
696
697 case colvarvalue::type_vector:
698 case colvarvalue::type_unitvector:
699 for (hill_iter h = h_first; h != h_last; h++) {
700 if (h->value() == 0.0) continue;
701 colvarvalue const ¢er = h->centers[i];
702 cvm::real const half_width = 0.5 * h->widths[i];
703 forces[i].rvector_value +=
704 ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
705 (colvars[i]->dist2_lgrad (x, center)).rvector_value );
706 }
707 break;
708
709 case colvarvalue::type_quaternion:
710 for (hill_iter h = h_first; h != h_last; h++) {
711 if (h->value() == 0.0) continue;
712 colvarvalue const ¢er = h->centers[i];
713 cvm::real const half_width = 0.5 * h->widths[i];
714 forces[i].quaternion_value +=
715 ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
716 (colvars[i]->dist2_lgrad (x, center)).quaternion_value );
717 }
718 break;
719
720 case colvarvalue::type_notset:
721 break;
722 }
723 }
724
725
726 // **********************************************************************
727 // grid management functions
728 // **********************************************************************
729
project_hills(colvarbias_meta::hill_iter h_first,colvarbias_meta::hill_iter h_last,colvar_grid_scalar * he,colvar_grid_gradient * hg,bool print_progress)730 void colvarbias_meta::project_hills (colvarbias_meta::hill_iter h_first,
731 colvarbias_meta::hill_iter h_last,
732 colvar_grid_scalar *he,
733 colvar_grid_gradient *hg,
734 bool print_progress)
735 {
736 if (cvm::debug())
737 cvm::log ("Metadynamics bias \""+this->name+"\""+
738 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
739 ": projecting hills.\n");
740
741 // TODO: improve it by looping over a small subgrid instead of the whole grid
742
743 std::vector<colvarvalue> colvar_values (colvars.size());
744 std::vector<cvm::real> colvar_forces_scalar (colvars.size());
745
746 std::vector<int> he_ix = he->new_index();
747 std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
748 cvm::real hills_energy_here = 0.0;
749 std::vector<colvarvalue> hills_forces_here (colvars.size(), 0.0);
750
751 size_t count = 0;
752 size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1)));
753
754 if (hg != NULL) {
755
756 // loop over the points of the grid
757 for ( ;
758 (he->index_ok (he_ix)) && (hg->index_ok (hg_ix));
759 count++) {
760
761 for (size_t i = 0; i < colvars.size(); i++) {
762 colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
763 }
764
765 // loop over the hills and increment the energy grid locally
766 hills_energy_here = 0.0;
767 calc_hills (h_first, h_last, hills_energy_here, colvar_values);
768 he->acc_value (he_ix, hills_energy_here);
769
770 for (size_t i = 0; i < colvars.size(); i++) {
771 hills_forces_here[i].reset();
772 calc_hills_force (i, h_first, h_last, hills_forces_here, colvar_values);
773 colvar_forces_scalar[i] = hills_forces_here[i].real_value;
774 }
775 hg->acc_force (hg_ix, &(colvar_forces_scalar.front()));
776
777 he->incr (he_ix);
778 hg->incr (hg_ix);
779
780 if ((count % print_frequency) == 0) {
781 if (print_progress) {
782 cvm::real const progress = cvm::real (count) / cvm::real (hg->number_of_points());
783 std::ostringstream os;
784 os.setf (std::ios::fixed, std::ios::floatfield);
785 os << std::setw (6) << std::setprecision (2)
786 << 100.0 * progress
787 << "% done.";
788 cvm::log (os.str());
789 }
790 }
791 }
792
793 } else {
794
795 // simpler version, with just the energy
796
797 for ( ; (he->index_ok (he_ix)); ) {
798
799 for (size_t i = 0; i < colvars.size(); i++) {
800 colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
801 }
802
803 hills_energy_here = 0.0;
804 calc_hills (h_first, h_last, hills_energy_here, colvar_values);
805 he->acc_value (he_ix, hills_energy_here);
806
807 he->incr (he_ix);
808
809 count++;
810 if ((count % print_frequency) == 0) {
811 if (print_progress) {
812 cvm::real const progress = cvm::real (count) / cvm::real (he->number_of_points());
813 std::ostringstream os;
814 os.setf (std::ios::fixed, std::ios::floatfield);
815 os << std::setw (6) << std::setprecision (2)
816 << 100.0 * progress
817 << "% done.";
818 cvm::log (os.str());
819 }
820 }
821 }
822 }
823
824 if (print_progress) {
825 cvm::log ("100.00% done.");
826 }
827
828 if (! keep_hills) {
829 hills.erase (hills.begin(), hills.end());
830 }
831 }
832
833
recount_hills_off_grid(colvarbias_meta::hill_iter h_first,colvarbias_meta::hill_iter h_last,colvar_grid_scalar * he)834 void colvarbias_meta::recount_hills_off_grid (colvarbias_meta::hill_iter h_first,
835 colvarbias_meta::hill_iter h_last,
836 colvar_grid_scalar *he)
837 {
838 hills_off_grid.clear();
839
840 for (hill_iter h = h_first; h != h_last; h++) {
841 cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers, true);
842 if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
843 hills_off_grid.push_back (*h);
844 }
845 }
846 }
847
848
849
850 // **********************************************************************
851 // multiple replicas functions
852 // **********************************************************************
853
854
update_replicas_registry()855 void colvarbias_meta::update_replicas_registry()
856 {
857 if (cvm::debug())
858 cvm::log ("Metadynamics bias \""+this->name+"\""+
859 ": updating the list of replicas, currently containing "+
860 cvm::to_str (replicas.size())+" elements.\n");
861
862 {
863 // copy the whole file into a string for convenience
864 std::string line ("");
865 std::ifstream reg_file (replicas_registry_file.c_str());
866 if (reg_file.good()) {
867 replicas_registry.clear();
868 while (colvarparse::getline_nocomments (reg_file, line))
869 replicas_registry.append (line+"\n");
870 } else {
871 cvm::fatal_error ("Error: failed to open file \""+replicas_registry_file+
872 "\" for reading.\n");
873 }
874 }
875
876 // now parse it
877 std::istringstream reg_is (replicas_registry);
878 if (reg_is.good()) {
879
880 std::string new_replica ("");
881 std::string new_replica_file ("");
882 while ((reg_is >> new_replica) && new_replica.size() &&
883 (reg_is >> new_replica_file) && new_replica_file.size()) {
884
885 if (new_replica == this->replica_id) {
886 // this is the record for this same replica, skip it
887 if (cvm::debug())
888 cvm::log ("Metadynamics bias \""+this->name+"\""+
889 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
890 ": skipping this replica's own record: \""+
891 new_replica+"\", \""+new_replica_file+"\"\n");
892 new_replica_file.clear();
893 new_replica.clear();
894 continue;
895 }
896
897 bool already_loaded = false;
898 for (size_t ir = 0; ir < replicas.size(); ir++) {
899 if (new_replica == (replicas[ir])->replica_id) {
900 // this replica was already added
901 if (cvm::debug())
902 cvm::log ("Metadynamics bias \""+this->name+"\""+
903 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
904 ": skipping a replica already loaded, \""+
905 (replicas[ir])->replica_id+"\".\n");
906 already_loaded = true;
907 break;
908 }
909 }
910
911 if (!already_loaded) {
912 // add this replica to the registry
913 cvm::log ("Metadynamics bias \""+this->name+"\""+
914 ": accessing replica \""+new_replica+"\".\n");
915 replicas.push_back (new colvarbias_meta());
916 (replicas.back())->replica_id = new_replica;
917 (replicas.back())->replica_list_file = new_replica_file;
918 (replicas.back())->replica_state_file = "";
919 (replicas.back())->replica_state_file_in_sync = false;
920
921 // Note: the following could become a copy constructor?
922 (replicas.back())->name = this->name;
923 (replicas.back())->colvars = colvars;
924 (replicas.back())->use_grids = use_grids;
925 (replicas.back())->dump_fes = false;
926 (replicas.back())->expand_grids = false;
927 (replicas.back())->rebin_grids = false;
928 (replicas.back())->keep_hills = false;
929 (replicas.back())->colvar_forces = colvar_forces;
930
931 (replicas.back())->comm = multiple_replicas;
932
933 if (use_grids) {
934 (replicas.back())->hills_energy = new colvar_grid_scalar (colvars);
935 (replicas.back())->hills_energy_gradients = new colvar_grid_gradient (colvars);
936 }
937 }
938 }
939
940 // continue the cycle
941 new_replica_file = "";
942 new_replica = "";
943 } else {
944 cvm::fatal_error ("Error: cannot read the replicas registry file \""+
945 replicas_registry+"\".\n");
946 }
947
948 // now (re)read the list file of each replica
949 for (size_t ir = 0; ir < replicas.size(); ir++) {
950 if (cvm::debug())
951 cvm::log ("Metadynamics bias \""+this->name+"\""+
952 ": reading the list file for replica \""+
953 (replicas[ir])->replica_id+"\".\n");
954
955 std::ifstream list_is ((replicas[ir])->replica_list_file.c_str());
956 std::string key;
957 std::string new_state_file, new_hills_file;
958 if (!(list_is >> key) ||
959 !(key == std::string ("stateFile")) ||
960 !(list_is >> new_state_file) ||
961 !(list_is >> key) ||
962 !(key == std::string ("hillsFile")) ||
963 !(list_is >> new_hills_file)) {
964 cvm::log ("Metadynamics bias \""+this->name+"\""+
965 ": failed to read the file \""+
966 (replicas[ir])->replica_list_file+"\": will try again after "+
967 cvm::to_str (replica_update_freq)+" steps.\n");
968 (replicas[ir])->update_status++;
969 } else {
970 (replicas[ir])->update_status = 0;
971 if (new_state_file != (replicas[ir])->replica_state_file) {
972 cvm::log ("Metadynamics bias \""+this->name+"\""+
973 ": replica \""+(replicas[ir])->replica_id+
974 "\" has supplied a new state file, \""+new_state_file+
975 "\".\n");
976 (replicas[ir])->replica_state_file_in_sync = false;
977 (replicas[ir])->replica_state_file = new_state_file;
978 (replicas[ir])->replica_hills_file = new_hills_file;
979 }
980 }
981 }
982
983
984 if (cvm::debug())
985 cvm::log ("Metadynamics bias \""+this->name+"\": the list of replicas contains "+
986 cvm::to_str (replicas.size())+" elements.\n");
987 }
988
989
read_replica_files()990 void colvarbias_meta::read_replica_files()
991 {
992 for (size_t ir = 0; ir < replicas.size(); ir++) {
993
994 if (! (replicas[ir])->replica_state_file_in_sync) {
995 // if a new state file is being read, the hills file is also new
996 (replicas[ir])->replica_hills_file_pos = 0;
997 }
998
999 // (re)read the state file if necessary
1000 if ( (! (replicas[ir])->has_data) ||
1001 (! (replicas[ir])->replica_state_file_in_sync) ) {
1002
1003 cvm::log ("Metadynamics bias \""+this->name+"\""+
1004 ": reading the state of replica \""+
1005 (replicas[ir])->replica_id+"\" from file \""+
1006 (replicas[ir])->replica_state_file+"\".\n");
1007
1008 std::ifstream is ((replicas[ir])->replica_state_file.c_str());
1009 if (! (replicas[ir])->read_restart (is)) {
1010 cvm::log ("Reading from file \""+(replicas[ir])->replica_state_file+
1011 "\" failed or incomplete: will try again in "+
1012 cvm::to_str (replica_update_freq)+" steps.\n");
1013 } else {
1014 // state file has been read successfully
1015 (replicas[ir])->replica_state_file_in_sync = true;
1016 (replicas[ir])->update_status = 0;
1017 }
1018 is.close();
1019 }
1020
1021 // now read the hills added after writing the state file
1022 if ((replicas[ir])->replica_hills_file.size()) {
1023
1024 if (cvm::debug())
1025 cvm::log ("Metadynamics bias \""+this->name+"\""+
1026 ": checking for new hills from replica \""+
1027 (replicas[ir])->replica_id+"\" in the file \""+
1028 (replicas[ir])->replica_hills_file+"\".\n");
1029
1030 // read hills from the other replicas' files; for each file, resume
1031 // the position recorded previously
1032
1033 std::ifstream is ((replicas[ir])->replica_hills_file.c_str());
1034 if (is.good()) {
1035
1036 // try to resume the previous position
1037 is.seekg ((replicas[ir])->replica_hills_file_pos, std::ios::beg);
1038 if (!is.good()){
1039 // if fail (the file may have been overwritten), reset this
1040 // position
1041 is.clear();
1042 is.seekg (0, std::ios::beg);
1043 // reset the counter
1044 (replicas[ir])->replica_hills_file_pos = 0;
1045 // schedule to reread the state file
1046 (replicas[ir])->replica_state_file_in_sync = false;
1047 // and record the failure
1048 (replicas[ir])->update_status++;
1049 cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+
1050 "\" at the previous position: will try again in "+
1051 cvm::to_str (replica_update_freq)+" steps.\n");
1052 } else {
1053
1054 while ((replicas[ir])->read_hill (is)) {
1055 // if (cvm::debug())
1056 cvm::log ("Metadynamics bias \""+this->name+"\""+
1057 ": received a hill from replica \""+
1058 (replicas[ir])->replica_id+
1059 "\" at step "+
1060 cvm::to_str (((replicas[ir])->hills.back()).it)+".\n");
1061 }
1062 is.clear();
1063 // store the position for the next read
1064 (replicas[ir])->replica_hills_file_pos = is.tellg();
1065 if (cvm::debug())
1066 cvm::log ("Metadynamics bias \""+this->name+"\""+
1067 ": stopped reading file \""+(replicas[ir])->replica_hills_file+
1068 "\" at position "+
1069 cvm::to_str ((replicas[ir])->replica_hills_file_pos)+".\n");
1070
1071 // test whether this is the end of the file
1072 is.seekg (0, std::ios::end);
1073 if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) {
1074 (replicas[ir])->update_status++;
1075 } else {
1076 (replicas[ir])->update_status = 0;
1077 }
1078 }
1079
1080 } else {
1081 cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+
1082 "\": will try again in "+
1083 cvm::to_str (replica_update_freq)+" steps.\n");
1084 (replicas[ir])->update_status++;
1085 // cvm::fatal_error ("Error: cannot read from file \""+
1086 // (replicas[ir])->replica_hills_file+"\".\n");
1087 }
1088 is.close();
1089 }
1090
1091 size_t const n_flush = (replica_update_freq/new_hill_freq + 1);
1092 if ((replicas[ir])->update_status > 3*n_flush) {
1093 // TODO: suspend the calculation?
1094 cvm::log ("WARNING: in metadynamics bias \""+this->name+"\""+
1095 " failed to read completely the output of replica \""+
1096 (replicas[ir])->replica_id+
1097 "\" after more than "+
1098 cvm::to_str ((replicas[ir])->update_status * replica_update_freq)+
1099 " steps. Ensure that it is still running.\n");
1100 }
1101 }
1102 }
1103
1104
1105 // **********************************************************************
1106 // input functions
1107 // **********************************************************************
1108
1109
read_restart(std::istream & is)1110 std::istream & colvarbias_meta::read_restart (std::istream& is)
1111 {
1112 size_t const start_pos = is.tellg();
1113
1114 if (comm == single_replica) {
1115 // if using a multiple replicas scheme, output messages
1116 // are printed before and after calling this function
1117 cvm::log ("Restarting metadynamics bias \""+this->name+"\""+
1118 ".\n");
1119 }
1120 std::string key, brace, conf;
1121 if ( !(is >> key) || !(key == "metadynamics") ||
1122 !(is >> brace) || !(brace == "{") ||
1123 !(is >> colvarparse::read_block ("configuration", conf)) ) {
1124
1125 if (comm == single_replica)
1126 cvm::log ("Error: in reading restart configuration for metadynamics bias \""+
1127 this->name+"\""+
1128 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
1129 (replica_state_file_in_sync ? ("at position "+
1130 cvm::to_str (start_pos)+
1131 " in the state file") : "")+".\n");
1132 is.clear();
1133 is.seekg (start_pos, std::ios::beg);
1134 is.setstate (std::ios::failbit);
1135 return is;
1136 }
1137
1138 std::string name = "";
1139 if ( colvarparse::get_keyval (conf, "name", name,
1140 std::string (""), colvarparse::parse_silent) &&
1141 (name != this->name) )
1142 cvm::fatal_error ("Error: in the restart file, the "
1143 "\"metadynamics\" block has a different name: different system?\n");
1144
1145 if (name.size() == 0) {
1146 cvm::fatal_error ("Error: \"metadynamics\" block within the restart file "
1147 "has no identifiers.\n");
1148 }
1149
1150 if (comm != single_replica) {
1151 std::string replica = "";
1152 if (colvarparse::get_keyval (conf, "replicaID", replica,
1153 std::string (""), colvarparse::parse_silent) &&
1154 (replica != this->replica_id))
1155 cvm::fatal_error ("Error: in the restart file, the "
1156 "\"metadynamics\" block has a different replicaID: different system?\n");
1157
1158 colvarparse::get_keyval (conf, "step", state_file_step,
1159 cvm::step_absolute(), colvarparse::parse_silent);
1160 }
1161
1162 bool grids_from_restart_file = use_grids;
1163
1164 if (use_grids) {
1165
1166 if (expand_grids) {
1167 // the boundaries of the colvars may have been changed; TODO:
1168 // this reallocation is only for backward-compatibility, and may
1169 // be deleted when grid_parameters (i.e. colvargrid's own
1170 // internal reallocation) has kicked in
1171 delete hills_energy;
1172 delete hills_energy_gradients;
1173 hills_energy = new colvar_grid_scalar (colvars);
1174 hills_energy_gradients = new colvar_grid_gradient (colvars);
1175 }
1176
1177 colvar_grid_scalar *hills_energy_backup = NULL;
1178 colvar_grid_gradient *hills_energy_gradients_backup = NULL;
1179
1180 if (has_data) {
1181 if (cvm::debug())
1182 cvm::log ("Backupping grids for metadynamics bias \""+
1183 this->name+"\""+
1184 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
1185 hills_energy_backup = hills_energy;
1186 hills_energy_gradients_backup = hills_energy_gradients;
1187 hills_energy = new colvar_grid_scalar (colvars);
1188 hills_energy_gradients = new colvar_grid_gradient (colvars);
1189 }
1190
1191 size_t const hills_energy_pos = is.tellg();
1192 if (!(is >> key)) {
1193 if (hills_energy_backup != NULL) {
1194 delete hills_energy;
1195 delete hills_energy_gradients;
1196 hills_energy = hills_energy_backup;
1197 hills_energy_gradients = hills_energy_gradients_backup;
1198 }
1199 is.clear();
1200 is.seekg (hills_energy_pos, std::ios::beg);
1201 is.setstate (std::ios::failbit);
1202 return is;
1203 } else if (!(key == std::string ("hills_energy")) ||
1204 !(hills_energy->read_restart (is))) {
1205 is.clear();
1206 is.seekg (hills_energy_pos, std::ios::beg);
1207 grids_from_restart_file = false;
1208 if (!rebin_grids) {
1209 if (hills_energy_backup == NULL)
1210 cvm::fatal_error ("Error: couldn't read the free energy grid for metadynamics bias \""+
1211 this->name+"\""+
1212 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
1213 "; if useGrids was off when the state file was written, "
1214 "enable rebinGrids now to regenerate the grids.\n");
1215 else {
1216 if (comm == single_replica)
1217 cvm::log ("Error: couldn't read the free energy grid for metadynamics bias \""+
1218 this->name+"\".\n");
1219 delete hills_energy;
1220 delete hills_energy_gradients;
1221 hills_energy = hills_energy_backup;
1222 hills_energy_gradients = hills_energy_gradients_backup;
1223 is.setstate (std::ios::failbit);
1224 return is;
1225 }
1226 }
1227 }
1228
1229 size_t const hills_energy_gradients_pos = is.tellg();
1230 if (!(is >> key)) {
1231 if (hills_energy_backup != NULL) {
1232 delete hills_energy;
1233 delete hills_energy_gradients;
1234 hills_energy = hills_energy_backup;
1235 hills_energy_gradients = hills_energy_gradients_backup;
1236 }
1237 is.clear();
1238 is.seekg (hills_energy_gradients_pos, std::ios::beg);
1239 is.setstate (std::ios::failbit);
1240 return is;
1241 } else if (!(key == std::string ("hills_energy_gradients")) ||
1242 !(hills_energy_gradients->read_restart (is))) {
1243 is.clear();
1244 is.seekg (hills_energy_gradients_pos, std::ios::beg);
1245 grids_from_restart_file = false;
1246 if (!rebin_grids) {
1247 if (hills_energy_backup == NULL)
1248 cvm::fatal_error ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
1249 this->name+"\""+
1250 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
1251 "; if useGrids was off when the state file was written, "
1252 "enable rebinGrids now to regenerate the grids.\n");
1253 else {
1254 if (comm == single_replica)
1255 cvm::log ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
1256 this->name+"\".\n");
1257 delete hills_energy;
1258 delete hills_energy_gradients;
1259 hills_energy = hills_energy_backup;
1260 hills_energy_gradients = hills_energy_gradients_backup;
1261 is.setstate (std::ios::failbit);
1262 return is;
1263 }
1264 }
1265 }
1266
1267 if (cvm::debug())
1268 cvm::log ("Successfully read new grids for bias \""+
1269 this->name+"\""+
1270 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
1271
1272 if (hills_energy_backup != NULL) {
1273 // now that we have successfully updated the grids, delete the
1274 // backup copies
1275 if (cvm::debug())
1276 cvm::log ("Deallocating the older grids.\n");
1277
1278 delete hills_energy_backup;
1279 delete hills_energy_gradients_backup;
1280 }
1281 }
1282
1283 bool const existing_hills = (hills.size() > 0);
1284 size_t const old_hills_size = hills.size();
1285 hill_iter old_hills_end = hills.end();
1286 hill_iter old_hills_off_grid_end = hills_off_grid.end();
1287
1288 // read the hills explicitly written (if there are any)
1289 while (read_hill (is)) {
1290 if (cvm::debug())
1291 cvm::log ("Read a previously saved hill under the "
1292 "metadynamics bias \""+
1293 this->name+"\", created at step "+
1294 cvm::to_str ((hills.back()).it)+".\n");
1295 }
1296 is.clear();
1297 new_hills_begin = hills.end();
1298 if (grids_from_restart_file) {
1299 if (hills.size() > old_hills_size)
1300 cvm::log ("Read "+cvm::to_str (hills.size())+
1301 " hills in addition to the grids.\n");
1302 } else {
1303 if (hills.size())
1304 cvm::log ("Read "+cvm::to_str (hills.size())+" hills.\n");
1305 }
1306
1307 if (rebin_grids) {
1308
1309 // allocate new grids (based on the new boundaries and widths just
1310 // read from the configuration file), and project onto them the
1311 // grids just read from the restart file
1312
1313 colvar_grid_scalar *new_hills_energy =
1314 new colvar_grid_scalar (colvars);
1315 colvar_grid_gradient *new_hills_energy_gradients =
1316 new colvar_grid_gradient (colvars);
1317
1318 if (!grids_from_restart_file || (keep_hills && hills.size())) {
1319 // if there are hills, recompute the new grids from them
1320 cvm::log ("Rebinning the energy and forces grids from "+
1321 cvm::to_str (hills.size())+" hills (this may take a while)...\n");
1322 project_hills (hills.begin(), hills.end(),
1323 new_hills_energy, new_hills_energy_gradients, true);
1324 cvm::log ("rebinning done.\n");
1325
1326 } else {
1327 // otherwise, use the grids in the restart file
1328 cvm::log ("Rebinning the energy and forces grids "
1329 "from the grids in the restart file.\n");
1330 new_hills_energy->map_grid (*hills_energy);
1331 new_hills_energy_gradients->map_grid (*hills_energy_gradients);
1332 }
1333
1334 delete hills_energy;
1335 delete hills_energy_gradients;
1336 hills_energy = new_hills_energy;
1337 hills_energy_gradients = new_hills_energy_gradients;
1338
1339 // assuming that some boundaries have expanded, eliminate those
1340 // off-grid hills that aren't necessary any more
1341 if (hills.size())
1342 recount_hills_off_grid (hills.begin(), hills.end(), hills_energy);
1343 }
1344
1345 if (use_grids) {
1346 if (hills_off_grid.size()) {
1347 cvm::log (cvm::to_str (hills_off_grid.size())+" hills are near the "
1348 "grid boundaries: they will be computed analytically "
1349 "and saved to the state files.\n");
1350 }
1351 }
1352
1353 is >> brace;
1354 if (brace != "}") {
1355 cvm::log ("Incomplete restart information for metadynamics bias \""+
1356 this->name+"\""+
1357 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
1358 ": no closing brace at position "+
1359 cvm::to_str (is.tellg())+" in the file.\n");
1360 is.setstate (std::ios::failbit);
1361 return is;
1362 }
1363
1364 if (cvm::debug())
1365 cvm::log ("colvarbias_meta::read_restart() done\n");
1366
1367 if (existing_hills) {
1368 hills.erase (hills.begin(), old_hills_end);
1369 hills_off_grid.erase (hills_off_grid.begin(), old_hills_off_grid_end);
1370 }
1371
1372 has_data = true;
1373
1374 if (comm != single_replica) {
1375 read_replica_files();
1376 }
1377
1378 return is;
1379 }
1380
1381
read_hill(std::istream & is)1382 std::istream & colvarbias_meta::read_hill (std::istream &is)
1383 {
1384 if (!is) return is; // do nothing if failbit is set
1385
1386 size_t const start_pos = is.tellg();
1387
1388 std::string data;
1389 if ( !(is >> read_block ("hill", data)) ) {
1390 is.clear();
1391 is.seekg (start_pos, std::ios::beg);
1392 is.setstate (std::ios::failbit);
1393 return is;
1394 }
1395
1396 size_t h_it;
1397 get_keyval (data, "step", h_it, 0, parse_silent);
1398 if (h_it <= state_file_step) {
1399 if (cvm::debug())
1400 cvm::log ("Skipping a hill older than the state file for metadynamics bias \""+
1401 this->name+"\""+
1402 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
1403 return is;
1404 }
1405
1406 cvm::real h_weight;
1407 get_keyval (data, "weight", h_weight, hill_weight, parse_silent);
1408
1409 std::vector<colvarvalue> h_centers (colvars.size());
1410 for (size_t i = 0; i < colvars.size(); i++) {
1411 h_centers[i].type ((colvars[i]->value()).type());
1412 }
1413 {
1414 // it is safer to read colvarvalue objects one at a time;
1415 // TODO: change this it later
1416 std::string centers_input;
1417 key_lookup (data, "centers", centers_input);
1418 std::istringstream centers_is (centers_input);
1419 for (size_t i = 0; i < colvars.size(); i++) {
1420 centers_is >> h_centers[i];
1421 }
1422 }
1423
1424 std::vector<cvm::real> h_widths (colvars.size());
1425 get_keyval (data, "widths", h_widths,
1426 std::vector<cvm::real> (colvars.size(), (std::sqrt (2.0 * PI) / 2.0)),
1427 parse_silent);
1428
1429 std::string h_replica = "";
1430 if (comm != single_replica) {
1431 get_keyval (data, "replicaID", h_replica, replica_id, parse_silent);
1432 if (h_replica != replica_id)
1433 cvm::fatal_error ("Error: trying to read a hill created by replica \""+h_replica+
1434 "\" for replica \""+replica_id+
1435 "\"; did you swap output files?\n");
1436 }
1437
1438 hill_iter const hills_end = hills.end();
1439 hills.push_back (hill (h_it, h_weight, h_centers, h_widths, h_replica));
1440 if (new_hills_begin == hills_end) {
1441 // if new_hills_begin is unset, set it for the first time
1442 new_hills_begin = hills.end();
1443 new_hills_begin--;
1444 }
1445
1446 if (use_grids) {
1447 // add this also to the list of hills that are off-grid, which will
1448 // be computed analytically
1449 cvm::real const min_dist =
1450 hills_energy->bin_distance_from_boundaries ((hills.back()).centers, true);
1451 if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
1452 hills_off_grid.push_back (hills.back());
1453 }
1454 }
1455
1456 has_data = true;
1457 return is;
1458 }
1459
1460
1461
1462
1463 // **********************************************************************
1464 // output functions
1465 // **********************************************************************
1466
write_restart(std::ostream & os)1467 std::ostream & colvarbias_meta::write_restart (std::ostream& os)
1468 {
1469 os << "metadynamics {\n"
1470 << " configuration {\n"
1471 << " step " << cvm::step_absolute() << "\n"
1472 << " name " << this->name << "\n";
1473 if (this->comm != single_replica)
1474 os << " replicaID " << this->replica_id << "\n";
1475 os << " }\n\n";
1476
1477 if (use_grids) {
1478
1479 // this is a very good time to project hills, if you haven't done
1480 // it already!
1481 project_hills (new_hills_begin, hills.end(),
1482 hills_energy, hills_energy_gradients);
1483 new_hills_begin = hills.end();
1484
1485 // write down the grids to the restart file
1486 os << " hills_energy\n";
1487 hills_energy->write_restart (os);
1488 os << " hills_energy_gradients\n";
1489 hills_energy_gradients->write_restart (os);
1490 }
1491
1492 if ( (!use_grids) || keep_hills ) {
1493 // write all hills currently in memory
1494 for (std::list<hill>::const_iterator h = this->hills.begin();
1495 h != this->hills.end();
1496 h++) {
1497 os << *h;
1498 }
1499 } else {
1500 // write just those that are near the grid boundaries
1501 for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
1502 h != this->hills_off_grid.end();
1503 h++) {
1504 os << *h;
1505 }
1506 }
1507
1508 os << "}\n\n";
1509
1510 if (comm != single_replica) {
1511 write_replica_state_file();
1512 // schedule to reread the state files of the other replicas (they
1513 // have also rewritten them)
1514 for (size_t ir = 0; ir < replicas.size(); ir++) {
1515 (replicas[ir])->replica_state_file_in_sync = false;
1516 }
1517 }
1518
1519 if (dump_fes) {
1520 write_pmf();
1521 }
1522
1523 return os;
1524 }
1525
1526
write_pmf()1527 void colvarbias_meta::write_pmf()
1528 {
1529 // allocate a new grid to store the pmf
1530 colvar_grid_scalar *pmf = new colvar_grid_scalar (*hills_energy);
1531 pmf->create();
1532
1533 std::string fes_file_name_prefix (cvm::output_prefix);
1534
1535 if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) {
1536 // if this is not the only free energy integrator, append
1537 // this bias's name, to distinguish it from the output of the other
1538 // biases producing a .pmf file
1539 // TODO: fix for ABF with updateBias == no
1540 fes_file_name_prefix += ("."+this->name);
1541 }
1542
1543 if ((comm == single_replica) || (dump_replica_fes)) {
1544 // output the PMF from this instance or replica
1545 pmf->reset();
1546 pmf->add_grid (*hills_energy);
1547 cvm::real const max = pmf->maximum_value();
1548 pmf->add_constant (-1.0 * max);
1549 pmf->multiply_constant (-1.0);
1550 if (well_tempered) {
1551 cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
1552 pmf->multiply_constant (well_temper_scale);
1553 }
1554 {
1555 std::string const fes_file_name (fes_file_name_prefix +
1556 ((comm != single_replica) ? ".partial" : "") +
1557 (dump_fes_save ?
1558 "."+cvm::to_str (cvm::step_absolute()) : "") +
1559 ".pmf");
1560 cvm::backup_file (fes_file_name.c_str());
1561 std::ofstream fes_os (fes_file_name.c_str());
1562 pmf->write_multicol (fes_os);
1563 fes_os.close();
1564 }
1565 }
1566 if (comm != single_replica) {
1567 // output the combined PMF from all replicas
1568 pmf->reset();
1569 pmf->add_grid (*hills_energy);
1570 for (size_t ir = 0; ir < replicas.size(); ir++) {
1571 pmf->add_grid (*(replicas[ir]->hills_energy));
1572 }
1573 cvm::real const max = pmf->maximum_value();
1574 pmf->add_constant (-1.0 * max);
1575 pmf->multiply_constant (-1.0);
1576 if (well_tempered) {
1577 cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
1578 pmf->multiply_constant (well_temper_scale);
1579 }
1580 std::string const fes_file_name (fes_file_name_prefix +
1581 (dump_fes_save ?
1582 "."+cvm::to_str (cvm::step_absolute()) : "") +
1583 ".pmf");
1584 cvm::backup_file (fes_file_name.c_str());
1585 std::ofstream fes_os (fes_file_name.c_str());
1586 pmf->write_multicol (fes_os);
1587 fes_os.close();
1588 }
1589
1590 delete pmf;
1591 }
1592
1593
1594
write_replica_state_file()1595 void colvarbias_meta::write_replica_state_file()
1596 {
1597 // write down also the restart for the other replicas: TODO: this
1598 // is duplicated code, that could be cleaned up later
1599 cvm::backup_file (replica_state_file.c_str());
1600 std::ofstream rep_state_os (replica_state_file.c_str());
1601 if (!rep_state_os.good())
1602 cvm::fatal_error ("Error: in opening file \""+
1603 replica_state_file+"\" for writing.\n");
1604
1605 rep_state_os.setf (std::ios::scientific, std::ios::floatfield);
1606 rep_state_os << "\n"
1607 << "metadynamics {\n"
1608 << " configuration {\n"
1609 << " name " << this->name << "\n"
1610 << " step " << cvm::step_absolute() << "\n";
1611 if (this->comm != single_replica) {
1612 rep_state_os << " replicaID " << this->replica_id << "\n";
1613 }
1614 rep_state_os << " }\n\n";
1615 rep_state_os << " hills_energy\n";
1616 rep_state_os << std::setprecision (cvm::cv_prec)
1617 << std::setw (cvm::cv_width);
1618 hills_energy->write_restart (rep_state_os);
1619 rep_state_os << " hills_energy_gradients\n";
1620 rep_state_os << std::setprecision (cvm::cv_prec)
1621 << std::setw (cvm::cv_width);
1622 hills_energy_gradients->write_restart (rep_state_os);
1623
1624 if ( (!use_grids) || keep_hills ) {
1625 // write all hills currently in memory
1626 for (std::list<hill>::const_iterator h = this->hills.begin();
1627 h != this->hills.end();
1628 h++) {
1629 rep_state_os << *h;
1630 }
1631 } else {
1632 // write just those that are near the grid boundaries
1633 for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
1634 h != this->hills_off_grid.end();
1635 h++) {
1636 rep_state_os << *h;
1637 }
1638 }
1639
1640 rep_state_os << "}\n\n";
1641 rep_state_os.close();
1642
1643 // reopen the hills file
1644 replica_hills_os.close();
1645 replica_hills_os.open (replica_hills_file.c_str());
1646 if (!replica_hills_os.good())
1647 cvm::fatal_error ("Error: in opening file \""+
1648 replica_hills_file+"\" for writing.\n");
1649 replica_hills_os.setf (std::ios::scientific, std::ios::floatfield);
1650 }
1651
output_traj()1652 std::string colvarbias_meta::hill::output_traj()
1653 {
1654 std::ostringstream os;
1655 os.setf (std::ios::fixed, std::ios::floatfield);
1656 os << std::setw (cvm::it_width) << it << " ";
1657
1658 os.setf (std::ios::scientific, std::ios::floatfield);
1659
1660 os << " ";
1661 for (size_t i = 0; i < centers.size(); i++) {
1662 os << " ";
1663 os << std::setprecision (cvm::cv_prec)
1664 << std::setw (cvm::cv_width) << centers[i];
1665 }
1666
1667 os << " ";
1668 for (size_t i = 0; i < widths.size(); i++) {
1669 os << " ";
1670 os << std::setprecision (cvm::cv_prec)
1671 << std::setw (cvm::cv_width) << widths[i];
1672 }
1673
1674 os << " ";
1675 os << std::setprecision (cvm::en_prec)
1676 << std::setw (cvm::en_width) << W << "\n";
1677
1678 return os.str();
1679 }
1680
1681
operator <<(std::ostream & os,colvarbias_meta::hill const & h)1682 std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
1683 {
1684 os.setf (std::ios::scientific, std::ios::floatfield);
1685
1686 os << "hill {\n";
1687 os << " step " << std::setw (cvm::it_width) << h.it << "\n";
1688 os << " weight "
1689 << std::setprecision (cvm::en_prec)
1690 << std::setw (cvm::en_width)
1691 << h.W << "\n";
1692
1693 if (h.replica.size())
1694 os << " replicaID " << h.replica << "\n";
1695
1696 os << " centers ";
1697 for (size_t i = 0; i < (h.centers).size(); i++) {
1698 os << " "
1699 << std::setprecision (cvm::cv_prec)
1700 << std::setw (cvm::cv_width)
1701 << h.centers[i];
1702 }
1703 os << "\n";
1704
1705 os << " widths ";
1706 for (size_t i = 0; i < (h.widths).size(); i++) {
1707 os << " "
1708 << std::setprecision (cvm::cv_prec)
1709 << std::setw (cvm::cv_width)
1710 << h.widths[i];
1711 }
1712 os << "\n";
1713
1714 os << "}\n";
1715
1716 return os;
1717 }
1718