Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

colvarbias_meta.C

Go to the documentation of this file.
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 }

Generated on Fri Oct 4 02:43:36 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002