colvarbias_meta.C

00001 // -*- c++ -*-
00002 
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/Colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 #include <iostream>
00011 #include <sstream>
00012 #include <fstream>
00013 #include <iomanip>
00014 #include <algorithm>
00015 
00016 // used to set the absolute path of a replica file
00017 #if defined(_WIN32) && !defined(__CYGWIN__)
00018 #include <direct.h>
00019 #define CHDIR ::_chdir
00020 #define GETCWD ::_getcwd
00021 #define PATHSEP "\\"
00022 #else
00023 #include <unistd.h>
00024 #define CHDIR ::chdir
00025 #define GETCWD ::getcwd
00026 #define PATHSEP "/"
00027 #endif
00028 
00029 #include "colvarmodule.h"
00030 #include "colvarproxy.h"
00031 #include "colvar.h"
00032 #include "colvarbias_meta.h"
00033 
00034 
00035 colvarbias_meta::colvarbias_meta(char const *key)
00036 : colvarbias(key), colvarbias_ti(key)
00037 {
00038 new_hills_begin = hills.end();
00039 
00040 hill_weight = 0.0;
00041 hill_width = 0.0;
00042 
00043 new_hill_freq = 1000;
00044 
00045 use_grids = true;
00046 grids_freq = 0;
00047 rebin_grids = false;
00048 hills_energy = NULL;
00049 hills_energy_gradients = NULL;
00050 
00051 dump_fes = true;
00052 keep_hills = false;
00053 restart_keep_hills = false;
00054 dump_fes_save = false;
00055 dump_replica_fes = false;
00056 
00057 b_hills_traj = false;
00058 
00059 ebmeta_equil_steps = 0L;
00060 
00061 replica_update_freq = 0;
00062 replica_id.clear();
00063 }
00064 
00065 
00066 int colvarbias_meta::init(std::string const &conf)
00067 {
00068 int error_code = COLVARS_OK;
00069 size_t i = 0;
00070 
00071 error_code |= colvarbias::init(conf);
00072 error_code |= colvarbias_ti::init(conf);
00073 
00074 cvm::main()->cite_feature("Metadynamics colvar bias implementation");
00075 
00076 enable(f_cvb_calc_pmf);
00077 
00078 get_keyval(conf, "hillWeight", hill_weight, hill_weight);
00079 if (hill_weight > 0.0) {
00080 enable(f_cvb_apply_force);
00081 } else {
00082 cvm::error("Error: hillWeight must be provided, and a positive number.\n", COLVARS_INPUT_ERROR);
00083 }
00084 
00085 get_keyval(conf, "newHillFrequency", new_hill_freq, new_hill_freq);
00086 if (new_hill_freq > 0) {
00087 enable(f_cvb_history_dependent);
00088 if (grids_freq == 0) {
00089 grids_freq = new_hill_freq;
00090 }
00091 }
00092 
00093 get_keyval(conf, "gaussianSigmas", colvar_sigmas, colvar_sigmas);
00094 
00095 get_keyval(conf, "hillWidth", hill_width, hill_width);
00096 
00097 if ((colvar_sigmas.size() > 0) && (hill_width > 0.0)) {
00098 error_code |= cvm::error("Error: hillWidth and gaussianSigmas are "
00099 "mutually exclusive.", COLVARS_INPUT_ERROR);
00100 }
00101 
00102 if (hill_width > 0.0) {
00103 colvar_sigmas.resize(num_variables());
00104 // Print the calculated sigma parameters
00105 cvm::log("Half-widths of the Gaussian hills (sigma's):\n");
00106 for (i = 0; i < num_variables(); i++) {
00107 colvar_sigmas[i] = variables(i)->width * hill_width / 2.0;
00108 cvm::log(variables(i)->name+std::string(": ")+
00109 cvm::to_str(colvar_sigmas[i]));
00110 }
00111 }
00112 
00113 if (colvar_sigmas.size() == 0) {
00114 error_code |= cvm::error("Error: positive values are required for "
00115 "either hillWidth or gaussianSigmas.",
00116 COLVARS_INPUT_ERROR);
00117 }
00118 
00119 {
00120 bool b_replicas = false;
00121 get_keyval(conf, "multipleReplicas", b_replicas, false);
00122 if (b_replicas) {
00123 cvm::main()->cite_feature("Multiple-walker metadynamics colvar bias implementation");
00124 comm = multiple_replicas;
00125 } else {
00126 comm = single_replica;
00127 }
00128 }
00129 
00130 get_keyval(conf, "useGrids", use_grids, use_grids);
00131 
00132 if (use_grids) {
00133 
00134 for (i = 0; i < num_variables(); i++) {
00135 if (2.0*colvar_sigmas[i] < variables(i)->width) {
00136 cvm::log("Warning: gaussianSigmas is too narrow for the grid "
00137 "spacing along "+variables(i)->name+".");
00138 }
00139 }
00140 
00141 get_keyval(conf, "gridsUpdateFrequency", grids_freq, grids_freq);
00142 get_keyval(conf, "rebinGrids", rebin_grids, rebin_grids);
00143 
00144 expand_grids = false;
00145 for (i = 0; i < num_variables(); i++) {
00146 variables(i)->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature
00147 if (variables(i)->expand_boundaries) {
00148 expand_grids = true;
00149 cvm::log("Metadynamics bias \""+this->name+"\""+
00150 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00151 ": Will expand grids when the colvar \""+
00152 variables(i)->name+"\" approaches its boundaries.\n");
00153 }
00154 }
00155 
00156 get_keyval(conf, "writeFreeEnergyFile", dump_fes, dump_fes);
00157 
00158 get_keyval(conf, "keepHills", keep_hills, keep_hills);
00159 get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
00160 
00161 if (hills_energy == NULL) {
00162 hills_energy = new colvar_grid_scalar(colvars);
00163 hills_energy_gradients = new colvar_grid_gradient(colvars);
00164 }
00165 
00166 } else {
00167 
00168 dump_fes = false;
00169 }
00170 
00171 get_keyval(conf, "writeHillsTrajectory", b_hills_traj, b_hills_traj);
00172 
00173 error_code |= init_replicas_params(conf);
00174 error_code |= init_well_tempered_params(conf);
00175 error_code |= init_ebmeta_params(conf);
00176 
00177 if (cvm::debug())
00178 cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+
00179 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
00180 
00181 return error_code;
00182 }
00183 
00184 
00185 int colvarbias_meta::init_replicas_params(std::string const &conf)
00186 {
00187 colvarproxy *proxy = cvm::main()->proxy;
00188 
00189 // in all cases, the first replica is this bias itself
00190 if (replicas.size() == 0) {
00191 replicas.push_back(this);
00192 }
00193 
00194 if (comm != single_replica) {
00195 
00196 if (!get_keyval(conf, "writePartialFreeEnergyFile",
00197 dump_replica_fes, dump_replica_fes)) {
00198 get_keyval(conf, "dumpPartialFreeEnergyFile", dump_replica_fes,
00199 dump_replica_fes, colvarparse::parse_silent);
00200 }
00201 
00202 if (dump_replica_fes && (! dump_fes)) {
00203 dump_fes = true;
00204 cvm::log("Enabling \"writeFreeEnergyFile\".\n");
00205 }
00206 
00207 get_keyval(conf, "replicaID", replica_id, replica_id);
00208 if (!replica_id.size()) {
00209 if (proxy->replica_enabled() == COLVARS_OK) {
00210 // Obtain replicaID from the communicator
00211 replica_id = cvm::to_str(proxy->replica_index());
00212 cvm::log("Setting replicaID from communication layer: replicaID = "+
00213 replica_id+".\n");
00214 } else {
00215 return cvm::error("Error: using more than one replica, but replicaID "
00216 "could not be obtained.\n", COLVARS_INPUT_ERROR);
00217 }
00218 }
00219 
00220 get_keyval(conf, "replicasRegistry", replicas_registry_file,
00221 replicas_registry_file);
00222 if (!replicas_registry_file.size()) {
00223 return cvm::error("Error: the name of the \"replicasRegistry\" file "
00224 "must be provided.\n", COLVARS_INPUT_ERROR);
00225 }
00226 
00227 get_keyval(conf, "replicaUpdateFrequency",
00228 replica_update_freq, replica_update_freq);
00229 if (replica_update_freq == 0) {
00230 return cvm::error("Error: replicaUpdateFrequency must be positive.\n",
00231 COLVARS_INPUT_ERROR);
00232 }
00233 
00234 if (expand_grids) {
00235 return cvm::error("Error: expandBoundaries is not supported when "
00236 "using more than one replicas; please allocate "
00237 "wide enough boundaries for each colvar"
00238 "ahead of time.\n", COLVARS_INPUT_ERROR);
00239 }
00240 
00241 if (keep_hills) {
00242 return cvm::error("Error: multipleReplicas and keepHills are not "
00243 "supported together.\n", COLVARS_INPUT_ERROR);
00244 }
00245 }
00246 
00247 return COLVARS_OK;
00248 }
00249 
00250 
00251 int colvarbias_meta::init_well_tempered_params(std::string const &conf)
00252 {
00253 // for well-tempered metadynamics
00254 get_keyval(conf, "wellTempered", well_tempered, false);
00255 get_keyval(conf, "biasTemperature", bias_temperature, -1.0);
00256 if ((bias_temperature == -1.0) && well_tempered) {
00257 cvm::error("Error: biasTemperature must be set to a positive value.\n",
00258 COLVARS_INPUT_ERROR);
00259 }
00260 if (well_tempered) {
00261 cvm::log("Well-tempered metadynamics is used.\n");
00262 cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
00263 }
00264 return COLVARS_OK;
00265 }
00266 
00267 
00268 int colvarbias_meta::init_ebmeta_params(std::string const &conf)
00269 {
00270 int error_code = COLVARS_OK;
00271 // for ebmeta
00272 target_dist = NULL;
00273 get_keyval(conf, "ebMeta", ebmeta, false);
00274 if(ebmeta){
00275 cvm::main()->cite_feature("Ensemble-biased metadynamics (ebMetaD)");
00276 if (use_grids && expand_grids) {
00277 error_code |= cvm::error("Error: expandBoundaries is not supported with "
00278 "ebMeta; please allocate wide enough boundaries "
00279 "for each colvar ahead of time and set "
00280 "targetDistFile accordingly.\n",
00281 COLVARS_INPUT_ERROR);
00282 }
00283 target_dist = new colvar_grid_scalar();
00284 error_code |= target_dist->init_from_colvars(colvars);
00285 std::string target_dist_file;
00286 get_keyval(conf, "targetDistFile", target_dist_file);
00287 error_code |= target_dist->read_multicol(target_dist_file,
00288 "ebMeta target histogram");
00289 cvm::real min_val = target_dist->minimum_value();
00290 cvm::real max_val = target_dist->maximum_value();
00291 if (min_val < 0.0) {
00292 error_code |= cvm::error("Error: Target distribution of EBMetaD "
00293 "has negative values!.\n",
00294 COLVARS_INPUT_ERROR);
00295 }
00296 cvm::real target_dist_min_val;
00297 get_keyval(conf, "targetDistMinVal", target_dist_min_val, 1/1000000.0);
00298 if(target_dist_min_val>0 && target_dist_min_val<1){
00299 target_dist_min_val=max_val*target_dist_min_val;
00300 target_dist->remove_small_values(target_dist_min_val);
00301 } else {
00302 if (target_dist_min_val==0) {
00303 cvm::log("NOTE: targetDistMinVal is set to zero, the minimum value of the target \n");
00304 cvm::log(" distribution will be set as the minimum positive value.\n");
00305 cvm::real min_pos_val = target_dist->minimum_pos_value();
00306 if (min_pos_val <= 0.0){
00307 error_code |= cvm::error("Error: Target distribution of EBMetaD has "
00308 "negative or zero minimum positive value.\n",
00309 COLVARS_INPUT_ERROR);
00310 }
00311 if (min_val == 0.0){
00312 cvm::log("WARNING: Target distribution has zero values.\n");
00313 cvm::log("Zeros will be converted to the minimum positive value.\n");
00314 target_dist->remove_small_values(min_pos_val);
00315 }
00316 } else {
00317 error_code |= cvm::error("Error: targetDistMinVal must be a value "
00318 "between 0 and 1.\n", COLVARS_INPUT_ERROR);
00319 }
00320 }
00321 // normalize target distribution and multiply by effective volume = exp(differential entropy)
00322 target_dist->multiply_constant(1.0/target_dist->integral());
00323 cvm::real volume = cvm::exp(target_dist->entropy());
00324 target_dist->multiply_constant(volume);
00325 get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, ebmeta_equil_steps);
00326 }
00327 
00328 return error_code;
00329 }
00330 
00331 
00332 colvarbias_meta::~colvarbias_meta()
00333 {
00334 colvarbias_meta::clear_state_data();
00335 colvarproxy *proxy = cvm::main()->proxy;
00336 
00337 proxy->close_output_stream(replica_hills_file);
00338 
00339 proxy->close_output_stream(hills_traj_file_name());
00340 
00341 if (target_dist) {
00342 delete target_dist;
00343 target_dist = NULL;
00344 }
00345 }
00346 
00347 
00348 int colvarbias_meta::clear_state_data()
00349 {
00350 if (hills_energy) {
00351 delete hills_energy;
00352 hills_energy = NULL;
00353 }
00354 
00355 if (hills_energy_gradients) {
00356 delete hills_energy_gradients;
00357 hills_energy_gradients = NULL;
00358 }
00359 
00360 hills.clear();
00361 hills_off_grid.clear();
00362 
00363 return COLVARS_OK;
00364 }
00365 
00366 
00367 // **********************************************************************
00368 // Hill management member functions
00369 // **********************************************************************
00370 
00371 std::list<colvarbias_meta::hill>::const_iterator
00372 colvarbias_meta::add_hill(colvarbias_meta::hill const &h)
00373 {
00374 hill_iter const hills_end = hills.end();
00375 hills.push_back(h);
00376 if (new_hills_begin == hills_end) {
00377 // if new_hills_begin is unset, set it for the first time
00378 new_hills_begin = hills.end();
00379 new_hills_begin--;
00380 }
00381 
00382 if (use_grids) {
00383 
00384 // also add it to the list of hills that are off-grid, which may
00385 // need to be computed analytically when the colvar returns
00386 // off-grid
00387 cvm::real const min_dist = hills_energy->bin_distance_from_boundaries(h.centers, true);
00388 if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) {
00389 hills_off_grid.push_back(h);
00390 }
00391 }
00392 
00393 // output to trajectory (if specified)
00394 if (b_hills_traj) {
00395 // Open trajectory file or access the one already open
00396 std::ostream &hills_traj_os =
00397 cvm::proxy->output_stream(hills_traj_file_name());
00398 hills_traj_os << (hills.back()).output_traj();
00399 cvm::proxy->flush_output_stream(hills_traj_file_name());
00400 }
00401 
00402 has_data = true;
00403 return hills.end();
00404 }
00405 
00406 
00407 std::list<colvarbias_meta::hill>::const_iterator
00408 colvarbias_meta::delete_hill(hill_iter &h)
00409 {
00410 if (cvm::debug()) {
00411 cvm::log("Deleting hill from the metadynamics bias \""+this->name+"\""+
00412 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00413 ", with step number "+
00414 cvm::to_str(h->it)+(h->replica.size() ?
00415 ", replica id \""+h->replica :
00416 "")+".\n");
00417 }
00418 
00419 if (use_grids && !hills_off_grid.empty()) {
00420 for (hill_iter hoff = hills_off_grid.begin();
00421 hoff != hills_off_grid.end(); hoff++) {
00422 if (*h == *hoff) {
00423 hills_off_grid.erase(hoff);
00424 break;
00425 }
00426 }
00427 }
00428 
00429 if (b_hills_traj) {
00430 // output to the trajectory
00431 std::ostream &hills_traj_os =
00432 cvm::proxy->output_stream(hills_traj_file_name());
00433 hills_traj_os << "# DELETED this hill: "
00434 << (hills.back()).output_traj()
00435 << "\n";
00436 cvm::proxy->flush_output_stream(hills_traj_file_name());
00437 }
00438 
00439 return hills.erase(h);
00440 }
00441 
00442 
00443 int colvarbias_meta::update()
00444 {
00445 int error_code = COLVARS_OK;
00446 
00447 // update base class
00448 error_code |= colvarbias::update();
00449 
00450 // update the TI estimator (if defined)
00451 error_code |= colvarbias_ti::update();
00452 
00453 // update grid definition, if needed
00454 error_code |= update_grid_params();
00455 // add new biasing energy/forces
00456 error_code |= update_bias();
00457 // update grid content to reflect new bias
00458 error_code |= update_grid_data();
00459 
00460 if (comm != single_replica &&
00461 (cvm::step_absolute() % replica_update_freq) == 0) {
00462 // sync with the other replicas (if needed)
00463 error_code |= replica_share();
00464 }
00465 
00466 error_code |= calc_energy(NULL);
00467 error_code |= calc_forces(NULL);
00468 
00469 return error_code;
00470 }
00471 
00472 
00473 int colvarbias_meta::update_grid_params()
00474 {
00475 if (use_grids) {
00476 
00477 std::vector<int> curr_bin = hills_energy->get_colvars_index();
00478 if (cvm::debug()) {
00479 cvm::log("Metadynamics bias \""+this->name+"\""+
00480 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00481 ": current coordinates on the grid: "+
00482 cvm::to_str(curr_bin)+".\n");
00483 }
00484 
00485 if (expand_grids) {
00486 // first of all, expand the grids, if specified
00487 bool changed_grids = false;
00488 int const min_buffer =
00489 (3 * (size_t) cvm::floor(hill_width)) + 1;
00490 
00491 std::vector<int> new_sizes(hills_energy->sizes());
00492 std::vector<colvarvalue> new_lower_boundaries(hills_energy->lower_boundaries);
00493 std::vector<colvarvalue> new_upper_boundaries(hills_energy->upper_boundaries);
00494 
00495 for (size_t i = 0; i < num_variables(); i++) {
00496 
00497 if (! variables(i)->expand_boundaries)
00498 continue;
00499 
00500 cvm::real &new_lb = new_lower_boundaries[i].real_value;
00501 cvm::real &new_ub = new_upper_boundaries[i].real_value;
00502 int &new_size = new_sizes[i];
00503 bool changed_lb = false, changed_ub = false;
00504 
00505 if (!variables(i)->is_enabled(f_cv_hard_lower_boundary))
00506 if (curr_bin[i] < min_buffer) {
00507 int const extra_points = (min_buffer - curr_bin[i]);
00508 new_lb -= extra_points * variables(i)->width;
00509 new_size += extra_points;
00510 // changed offset in this direction => the pointer needs to
00511 // be changed, too
00512 curr_bin[i] += extra_points;
00513 
00514 changed_lb = true;
00515 cvm::log("Metadynamics bias \""+this->name+"\""+
00516 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00517 ": new lower boundary for colvar \""+
00518 variables(i)->name+"\", at "+
00519 cvm::to_str(new_lower_boundaries[i])+".\n");
00520 }
00521 
00522 if (!variables(i)->is_enabled(f_cv_hard_upper_boundary))
00523 if (curr_bin[i] > new_size - min_buffer - 1) {
00524 int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
00525 new_ub += extra_points * variables(i)->width;
00526 new_size += extra_points;
00527 
00528 changed_ub = true;
00529 cvm::log("Metadynamics bias \""+this->name+"\""+
00530 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00531 ": new upper boundary for colvar \""+
00532 variables(i)->name+"\", at "+
00533 cvm::to_str(new_upper_boundaries[i])+".\n");
00534 }
00535 
00536 if (changed_lb || changed_ub)
00537 changed_grids = true;
00538 }
00539 
00540 if (changed_grids) {
00541 
00542 // map everything into new grids
00543 
00544 colvar_grid_scalar *new_hills_energy =
00545 new colvar_grid_scalar(*hills_energy);
00546 colvar_grid_gradient *new_hills_energy_gradients =
00547 new colvar_grid_gradient(*hills_energy_gradients);
00548 
00549 // supply new boundaries to the new grids
00550 
00551 new_hills_energy->lower_boundaries = new_lower_boundaries;
00552 new_hills_energy->upper_boundaries = new_upper_boundaries;
00553 new_hills_energy->setup(new_sizes, 0.0, 1);
00554 
00555 new_hills_energy_gradients->lower_boundaries = new_lower_boundaries;
00556 new_hills_energy_gradients->upper_boundaries = new_upper_boundaries;
00557 new_hills_energy_gradients->setup(new_sizes, 0.0, num_variables());
00558 
00559 new_hills_energy->map_grid(*hills_energy);
00560 new_hills_energy_gradients->map_grid(*hills_energy_gradients);
00561 
00562 delete hills_energy;
00563 delete hills_energy_gradients;
00564 hills_energy = new_hills_energy;
00565 hills_energy_gradients = new_hills_energy_gradients;
00566 
00567 curr_bin = hills_energy->get_colvars_index();
00568 if (cvm::debug())
00569 cvm::log("Coordinates on the new grid: "+
00570 cvm::to_str(curr_bin)+".\n");
00571 }
00572 }
00573 }
00574 return COLVARS_OK;
00575 }
00576 
00577 
00578 int colvarbias_meta::update_bias()
00579 {
00580 colvarproxy *proxy = cvm::main()->proxy;
00581 // add a new hill if the required time interval has passed
00582 if (((cvm::step_absolute() % new_hill_freq) == 0) &&
00583 can_accumulate_data() && is_enabled(f_cvb_history_dependent)) {
00584 
00585 if (cvm::debug()) {
00586 cvm::log("Metadynamics bias \""+this->name+"\""+
00587 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00588 ": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n");
00589 }
00590 
00591 cvm::real hills_scale=1.0;
00592 
00593 if (ebmeta) {
00594 hills_scale *= 1.0/target_dist->value(target_dist->get_colvars_index());
00595 if(cvm::step_absolute() <= ebmeta_equil_steps) {
00596 cvm::real const hills_lambda =
00597 (cvm::real(ebmeta_equil_steps - cvm::step_absolute())) /
00598 (cvm::real(ebmeta_equil_steps));
00599 hills_scale = hills_lambda + (1-hills_lambda)*hills_scale;
00600 }
00601 }
00602 
00603 if (well_tempered) {
00604 cvm::real hills_energy_sum_here = 0.0;
00605 if (use_grids) {
00606 std::vector<int> curr_bin = hills_energy->get_colvars_index();
00607 hills_energy_sum_here = hills_energy->value(curr_bin);
00608 } else {
00609 calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here, NULL);
00610 }
00611 hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*proxy->boltzmann()));
00612 }
00613 
00614 switch (comm) {
00615 
00616 case single_replica:
00617 
00618 add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
00619 colvar_values, colvar_sigmas));
00620 
00621 break;
00622 
00623 case multiple_replicas:
00624 add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale,
00625 colvar_values, colvar_sigmas, replica_id));
00626 std::ostream &replica_hills_os =
00627 cvm::proxy->output_stream(replica_hills_file);
00628 if (replica_hills_os) {
00629 replica_hills_os << hills.back();
00630 } else {
00631 return cvm::error("Error: in metadynamics bias \""+this->name+"\""+
00632 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00633 " while writing hills for the other replicas.\n", COLVARS_FILE_ERROR);
00634 }
00635 break;
00636 }
00637 }
00638 
00639 return COLVARS_OK;
00640 }
00641 
00642 
00643 int colvarbias_meta::update_grid_data()
00644 {
00645 if ((cvm::step_absolute() % grids_freq) == 0) {
00646 // map the most recent gaussians to the grids
00647 project_hills(new_hills_begin, hills.end(),
00648 hills_energy, hills_energy_gradients);
00649 new_hills_begin = hills.end();
00650 
00651 // TODO: we may want to condense all into one replicas array,
00652 // including "this" as the first element
00653 if (comm == multiple_replicas) {
00654 for (size_t ir = 0; ir < replicas.size(); ir++) {
00655 replicas[ir]->project_hills(replicas[ir]->new_hills_begin,
00656 replicas[ir]->hills.end(),
00657 replicas[ir]->hills_energy,
00658 replicas[ir]->hills_energy_gradients);
00659 replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
00660 }
00661 }
00662 }
00663 
00664 return COLVARS_OK;
00665 }
00666 
00667 
00668 int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
00669 {
00670 size_t ir = 0;
00671 
00672 for (ir = 0; ir < replicas.size(); ir++) {
00673 replicas[ir]->bias_energy = 0.0;
00674 }
00675 
00676 std::vector<int> const curr_bin = values ?
00677 hills_energy->get_colvars_index(*values) :
00678 hills_energy->get_colvars_index();
00679 
00680 if (hills_energy->index_ok(curr_bin)) {
00681 // index is within the grid: get the energy from there
00682 for (ir = 0; ir < replicas.size(); ir++) {
00683 
00684 bias_energy += replicas[ir]->hills_energy->value(curr_bin);
00685 if (cvm::debug()) {
00686 cvm::log("Metadynamics bias \""+this->name+"\""+
00687 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00688 ": current coordinates on the grid: "+
00689 cvm::to_str(curr_bin)+".\n");
00690 cvm::log("Grid energy = "+cvm::to_str(bias_energy)+".\n");
00691 }
00692 }
00693 } else {
00694 // off the grid: compute analytically only the hills at the grid's edges
00695 for (ir = 0; ir < replicas.size(); ir++) {
00696 calc_hills(replicas[ir]->hills_off_grid.begin(),
00697 replicas[ir]->hills_off_grid.end(),
00698 bias_energy,
00699 values);
00700 }
00701 }
00702 
00703 // now include the hills that have not been binned yet (starting
00704 // from new_hills_begin)
00705 
00706 for (ir = 0; ir < replicas.size(); ir++) {
00707 calc_hills(replicas[ir]->new_hills_begin,
00708 replicas[ir]->hills.end(),
00709 bias_energy,
00710 values);
00711 if (cvm::debug()) {
00712 cvm::log("Hills energy = "+cvm::to_str(bias_energy)+".\n");
00713 }
00714 }
00715 
00716 return COLVARS_OK;
00717 }
00718 
00719 
00720 int colvarbias_meta::calc_forces(std::vector<colvarvalue> const *values)
00721 {
00722 size_t ir = 0, ic = 0;
00723 for (ir = 0; ir < replicas.size(); ir++) {
00724 for (ic = 0; ic < num_variables(); ic++) {
00725 replicas[ir]->colvar_forces[ic].reset();
00726 }
00727 }
00728 
00729 std::vector<int> const curr_bin = values ?
00730 hills_energy->get_colvars_index(*values) :
00731 hills_energy->get_colvars_index();
00732 
00733 if (hills_energy->index_ok(curr_bin)) {
00734 for (ir = 0; ir < replicas.size(); ir++) {
00735 cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value(curr_bin));
00736 for (ic = 0; ic < num_variables(); ic++) {
00737 // the gradients are stored, not the forces
00738 colvar_forces[ic].real_value += -1.0 * f[ic];
00739 }
00740 }
00741 } else {
00742 // off the grid: compute analytically only the hills at the grid's edges
00743 for (ir = 0; ir < replicas.size(); ir++) {
00744 for (ic = 0; ic < num_variables(); ic++) {
00745 calc_hills_force(ic,
00746 replicas[ir]->hills_off_grid.begin(),
00747 replicas[ir]->hills_off_grid.end(),
00748 colvar_forces,
00749 values);
00750 }
00751 }
00752 }
00753 
00754 // now include the hills that have not been binned yet (starting
00755 // from new_hills_begin)
00756 
00757 if (cvm::debug()) {
00758 cvm::log("Metadynamics bias \""+this->name+"\""+
00759 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00760 ": adding the forces from the other replicas.\n");
00761 }
00762 
00763 for (ir = 0; ir < replicas.size(); ir++) {
00764 for (ic = 0; ic < num_variables(); ic++) {
00765 calc_hills_force(ic,
00766 replicas[ir]->new_hills_begin,
00767 replicas[ir]->hills.end(),
00768 colvar_forces,
00769 values);
00770 if (cvm::debug()) {
00771 cvm::log("Hills forces = "+cvm::to_str(colvar_forces)+".\n");
00772 }
00773 }
00774 }
00775 
00776 return COLVARS_OK;
00777 }
00778 
00779 
00780 
00781 void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first,
00782 colvarbias_meta::hill_iter h_last,
00783 cvm::real &energy,
00784 std::vector<colvarvalue> const *values)
00785 {
00786 size_t i = 0;
00787 
00788 for (hill_iter h = h_first; h != h_last; h++) {
00789 
00790 // compute the gaussian exponent
00791 cvm::real cv_sqdev = 0.0;
00792 for (i = 0; i < num_variables(); i++) {
00793 colvarvalue const &x = values ? (*values)[i] : colvar_values[i];
00794 colvarvalue const &center = h->centers[i];
00795 cvm::real const sigma = h->sigmas[i];
00796 cv_sqdev += (variables(i)->dist2(x, center)) / (sigma*sigma);
00797 }
00798 
00799 // compute the gaussian
00800 if (cv_sqdev > 23.0) {
00801 // set it to zero if the exponent is more negative than log(1.0E-06)
00802 h->value(0.0);
00803 } else {
00804 h->value(cvm::exp(-0.5*cv_sqdev));
00805 }
00806 energy += h->energy();
00807 }
00808 }
00809 
00810 
00811 void colvarbias_meta::calc_hills_force(size_t const &i,
00812 colvarbias_meta::hill_iter h_first,
00813 colvarbias_meta::hill_iter h_last,
00814 std::vector<colvarvalue> &forces,
00815 std::vector<colvarvalue> const *values)
00816 {
00817 // Retrieve the value of the colvar
00818 colvarvalue const x(values ? (*values)[i] : colvar_values[i]);
00819 
00820 // do the type check only once (all colvarvalues in the hills series
00821 // were already saved with their types matching those in the
00822 // colvars)
00823 
00824 hill_iter h;
00825 switch (x.type()) {
00826 
00827 case colvarvalue::type_scalar:
00828 for (h = h_first; h != h_last; h++) {
00829 if (h->value() == 0.0) continue;
00830 colvarvalue const &center = h->centers[i];
00831 cvm::real const sigma = h->sigmas[i];
00832 forces[i].real_value +=
00833 ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
00834 (variables(i)->dist2_lgrad(x, center)).real_value );
00835 }
00836 break;
00837 
00838 case colvarvalue::type_3vector:
00839 case colvarvalue::type_unit3vector:
00840 case colvarvalue::type_unit3vectorderiv:
00841 for (h = h_first; h != h_last; h++) {
00842 if (h->value() == 0.0) continue;
00843 colvarvalue const &center = h->centers[i];
00844 cvm::real const sigma = h->sigmas[i];
00845 forces[i].rvector_value +=
00846 ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
00847 (variables(i)->dist2_lgrad(x, center)).rvector_value );
00848 }
00849 break;
00850 
00851 case colvarvalue::type_quaternion:
00852 case colvarvalue::type_quaternionderiv:
00853 for (h = h_first; h != h_last; h++) {
00854 if (h->value() == 0.0) continue;
00855 colvarvalue const &center = h->centers[i];
00856 cvm::real const sigma = h->sigmas[i];
00857 forces[i].quaternion_value +=
00858 ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
00859 (variables(i)->dist2_lgrad(x, center)).quaternion_value );
00860 }
00861 break;
00862 
00863 case colvarvalue::type_vector:
00864 for (h = h_first; h != h_last; h++) {
00865 if (h->value() == 0.0) continue;
00866 colvarvalue const &center = h->centers[i];
00867 cvm::real const sigma = h->sigmas[i];
00868 forces[i].vector1d_value +=
00869 ( h->weight() * h->value() * (0.5 / (sigma*sigma)) *
00870 (variables(i)->dist2_lgrad(x, center)).vector1d_value );
00871 }
00872 break;
00873 
00874 case colvarvalue::type_notset:
00875 case colvarvalue::type_all:
00876 default:
00877 break;
00878 }
00879 }
00880 
00881 
00882 // **********************************************************************
00883 // grid management functions
00884 // **********************************************************************
00885 
00886 void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
00887 colvarbias_meta::hill_iter h_last,
00888 colvar_grid_scalar *he,
00889 colvar_grid_gradient *hg,
00890 bool print_progress)
00891 {
00892 if (cvm::debug())
00893 cvm::log("Metadynamics bias \""+this->name+"\""+
00894 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00895 ": projecting hills.\n");
00896 
00897 // TODO: improve it by looping over a small subgrid instead of the whole grid
00898 
00899 std::vector<colvarvalue> new_colvar_values(num_variables());
00900 std::vector<cvm::real> colvar_forces_scalar(num_variables());
00901 
00902 std::vector<int> he_ix = he->new_index();
00903 std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
00904 cvm::real hills_energy_here = 0.0;
00905 std::vector<colvarvalue> hills_forces_here(num_variables(), 0.0);
00906 
00907 size_t count = 0;
00908 size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1)));
00909 
00910 if (hg != NULL) {
00911 
00912 // loop over the points of the grid
00913 for ( ;
00914 (he->index_ok(he_ix)) && (hg->index_ok(hg_ix));
00915 count++) {
00916 size_t i;
00917 for (i = 0; i < num_variables(); i++) {
00918 new_colvar_values[i] = he->bin_to_value_scalar(he_ix[i], i);
00919 }
00920 
00921 // loop over the hills and increment the energy grid locally
00922 hills_energy_here = 0.0;
00923 calc_hills(h_first, h_last, hills_energy_here, &new_colvar_values);
00924 he->acc_value(he_ix, hills_energy_here);
00925 
00926 for (i = 0; i < num_variables(); i++) {
00927 hills_forces_here[i].reset();
00928 calc_hills_force(i, h_first, h_last, hills_forces_here, &new_colvar_values);
00929 colvar_forces_scalar[i] = hills_forces_here[i].real_value;
00930 }
00931 hg->acc_force(hg_ix, &(colvar_forces_scalar.front()));
00932 
00933 he->incr(he_ix);
00934 hg->incr(hg_ix);
00935 
00936 if ((count % print_frequency) == 0) {
00937 if (print_progress) {
00938 cvm::real const progress = cvm::real(count) / cvm::real(hg->number_of_points());
00939 std::ostringstream os;
00940 os.setf(std::ios::fixed, std::ios::floatfield);
00941 os << std::setw(6) << std::setprecision(2)
00942 << 100.0 * progress
00943 << "% done.";
00944 cvm::log(os.str());
00945 }
00946 }
00947 }
00948 
00949 } else {
00950 cvm::error("No grid object provided in metadynamics::project_hills()\n",
00951 COLVARS_BUG_ERROR);
00952 }
00953 
00954 if (print_progress) {
00955 cvm::log("100.00% done.\n");
00956 }
00957 
00958 if (! keep_hills) {
00959 hills.erase(hills.begin(), hills.end());
00960 }
00961 }
00962 
00963 
00964 void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first,
00965 colvarbias_meta::hill_iter h_last,
00966 colvar_grid_scalar * /* he */)
00967 {
00968 hills_off_grid.clear();
00969 
00970 for (hill_iter h = h_first; h != h_last; h++) {
00971 cvm::real const min_dist = hills_energy->bin_distance_from_boundaries(h->centers, true);
00972 if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) {
00973 hills_off_grid.push_back(*h);
00974 }
00975 }
00976 }
00977 
00978 
00979 
00980 // **********************************************************************
00981 // multiple replicas functions
00982 // **********************************************************************
00983 
00984 
00985 int colvarbias_meta::replica_share()
00986 {
00987 int error_code = COLVARS_OK;
00988 colvarproxy *proxy = cvm::proxy;
00989 // sync with the other replicas (if needed)
00990 if (comm == multiple_replicas) {
00991 // reread the replicas registry
00992 error_code |= update_replicas_registry();
00993 // empty the output buffer
00994 error_code |= proxy->flush_output_stream(replica_hills_file);
00995 error_code |= read_replica_files();
00996 }
00997 return error_code;
00998 }
00999 
01000 
01001 int colvarbias_meta::update_replicas_registry()
01002 {
01003 int error_code = COLVARS_OK;
01004 
01005 if (cvm::debug())
01006 cvm::log("Metadynamics bias \""+this->name+"\""+
01007 ": updating the list of replicas, currently containing "+
01008 cvm::to_str(replicas.size())+" elements.\n");
01009 
01010 {
01011 // copy the whole file into a string for convenience
01012 std::string line("");
01013 std::ifstream reg_file(replicas_registry_file.c_str());
01014 if (reg_file.is_open()) {
01015 replicas_registry.clear();
01016 while (colvarparse::getline_nocomments(reg_file, line))
01017 replicas_registry.append(line+"\n");
01018 } else {
01019 error_code |= cvm::error("Error: failed to open file \""+
01020 replicas_registry_file+"\" for reading.\n",
01021 COLVARS_FILE_ERROR);
01022 }
01023 }
01024 
01025 // now parse it
01026 std::istringstream reg_is(replicas_registry);
01027 if (reg_is.good()) {
01028 
01029 std::string new_replica("");
01030 std::string new_replica_file("");
01031 while ((reg_is >> new_replica) && new_replica.size() &&
01032 (reg_is >> new_replica_file) && new_replica_file.size()) {
01033 
01034 if (new_replica == this->replica_id) {
01035 // this is the record for this same replica, skip it
01036 new_replica_file.clear();
01037 new_replica.clear();
01038 continue;
01039 }
01040 
01041 bool already_loaded = false;
01042 for (size_t ir = 0; ir < replicas.size(); ir++) {
01043 if (new_replica == (replicas[ir])->replica_id) {
01044 // this replica was already added
01045 if (cvm::debug())
01046 cvm::log("Metadynamics bias \""+this->name+"\""+
01047 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01048 ": skipping a replica already loaded, \""+
01049 (replicas[ir])->replica_id+"\".\n");
01050 already_loaded = true;
01051 break;
01052 }
01053 }
01054 
01055 if (!already_loaded) {
01056 // add this replica to the registry
01057 cvm::log("Metadynamics bias \""+this->name+"\""+
01058 ": accessing replica \""+new_replica+"\".\n");
01059 replicas.push_back(new colvarbias_meta("metadynamics"));
01060 (replicas.back())->replica_id = new_replica;
01061 (replicas.back())->replica_list_file = new_replica_file;
01062 (replicas.back())->replica_state_file = "";
01063 (replicas.back())->replica_state_file_in_sync = false;
01064 
01065 // Note: the following could become a copy constructor?
01066 (replicas.back())->name = this->name;
01067 (replicas.back())->colvars = colvars;
01068 (replicas.back())->use_grids = use_grids;
01069 (replicas.back())->dump_fes = false;
01070 (replicas.back())->expand_grids = false;
01071 (replicas.back())->rebin_grids = false;
01072 (replicas.back())->keep_hills = false;
01073 (replicas.back())->colvar_forces = colvar_forces;
01074 
01075 (replicas.back())->comm = multiple_replicas;
01076 
01077 if (use_grids) {
01078 (replicas.back())->hills_energy = new colvar_grid_scalar(colvars);
01079 (replicas.back())->hills_energy_gradients = new colvar_grid_gradient(colvars);
01080 }
01081 if (is_enabled(f_cvb_calc_ti_samples)) {
01082 (replicas.back())->enable(f_cvb_calc_ti_samples);
01083 (replicas.back())->colvarbias_ti::init_grids();
01084 }
01085 (replicas.back())->update_status = 1;
01086 }
01087 }
01088 } else {
01089 error_code |= cvm::error("Error: cannot read the replicas registry file \""+
01090 replicas_registry+"\".\n", COLVARS_FILE_ERROR);
01091 }
01092 
01093 // now (re)read the list file of each replica
01094 for (size_t ir = 0; ir < replicas.size(); ir++) {
01095 if (cvm::debug())
01096 cvm::log("Metadynamics bias \""+this->name+"\""+
01097 ": reading the list file for replica \""+
01098 (replicas[ir])->replica_id+"\".\n");
01099 
01100 std::ifstream list_is((replicas[ir])->replica_list_file.c_str());
01101 std::string key;
01102 std::string new_state_file, new_hills_file;
01103 if (!(list_is >> key) ||
01104 !(key == std::string("stateFile")) ||
01105 !(list_is >> new_state_file) ||
01106 !(list_is >> key) ||
01107 !(key == std::string("hillsFile")) ||
01108 !(list_is >> new_hills_file)) {
01109 cvm::log("Metadynamics bias \""+this->name+"\""+
01110 ": failed to read the file \""+
01111 (replicas[ir])->replica_list_file+"\": will try again after "+
01112 cvm::to_str(replica_update_freq)+" steps.\n");
01113 (replicas[ir])->update_status++;
01114 } else {
01115 if (new_state_file != (replicas[ir])->replica_state_file) {
01116 cvm::log("Metadynamics bias \""+this->name+"\""+
01117 ": replica \""+(replicas[ir])->replica_id+
01118 "\" has supplied a new state file, \""+new_state_file+
01119 "\".\n");
01120 (replicas[ir])->replica_state_file_in_sync = false;
01121 (replicas[ir])->replica_state_file = new_state_file;
01122 (replicas[ir])->replica_hills_file = new_hills_file;
01123 }
01124 }
01125 }
01126 
01127 if (cvm::debug())
01128 cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+
01129 cvm::to_str(replicas.size())+" elements.\n");
01130 
01131 return error_code;
01132 }
01133 
01134 
01135 int colvarbias_meta::read_replica_files()
01136 {
01137 // Note: we start from the 2nd replica.
01138 for (size_t ir = 1; ir < replicas.size(); ir++) {
01139 
01140 // (re)read the state file if necessary
01141 if ( (! (replicas[ir])->has_data) ||
01142 (! (replicas[ir])->replica_state_file_in_sync) ) {
01143 if ((replicas[ir])->replica_state_file.size()) {
01144 cvm::log("Metadynamics bias \""+this->name+"\""+
01145 ": reading the state of replica \""+
01146 (replicas[ir])->replica_id+"\" from file \""+
01147 (replicas[ir])->replica_state_file+"\".\n");
01148 std::ifstream is((replicas[ir])->replica_state_file.c_str());
01149 if ((replicas[ir])->read_state(is)) {
01150 // state file has been read successfully
01151 (replicas[ir])->replica_state_file_in_sync = true;
01152 (replicas[ir])->update_status = 0;
01153 } else {
01154 cvm::log("Failed to read the file \""+
01155 (replicas[ir])->replica_state_file+
01156 "\": will try again in "+
01157 cvm::to_str(replica_update_freq)+" steps.\n");
01158 (replicas[ir])->replica_state_file_in_sync = false;
01159 (replicas[ir])->update_status++;
01160 }
01161 is.close();
01162 } else {
01163 cvm::log("Metadynamics bias \""+this->name+"\""+
01164 ": the state file of replica \""+
01165 (replicas[ir])->replica_id+"\" is currently undefined: "
01166 "will try again after "+
01167 cvm::to_str(replica_update_freq)+" steps.\n");
01168 (replicas[ir])->update_status++;
01169 }
01170 }
01171 
01172 if (! (replicas[ir])->replica_state_file_in_sync) {
01173 // if a new state file is being read, the hills file is also new
01174 (replicas[ir])->replica_hills_file_pos = 0;
01175 }
01176 
01177 // now read the hills added after writing the state file
01178 if ((replicas[ir])->replica_hills_file.size()) {
01179 
01180 if (cvm::debug())
01181 cvm::log("Metadynamics bias \""+this->name+"\""+
01182 ": checking for new hills from replica \""+
01183 (replicas[ir])->replica_id+"\" in the file \""+
01184 (replicas[ir])->replica_hills_file+"\".\n");
01185 
01186 // read hills from the other replicas' files
01187 
01188 std::ifstream is((replicas[ir])->replica_hills_file.c_str());
01189 if (is.is_open()) {
01190 
01191 // try to resume the previous position (if not the beginning)
01192 if ((replicas[ir])->replica_hills_file_pos > 0) {
01193 is.seekg((replicas[ir])->replica_hills_file_pos, std::ios::beg);
01194 }
01195 
01196 if (!is.is_open()){
01197 // if fail (the file may have been overwritten), reset this
01198 // position
01199 is.clear();
01200 is.seekg(0, std::ios::beg);
01201 // reset the counter
01202 (replicas[ir])->replica_hills_file_pos = 0;
01203 // schedule to reread the state file
01204 (replicas[ir])->replica_state_file_in_sync = false;
01205 // and record the failure
01206 (replicas[ir])->update_status++;
01207 cvm::log("Failed to read the file \""+(replicas[ir])->replica_hills_file+
01208 "\" at the previous position: will try again in "+
01209 cvm::to_str(replica_update_freq)+" steps.\n");
01210 } else {
01211 
01212 while ((replicas[ir])->read_hill(is)) {
01213 cvm::log("Metadynamics bias \""+this->name+"\""+
01214 ": received a hill from replica \""+
01215 (replicas[ir])->replica_id+
01216 "\" at step "+
01217 cvm::to_str(((replicas[ir])->hills.back()).it)+".\n");
01218 }
01219 is.clear();
01220 // store the position for the next read
01221 (replicas[ir])->replica_hills_file_pos = is.tellg();
01222 if (cvm::debug()) {
01223 cvm::log("Metadynamics bias \""+this->name+"\""+
01224 ": stopped reading file \""+
01225 (replicas[ir])->replica_hills_file+
01226 "\" at position "+
01227 cvm::to_str((replicas[ir])->replica_hills_file_pos)+".\n");
01228 }
01229 
01230 // test whether this is the end of the file
01231 is.seekg(0, std::ios::end);
01232 if (is.tellg() > (replicas[ir])->replica_hills_file_pos + ((std::streampos) 1)) {
01233 (replicas[ir])->update_status++;
01234 } else {
01235 (replicas[ir])->update_status = 0;
01236 }
01237 }
01238 
01239 } else {
01240 cvm::log("Failed to read the file \""+
01241 (replicas[ir])->replica_hills_file+
01242 "\": will try again in "+
01243 cvm::to_str(replica_update_freq)+" steps.\n");
01244 (replicas[ir])->update_status++;
01245 }
01246 is.close();
01247 }
01248 
01249 size_t const n_flush = (replica_update_freq/new_hill_freq + 1);
01250 if ((replicas[ir])->update_status > 3*n_flush) {
01251 // TODO: suspend the calculation?
01252 cvm::log("WARNING: metadynamics bias \""+this->name+"\""+
01253 " could not read information from replica \""+
01254 (replicas[ir])->replica_id+
01255 "\" after more than "+
01256 cvm::to_str((replicas[ir])->update_status * replica_update_freq)+
01257 " steps. Ensure that it is still running.\n");
01258 }
01259 }
01260 return COLVARS_OK;
01261 }
01262 
01263 
01264 int colvarbias_meta::set_state_params(std::string const &state_conf)
01265 {
01266 int error_code = colvarbias::set_state_params(state_conf);
01267 
01268 if (error_code != COLVARS_OK) {
01269 return error_code;
01270 }
01271 
01272 colvarparse::get_keyval(state_conf, "keepHills", restart_keep_hills, false,
01273 colvarparse::parse_restart);
01274 
01275 if ((!restart_keep_hills) && (cvm::main()->restart_version_number() < 20210604)) {
01276 if (keep_hills) {
01277 cvm::log("Warning: could not ensure that keepHills was enabled when "
01278 "this state file was written; because it is enabled now, "
01279 "it is assumed that it was also then, but please verify.\n");
01280 restart_keep_hills = true;
01281 }
01282 } else {
01283 if (restart_keep_hills) {
01284 cvm::log("This state file/stream contains explicit hills.\n");
01285 }
01286 }
01287 
01288 std::string check_replica = "";
01289 if (colvarparse::get_keyval(state_conf, "replicaID", check_replica,
01290 std::string(""), colvarparse::parse_restart) &&
01291 (check_replica != this->replica_id)) {
01292 return cvm::error("Error: in the state file , the "
01293 "\"metadynamics\" block has a different replicaID ("+
01294 check_replica+" instead of "+replica_id+").\n",
01295 COLVARS_INPUT_ERROR);
01296 }
01297 
01298 return COLVARS_OK;
01299 }
01300 
01301 
01302 std::istream & colvarbias_meta::read_state_data(std::istream& is)
01303 {
01304 if (use_grids) {
01305 
01306 if (expand_grids) {
01307 // the boundaries of the colvars may have been changed; TODO:
01308 // this reallocation is only for backward-compatibility, and may
01309 // be deleted when grid_parameters (i.e. colvargrid's own
01310 // internal reallocation) has kicked in
01311 delete hills_energy;
01312 delete hills_energy_gradients;
01313 hills_energy = new colvar_grid_scalar(colvars);
01314 hills_energy_gradients = new colvar_grid_gradient(colvars);
01315 }
01316 
01317 colvar_grid_scalar *hills_energy_backup = NULL;
01318 colvar_grid_gradient *hills_energy_gradients_backup = NULL;
01319 
01320 if (has_data) {
01321 if (cvm::debug())
01322 cvm::log("Backupping grids for metadynamics bias \""+
01323 this->name+"\""+
01324 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
01325 hills_energy_backup = hills_energy;
01326 hills_energy_gradients_backup = hills_energy_gradients;
01327 hills_energy = new colvar_grid_scalar(colvars);
01328 hills_energy_gradients = new colvar_grid_gradient(colvars);
01329 }
01330 
01331 std::streampos const hills_energy_pos = is.tellg();
01332 std::string key;
01333 if (!(is >> key)) {
01334 if (hills_energy_backup != NULL) {
01335 delete hills_energy;
01336 delete hills_energy_gradients;
01337 hills_energy = hills_energy_backup;
01338 hills_energy_gradients = hills_energy_gradients_backup;
01339 }
01340 is.clear();
01341 is.seekg(hills_energy_pos, std::ios::beg);
01342 is.setstate(std::ios::failbit);
01343 return is;
01344 } else if (!(key == std::string("hills_energy")) ||
01345 !(hills_energy->read_restart(is))) {
01346 is.clear();
01347 is.seekg(hills_energy_pos, std::ios::beg);
01348 if (!rebin_grids) {
01349 if ((hills_energy_backup == NULL) || (comm == single_replica)) {
01350 cvm::error("Error: couldn't read the energy grid for metadynamics bias \""+
01351 this->name+"\""+
01352 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01353 "; if useGrids was off when the state file was written, "
01354 "enable rebinGrids now to regenerate the grids.\n");
01355 } else {
01356 delete hills_energy;
01357 delete hills_energy_gradients;
01358 hills_energy = hills_energy_backup;
01359 hills_energy_gradients = hills_energy_gradients_backup;
01360 is.setstate(std::ios::failbit);
01361 return is;
01362 }
01363 }
01364 }
01365 
01366 std::streampos const hills_energy_gradients_pos = is.tellg();
01367 if (!(is >> key)) {
01368 if (hills_energy_backup != NULL) {
01369 delete hills_energy;
01370 delete hills_energy_gradients;
01371 hills_energy = hills_energy_backup;
01372 hills_energy_gradients = hills_energy_gradients_backup;
01373 }
01374 is.clear();
01375 is.seekg(hills_energy_gradients_pos, std::ios::beg);
01376 is.setstate(std::ios::failbit);
01377 return is;
01378 } else if (!(key == std::string("hills_energy_gradients")) ||
01379 !(hills_energy_gradients->read_restart(is))) {
01380 is.clear();
01381 is.seekg(hills_energy_gradients_pos, std::ios::beg);
01382 if (!rebin_grids) {
01383 if ((hills_energy_backup == NULL) || (comm == single_replica)) {
01384 cvm::error("Error: couldn't read the gradients grid for metadynamics bias \""+
01385 this->name+"\""+
01386 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01387 "; if useGrids was off when the state file was written, "
01388 "enable rebinGrids now to regenerate the grids.\n");
01389 } else {
01390 delete hills_energy;
01391 delete hills_energy_gradients;
01392 hills_energy = hills_energy_backup;
01393 hills_energy_gradients = hills_energy_gradients_backup;
01394 is.setstate(std::ios::failbit);
01395 return is;
01396 }
01397 }
01398 }
01399 
01400 if (cvm::debug())
01401 cvm::log("Successfully read new grids for bias \""+
01402 this->name+"\""+
01403 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
01404 
01405 cvm::log(" read biasing energy and forces from grids.\n");
01406 
01407 if (hills_energy_backup != NULL) {
01408 // now that we have successfully updated the grids, delete the
01409 // backup copies
01410 if (cvm::debug())
01411 cvm::log("Deallocating the older grids.\n");
01412 
01413 delete hills_energy_backup;
01414 delete hills_energy_gradients_backup;
01415 }
01416 }
01417 
01418 // Save references to the end of the list of existing hills, so that it can
01419 // be cleared if hills are read successfully state
01420 bool const existing_hills = !hills.empty();
01421 size_t const old_hills_size = hills.size();
01422 hill_iter old_hills_end = hills.end();
01423 hill_iter old_hills_off_grid_end = hills_off_grid.end();
01424 if (cvm::debug()) {
01425 cvm::log("Before reading hills from the state file, there are "+
01426 cvm::to_str(hills.size())+" hills in memory.\n");
01427 }
01428 
01429 // Read any hills following the grid data (if any)
01430 while (read_hill(is)) {
01431 if (cvm::debug()) {
01432 cvm::log("Read a previously saved hill under the "
01433 "metadynamics bias \""+
01434 this->name+"\", created at step "+
01435 cvm::to_str((hills.back()).it)+".\n");
01436 }
01437 }
01438 is.clear();
01439 new_hills_begin = hills.end();
01440 cvm::log(" read "+cvm::to_str(hills.size() - old_hills_size)+
01441 " additional explicit hills.\n");
01442 
01443 if (existing_hills) {
01444 hills.erase(hills.begin(), old_hills_end);
01445 hills_off_grid.erase(hills_off_grid.begin(), old_hills_off_grid_end);
01446 if (cvm::debug()) {
01447 cvm::log("After pruning the old hills, there are now "+
01448 cvm::to_str(hills.size())+" hills in memory.\n");
01449 }
01450 }
01451 
01452 if (rebin_grids) {
01453 
01454 // allocate new grids (based on the new boundaries and widths just
01455 // read from the configuration file), and project onto them the
01456 // grids just read from the restart file
01457 
01458 colvar_grid_scalar *new_hills_energy =
01459 new colvar_grid_scalar(colvars);
01460 colvar_grid_gradient *new_hills_energy_gradients =
01461 new colvar_grid_gradient(colvars);
01462 
01463 if (cvm::debug()) {
01464 std::ostringstream tmp_os;
01465 tmp_os << "hills_energy parameters:\n";
01466 hills_energy->write_params(tmp_os);
01467 tmp_os << "new_hills_energy parameters:\n";
01468 new_hills_energy->write_params(tmp_os);
01469 cvm::log(tmp_os.str());
01470 }
01471 
01472 if (restart_keep_hills && !hills.empty()) {
01473 // if there are hills, recompute the new grids from them
01474 cvm::log("Rebinning the energy and forces grids from "+
01475 cvm::to_str(hills.size())+" hills (this may take a while)...\n");
01476 project_hills(hills.begin(), hills.end(),
01477 new_hills_energy, new_hills_energy_gradients, true);
01478 cvm::log("rebinning done.\n");
01479 
01480 } else {
01481 // otherwise, use the grids in the restart file
01482 cvm::log("Rebinning the energy and forces grids "
01483 "from the grids in the restart file.\n");
01484 new_hills_energy->map_grid(*hills_energy);
01485 new_hills_energy_gradients->map_grid(*hills_energy_gradients);
01486 }
01487 
01488 delete hills_energy;
01489 delete hills_energy_gradients;
01490 hills_energy = new_hills_energy;
01491 hills_energy_gradients = new_hills_energy_gradients;
01492 
01493 // assuming that some boundaries have expanded, eliminate those
01494 // off-grid hills that aren't necessary any more
01495 if (!hills.empty())
01496 recount_hills_off_grid(hills.begin(), hills.end(), hills_energy);
01497 }
01498 
01499 if (use_grids) {
01500 if (!hills_off_grid.empty()) {
01501 cvm::log(cvm::to_str(hills_off_grid.size())+" hills are near the "
01502 "grid boundaries: they will be computed analytically "
01503 "and saved to the state files.\n");
01504 }
01505 }
01506 
01507 colvarbias_ti::read_state_data(is);
01508 
01509 if (cvm::debug())
01510 cvm::log("colvarbias_meta::read_restart() done\n");
01511 
01512 has_data = true;
01513 
01514 if (comm != single_replica) {
01515 read_replica_files();
01516 }
01517 
01518 return is;
01519 }
01520 
01521 
01522 inline std::istream & reset_istream(std::istream &is, size_t start_pos)
01523 {
01524 is.clear();
01525 is.seekg(start_pos, std::ios::beg);
01526 is.setstate(std::ios::failbit);
01527 return is;
01528 }
01529 
01530 
01531 std::istream & colvarbias_meta::read_hill(std::istream &is)
01532 {
01533 if (!is) return is; // do nothing if failbit is set
01534 
01535 std::streampos const start_pos = is.tellg();
01536 size_t i = 0;
01537 
01538 std::string data;
01539 if ( !(is >> read_block("hill", &data)) ) {
01540 return reset_istream(is, start_pos);
01541 }
01542 
01543 std::istringstream data_is(data);
01544 
01545 cvm::step_number h_it = 0L;
01546 cvm::real h_weight;
01547 std::vector<colvarvalue> h_centers(num_variables());
01548 for (i = 0; i < num_variables(); i++) {
01549 h_centers[i].type(variables(i)->value());
01550 }
01551 std::vector<cvm::real> h_sigmas(num_variables());
01552 std::string h_replica;
01553 
01554 std::string keyword;
01555 while (data_is >> keyword) {
01556 
01557 if (keyword == "step") {
01558 if ( !(data_is >> h_it)) {
01559 return reset_istream(is, start_pos);
01560 }
01561 if ((h_it <= state_file_step) && !restart_keep_hills) {
01562 if (cvm::debug())
01563 cvm::log("Skipping a hill older than the state file for metadynamics bias \""+
01564 this->name+"\""+
01565 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
01566 return is;
01567 }
01568 }
01569 
01570 if (keyword == "weight") {
01571 if ( !(data_is >> h_weight)) {
01572 return reset_istream(is, start_pos);
01573 }
01574 }
01575 
01576 if (keyword == "centers") {
01577 for (i = 0; i < num_variables(); i++) {
01578 if ( !(data_is >> h_centers[i])) {
01579 return reset_istream(is, start_pos);
01580 }
01581 }
01582 }
01583 
01584 if (keyword == "widths") {
01585 for (i = 0; i < num_variables(); i++) {
01586 if ( !(data_is >> h_sigmas[i])) {
01587 return reset_istream(is, start_pos);
01588 }
01589 // For backward compatibility, read the widths instead of the sigmas
01590 h_sigmas[i] /= 2.0;
01591 }
01592 }
01593 
01594 if (comm != single_replica) {
01595 if (keyword == "replicaID") {
01596 if ( !(data_is >> h_replica)) {
01597 return reset_istream(is, start_pos);
01598 }
01599 if (h_replica != replica_id) {
01600 cvm::error("Error: trying to read a hill created by replica \""+
01601 h_replica+"\" for replica \""+replica_id+
01602 "\"; did you swap output files?\n", COLVARS_INPUT_ERROR);
01603 }
01604 }
01605 }
01606 }
01607 
01608 hill_iter const hills_end = hills.end();
01609 hills.push_back(hill(h_it, h_weight, h_centers, h_sigmas, h_replica));
01610 if (new_hills_begin == hills_end) {
01611 // if new_hills_begin is unset, set it for the first time
01612 new_hills_begin = hills.end();
01613 new_hills_begin--;
01614 }
01615 
01616 if (use_grids) {
01617 // add this also to the list of hills that are off-grid, which will
01618 // be computed analytically
01619 cvm::real const min_dist =
01620 hills_energy->bin_distance_from_boundaries((hills.back()).centers, true);
01621 if (min_dist < (3.0 * cvm::floor(hill_width)) + 1.0) {
01622 hills_off_grid.push_back(hills.back());
01623 }
01624 }
01625 
01626 has_data = true;
01627 return is;
01628 }
01629 
01630 
01631 int colvarbias_meta::setup_output()
01632 {
01633 int error_code = COLVARS_OK;
01634 
01635 output_prefix = cvm::output_prefix();
01636 if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) > 1) {
01637 // if this is not the only free energy integrator, append
01638 // this bias's name, to distinguish it from the output of the other
01639 // biases producing a .pmf file
01640 output_prefix += ("."+this->name);
01641 }
01642 
01643 if (comm == multiple_replicas) {
01644 
01645 // TODO: one may want to specify the path manually for intricated filesystems?
01646 char *pwd = new char[3001];
01647 if (GETCWD(pwd, 3000) == NULL) {
01648 return cvm::error("Error: cannot get the path of the current working directory.\n",
01649 COLVARS_BUG_ERROR);
01650 }
01651 
01652 replica_list_file =
01653 (std::string(pwd)+std::string(PATHSEP)+
01654 this->name+"."+replica_id+".files.txt");
01655 // replica_hills_file and replica_state_file are those written
01656 // by the current replica; within the mirror biases, they are
01657 // those by another replica
01658 replica_hills_file =
01659 (std::string(pwd)+std::string(PATHSEP)+
01660 cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".hills");
01661 replica_state_file =
01662 (std::string(pwd)+std::string(PATHSEP)+
01663 cvm::output_prefix()+".colvars."+this->name+"."+replica_id+".state");
01664 delete[] pwd;
01665 
01666 // now register this replica
01667 
01668 // first check that it isn't already there
01669 bool registered_replica = false;
01670 std::ifstream reg_is(replicas_registry_file.c_str());
01671 if (reg_is.is_open()) { // the file may not be there yet
01672 std::string existing_replica("");
01673 std::string existing_replica_file("");
01674 while ((reg_is >> existing_replica) && existing_replica.size() &&
01675 (reg_is >> existing_replica_file) && existing_replica_file.size()) {
01676 if (existing_replica == replica_id) {
01677 // this replica was already registered
01678 replica_list_file = existing_replica_file;
01679 reg_is.close();
01680 registered_replica = true;
01681 break;
01682 }
01683 }
01684 reg_is.close();
01685 }
01686 
01687 // if this replica was not included yet, we should generate a
01688 // new record for it: but first, we write this replica's files,
01689 // for the others to read
01690 
01691 // open the "hills" buffer file
01692 reopen_replica_buffer_file();
01693 
01694 // write the state file (so that there is always one available)
01695 write_replica_state_file();
01696 
01697 // schedule to read the state files of the other replicas
01698 for (size_t ir = 0; ir < replicas.size(); ir++) {
01699 (replicas[ir])->replica_state_file_in_sync = false;
01700 }
01701 
01702 // if we're running without grids, use a growing list of "hills" files
01703 // otherwise, just one state file and one "hills" file as buffer
01704 std::ostream &list_os = cvm::proxy->output_stream(replica_list_file);
01705 if (list_os) {
01706 list_os << "stateFile " << replica_state_file << "\n";
01707 list_os << "hillsFile " << replica_hills_file << "\n";
01708 cvm::proxy->close_output_stream(replica_list_file);
01709 } else {
01710 error_code |= COLVARS_FILE_ERROR;
01711 }
01712 
01713 // finally, add a new record for this replica to the registry
01714 if (! registered_replica) {
01715 std::ofstream reg_os(replicas_registry_file.c_str(), std::ios::app);
01716 if (!reg_os) {
01717 return cvm::get_error();
01718 }
01719 reg_os << replica_id << " " << replica_list_file << "\n";
01720 cvm::proxy->close_output_stream(replicas_registry_file);
01721 }
01722 }
01723 
01724 if (b_hills_traj) {
01725 std::ostream &hills_traj_os =
01726 cvm::proxy->output_stream(hills_traj_file_name());
01727 if (!hills_traj_os) {
01728 error_code |= COLVARS_FILE_ERROR;
01729 }
01730 }
01731 
01732 return error_code;
01733 }
01734 
01735 
01736 std::string const colvarbias_meta::hills_traj_file_name() const
01737 {
01738 return std::string(cvm::output_prefix()+
01739 ".colvars."+this->name+
01740 ( (comm != single_replica) ?
01741 ("."+replica_id) :
01742 ("") )+
01743 ".hills.traj");
01744 }
01745 
01746 
01747 std::string const colvarbias_meta::get_state_params() const
01748 {
01749 std::ostringstream os;
01750 if (keep_hills) {
01751 os << "keepHills on" << "\n";
01752 }
01753 if (this->comm != single_replica) {
01754 os << "replicaID " << this->replica_id << "\n";
01755 }
01756 return (colvarbias::get_state_params() + os.str());
01757 }
01758 
01759 
01760 std::ostream & colvarbias_meta::write_state_data(std::ostream& os)
01761 {
01762 if (use_grids) {
01763 
01764 // this is a very good time to project hills, if you haven't done
01765 // it already!
01766 project_hills(new_hills_begin, hills.end(),
01767 hills_energy, hills_energy_gradients);
01768 new_hills_begin = hills.end();
01769 
01770 // write down the grids to the restart file
01771 os << " hills_energy\n";
01772 hills_energy->write_restart(os);
01773 os << " hills_energy_gradients\n";
01774 hills_energy_gradients->write_restart(os);
01775 }
01776 
01777 if ( (!use_grids) || keep_hills ) {
01778 // write all hills currently in memory
01779 for (std::list<hill>::const_iterator h = this->hills.begin();
01780 h != this->hills.end();
01781 h++) {
01782 os << *h;
01783 }
01784 } else {
01785 // write just those that are near the grid boundaries
01786 for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
01787 h != this->hills_off_grid.end();
01788 h++) {
01789 os << *h;
01790 }
01791 }
01792 
01793 colvarbias_ti::write_state_data(os);
01794 return os;
01795 }
01796 
01797 
01798 int colvarbias_meta::write_state_to_replicas()
01799 {
01800 int error_code = COLVARS_OK;
01801 if (comm != single_replica) {
01802 error_code |= write_replica_state_file();
01803 error_code |= reopen_replica_buffer_file();
01804 // schedule to reread the state files of the other replicas
01805 for (size_t ir = 0; ir < replicas.size(); ir++) {
01806 (replicas[ir])->replica_state_file_in_sync = false;
01807 }
01808 }
01809 return error_code;
01810 }
01811 
01812 
01813 int colvarbias_meta::write_output_files()
01814 {
01815 colvarbias_ti::write_output_files();
01816 if (dump_fes) {
01817 write_pmf();
01818 }
01819 return COLVARS_OK;
01820 }
01821 
01822 
01823 void colvarbias_meta::write_pmf()
01824 {
01825 colvarproxy *proxy = cvm::main()->proxy;
01826 // allocate a new grid to store the pmf
01827 colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy);
01828 pmf->setup();
01829 
01830 if ((comm == single_replica) || (dump_replica_fes)) {
01831 // output the PMF from this instance or replica
01832 pmf->reset();
01833 pmf->add_grid(*hills_energy);
01834 
01835 if (ebmeta) {
01836 int nt_points=pmf->number_of_points();
01837 for (int i = 0; i < nt_points; i++) {
01838 cvm::real pmf_val=0.0;
01839 cvm::real target_val=target_dist->value(i);
01840 if (target_val>0) {
01841 pmf_val=pmf->value(i);
01842 pmf_val=pmf_val + proxy->target_temperature() * proxy->boltzmann() * cvm::logn(target_val);
01843 }
01844 pmf->set_value(i,pmf_val);
01845 }
01846 }
01847 
01848 cvm::real const max = pmf->maximum_value();
01849 pmf->add_constant(-1.0 * max);
01850 pmf->multiply_constant(-1.0);
01851 if (well_tempered) {
01852 cvm::real const well_temper_scale = (bias_temperature + proxy->target_temperature()) / bias_temperature;
01853 pmf->multiply_constant(well_temper_scale);
01854 }
01855 {
01856 std::string const fes_file_name(this->output_prefix +
01857 ((comm != single_replica) ? ".partial" : "") +
01858 (dump_fes_save ?
01859 "."+cvm::to_str(cvm::step_absolute()) : "") +
01860 ".pmf");
01861 pmf->write_multicol(fes_file_name, "PMF file");
01862 }
01863 }
01864 
01865 if (comm != single_replica) {
01866 // output the combined PMF from all replicas
01867 pmf->reset();
01868 // current replica already included in the pools of replicas
01869 for (size_t ir = 0; ir < replicas.size(); ir++) {
01870 pmf->add_grid(*(replicas[ir]->hills_energy));
01871 }
01872 
01873 if (ebmeta) {
01874 int nt_points=pmf->number_of_points();
01875 for (int i = 0; i < nt_points; i++) {
01876 cvm::real pmf_val=0.0;
01877 cvm::real target_val=target_dist->value(i);
01878 if (target_val>0) {
01879 pmf_val=pmf->value(i);
01880 pmf_val=pmf_val + proxy->target_temperature() * proxy->boltzmann() * cvm::logn(target_val);
01881 }
01882 pmf->set_value(i,pmf_val);
01883 }
01884 }
01885 
01886 cvm::real const max = pmf->maximum_value();
01887 pmf->add_constant(-1.0 * max);
01888 pmf->multiply_constant(-1.0);
01889 if (well_tempered) {
01890 cvm::real const well_temper_scale = (bias_temperature + proxy->target_temperature()) / bias_temperature;
01891 pmf->multiply_constant(well_temper_scale);
01892 }
01893 std::string const fes_file_name(this->output_prefix +
01894 (dump_fes_save ?
01895 "."+cvm::to_str(cvm::step_absolute()) : "") +
01896 ".pmf");
01897 pmf->write_multicol(fes_file_name, "partial PMF file");
01898 }
01899 
01900 delete pmf;
01901 }
01902 
01903 
01904 
01905 int colvarbias_meta::write_replica_state_file()
01906 {
01907 colvarproxy *proxy = cvm::proxy;
01908 
01909 if (cvm::debug()) {
01910 cvm::log("Writing replica state file for bias \""+name+"\"\n");
01911 }
01912 
01913 int error_code = COLVARS_OK;
01914 
01915 // Write to temporary state file
01916 std::string const tmp_state_file(replica_state_file+".tmp");
01917 error_code |= proxy->remove_file(tmp_state_file);
01918 std::ostream &rep_state_os = cvm::proxy->output_stream(tmp_state_file);
01919 if (rep_state_os) {
01920 if (!write_state(rep_state_os)) {
01921 error_code |= cvm::error("Error: in writing to temporary file \""+
01922 tmp_state_file+"\".\n", COLVARS_FILE_ERROR);
01923 }
01924 }
01925 error_code |= proxy->close_output_stream(tmp_state_file);
01926 
01927 error_code |= proxy->rename_file(tmp_state_file, replica_state_file);
01928 
01929 return error_code;
01930 }
01931 
01932 
01933 int colvarbias_meta::reopen_replica_buffer_file()
01934 {
01935 int error_code = COLVARS_OK;
01936 colvarproxy *proxy = cvm::proxy;
01937 if (proxy->output_stream(replica_hills_file)) {
01938 error_code |= proxy->close_output_stream(replica_hills_file);
01939 }
01940 error_code |= proxy->remove_file(replica_hills_file);
01941 std::ostream &replica_hills_os = proxy->output_stream(replica_hills_file);
01942 if (replica_hills_os) {
01943 replica_hills_os.setf(std::ios::scientific, std::ios::floatfield);
01944 } else {
01945 error_code |= COLVARS_FILE_ERROR;
01946 }
01947 return error_code;
01948 }
01949 
01950 
01951 std::string colvarbias_meta::hill::output_traj()
01952 {
01953 std::ostringstream os;
01954 os.setf(std::ios::fixed, std::ios::floatfield);
01955 os << std::setw(cvm::it_width) << it << " ";
01956 
01957 os.setf(std::ios::scientific, std::ios::floatfield);
01958 
01959 size_t i;
01960 os << " ";
01961 for (i = 0; i < centers.size(); i++) {
01962 os << " ";
01963 os << std::setprecision(cvm::cv_prec)
01964 << std::setw(cvm::cv_width) << centers[i];
01965 }
01966 
01967 os << " ";
01968 for (i = 0; i < sigmas.size(); i++) {
01969 os << " ";
01970 os << std::setprecision(cvm::cv_prec)
01971 << std::setw(cvm::cv_width) << sigmas[i];
01972 }
01973 
01974 os << " ";
01975 os << std::setprecision(cvm::en_prec)
01976 << std::setw(cvm::en_width) << W << "\n";
01977 
01978 return os.str();
01979 }
01980 
01981 
01982 colvarbias_meta::hill::hill(cvm::step_number it_in,
01983 cvm::real W_in,
01984 std::vector<colvarvalue> const &cv_values,
01985 std::vector<cvm::real> const &cv_sigmas,
01986 std::string const &replica_in)
01987 : it(it_in),
01988 sW(1.0),
01989 W(W_in),
01990 centers(cv_values.size()),
01991 sigmas(cv_values.size()),
01992 replica(replica_in)
01993 {
01994 hill_value = 0.0;
01995 for (size_t i = 0; i < cv_values.size(); i++) {
01996 centers[i].type(cv_values[i]);
01997 centers[i] = cv_values[i];
01998 sigmas[i] = cv_sigmas[i];
01999 }
02000 if (cvm::debug()) {
02001 cvm::log("New hill, applied to "+cvm::to_str(cv_values.size())+
02002 " collective variables, with centers "+
02003 cvm::to_str(centers)+", sigmas "+
02004 cvm::to_str(sigmas)+" and weight "+
02005 cvm::to_str(W)+".\n");
02006 }
02007 }
02008 
02009 
02010 colvarbias_meta::hill::hill(colvarbias_meta::hill const &h)
02011 : it(h.it),
02012 hill_value(0.0),
02013 sW(1.0),
02014 W(h.W),
02015 centers(h.centers),
02016 sigmas(h.sigmas),
02017 replica(h.replica)
02018 {
02019 hill_value = 0.0;
02020 }
02021 
02022 
02023 colvarbias_meta::hill &
02024 colvarbias_meta::hill::operator = (colvarbias_meta::hill const &h)
02025 {
02026 it = h.it;
02027 hill_value = 0.0;
02028 sW = 1.0;
02029 W = h.W;
02030 centers = h.centers;
02031 sigmas = h.sigmas;
02032 replica = h.replica;
02033 hill_value = h.hill_value;
02034 return *this;
02035 }
02036 
02037 
02038 colvarbias_meta::hill::~hill()
02039 {}
02040 
02041 
02042 std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
02043 {
02044 os.setf(std::ios::scientific, std::ios::floatfield);
02045 
02046 os << "hill {\n";
02047 os << " step " << std::setw(cvm::it_width) << h.it << "\n";
02048 os << " weight "
02049 << std::setprecision(cvm::en_prec)
02050 << std::setw(cvm::en_width)
02051 << h.W << "\n";
02052 
02053 if (h.replica.size())
02054 os << " replicaID " << h.replica << "\n";
02055 
02056 size_t i;
02057 os << " centers ";
02058 for (i = 0; i < (h.centers).size(); i++) {
02059 os << " "
02060 << std::setprecision(cvm::cv_prec)
02061 << std::setw(cvm::cv_width)
02062 << h.centers[i];
02063 }
02064 os << "\n";
02065 
02066 // For backward compatibility, write the widths instead of the sigmas
02067 os << " widths ";
02068 for (i = 0; i < (h.sigmas).size(); i++) {
02069 os << " "
02070 << std::setprecision(cvm::cv_prec)
02071 << std::setw(cvm::cv_width)
02072 << 2.0 * h.sigmas[i];
02073 }
02074 os << "\n";
02075 
02076 os << "}\n";
02077 
02078 return os;
02079 }