00001 #include <iostream>
00002 #include <sstream>
00003 #include <fstream>
00004 #include <iomanip>
00005 #include <cmath>
00006 #include <algorithm>
00007
00008
00009 #if defined (WIN32) && !defined(__CYGWIN__)
00010 #include <direct.h>
00011 #define CHDIR ::_chdir
00012 #define GETCWD ::_getcwd
00013 #define PATHSEP "\\"
00014 #else
00015 #include <unistd.h>
00016 #define CHDIR ::chdir
00017 #define GETCWD ::getcwd
00018 #define PATHSEP "/"
00019 #endif
00020
00021
00022 #include "colvar.h"
00023 #include "colvarbias_meta.h"
00024
00025
00026 colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
00027 : colvarbias (conf, key),
00028 comm (single_replica),
00029 new_hills_begin (hills.end())
00030
00031
00032 {
00033 if (cvm::n_abf_biases > 0)
00034 cvm::log ("Warning: ABF and metadynamics running at the "
00035 "same time can give inconsistent results.\n");
00036
00037 get_keyval (conf, "hillWeight", hill_weight, 0.1);
00038 if (hill_weight == 0.0)
00039 cvm::log ("Warning: zero weight specified, "
00040 "this bias will have no effect.\n");
00041
00042 get_keyval (conf, "newHillFrequency", new_hill_freq, 1000);
00043
00044 get_keyval (conf, "hillWidth", hill_width, ::sqrt (2.0 * PI) / 2.0);
00045
00046 {
00047 bool b_replicas = false;
00048 get_keyval (conf, "multipleReplicas", b_replicas, false);
00049 if (b_replicas)
00050 comm = multiple_replicas;
00051 else
00052 comm = single_replica;
00053 }
00054
00055 get_keyval (conf, "useGrids", use_grids, true);
00056 if (use_grids && (comm != single_replica)) {
00057 cvm::log ("Warning: calculations with a grid and "
00058 "multiple replicas is currently not supported: "
00059 "setting useGrids to \"no\".\n");
00060 use_grids = false;
00061 }
00062
00063 if (use_grids) {
00064 get_keyval (conf, "gridsUpdateFrequency", grids_freq, new_hill_freq);
00065
00066 get_keyval (conf, "rebinGrids", rebin_grids, false);
00067
00068 expand_grids = false;
00069 for (size_t i = 0; i < colvars.size(); i++) {
00070 if (colvars[i]->expand_boundaries) {
00071 expand_grids = true;
00072 cvm::log ("Will expand the metadynamics grid when the colvar \""+
00073 colvars[i]->name+"\" approaches the boundaries.\n");
00074 }
00075 }
00076
00077
00078
00079 get_keyval (conf, "keepHills", keep_hills, false);
00080 get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true);
00081 get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false);
00082
00083 for (size_t i = 0; i < colvars.size(); i++) {
00084 colvars[i]->enable (colvar::task_grid);
00085 }
00086
00087 hills_energy = new colvar_grid_scalar (colvars);
00088 hills_energy_gradients = new colvar_grid_gradient (colvars);
00089 }
00090
00091 if (comm != single_replica) {
00092
00093 get_keyval (conf, "replicaID", replica, std::string (""));
00094 if (!replica.size())
00095 cvm::fatal_error ("Error: you must define an id for this replica "
00096 "when using more than one.\n");
00097
00098 get_keyval (conf, "replicaFilesRegistry",
00099 replica_files_registry,
00100 (this->name+".replica_files.txt"));
00101
00102 get_keyval (conf, "replicaUpdateFrequency",
00103 replica_update_freq, new_hill_freq);
00104
00105 char *pwd = new char[321];
00106 if (GETCWD (pwd, 320) == NULL)
00107 cvm::fatal_error ("Error: cannot read the current working directory.\n");
00108 replica_out_file_name =
00109 (std::string (pwd)+std::string (PATHSEP)+
00110 cvm::output_prefix+".colvars."+this->name+
00111 "."+replica+".hills");
00112 delete pwd;
00113
00114 replica_out_file.open (replica_out_file_name.c_str());
00115 if (!replica_out_file.good())
00116 cvm::fatal_error ("Error: in opening hills output file \""+
00117 replica_out_file_name+"\".\n");
00118 register_replica_file (replica_out_file_name);
00119 }
00120
00121 get_keyval (conf, "writeHillsTrajectory", b_hills_traj, false);
00122 if (b_hills_traj) {
00123
00124 std::string const traj_file_name (cvm::output_prefix+
00125 ".colvars."+this->name+
00126 ( (comm != single_replica) ?
00127 ("."+replica) :
00128 ("") )+
00129 ".hills.traj");
00130 hills_traj_os.open (traj_file_name.c_str());
00131 if (!hills_traj_os.good())
00132 cvm::fatal_error ("Error: in opening hills output file \"" +
00133 traj_file_name + "\".\n");
00134 }
00135
00136 if (cvm::debug())
00137 cvm::log ("Done initializing the metadynamics bias \""+this->name+"\".\n");
00138
00139 save_delimiters = false;
00140 }
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 colvarbias_meta::~colvarbias_meta()
00201 {
00202 if (hills_energy) {
00203 delete hills_energy;
00204 hills_energy = NULL;
00205 }
00206
00207 if (hills_energy_gradients) {
00208 delete hills_energy_gradients;
00209 hills_energy_gradients = NULL;
00210 }
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232 if (replica_out_file.good())
00233 replica_out_file.close();
00234
00235 if (hills_traj_os.good())
00236 hills_traj_os.close();
00237 }
00238
00239
00240
00241
00242
00243
00244
00245 std::list<colvarbias_meta::hill>::const_iterator
00246 colvarbias_meta::create_hill (colvarbias_meta::hill const &h)
00247 {
00248 hill_iter const hills_end = hills.end();
00249 hills.push_back (h);
00250 if (new_hills_begin == hills_end) {
00251
00252 new_hills_begin = hills.end();
00253 new_hills_begin--;
00254 }
00255
00256 if (use_grids) {
00257
00258
00259
00260
00261 cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers);
00262 if (min_dist < (3.0 * ::floor (hill_width)) + 1.0) {
00263 hills_off_grid.push_back (hills.back());
00264 }
00265 }
00266
00267
00268 if (hills_traj_os.good()) {
00269 hills_traj_os << (hills.back()).output_traj();
00270 if (cvm::debug()) {
00271 hills_traj_os.flush();
00272 }
00273 }
00274
00275
00276 if (comm != single_replica) {
00277 if (replica_out_file.good()) {
00278 replica_out_file << hills.back();
00279 } else {
00280 cvm::fatal_error ("Error in writing to the hills database for replica \""+
00281 replica+"\".\n");
00282 }
00283 }
00284
00285 return hills.end();
00286 }
00287
00288
00289 std::list<colvarbias_meta::hill>::const_iterator
00290 colvarbias_meta::delete_hill (hill_iter &h)
00291 {
00292 if (cvm::debug()) {
00293 cvm::log ("Deleting hill from the metadynamics bias \""+this->name+
00294 "\", with step number "+
00295 cvm::to_str (h->it)+(h->replica.size() ?
00296 ", replica id \""+h->replica :
00297 "")+".\n");
00298 }
00299
00300 if (use_grids && hills_off_grid.size()) {
00301 for (hill_iter hoff = hills_off_grid.begin();
00302 hoff != hills_off_grid.end(); hoff++) {
00303 if (*h == *hoff) {
00304 hills_off_grid.erase (hoff);
00305 break;
00306 }
00307 }
00308 }
00309
00310 if (hills_traj_os.good()) {
00311
00312 hills_traj_os << "# DELETED this hill: "
00313 << (hills.back()).output_traj()
00314 << "\n";
00315 if (cvm::debug())
00316 hills_traj_os.flush();
00317 }
00318
00319 return hills.erase (h);
00320 }
00321
00322
00323 void colvarbias_meta::update()
00324 {
00325 if (cvm::debug())
00326 cvm::log ("Updating the metadynamics bias \""+this->name+"\".\n");
00327
00328 if (use_grids) {
00329
00330 std::vector<int> curr_bin = hills_energy->get_colvars_index();
00331
00332 if (expand_grids) {
00333
00334
00335 if (cvm::debug())
00336 cvm::log ("Current coordinates on the grid: "+
00337 cvm::to_str (curr_bin)+".\n");
00338
00339 if (cvm::debug())
00340 cvm::log ("Checking if the grid is big enough.\n");
00341
00342 bool changed_grids = false;
00343 int const min_buffer =
00344 (3 * (size_t) ::floor (hill_width)) + 1;
00345
00346 std::vector<int> new_sizes (hills_energy->sizes());
00347 std::vector<colvarvalue> new_lower_boundaries (hills_energy->lower_boundaries);
00348 std::vector<colvarvalue> new_upper_boundaries (hills_energy->upper_boundaries);
00349
00350 for (size_t i = 0; i < colvars.size(); i++) {
00351
00352 if (! colvars[i]->expand_boundaries)
00353 continue;
00354
00355 cvm::real &new_lb = new_lower_boundaries[i].real_value;
00356 cvm::real &new_ub = new_upper_boundaries[i].real_value;
00357 int &new_size = new_sizes[i];
00358 bool changed_lb = false, changed_ub = false;
00359
00360 if (curr_bin[i] < min_buffer) {
00361 int const extra_points = (min_buffer - curr_bin[i]);
00362 new_lb -= extra_points * colvars[i]->width;
00363 new_size += extra_points;
00364
00365
00366 curr_bin[i] += extra_points;
00367
00368 changed_lb = true;
00369 cvm::log ("Metadynamics"+((cvm::n_meta_biases > 1) ? " \""+name+"\"" : "")+
00370 ": new lower boundary for colvar \""+
00371 colvars[i]->name+"\", at "+
00372 cvm::to_str (new_lower_boundaries[i])+".\n");
00373 }
00374
00375 if (curr_bin[i] > new_size - min_buffer - 1) {
00376 int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
00377 new_ub += extra_points * colvars[i]->width;
00378 new_size += extra_points;
00379
00380 changed_ub = true;
00381 cvm::log ("Metadynamics"+((cvm::n_meta_biases > 1) ? " \""+name+"\"" : "")+
00382 ": new upper boundary for colvar \""+
00383 colvars[i]->name+"\", at "+
00384 cvm::to_str (new_upper_boundaries[i])+".\n");
00385 }
00386
00387 if (changed_lb || changed_ub)
00388 changed_grids = true;
00389 }
00390
00391 if (changed_grids) {
00392
00393
00394 if (cvm::debug())
00395 cvm::log ("Expanding the grids for the metadynamics bias \""+
00396 name+"\".\n");
00397
00398 colvar_grid_scalar *new_hills_energy =
00399 new colvar_grid_scalar (*hills_energy);
00400 colvar_grid_gradient *new_hills_energy_gradients =
00401 new colvar_grid_gradient (*hills_energy_gradients);
00402
00403
00404
00405 new_hills_energy->lower_boundaries = new_lower_boundaries;
00406 new_hills_energy->upper_boundaries = new_upper_boundaries;
00407 new_hills_energy->create (new_sizes, 0.0, 1);
00408
00409 new_hills_energy_gradients->lower_boundaries = new_lower_boundaries;
00410 new_hills_energy_gradients->upper_boundaries = new_upper_boundaries;
00411 new_hills_energy_gradients->create (new_sizes, 0.0, colvars.size());
00412
00413 new_hills_energy->map_grid (*hills_energy);
00414 new_hills_energy_gradients->map_grid (*hills_energy_gradients);
00415
00416 delete hills_energy;
00417 delete hills_energy_gradients;
00418 hills_energy = new_hills_energy;
00419 hills_energy_gradients = new_hills_energy_gradients;
00420
00421 curr_bin = hills_energy->get_colvars_index();
00422 if (cvm::debug())
00423 cvm::log ("Coordinates on the new grid: "+
00424 cvm::to_str (curr_bin)+".\n");
00425 }
00426 }
00427 }
00428
00429
00430 if ((cvm::it % new_hill_freq) == 0) {
00431
00432 if (cvm::debug())
00433 cvm::log ("Adding a new hill under the bias \""+
00434 this->name+"\", at step "+cvm::to_str (cvm::it)+".\n");
00435
00436 switch (comm) {
00437
00438 case single_replica:
00439 create_hill (hill (hill_weight, colvars, hill_width));
00440 break;
00441
00442 case multiple_replicas:
00443 create_hill (hill (hill_weight, colvars, hill_width, replica));
00444 break;
00445 }
00446 }
00447
00448
00449 if (comm != single_replica) {
00450 if ((cvm::it % replica_update_freq) == 0) {
00451
00452
00453 replica_out_file.flush();
00454
00455
00456 update_replica_files_registry();
00457 read_replica_files();
00458 }
00459 }
00460
00461
00462 colvar_energy = 0.0;
00463 for (size_t i = 0; i < colvars.size(); i++) {
00464 colvar_forces[i].reset();
00465 }
00466
00467
00468 if (use_grids) {
00469
00470
00471
00472 if ((cvm::step_relative() % grids_freq) == 0) {
00473
00474 project_hills (new_hills_begin, hills.end(),
00475 hills_energy, hills_energy_gradients);
00476 new_hills_begin = hills.end();
00477 }
00478
00479 std::vector<int> curr_bin = hills_energy->get_colvars_index();
00480 if (cvm::debug())
00481 cvm::log ("Current coordinates on the grid: "+
00482 cvm::to_str (curr_bin)+".\n");
00483
00484 if (hills_energy->index_ok (curr_bin)) {
00485
00486
00487
00488 colvar_energy += hills_energy->value (curr_bin);
00489 for (size_t i = 0; i < colvars.size(); i++) {
00490 cvm::real const *f = &(hills_energy_gradients->value (curr_bin));
00491 colvar_forces[i].real_value += -1.0 * f[i];
00492
00493
00494 }
00495
00496 } else {
00497
00498
00499
00500
00501 calc_hills (hills_off_grid.begin(), hills_off_grid.end(), colvar_energy);
00502 for (size_t i = 0; i < colvars.size(); i++) {
00503 calc_hills_force (i, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces);
00504 }
00505 }
00506
00507 } else {
00508
00509
00510 calc_hills (hills.begin(), hills.end(), colvar_energy);
00511
00512
00513
00514 for (size_t i = 0; i < colvars.size(); i++) {
00515 calc_hills_force (i, hills.begin(), hills.end(), colvar_forces);
00516 }
00517 }
00518
00519 if (cvm::debug())
00520 cvm::log ("Hills energy = "+cvm::to_str (colvar_energy)+
00521 ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
00522 }
00523
00524
00525 void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter h_first,
00526 colvarbias_meta::hill_iter h_last,
00527 cvm::real &energy,
00528 std::vector<colvarvalue> const &colvar_values)
00529 {
00530 std::vector<colvarvalue> curr_values (colvars.size());
00531 for (size_t i = 0; i < colvars.size(); i++) {
00532 curr_values[i].type (colvars[i]->type());
00533 }
00534
00535 if (colvar_values.size()) {
00536 for (size_t i = 0; i < colvars.size(); i++) {
00537 curr_values[i] = colvar_values[i];
00538 }
00539 } else {
00540 for (size_t i = 0; i < colvars.size(); i++) {
00541 curr_values[i] = colvars[i]->value();
00542 }
00543 }
00544
00545 for (hill_iter h = h_first; h != h_last; h++) {
00546
00547
00548 cvm::real cv_sqdev = 0.0;
00549 for (size_t i = 0; i < colvars.size(); i++) {
00550 colvarvalue const &x = curr_values[i];
00551 colvarvalue const ¢er = h->centers[i];
00552 cvm::real const half_width = 0.5 * h->widths[i];
00553 cv_sqdev += (colvars[i]->dist2 (x, center)) / (half_width*half_width);
00554 }
00555
00556
00557 if (cv_sqdev > 23.0) {
00558
00559 h->value (0.0);
00560 } else {
00561 h->value (::exp (-0.5*cv_sqdev));
00562 }
00563 energy += h->energy();
00564 }
00565 }
00566
00567
00568 void colvarbias_meta::calc_hills_force (size_t const &i,
00569 colvarbias_meta::hill_iter h_first,
00570 colvarbias_meta::hill_iter h_last,
00571 std::vector<colvarvalue> &forces,
00572 std::vector<colvarvalue> const &values)
00573 {
00574
00575 colvarvalue x (values.size() ? values[i].type() : colvars[i]->type());
00576 x = (values.size() ? values[i] : colvars[i]->value());
00577
00578
00579
00580
00581
00582 switch (colvars[i]->type()) {
00583
00584 case colvarvalue::type_scalar:
00585 for (hill_iter h = h_first; h != h_last; h++) {
00586 if (h->value() == 0.0) continue;
00587 colvarvalue const ¢er = h->centers[i];
00588 cvm::real const half_width = 0.5 * h->widths[i];
00589 forces[i].real_value +=
00590 ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
00591 (colvars[i]->dist2_lgrad (x, center)).real_value );
00592 }
00593 break;
00594
00595 case colvarvalue::type_vector:
00596 case colvarvalue::type_unitvector:
00597 for (hill_iter h = h_first; h != h_last; h++) {
00598 if (h->value() == 0.0) continue;
00599 colvarvalue const ¢er = h->centers[i];
00600 cvm::real const half_width = 0.5 * h->widths[i];
00601 forces[i].rvector_value +=
00602 ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
00603 (colvars[i]->dist2_lgrad (x, center)).rvector_value );
00604 }
00605 break;
00606
00607 case colvarvalue::type_quaternion:
00608 for (hill_iter h = h_first; h != h_last; h++) {
00609 if (h->value() == 0.0) continue;
00610 colvarvalue const ¢er = h->centers[i];
00611 cvm::real const half_width = 0.5 * h->widths[i];
00612 forces[i].quaternion_value +=
00613 ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
00614 (colvars[i]->dist2_lgrad (x, center)).quaternion_value );
00615 }
00616 break;
00617
00618 case colvarvalue::type_notset:
00619 break;
00620 }
00621 }
00622
00623
00624
00625
00626
00627
00628
00629 void colvarbias_meta::register_replica_file (std::string const &new_file)
00630 {
00631 std::ifstream reg_is (replica_files_registry.c_str());
00632 if (reg_is.good()) {
00633 std::string s;
00634 while (reg_is >> s)
00635 if (new_file == s) {
00636
00637 return;
00638 }
00639 }
00640 reg_is.close();
00641
00642 std::ofstream reg_os (replica_files_registry.c_str(), std::ios::app);
00643 reg_os << new_file << "\n";
00644 reg_os.close();
00645 }
00646
00647
00648 void colvarbias_meta::update_replica_files_registry()
00649 {
00650 if (cvm::debug())
00651 cvm::log ("Updating the list of replica files, currently containing "+
00652 cvm::to_str (replica_files.size())+" elements.\n");
00653
00654
00655 std::ifstream reg_is (replica_files_registry.c_str());
00656
00657 if (reg_is.good()) {
00658
00659 std::string s ("");
00660 while ((reg_is >> s) && s.size()) {
00661
00662 if (s == replica_out_file_name) {
00663
00664 if (cvm::debug())
00665 cvm::log ("Skipping this job's replica file, \""+s+"\".\n");
00666 s = "";
00667 continue;
00668 }
00669
00670 bool already_loaded = false;
00671 std::list<std::string>::iterator rsi = replica_files.begin();
00672 for ( ; rsi != replica_files.end(); rsi++) {
00673 if (s == *rsi) {
00674
00675 if (cvm::debug())
00676 cvm::log ("Skipping a replica file already loaded, \""+(*rsi)+"\".\n");
00677 already_loaded = true;;
00678 }
00679 }
00680
00681 if ( (replica_files.size() == 0) || (!already_loaded) ) {
00682
00683 if (cvm::debug())
00684 cvm::log ("New replica file found: \""+s+"\".\n");
00685 replica_files.push_back (s);
00686 replica_files_pos.push_back (0);
00687 }
00688
00689 }
00690
00691 s = "";
00692 } else {
00693 cvm::fatal_error ("Error: cannot read the replica files registry, \""+
00694 replica_files_registry+"\".\n");
00695 }
00696
00697 reg_is.close();
00698
00699 if (cvm::debug())
00700 cvm::log ("The list of replica files now contains "+
00701 cvm::to_str (replica_files.size())+" elements.\n");
00702 }
00703
00704
00705 void colvarbias_meta::read_replica_files()
00706 {
00707
00708
00709
00710 std::list<std::string>::iterator rsi = replica_files.begin();
00711 std::list<size_t>::iterator pi = replica_files_pos.begin();
00712 for ( ; rsi != replica_files.end(); pi++, rsi++) {
00713
00714 if (cvm::debug())
00715 cvm::log ("Checking for new hills in the replica file \""+
00716 (*rsi)+"\".\n");
00717
00718 std::ifstream is ((*rsi).c_str());
00719 if (is.good()) {
00720
00721 is.seekg ((*pi), std::ios::beg);
00722 while (read_hill (is)) {
00723
00724 if (cvm::debug())
00725 cvm::log ("Read a previously saved hill under the "
00726 "metadynamics bias \""+
00727 this->name+"\", created by replica "+
00728 (hills.back()).replica+
00729 " at step "+
00730 cvm::to_str ((hills.back()).it)+".\n");
00731
00732
00733 if (hills_traj_os.good()) {
00734 hills_traj_os << (hills.back()).output_traj();
00735 if (cvm::debug()) {
00736 hills_traj_os.flush();
00737 }
00738 }
00739 }
00740 is.clear();
00741
00742 *pi = is.tellg();
00743 if (cvm::debug())
00744 cvm::log ("Stopped reading at position "+
00745 cvm::to_str (*pi)+".\n");
00746 } else {
00747 cvm::fatal_error ("Error: cannot read from file \""+(*rsi)+"\".\n");
00748 }
00749
00750 is.close();
00751 }
00752
00753 }
00754
00755
00756
00757
00758
00759
00760 std::istream & colvarbias_meta::read_restart (std::istream& is)
00761 {
00762 size_t const start_pos = is.tellg();
00763
00764 cvm::log ("Restarting metadynamics bias \""+
00765 this->name+"\".\n");
00766 std::string key, brace, conf;
00767 if ( !(is >> key) || !(key == "metadynamics") ||
00768 !(is >> brace) || !(brace == "{") ||
00769 !(is >> colvarparse::read_block ("configuration", conf)) ) {
00770
00771 cvm::log ("Error: in reading restart configuration for metadynamics bias \""+
00772 this->name+"\" at position "+
00773 cvm::to_str (start_pos)+" in stream.\n");
00774 is.clear();
00775 is.seekg (start_pos, std::ios::beg);
00776 is.setstate (std::ios::failbit);
00777 return is;
00778 }
00779
00780 std::string name = "";
00781 if ( colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent) &&
00782 (name != this->name) )
00783 cvm::fatal_error ("Error: in the restart file, the "
00784 "\"metadynamics\" block has a wrong name: different system?\n");
00785
00786 if (name.size() == 0) {
00787 cvm::fatal_error ("Error: \"metadynamics\" block in the restart file "
00788 "has no identifiers.\n");
00789 }
00790
00791 if (this->comm == single_replica) {
00792
00793 if (use_grids) {
00794
00795 if (expand_grids) {
00796
00797
00798
00799 delete hills_energy;
00800 delete hills_energy_gradients;
00801 hills_energy = new colvar_grid_scalar (colvars);
00802 hills_energy_gradients = new colvar_grid_gradient (colvars);
00803 }
00804
00805 if ( !(is >> key) ||
00806 !(key == std::string ("hills_energy")) ||
00807 !(hills_energy->read_restart (is)) ) {
00808 cvm::log ("Error: in reading restart information for metadynamics bias \""+
00809 this->name+"\".\n");
00810 is.clear();
00811 is.seekg (start_pos, std::ios::beg);
00812 is.setstate (std::ios::failbit);
00813 return is;
00814 }
00815
00816 if ( !(is >> key) ||
00817 !(key == std::string ("hills_energy_gradients")) ||
00818 !(hills_energy_gradients->read_restart (is)) ) {
00819 cvm::log ("Error: in reading restart information for metadynamics bias \""+
00820 this->name+"\".\n");
00821 is.clear();
00822 is.seekg (start_pos, std::ios::beg);
00823 is.setstate (std::ios::failbit);
00824 return is;
00825 }
00826
00827 }
00828
00829
00830 while (this->read_hill (is)) {
00831 if (cvm::debug())
00832 cvm::log ("Read a previously saved hill under the "
00833 "metadynamics bias \""+
00834 this->name+"\", created at step "+
00835 cvm::to_str ((hills.back()).it)+".\n");
00836 }
00837 is.clear();
00838 new_hills_begin = hills.end();
00839
00840 if (rebin_grids) {
00841
00842
00843
00844
00845 colvar_grid_scalar *new_hills_energy =
00846 new colvar_grid_scalar (colvars);
00847 colvar_grid_gradient *new_hills_energy_gradients =
00848 new colvar_grid_gradient (colvars);
00849
00850 if (keep_hills && (hills.size() > 0)) {
00851
00852
00853 cvm::log ("rebinGrids and keepHills defined, recomputing "
00854 "analytically the energy and force grids; this may take a while...\n");
00855 project_hills (hills.begin(), hills.end(),
00856 new_hills_energy, new_hills_energy_gradients);
00857 cvm::log ("rebinning done.\n");
00858
00859 } else {
00860
00861
00862
00863 cvm::log ("rebinGrids defined, mapping energy and forces grid to new grids.\n");
00864 new_hills_energy->map_grid (*hills_energy);
00865 new_hills_energy_gradients->map_grid (*hills_energy_gradients);
00866 }
00867
00868 delete hills_energy;
00869 delete hills_energy_gradients;
00870 hills_energy = new_hills_energy;
00871 hills_energy_gradients = new_hills_energy_gradients;
00872
00873 if (hills.size())
00874 recount_hills_off_grid (hills.begin(), hills.end(), hills_energy);
00875 }
00876
00877 } else {
00878
00879
00880 cvm::log ("Multiple replicas have been defined for the "
00881 "metadynamics bias \""+this->name+"\", skipping "
00882 "the restart file and reading hills from the replica "
00883 "files, whose list is contained in \""+
00884 this->replica_files_registry+"\".\n");
00885 this->update_replica_files_registry();
00886 this->read_replica_files();
00887 }
00888
00889 is >> brace;
00890 if (brace != "}") {
00891 cvm::fatal_error ("Error: corrupted restart information for metadynamics bias \""+
00892 this->name+"\": no matching brace at position "+
00893 cvm::to_str (is.tellg())+" in the restart file.\n");
00894 is.setstate (std::ios::failbit);
00895 }
00896
00897 if (cvm::debug())
00898 cvm::log ("colvarbias_meta::read_restart() done\n");
00899
00900 return is;
00901 }
00902
00903
00904 std::istream & colvarbias_meta::read_hill (std::istream &is)
00905 {
00906 if (!is) return is;
00907
00908 size_t const start_pos = is.tellg();
00909
00910 std::string data;
00911 if ( !(is >> read_block ("hill", data)) ) {
00912 is.clear();
00913 is.seekg (start_pos, std::ios::beg);
00914 is.setstate (std::ios::failbit);
00915 return is;
00916 }
00917
00918 std::string h_replica = "";
00919
00920 int h_it;
00921 get_keyval (data, "step", h_it, 0, parse_silent);
00922 cvm::real h_weight;
00923 get_keyval (data, "weight", h_weight, hill_weight, parse_silent);
00924
00925 std::vector<colvarvalue> h_centers (colvars.size());
00926 for (size_t i = 0; i < colvars.size(); i++) {
00927 h_centers[i].type ((colvars[i]->value()).type());
00928 }
00929 {
00930
00931
00932 std::string centers_input;
00933 key_lookup (data, "centers", centers_input);
00934 std::istringstream centers_is (centers_input);
00935 for (size_t i = 0; i < colvars.size(); i++) {
00936 centers_is >> h_centers[i];
00937 }
00938 }
00939
00940 std::vector<cvm::real> h_widths (colvars.size());
00941 get_keyval (data, "widths", h_widths, std::vector<cvm::real> (colvars.size(), 1.0), parse_silent);
00942
00943 if (comm != single_replica)
00944 get_keyval (data, "replica", h_replica, replica, parse_silent);
00945
00946 hill_iter const hills_end = hills.end();
00947 hills.push_back (hill (h_it, h_weight, h_centers, h_widths, h_replica));
00948 if (new_hills_begin == hills_end) {
00949
00950 new_hills_begin = hills.end();
00951 new_hills_begin--;
00952 }
00953
00954 if (use_grids) {
00955
00956
00957
00958 cvm::real const min_dist =
00959 hills_energy->bin_distance_from_boundaries ((hills.back()).centers);
00960 if (min_dist < (3.0 * ::floor (hill_width)) + 1.0) {
00961 hills_off_grid.push_back (hills.back());
00962 }
00963 }
00964
00965 return is;
00966 }
00967
00968
00969
00970
00971
00972
00973 std::ostream & colvarbias_meta::write_restart (std::ostream& os)
00974 {
00975 os << "metadynamics {\n"
00976 << " configuration {\n"
00977
00978 << " name " << this->name << "\n";
00979 if (this->comm != single_replica)
00980 os << " replica " << this->replica << "\n";
00981 os << " }\n";
00982
00983 if (this->comm == single_replica) {
00984
00985 if (use_grids) {
00986
00987 os << " hills_energy\n";
00988 hills_energy->write_restart (os);
00989 os << " hills_energy_gradients\n";
00990 hills_energy_gradients->write_restart (os);
00991
00992 if (dump_fes) {
00993 cvm::real const max = hills_energy->maximum_value();
00994 hills_energy->add_constant (-1.0 * max);
00995 hills_energy->multiply_constant (-1.0);
00996
00997
00998
00999 std::string const fes_file_name =
01000 ((cvm::n_meta_biases == 1) && (cvm::n_abf_biases == 0)) ?
01001 std::string (cvm::output_prefix+
01002 (dump_fes_save ? "."+cvm::to_str (cvm::step_absolute()) : "")+
01003 ".pmf") :
01004 std::string (cvm::output_prefix+"."+this->name+
01005 (dump_fes_save ? "."+cvm::to_str (cvm::step_absolute()) : "")+
01006 ".pmf");
01007 std::ofstream fes_os (fes_file_name.c_str());
01008 hills_energy->write_multicol (fes_os);
01009
01010
01011 hills_energy->multiply_constant (-1.0);
01012 hills_energy->add_constant (max);
01013 }
01014
01015 }
01016
01017 if ( (!use_grids) || keep_hills ) {
01018
01019
01020 for (std::list<hill>::const_iterator h = this->hills.begin();
01021 h != this->hills.end();
01022 h++) {
01023 os << *h;
01024 }
01025
01026 } else {
01027
01028
01029 for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
01030 h != this->hills_off_grid.end();
01031 h++) {
01032 os << *h;
01033 }
01034 }
01035
01036 } else {
01037 cvm::log ("Hills from the metadynamics bias \""+
01038 this->name+"\" have been previously saved in the file \""+
01039 replica_out_file_name+"\".\n");
01040 }
01041
01042 os << "}\n\n";
01043 return os;
01044 }
01045
01046
01047 std::string colvarbias_meta::hill::output_traj()
01048 {
01049 std::ostringstream os;
01050 os.setf (std::ios::fixed, std::ios::floatfield);
01051 os << std::setw (cvm::it_width) << it << " ";
01052
01053 os.setf (std::ios::scientific, std::ios::floatfield);
01054
01055 os << " ";
01056 for (size_t i = 0; i < centers.size(); i++) {
01057 os << " ";
01058 os << std::setprecision (cvm::cv_prec)
01059 << std::setw (cvm::cv_width) << centers[i];
01060 }
01061
01062 os << " ";
01063 for (size_t i = 0; i < widths.size(); i++) {
01064 os << " ";
01065 os << std::setprecision (cvm::cv_prec)
01066 << std::setw (cvm::cv_width) << widths[i];
01067 }
01068
01069 os << " ";
01070 os << std::setprecision (cvm::en_prec)
01071 << std::setw (cvm::en_width) << W << "\n";
01072
01073 return os.str();
01074 }
01075
01076
01077 std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
01078 {
01079 os.setf (std::ios::scientific, std::ios::floatfield);
01080
01081 os << "hill {\n";
01082 os << " step " << std::setw (cvm::it_width) << h.it << "\n";
01083 os << " weight "
01084 << std::setprecision (cvm::en_prec)
01085 << std::setw (cvm::en_width)
01086 << h.W << "\n";
01087
01088 if (h.replica.size())
01089 os << " replica " << h.replica << "\n";
01090
01091 os << " centers ";
01092 for (size_t i = 0; i < (h.centers).size(); i++) {
01093 os << " "
01094 << std::setprecision (cvm::cv_prec)
01095 << std::setw (cvm::cv_width)
01096 << h.centers[i];
01097 }
01098 os << "\n";
01099
01100 os << " widths ";
01101 for (size_t i = 0; i < (h.widths).size(); i++) {
01102 os << " "
01103 << std::setprecision (cvm::cv_prec)
01104 << std::setw (cvm::cv_width)
01105 << h.widths[i];
01106 }
01107 os << "\n";
01108
01109 os << "}\n";
01110
01111 return os;
01112 }
01113
01114
01115 void colvarbias_meta::project_hills (colvarbias_meta::hill_iter h_first,
01116 colvarbias_meta::hill_iter h_last,
01117 colvar_grid_scalar *he,
01118 colvar_grid_gradient *hg,
01119 cvm::real const scale_factor)
01120 {
01121 if (cvm::debug())
01122 cvm::log ("Projecting hills.\n");
01123
01124 std::vector<colvarvalue> colvar_values (colvars.size());
01125 std::vector<cvm::real> colvar_forces_scalar (colvars.size());
01126
01127 std::vector<int> he_ix = he->new_index();
01128 std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
01129 cvm::real hills_energy_here = 0.0;
01130 std::vector<colvarvalue> hills_forces_here (colvars.size(), 0.0);
01131
01132 for ( ;
01133 he->index_ok (he_ix) && ((hg != NULL) ? hg->index_ok (hg_ix) : true);
01134 ) {
01135
01136 for (size_t i = 0; i < colvars.size(); i++) {
01137 colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
01138 }
01139
01140
01141 hills_energy_here = 0.0;
01142 calc_hills (h_first, h_last, hills_energy_here, colvar_values);
01143 he->acc_value (he_ix, scale_factor * hills_energy_here);
01144
01145 if (hg != NULL) {
01146 for (size_t i = 0; i < colvars.size(); i++) {
01147 hills_forces_here[i].reset();
01148 calc_hills_force (i, h_first, h_last, hills_forces_here, colvar_values);
01149 colvar_forces_scalar[i] = scale_factor * hills_forces_here[i].real_value;
01150 }
01151 hg->acc_force (hg_ix, &(colvar_forces_scalar.front()));
01152 }
01153
01154 he->incr (he_ix);
01155 if (hg != NULL)
01156 hg->incr (hg_ix);
01157 }
01158
01159 if (! keep_hills) {
01160 hills.erase (hills.begin(), hills.end());
01161 }
01162 }
01163
01164
01165 void colvarbias_meta::recount_hills_off_grid (colvarbias_meta::hill_iter h_first,
01166 colvarbias_meta::hill_iter h_last,
01167 colvar_grid_scalar *he)
01168 {
01169 hills_off_grid.clear();
01170
01171 for (hill_iter h = h_first; h != h_last; h++) {
01172 cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers);
01173 if (min_dist < (3.0 * ::floor (hill_width)) + 1.0) {
01174 hills_off_grid.push_back (*h);
01175 }
01176 }
01177 }
01178
01179
01180 void colvarbias_meta::analyse()
01181 {}
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204
01205
01206
01207
01208
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221
01222
01223
01224
01225
01226
01227
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252
01253
01254
01255
01256
01257
01258
01259
01260
01261
01262
01263
01264
01265
01266
01267
01268
01269
01270
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
01285
01286
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297
01298
01299
01300
01301
01302
01303
01304
01305
01306
01307
01308
01309
01310
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363
01364
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376
01377
01378
01379
01380