Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Namespace Members | Class Members | File Members

colvarbias_meta.C

Go to the documentation of this file.
00001 #include <iostream>
00002 #include <sstream>
00003 #include <fstream>
00004 #include <iomanip>
00005 #include <cmath>
00006 #include <algorithm>
00007 
00008 // used to set the absolute path of a replica file
00009 #if defined (WIN32) && !defined(__CYGWIN__)
00010 #include <direct.h>
00011 #define CHDIR ::_chdir
00012 #define GETCWD ::_getcwd
00013 #define PATHSEP "\\"
00014 #else
00015 #include <unistd.h>
00016 #define CHDIR ::chdir
00017 #define GETCWD ::getcwd
00018 #define PATHSEP "/"
00019 #endif
00020 
00021 
00022 #include "colvar.h"
00023 #include "colvarbias_meta.h"
00024 
00025 
00026 colvarbias_meta::colvarbias_meta()
00027   : colvarbias(),
00028     state_file_step (0),
00029     new_hills_begin (hills.end())
00030 {
00031 }
00032 
00033 
00034 colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
00035   : colvarbias (conf, key),
00036     state_file_step (0),
00037     new_hills_begin (hills.end())
00038 {
00039   if (cvm::n_abf_biases > 0)
00040     cvm::log ("Warning: ABF and metadynamics running at the "
00041               "same time can lead to inconsistent results.\n");
00042 
00043   get_keyval (conf, "hillWeight", hill_weight, 0.01);
00044   if (hill_weight == 0.0)
00045     cvm::log ("Warning: hillWeight has been set to zero, "
00046               "this bias will have no effect.\n");
00047 
00048   get_keyval (conf, "newHillFrequency", new_hill_freq, 1000);
00049 
00050   get_keyval (conf, "hillWidth", hill_width, std::sqrt (2.0 * PI) / 2.0);
00051 
00052   {
00053     bool b_replicas = false;
00054     get_keyval (conf, "multipleReplicas", b_replicas, false);
00055     if (b_replicas)
00056       comm = multiple_replicas;
00057     else 
00058       comm = single_replica;
00059   }
00060 
00061   get_keyval (conf, "useGrids", use_grids, true);
00062 
00063   if (use_grids) {
00064     get_keyval (conf, "gridsUpdateFrequency", grids_freq, new_hill_freq);
00065     get_keyval (conf, "rebinGrids", rebin_grids, false);
00066 
00067     expand_grids = false;
00068     for (size_t i = 0; i < colvars.size(); i++) {
00069       if (colvars[i]->expand_boundaries) {
00070         expand_grids = true;
00071         cvm::log ("Metadynamics bias \""+this->name+"\""+
00072                   ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00073                   ": Will expand grids when the colvar \""+
00074                   colvars[i]->name+"\" approaches its boundaries.\n");
00075       }
00076     }
00077 
00078     get_keyval (conf, "keepHills", keep_hills, false);
00079     get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true);
00080     get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false);
00081 
00082     for (size_t i = 0; i < colvars.size(); i++) {
00083       colvars[i]->enable (colvar::task_grid);
00084     }
00085 
00086     hills_energy           = new colvar_grid_scalar   (colvars);
00087     hills_energy_gradients = new colvar_grid_gradient (colvars);
00088   } else {
00089     rebin_grids = false;
00090     keep_hills = false;
00091     dump_fes = false;
00092     dump_fes_save = false;
00093     dump_replica_fes = false;
00094   }
00095 
00096   if (comm != single_replica) {
00097 
00098     if (expand_grids)
00099       cvm::fatal_error ("Error: expandBoundaries is not supported when "
00100                         "using more than one replicas; please allocate "
00101                         "wide enough boundaries for each colvar"
00102                         "ahead of time.\n");
00103 
00104     if (get_keyval (conf, "dumpPartialFreeEnergyFile", dump_replica_fes, false)) {
00105       if (dump_replica_fes && (! dump_fes)) {
00106         cvm::log ("Enabling \"dumpFreeEnergyFile\".\n");
00107       }
00108     }
00109 
00110     get_keyval (conf, "replicaID", replica_id, std::string (""));
00111     if (!replica_id.size())
00112       cvm::fatal_error ("Error: replicaID must be defined "
00113                         "when using more than one replica.\n");
00114 
00115     get_keyval (conf, "replicasRegistry",
00116                 replicas_registry_file, 
00117                 (this->name+".replicas.registry.txt"));
00118 
00119     get_keyval (conf, "replicaUpdateFrequency",
00120                 replica_update_freq, new_hill_freq);
00121 
00122     if (keep_hills)
00123       cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+
00124                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00125                 ": keepHills with more than one replica can lead to a very "
00126                 "large amount input/output and slow down your calculations.  "
00127                 "Please consider disabling it.\n");
00128 
00129 
00130     {
00131       // TODO: one may want to specify the path manually for intricated filesystems?
00132       char *pwd = new char[3001];
00133       if (GETCWD (pwd, 3000) == NULL)
00134         cvm::fatal_error ("Error: cannot get the path of the current working directory.\n");
00135       replica_list_file =
00136         (std::string (pwd)+std::string (PATHSEP)+
00137          this->name+"."+replica_id+".files.txt");
00138       // replica_hills_file and replica_state_file are those written
00139       // by the current replica; within the mirror biases, they are
00140       // those by another replica
00141       replica_hills_file =
00142         (std::string (pwd)+std::string (PATHSEP)+
00143          cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills");
00144       replica_state_file =
00145         (std::string (pwd)+std::string (PATHSEP)+
00146          cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state");
00147       delete pwd;
00148     }
00149 
00150     // now register this replica
00151 
00152     // first check that it isn't already there
00153     bool registered_replica = false;
00154     std::ifstream reg_is (replicas_registry_file.c_str());
00155     if (reg_is.good()) {  // the file may not be there yet
00156       std::string existing_replica ("");
00157       std::string existing_replica_file ("");
00158       while ((reg_is >> existing_replica) && existing_replica.size() &&
00159              (reg_is >> existing_replica_file) && existing_replica_file.size()) {
00160         if (existing_replica == replica_id) {
00161           // this replica was already registered
00162           replica_list_file = existing_replica_file;
00163           reg_is.close();
00164           registered_replica = true;
00165           break;
00166         }
00167       }
00168       reg_is.close();
00169     }
00170 
00171     // if this replica was not included yet, we should generate a
00172     // new record for it: but first, we write this replica's files,
00173     // for the others to read
00174       
00175     // open the "hills" buffer file
00176     replica_hills_os.open (replica_hills_file.c_str());
00177     if (!replica_hills_os.good())
00178       cvm::fatal_error ("Error: in opening file \""+
00179                         replica_hills_file+"\" for writing.\n");
00180     replica_hills_os.setf (std::ios::scientific, std::ios::floatfield);
00181 
00182     // write the state file (so that there is always one available)
00183     write_replica_state_file();
00184     // schedule to read the state files of the other replicas
00185     for (size_t ir = 0; ir < replicas.size(); ir++) {
00186       (replicas[ir])->replica_state_file_in_sync = false;
00187     }
00188 
00189     // if we're running without grids, use a growing list of "hills" files
00190     // otherwise, just one state file and one "hills" file as buffer
00191     std::ofstream list_os (replica_list_file.c_str(), 
00192                            (use_grids ? std::ios::trunc : std::ios::app));
00193     if (! list_os.good())
00194       cvm::fatal_error ("Error: in opening file \""+
00195                         replica_list_file+"\" for writing.\n");
00196     list_os << "stateFile " << replica_state_file << "\n";
00197     list_os << "hillsFile " << replica_hills_file << "\n";
00198     list_os.close();
00199 
00200     // finally, if add a new record for this replica to the registry
00201     if (! registered_replica) {
00202       std::ofstream reg_os (replicas_registry_file.c_str(), std::ios::app);
00203       if (! reg_os.good())
00204         cvm::fatal_error ("Error: in opening file \""+
00205                           replicas_registry_file+"\" for writing.\n");
00206       reg_os << replica_id << " " << replica_list_file << "\n";
00207       reg_os.close();
00208     }
00209   }
00210 
00211   get_keyval (conf, "writeHillsTrajectory", b_hills_traj, false);
00212   if (b_hills_traj) {
00213     std::string const traj_file_name (cvm::output_prefix+
00214                                       ".colvars."+this->name+
00215                                       ( (comm != single_replica) ?
00216                                         ("."+replica_id) :
00217                                         ("") )+
00218                                       ".hills.traj");
00219     hills_traj_os.open (traj_file_name.c_str());
00220     if (!hills_traj_os.good())
00221       cvm::fatal_error ("Error: in opening hills output file \"" +
00222                         traj_file_name + "\".\n");
00223   }
00224 
00225   // for well-tempered metadynamics
00226   get_keyval (conf, "wellTempered", well_tempered, false);
00227   get_keyval (conf, "biasTemperature", bias_temperature, -1.0);
00228   if ((bias_temperature == -1.0) && well_tempered) {
00229     cvm::fatal_error ("Error: biasTemperature is not set.\n");
00230   }
00231   if (well_tempered) {
00232     cvm::log("Well-tempered metadynamics is used.\n");
00233     cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
00234   }
00235   
00236   if (cvm::debug())
00237     cvm::log ("Done initializing the metadynamics bias \""+this->name+"\""+
00238               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
00239 
00240   save_delimiters = false;
00241 }
00242 
00243 
00244 colvarbias_meta::~colvarbias_meta()
00245 {
00246   if (hills_energy) {
00247     delete hills_energy;
00248     hills_energy = NULL;
00249   }
00250 
00251   if (hills_energy_gradients) {
00252     delete hills_energy_gradients;
00253     hills_energy_gradients = NULL;
00254   }
00255 
00256   if (replica_hills_os.good())
00257     replica_hills_os.close();
00258 
00259   if (hills_traj_os.good())
00260     hills_traj_os.close();
00261 }
00262 
00263 
00264 
00265 // **********************************************************************
00266 // Hill management member functions
00267 // **********************************************************************
00268 
00269 std::list<colvarbias_meta::hill>::const_iterator 
00270 colvarbias_meta::create_hill (colvarbias_meta::hill const &h)
00271 {
00272   hill_iter const hills_end = hills.end();
00273   hills.push_back (h);
00274   if (new_hills_begin == hills_end) {
00275     // if new_hills_begin is unset, set it for the first time
00276     new_hills_begin = hills.end();
00277     new_hills_begin--;
00278   }
00279 
00280   if (use_grids) {
00281 
00282     // also add it to the list of hills that are off-grid, which may
00283     // need to be computed analytically when the colvar returns
00284     // off-grid
00285     cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers, true);
00286     if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
00287       hills_off_grid.push_back (h);
00288     }
00289   }
00290 
00291   // output to trajectory (if specified)
00292   if (hills_traj_os.good()) {
00293     hills_traj_os << (hills.back()).output_traj();
00294     if (cvm::debug()) {
00295       hills_traj_os.flush();
00296     }
00297   }
00298 
00299   has_data = true;
00300   return hills.end();
00301 }
00302 
00303 
00304 std::list<colvarbias_meta::hill>::const_iterator
00305 colvarbias_meta::delete_hill (hill_iter &h)
00306 {
00307   if (cvm::debug()) {
00308     cvm::log ("Deleting hill from the metadynamics bias \""+this->name+"\""+
00309               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00310               ", with step number "+
00311               cvm::to_str (h->it)+(h->replica.size() ?
00312                                    ", replica id \""+h->replica :
00313                                    "")+".\n");
00314   }
00315 
00316   if (use_grids && hills_off_grid.size()) {
00317     for (hill_iter hoff = hills_off_grid.begin();
00318          hoff != hills_off_grid.end(); hoff++) {
00319       if (*h == *hoff) {
00320         hills_off_grid.erase (hoff);
00321         break;
00322       }
00323     }
00324   }
00325 
00326   if (hills_traj_os.good()) {
00327     // output to the trajectory 
00328     hills_traj_os << "# DELETED this hill: "
00329                   << (hills.back()).output_traj()
00330                   << "\n";
00331     if (cvm::debug())
00332       hills_traj_os.flush();
00333   }
00334 
00335   return hills.erase (h);
00336 }
00337 
00338 
00339 cvm::real colvarbias_meta::update()
00340 {
00341   if (cvm::debug())
00342     cvm::log ("Updating the metadynamics bias \""+this->name+"\""+
00343               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
00344 
00345   if (use_grids) {
00346 
00347     std::vector<int> curr_bin = hills_energy->get_colvars_index();
00348 
00349     if (expand_grids) {
00350 
00351       // first of all, expand the grids, if specified
00352       if (cvm::debug())
00353         cvm::log ("Metadynamics bias \""+this->name+"\""+
00354                   ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00355                   ": current coordinates on the grid: "+
00356                   cvm::to_str (curr_bin)+".\n");
00357 
00358       bool changed_grids = false;
00359       int const min_buffer =
00360         (3 * (size_t) std::floor (hill_width)) + 1;
00361 
00362       std::vector<int>         new_sizes (hills_energy->sizes());
00363       std::vector<colvarvalue> new_lower_boundaries (hills_energy->lower_boundaries);
00364       std::vector<colvarvalue> new_upper_boundaries (hills_energy->upper_boundaries);
00365 
00366       for (size_t i = 0; i < colvars.size(); i++) {
00367 
00368         if (! colvars[i]->expand_boundaries)
00369           continue;
00370 
00371         cvm::real &new_lb   = new_lower_boundaries[i].real_value;
00372         cvm::real &new_ub   = new_upper_boundaries[i].real_value;
00373         int       &new_size = new_sizes[i];
00374         bool changed_lb = false, changed_ub = false;
00375 
00376         if (!colvars[i]->hard_lower_boundary)
00377           if (curr_bin[i] < min_buffer) {
00378             int const extra_points = (min_buffer - curr_bin[i]);
00379             new_lb -= extra_points * colvars[i]->width;
00380             new_size += extra_points;
00381             // changed offset in this direction => the pointer needs to
00382             // be changed, too
00383             curr_bin[i] += extra_points;
00384 
00385             changed_lb = true;
00386             cvm::log ("Metadynamics bias \""+this->name+"\""+
00387                       ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00388                       ": new lower boundary for colvar \""+
00389                       colvars[i]->name+"\", at "+
00390                       cvm::to_str (new_lower_boundaries[i])+".\n");
00391           }
00392 
00393         if (!colvars[i]->hard_upper_boundary)
00394           if (curr_bin[i] > new_size - min_buffer - 1) {
00395             int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
00396             new_ub += extra_points * colvars[i]->width;
00397             new_size += extra_points;
00398 
00399             changed_ub = true;
00400             cvm::log ("Metadynamics bias \""+this->name+"\""+
00401                       ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00402                       ": new upper boundary for colvar \""+
00403                       colvars[i]->name+"\", at "+
00404                       cvm::to_str (new_upper_boundaries[i])+".\n");
00405           }
00406 
00407         if (changed_lb || changed_ub)
00408           changed_grids = true;
00409       }
00410 
00411       if (changed_grids) {
00412 
00413         // map everything into new grids
00414 
00415         colvar_grid_scalar *new_hills_energy =
00416           new colvar_grid_scalar (*hills_energy);
00417         colvar_grid_gradient *new_hills_energy_gradients =
00418           new colvar_grid_gradient (*hills_energy_gradients);
00419 
00420         // supply new boundaries to the new grids
00421 
00422         new_hills_energy->lower_boundaries = new_lower_boundaries;
00423         new_hills_energy->upper_boundaries = new_upper_boundaries;
00424         new_hills_energy->create (new_sizes, 0.0, 1);
00425 
00426         new_hills_energy_gradients->lower_boundaries = new_lower_boundaries;
00427         new_hills_energy_gradients->upper_boundaries = new_upper_boundaries;
00428         new_hills_energy_gradients->create (new_sizes, 0.0, colvars.size());
00429 
00430         new_hills_energy->map_grid (*hills_energy);
00431         new_hills_energy_gradients->map_grid (*hills_energy_gradients);
00432 
00433         delete hills_energy;
00434         delete hills_energy_gradients;
00435         hills_energy = new_hills_energy;
00436         hills_energy_gradients = new_hills_energy_gradients;
00437 
00438         curr_bin = hills_energy->get_colvars_index();
00439         if (cvm::debug())
00440           cvm::log ("Coordinates on the new grid: "+
00441                     cvm::to_str (curr_bin)+".\n");
00442       }
00443     }
00444   }
00445 
00446   // add a new hill if the required time interval has passed
00447   if ((cvm::step_absolute() % new_hill_freq) == 0) {
00448 
00449     if (cvm::debug())
00450       cvm::log ("Metadynamics bias \""+this->name+"\""+
00451                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00452                 ": adding a new hill at step "+cvm::to_str (cvm::step_absolute())+".\n");
00453 
00454     switch (comm) {
00455 
00456     case single_replica:
00457       if (well_tempered) {
00458         std::vector<int> curr_bin = hills_energy->get_colvars_index();
00459         cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin);
00460         cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
00461         create_hill (hill ((hill_weight*exp_weight), colvars, hill_width));
00462       } else { 
00463         create_hill (hill (hill_weight, colvars, hill_width));
00464       } 
00465       break;
00466 
00467     case multiple_replicas:
00468       if (well_tempered) {
00469         std::vector<int> curr_bin = hills_energy->get_colvars_index();
00470         cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin);
00471         cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
00472         create_hill (hill ((hill_weight*exp_weight), colvars, hill_width, replica_id));
00473       } else { 
00474         create_hill (hill (hill_weight, colvars, hill_width, replica_id));
00475       } 
00476       if (replica_hills_os.good()) {
00477         replica_hills_os << hills.back();
00478       } else {
00479         cvm::fatal_error ("Error: in metadynamics bias \""+this->name+"\""+
00480                           ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00481                           " while writing hills for the other replicas.\n");
00482       }
00483       break;
00484     }
00485   }
00486 
00487   // sync with the other replicas (if needed)
00488   if (comm != single_replica) {
00489 
00490     // reread the replicas registry
00491     if ((cvm::step_absolute() % replica_update_freq) == 0) {
00492       update_replicas_registry();
00493       // empty the output buffer
00494       replica_hills_os.flush();
00495 
00496       read_replica_files();
00497     }
00498   }
00499 
00500   // calculate the biasing energy and forces
00501   bias_energy = 0.0;
00502   for (size_t ir = 0; ir < colvars.size(); ir++) {
00503     colvar_forces[ir].reset();
00504   }
00505   if (comm == multiple_replicas)
00506     for (size_t ir = 0; ir < replicas.size(); ir++) {
00507       replicas[ir]->bias_energy = 0.0;
00508       for (size_t ic = 0; ic < colvars.size(); ic++) {
00509         replicas[ir]->colvar_forces[ic].reset();
00510       }
00511     }
00512 
00513   if (use_grids) {
00514 
00515     // get the forces from the grid
00516 
00517     if ((cvm::step_absolute() % grids_freq) == 0) {
00518       // map the most recent gaussians to the grids
00519       project_hills (new_hills_begin, hills.end(),
00520                      hills_energy,    hills_energy_gradients);
00521       new_hills_begin = hills.end();
00522 
00523       // TODO: we may want to condense all into one replicas array,
00524       // including "this" as the first element
00525       if (comm == multiple_replicas) {
00526         for (size_t ir = 0; ir < replicas.size(); ir++) {
00527           replicas[ir]->project_hills (replicas[ir]->new_hills_begin,
00528                                        replicas[ir]->hills.end(),
00529                                        replicas[ir]->hills_energy,
00530                                        replicas[ir]->hills_energy_gradients);
00531           replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
00532         }
00533       }
00534     }
00535 
00536     std::vector<int> curr_bin = hills_energy->get_colvars_index();
00537     if (cvm::debug())
00538       cvm::log ("Metadynamics bias \""+this->name+"\""+
00539                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00540                 ": current coordinates on the grid: "+
00541                 cvm::to_str (curr_bin)+".\n");
00542 
00543     if (hills_energy->index_ok (curr_bin)) {
00544 
00545       // within the grid: add the energy and the forces from there
00546 
00547       bias_energy += hills_energy->value (curr_bin);
00548       for (size_t ic = 0; ic < colvars.size(); ic++) {
00549         cvm::real const *f = &(hills_energy_gradients->value (curr_bin));
00550         colvar_forces[ic].real_value += -1.0 * f[ic];
00551         // the gradients are stored, not the forces
00552       }
00553       if (comm == multiple_replicas)
00554         for (size_t ir = 0; ir < replicas.size(); ir++) {
00555           bias_energy += replicas[ir]->hills_energy->value (curr_bin);
00556           cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value (curr_bin));
00557           for (size_t ic = 0; ic < colvars.size(); ic++) {
00558             colvar_forces[ic].real_value += -1.0 * f[ic];
00559           }
00560         }
00561 
00562     } else {
00563 
00564       // off the grid: compute analytically only the hills at the grid's edges
00565 
00566       calc_hills (hills_off_grid.begin(), hills_off_grid.end(), bias_energy);
00567       for (size_t ic = 0; ic < colvars.size(); ic++) {
00568         calc_hills_force (ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces);
00569       }
00570 
00571       if (comm == multiple_replicas)
00572         for (size_t ir = 0; ir < replicas.size(); ir++) {
00573           calc_hills (replicas[ir]->hills_off_grid.begin(),
00574                       replicas[ir]->hills_off_grid.end(),
00575                       bias_energy);
00576           for (size_t ic = 0; ic < colvars.size(); ic++) {
00577             calc_hills_force (ic,
00578                               replicas[ir]->hills_off_grid.begin(),
00579                               replicas[ir]->hills_off_grid.end(),
00580                               colvar_forces);
00581           }
00582         }
00583     }
00584   }
00585 
00586   // now include the hills that have not been binned yet (starting
00587   // from new_hills_begin)
00588 
00589   calc_hills (new_hills_begin, hills.end(), bias_energy);
00590   for (size_t ic = 0; ic < colvars.size(); ic++) {
00591     calc_hills_force (ic, new_hills_begin, hills.end(), colvar_forces);
00592   }
00593 
00594   if (cvm::debug()) 
00595     cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
00596               ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
00597 
00598   if (cvm::debug()) 
00599     cvm::log ("Metadynamics bias \""+this->name+"\""+
00600               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00601               ": adding the forces from the other replicas.\n");
00602 
00603   if (comm == multiple_replicas)
00604     for (size_t ir = 0; ir < replicas.size(); ir++) {
00605       calc_hills (replicas[ir]->new_hills_begin,
00606                   replicas[ir]->hills.end(),
00607                   bias_energy);
00608       for (size_t ic = 0; ic < colvars.size(); ic++) {
00609         calc_hills_force (ic,
00610                           replicas[ir]->new_hills_begin,
00611                           replicas[ir]->hills.end(),
00612                           colvar_forces);
00613       }
00614       if (cvm::debug()) 
00615         cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
00616                   ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
00617     }
00618 
00619   return bias_energy;
00620 }
00621 
00622 
00623 void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter      h_first,
00624                                   colvarbias_meta::hill_iter      h_last,
00625                                   cvm::real                      &energy,
00626                                   std::vector<colvarvalue> const &colvar_values)
00627 {
00628   std::vector<colvarvalue> curr_values (colvars.size());
00629   for (size_t i = 0; i < colvars.size(); i++) {
00630     curr_values[i].type (colvars[i]->type());
00631   }
00632 
00633   if (colvar_values.size()) {
00634     for (size_t i = 0; i < colvars.size(); i++) {
00635       curr_values[i] = colvar_values[i];
00636     }
00637   } else {
00638     for (size_t i = 0; i < colvars.size(); i++) {
00639       curr_values[i] = colvars[i]->value();
00640     }
00641   }
00642 
00643   for (hill_iter h = h_first; h != h_last; h++) {
00644       
00645     // compute the gaussian exponent
00646     cvm::real cv_sqdev = 0.0;
00647     for (size_t i = 0; i < colvars.size(); i++) {
00648       colvarvalue const &x  = curr_values[i];
00649       colvarvalue const &center = h->centers[i];
00650       cvm::real const    half_width = 0.5 * h->widths[i];
00651       cv_sqdev += (colvars[i]->dist2 (x, center)) / (half_width*half_width);
00652     }
00653 
00654     // compute the gaussian
00655     if (cv_sqdev > 23.0) {
00656       // set it to zero if the exponent is more negative than log(1.0E-05)
00657       h->value (0.0);
00658     } else {
00659       h->value (std::exp (-0.5*cv_sqdev));
00660     }
00661     energy += h->energy();
00662   }
00663 }
00664 
00665 
00666 void colvarbias_meta::calc_hills_force (size_t const &i,
00667                                         colvarbias_meta::hill_iter      h_first,
00668                                         colvarbias_meta::hill_iter      h_last,
00669                                         std::vector<colvarvalue>       &forces,
00670                                         std::vector<colvarvalue> const &values)
00671 {
00672   // Retrieve the value of the colvar
00673   colvarvalue x (values.size() ? values[i].type() : colvars[i]->type());
00674   x = (values.size() ? values[i] : colvars[i]->value());
00675 
00676   // do the type check only once (all colvarvalues in the hills series
00677   // were already saved with their types matching those in the
00678   // colvars)
00679 
00680   switch (colvars[i]->type()) {
00681 
00682   case colvarvalue::type_scalar:
00683     for (hill_iter h = h_first; h != h_last; h++) {
00684       if (h->value() == 0.0) continue;
00685       colvarvalue const &center = h->centers[i];
00686       cvm::real const    half_width = 0.5 * h->widths[i];
00687       forces[i].real_value += 
00688         ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * 
00689           (colvars[i]->dist2_lgrad (x, center)).real_value );
00690     }
00691     break;
00692 
00693   case colvarvalue::type_vector:
00694   case colvarvalue::type_unitvector:
00695     for (hill_iter h = h_first; h != h_last; h++) {
00696       if (h->value() == 0.0) continue;
00697       colvarvalue const &center = h->centers[i];
00698       cvm::real const    half_width = 0.5 * h->widths[i];
00699       forces[i].rvector_value +=
00700         ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
00701           (colvars[i]->dist2_lgrad (x, center)).rvector_value );
00702     }
00703     break;
00704 
00705   case colvarvalue::type_quaternion:
00706     for (hill_iter h = h_first; h != h_last; h++) {
00707       if (h->value() == 0.0) continue;
00708       colvarvalue const &center = h->centers[i];
00709       cvm::real const    half_width = 0.5 * h->widths[i];
00710       forces[i].quaternion_value +=
00711         ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
00712           (colvars[i]->dist2_lgrad (x, center)).quaternion_value );
00713     }
00714     break;
00715 
00716   case colvarvalue::type_notset:
00717     break;
00718   }
00719 }
00720 
00721 
00722 // **********************************************************************
00723 // grid management functions
00724 // **********************************************************************
00725 
00726 void colvarbias_meta::project_hills (colvarbias_meta::hill_iter  h_first,
00727                                      colvarbias_meta::hill_iter  h_last,
00728                                      colvar_grid_scalar         *he,
00729                                      colvar_grid_gradient       *hg,
00730                                      bool print_progress)
00731 {
00732   if (cvm::debug())
00733     cvm::log ("Metadynamics bias \""+this->name+"\""+
00734               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00735               ": projecting hills.\n");
00736 
00737   // TODO: improve it by looping over a small subgrid instead of the whole grid
00738 
00739   std::vector<colvarvalue> colvar_values (colvars.size());
00740   std::vector<cvm::real> colvar_forces_scalar (colvars.size());
00741 
00742   std::vector<int> he_ix = he->new_index();
00743   std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
00744   cvm::real hills_energy_here = 0.0;
00745   std::vector<colvarvalue> hills_forces_here (colvars.size(), 0.0);
00746 
00747   size_t count = 0;
00748   size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1)));
00749 
00750   if (hg != NULL) {
00751 
00752     // loop over the points of the grid
00753     for ( ;
00754           (he->index_ok (he_ix)) && (hg->index_ok (hg_ix));
00755           count++) {
00756 
00757       for (size_t i = 0; i < colvars.size(); i++) {
00758         colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
00759       }
00760 
00761       // loop over the hills and increment the energy grid locally
00762       hills_energy_here = 0.0;
00763       calc_hills (h_first, h_last, hills_energy_here, colvar_values);
00764       he->acc_value (he_ix, hills_energy_here);
00765 
00766       for (size_t i = 0; i < colvars.size(); i++) {
00767         hills_forces_here[i].reset();
00768         calc_hills_force (i, h_first, h_last, hills_forces_here, colvar_values);
00769         colvar_forces_scalar[i] = hills_forces_here[i].real_value;
00770       }
00771       hg->acc_force (hg_ix, &(colvar_forces_scalar.front()));
00772 
00773       he->incr (he_ix);
00774       hg->incr (hg_ix);
00775 
00776       if ((count % print_frequency) == 0) {
00777         if (print_progress) {
00778           cvm::real const progress = cvm::real (count) / cvm::real (hg->number_of_points());
00779           std::ostringstream os;
00780           os.setf (std::ios::fixed, std::ios::floatfield);
00781           os << std::setw (6) << std::setprecision (2)
00782              << 100.0 * progress
00783              << "% done.";
00784           cvm::log (os.str());
00785         }
00786       }
00787     }
00788 
00789   } else {
00790 
00791     // simpler version, with just the energy
00792 
00793     for ( ; (he->index_ok (he_ix)); ) {
00794 
00795       for (size_t i = 0; i < colvars.size(); i++) {
00796         colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
00797       }
00798 
00799       hills_energy_here = 0.0;
00800       calc_hills (h_first, h_last, hills_energy_here, colvar_values);
00801       he->acc_value (he_ix, hills_energy_here);
00802 
00803       he->incr (he_ix);
00804 
00805       count++;
00806       if ((count % print_frequency) == 0) {
00807         if (print_progress) {
00808           cvm::real const progress = cvm::real (count) / cvm::real (he->number_of_points());
00809           std::ostringstream os;
00810           os.setf (std::ios::fixed, std::ios::floatfield);
00811           os << std::setw (6) << std::setprecision (2)
00812              << 100.0 * progress
00813              << "% done.";
00814           cvm::log (os.str());
00815         }
00816       }
00817     }
00818   }
00819 
00820   if (print_progress) {
00821     cvm::log ("100.00% done.");
00822   }
00823 
00824   if (! keep_hills) {
00825     hills.erase (hills.begin(), hills.end());
00826   }
00827 }
00828 
00829 
00830 void colvarbias_meta::recount_hills_off_grid (colvarbias_meta::hill_iter  h_first,
00831                                               colvarbias_meta::hill_iter  h_last,
00832                                               colvar_grid_scalar         *he)
00833 {
00834   hills_off_grid.clear();
00835 
00836   for (hill_iter h = h_first; h != h_last; h++) {
00837     cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers, true);
00838     if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
00839       hills_off_grid.push_back (*h);
00840     }
00841   }
00842 }
00843 
00844 
00845 
00846 // **********************************************************************
00847 // multiple replicas functions
00848 // **********************************************************************
00849 
00850 
00851 void colvarbias_meta::update_replicas_registry()
00852 {
00853   if (cvm::debug())
00854     cvm::log ("Metadynamics bias \""+this->name+"\""+
00855               ": updating the list of replicas, currently containing "+
00856               cvm::to_str (replicas.size())+" elements.\n");
00857 
00858   {
00859     // copy the whole file into a string for convenience
00860     std::string line ("");
00861     std::ifstream reg_file (replicas_registry_file.c_str());
00862     if (reg_file.good()) {
00863       replicas_registry.clear();
00864       while (colvarparse::getline_nocomments (reg_file, line))
00865         replicas_registry.append (line+"\n");
00866     } else {
00867       cvm::fatal_error ("Error: failed to open file \""+replicas_registry_file+
00868                         "\" for reading.\n");
00869     }
00870   }
00871 
00872   // now parse it
00873   std::istringstream reg_is (replicas_registry);
00874   if (reg_is.good()) {
00875 
00876     std::string new_replica ("");
00877     std::string new_replica_file ("");
00878     while ((reg_is >> new_replica) && new_replica.size() &&
00879            (reg_is >> new_replica_file) && new_replica_file.size()) {
00880 
00881       if (new_replica == this->replica_id) {
00882         // this is the record for this same replica, skip it
00883         if (cvm::debug())
00884           cvm::log ("Metadynamics bias \""+this->name+"\""+
00885                     ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00886                     ": skipping this replica's own record: \""+
00887                     new_replica+"\", \""+new_replica_file+"\"\n");
00888         new_replica_file.clear();
00889         new_replica.clear();
00890         continue;
00891       }
00892 
00893       bool already_loaded = false;
00894       for (size_t ir = 0; ir < replicas.size(); ir++) {
00895         if (new_replica == (replicas[ir])->replica_id) {
00896           // this replica was already added
00897           if (cvm::debug())
00898             cvm::log ("Metadynamics bias \""+this->name+"\""+
00899                       ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
00900                       ": skipping a replica already loaded, \""+
00901                       (replicas[ir])->replica_id+"\".\n");
00902           already_loaded = true;
00903           break;
00904         }
00905       }
00906 
00907       if (!already_loaded) {
00908         // add this replica to the registry
00909         cvm::log ("Metadynamics bias \""+this->name+"\""+
00910                   ": accessing replica \""+new_replica+"\".\n");
00911         replicas.push_back (new colvarbias_meta());
00912         (replicas.back())->replica_id = new_replica;
00913         (replicas.back())->replica_list_file = new_replica_file;
00914         (replicas.back())->replica_state_file = "";
00915         (replicas.back())->replica_state_file_in_sync = false;
00916 
00917         // Note: the following could become a copy constructor?
00918         (replicas.back())->name = this->name;
00919         (replicas.back())->colvars = colvars;
00920         (replicas.back())->use_grids = use_grids;
00921         (replicas.back())->dump_fes = false;
00922         (replicas.back())->expand_grids = false;
00923         (replicas.back())->rebin_grids = false;
00924         (replicas.back())->keep_hills = false;
00925         (replicas.back())->colvar_forces = colvar_forces;
00926 
00927         (replicas.back())->comm = multiple_replicas;
00928 
00929         if (use_grids) {
00930           (replicas.back())->hills_energy           = new colvar_grid_scalar   (colvars);
00931           (replicas.back())->hills_energy_gradients = new colvar_grid_gradient (colvars);
00932         }
00933       }
00934     }
00935 
00936     // continue the cycle
00937     new_replica_file = "";
00938     new_replica = "";
00939   } else {
00940     cvm::fatal_error ("Error: cannot read the replicas registry file \""+
00941                       replicas_registry+"\".\n");
00942   }
00943 
00944   // now (re)read the list file of each replica
00945   for (size_t ir = 0; ir < replicas.size(); ir++) {
00946     if (cvm::debug())
00947       cvm::log ("Metadynamics bias \""+this->name+"\""+
00948                 ": reading the list file for replica \""+
00949                 (replicas[ir])->replica_id+"\".\n");
00950 
00951     std::ifstream list_is ((replicas[ir])->replica_list_file.c_str());
00952     std::string key;
00953     std::string new_state_file, new_hills_file;
00954     if (!(list_is >> key) ||
00955         !(key == std::string ("stateFile")) ||
00956         !(list_is >> new_state_file) ||
00957         !(list_is >> key) ||
00958         !(key == std::string ("hillsFile")) ||
00959         !(list_is >> new_hills_file)) {
00960       cvm::log ("Metadynamics bias \""+this->name+"\""+
00961                 ": failed to read the file \""+
00962                 (replicas[ir])->replica_list_file+"\": will try again after "+
00963                 cvm::to_str (replica_update_freq)+" steps.\n");
00964       (replicas[ir])->update_status++;
00965     } else {
00966       (replicas[ir])->update_status = 0;
00967       if (new_state_file != (replicas[ir])->replica_state_file) {
00968         cvm::log ("Metadynamics bias \""+this->name+"\""+
00969                   ": replica \""+(replicas[ir])->replica_id+
00970                   "\" has supplied a new state file, \""+new_state_file+
00971                   "\".\n");
00972         (replicas[ir])->replica_state_file_in_sync = false;
00973         (replicas[ir])->replica_state_file = new_state_file;
00974         (replicas[ir])->replica_hills_file = new_hills_file;
00975       }
00976     }
00977   }
00978 
00979 
00980   if (cvm::debug())
00981     cvm::log ("Metadynamics bias \""+this->name+"\": the list of replicas contains "+
00982               cvm::to_str (replicas.size())+" elements.\n");
00983 }
00984 
00985 
00986 void colvarbias_meta::read_replica_files()
00987 {
00988   for (size_t ir = 0; ir < replicas.size(); ir++) {
00989 
00990     if (! (replicas[ir])->replica_state_file_in_sync) {
00991       // if a new state file is being read, the hills file is also new
00992       (replicas[ir])->replica_hills_file_pos = 0;
00993     }
00994 
00995     // (re)read the state file if necessary
00996     if ( (! (replicas[ir])->has_data) || 
00997          (! (replicas[ir])->replica_state_file_in_sync) ) {
00998 
00999       cvm::log ("Metadynamics bias \""+this->name+"\""+
01000                 ": reading the state of replica \""+
01001                 (replicas[ir])->replica_id+"\" from file \""+
01002                 (replicas[ir])->replica_state_file+"\".\n");
01003 
01004       std::ifstream is ((replicas[ir])->replica_state_file.c_str());
01005       if (! (replicas[ir])->read_restart (is)) {
01006         cvm::log ("Reading from file \""+(replicas[ir])->replica_state_file+
01007                   "\" failed or incomplete: will try again in "+
01008                   cvm::to_str (replica_update_freq)+" steps.\n");
01009       } else {
01010         // state file has been read successfully
01011         (replicas[ir])->replica_state_file_in_sync = true;
01012         (replicas[ir])->update_status = 0;
01013       }
01014       is.close();
01015     }
01016     
01017     // now read the hills added after writing the state file
01018     if ((replicas[ir])->replica_hills_file.size()) {
01019 
01020       if (cvm::debug()) 
01021         cvm::log ("Metadynamics bias \""+this->name+"\""+
01022                   ": checking for new hills from replica \""+
01023                   (replicas[ir])->replica_id+"\" in the file \""+
01024                   (replicas[ir])->replica_hills_file+"\".\n");
01025 
01026       // read hills from the other replicas' files; for each file, resume
01027       // the position recorded previously
01028 
01029       std::ifstream is ((replicas[ir])->replica_hills_file.c_str());
01030       if (is.good()) {
01031 
01032         // try to resume the previous position
01033         is.seekg ((replicas[ir])->replica_hills_file_pos, std::ios::beg);
01034         if (!is.good()){
01035           // if fail (the file may have been overwritten), reset this
01036           // position
01037           is.clear();
01038           is.seekg (0, std::ios::beg);
01039           // reset the counter
01040           (replicas[ir])->replica_hills_file_pos = 0;
01041           // schedule to reread the state file
01042           (replicas[ir])->replica_state_file_in_sync = false;
01043           // and record the failure
01044           (replicas[ir])->update_status++;
01045           cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+
01046                     "\" at the previous position: will try again in "+
01047                     cvm::to_str (replica_update_freq)+" steps.\n");
01048         } else {
01049 
01050           while ((replicas[ir])->read_hill (is)) {
01051             //           if (cvm::debug())
01052               cvm::log ("Metadynamics bias \""+this->name+"\""+
01053                         ": received a hill from replica \""+
01054                         (replicas[ir])->replica_id+
01055                         "\" at step "+
01056                         cvm::to_str (((replicas[ir])->hills.back()).it)+".\n");
01057           }
01058           is.clear();
01059           // store the position for the next read
01060           (replicas[ir])->replica_hills_file_pos = is.tellg();
01061           if (cvm::debug())
01062             cvm::log ("Metadynamics bias \""+this->name+"\""+
01063                       ": stopped reading file \""+(replicas[ir])->replica_hills_file+
01064                       "\" at position "+
01065                       cvm::to_str ((replicas[ir])->replica_hills_file_pos)+".\n");
01066 
01067           // test whether this is the end of the file 
01068           is.seekg (0, std::ios::end);
01069           if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) {
01070             (replicas[ir])->update_status++;
01071           } else {
01072             (replicas[ir])->update_status = 0;
01073           }
01074         }
01075 
01076       } else {
01077         cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+
01078                   "\": will try again in "+
01079                   cvm::to_str (replica_update_freq)+" steps.\n");
01080         (replicas[ir])->update_status++;
01081         // cvm::fatal_error ("Error: cannot read from file \""+
01082         //                   (replicas[ir])->replica_hills_file+"\".\n");
01083       }
01084       is.close();
01085     }
01086 
01087     size_t const n_flush = (replica_update_freq/new_hill_freq + 1);
01088     if ((replicas[ir])->update_status > 3*n_flush) {
01089       // TODO: suspend the calculation?
01090       cvm::log ("WARNING: in metadynamics bias \""+this->name+"\""+
01091                 " failed to read completely the output of replica \""+
01092                 (replicas[ir])->replica_id+
01093                 "\" after more than "+
01094                 cvm::to_str ((replicas[ir])->update_status * replica_update_freq)+
01095                 " steps.  Ensure that it is still running.\n");
01096     }
01097   }
01098 }
01099 
01100 
01101 // **********************************************************************
01102 // input functions
01103 // **********************************************************************
01104 
01105 
01106 std::istream & colvarbias_meta::read_restart (std::istream& is)
01107 {
01108   size_t const start_pos = is.tellg();
01109 
01110   if (comm == single_replica) {
01111     // if using a multiple replicas scheme, output messages 
01112     // are printed before and after calling this function
01113     cvm::log ("Restarting metadynamics bias \""+this->name+"\""+
01114               ".\n");
01115   }
01116   std::string key, brace, conf;
01117   if ( !(is >> key)   || !(key == "metadynamics") ||
01118        !(is >> brace) || !(brace == "{") ||
01119        !(is >> colvarparse::read_block ("configuration", conf)) ) {
01120 
01121     if (comm == single_replica) 
01122       cvm::log ("Error: in reading restart configuration for metadynamics bias \""+
01123                 this->name+"\""+
01124                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01125                 (replica_state_file_in_sync ? ("at position "+
01126                                                cvm::to_str (start_pos)+
01127                                                " in the state file") : "")+".\n");
01128     is.clear();
01129     is.seekg (start_pos, std::ios::beg);
01130     is.setstate (std::ios::failbit);
01131     return is;
01132   }
01133 
01134   std::string name = "";
01135   if ( colvarparse::get_keyval (conf, "name", name,
01136                                 std::string (""), colvarparse::parse_silent) &&
01137        (name != this->name) )
01138     cvm::fatal_error ("Error: in the restart file, the "
01139                       "\"metadynamics\" block has a different name: different system?\n");
01140 
01141   if (name.size() == 0) {
01142     cvm::fatal_error ("Error: \"metadynamics\" block within the restart file "
01143                       "has no identifiers.\n");
01144   }
01145 
01146   if (comm != single_replica) {
01147     std::string replica = "";
01148     if (colvarparse::get_keyval (conf, "replicaID", replica,
01149                                  std::string (""), colvarparse::parse_silent) &&
01150         (replica != this->replica_id))
01151       cvm::fatal_error ("Error: in the restart file, the "
01152                         "\"metadynamics\" block has a different replicaID: different system?\n");
01153 
01154     colvarparse::get_keyval (conf, "step", state_file_step,
01155                              cvm::step_absolute(), colvarparse::parse_silent);
01156   }
01157 
01158   bool grids_from_restart_file = use_grids;
01159 
01160   if (use_grids) {
01161 
01162     if (expand_grids) {
01163       // the boundaries of the colvars may have been changed; TODO:
01164       // this reallocation is only for backward-compatibility, and may
01165       // be deleted when grid_parameters (i.e. colvargrid's own
01166       // internal reallocation) has kicked in
01167       delete hills_energy;
01168       delete hills_energy_gradients;
01169       hills_energy = new colvar_grid_scalar (colvars);
01170       hills_energy_gradients = new colvar_grid_gradient (colvars);
01171     }
01172 
01173     colvar_grid_scalar   *hills_energy_backup = NULL;
01174     colvar_grid_gradient *hills_energy_gradients_backup = NULL;
01175 
01176     if (has_data) {
01177       if (cvm::debug())
01178         cvm::log ("Backupping grids for metadynamics bias \""+
01179                   this->name+"\""+
01180                   ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
01181       hills_energy_backup           = hills_energy;
01182       hills_energy_gradients_backup = hills_energy_gradients;
01183       hills_energy                  = new colvar_grid_scalar (colvars);
01184       hills_energy_gradients        = new colvar_grid_gradient (colvars);
01185     }
01186 
01187     size_t const hills_energy_pos = is.tellg();
01188     if (!(is >> key)) {
01189       if (hills_energy_backup != NULL) {
01190         delete hills_energy;
01191         delete hills_energy_gradients;
01192         hills_energy           = hills_energy_backup;
01193         hills_energy_gradients = hills_energy_gradients_backup;
01194       }
01195       is.clear();
01196       is.seekg (hills_energy_pos, std::ios::beg);
01197       is.setstate (std::ios::failbit);
01198       return is;
01199     } else if (!(key == std::string ("hills_energy")) ||
01200                !(hills_energy->read_restart (is))) {
01201       is.clear();
01202       is.seekg (hills_energy_pos, std::ios::beg);
01203       grids_from_restart_file = false;
01204       if (!rebin_grids) {
01205         if (hills_energy_backup == NULL)
01206           cvm::fatal_error ("Error: couldn't read the free energy grid for metadynamics bias \""+
01207                             this->name+"\""+
01208                             ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01209                             "; if useGrids was off when the state file was written, "
01210                             "enable rebinGrids now to regenerate the grids.\n");
01211         else {
01212           if (comm == single_replica)
01213             cvm::log ("Error: couldn't read the free energy grid for metadynamics bias \""+
01214                       this->name+"\".\n");
01215           delete hills_energy;
01216           delete hills_energy_gradients;
01217           hills_energy           = hills_energy_backup;
01218           hills_energy_gradients = hills_energy_gradients_backup;
01219           is.setstate (std::ios::failbit);
01220           return is;
01221         }
01222       }
01223     }
01224 
01225     size_t const hills_energy_gradients_pos = is.tellg();
01226     if (!(is >> key)) {
01227       if (hills_energy_backup != NULL)  {
01228         delete hills_energy;
01229         delete hills_energy_gradients;
01230         hills_energy           = hills_energy_backup;
01231         hills_energy_gradients = hills_energy_gradients_backup;
01232       }
01233       is.clear();
01234       is.seekg (hills_energy_gradients_pos, std::ios::beg);
01235       is.setstate (std::ios::failbit);
01236       return is;
01237     } else if (!(key == std::string ("hills_energy_gradients")) ||
01238                !(hills_energy_gradients->read_restart (is))) {
01239       is.clear();
01240       is.seekg (hills_energy_gradients_pos, std::ios::beg);
01241       grids_from_restart_file = false;
01242       if (!rebin_grids) { 
01243         if (hills_energy_backup == NULL)
01244           cvm::fatal_error ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
01245                             this->name+"\""+
01246                             ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01247                             "; if useGrids was off when the state file was written, "
01248                             "enable rebinGrids now to regenerate the grids.\n");
01249         else {
01250           if (comm == single_replica)
01251             cvm::log ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
01252                       this->name+"\".\n");
01253           delete hills_energy;
01254           delete hills_energy_gradients;
01255           hills_energy           = hills_energy_backup;
01256           hills_energy_gradients = hills_energy_gradients_backup;
01257           is.setstate (std::ios::failbit);
01258           return is;
01259         }
01260       }
01261     }
01262 
01263     if (cvm::debug())
01264       cvm::log ("Successfully read new grids for bias \""+
01265                 this->name+"\""+
01266                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
01267 
01268     if (hills_energy_backup != NULL) {
01269       // now that we have successfully updated the grids, delete the
01270       // backup copies
01271       if (cvm::debug())
01272         cvm::log ("Deallocating the older grids.\n");
01273 
01274       delete hills_energy_backup;
01275       delete hills_energy_gradients_backup;
01276     }
01277   }
01278 
01279   bool const existing_hills = (hills.size() > 0);
01280   size_t const old_hills_size = hills.size();
01281   hill_iter old_hills_end = hills.end();
01282   hill_iter old_hills_off_grid_end = hills_off_grid.end();
01283 
01284   // read the hills explicitly written (if there are any)
01285   while (read_hill (is)) {
01286     if (cvm::debug()) 
01287       cvm::log ("Read a previously saved hill under the "
01288                 "metadynamics bias \""+
01289                 this->name+"\", created at step "+
01290                 cvm::to_str ((hills.back()).it)+".\n");
01291   }
01292   is.clear();
01293   new_hills_begin = hills.end();
01294   if (grids_from_restart_file) {
01295     if (hills.size() > old_hills_size)
01296       cvm::log ("Read "+cvm::to_str (hills.size())+
01297                 " hills in addition to the grids.\n");
01298   } else {
01299     if (hills.size())
01300       cvm::log ("Read "+cvm::to_str (hills.size())+" hills.\n");
01301   }
01302 
01303   if (rebin_grids) {
01304 
01305     // allocate new grids (based on the new boundaries and widths just
01306     // read from the configuration file), and project onto them the
01307     // grids just read from the restart file
01308 
01309     colvar_grid_scalar   *new_hills_energy =
01310       new colvar_grid_scalar (colvars);
01311     colvar_grid_gradient *new_hills_energy_gradients =
01312       new colvar_grid_gradient (colvars);
01313 
01314     if (!grids_from_restart_file || (keep_hills && hills.size())) {
01315       // if there are hills, recompute the new grids from them
01316       cvm::log ("Rebinning the energy and forces grids from "+
01317                 cvm::to_str (hills.size())+" hills (this may take a while)...\n");
01318       project_hills (hills.begin(), hills.end(),
01319                      new_hills_energy, new_hills_energy_gradients, true);
01320       cvm::log ("rebinning done.\n");
01321 
01322     } else {
01323       // otherwise, use the grids in the restart file
01324       cvm::log ("Rebinning the energy and forces grids "
01325                 "from the grids in the restart file.\n");
01326       new_hills_energy->map_grid (*hills_energy);
01327       new_hills_energy_gradients->map_grid (*hills_energy_gradients);
01328     }
01329 
01330     delete hills_energy;
01331     delete hills_energy_gradients;
01332     hills_energy = new_hills_energy;
01333     hills_energy_gradients = new_hills_energy_gradients;
01334 
01335     // assuming that some boundaries have expanded, eliminate those
01336     // off-grid hills that aren't necessary any more
01337     if (hills.size())
01338       recount_hills_off_grid (hills.begin(), hills.end(), hills_energy);
01339   }
01340 
01341   if (use_grids) {
01342     if (hills_off_grid.size()) {
01343       cvm::log (cvm::to_str (hills_off_grid.size())+" hills are near the "
01344                 "grid boundaries: they will be computed analytically "
01345                 "and saved to the state files.\n");
01346     }
01347   }
01348 
01349   is >> brace;
01350   if (brace != "}") {
01351     cvm::log ("Incomplete restart information for metadynamics bias \""+
01352               this->name+"\""+
01353               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
01354               ": no closing brace at position "+
01355               cvm::to_str (is.tellg())+" in the file.\n");
01356     is.setstate (std::ios::failbit);
01357     return is;
01358   }
01359 
01360   if (cvm::debug())
01361     cvm::log ("colvarbias_meta::read_restart() done\n");
01362 
01363   if (existing_hills) {
01364     hills.erase (hills.begin(), old_hills_end);
01365     hills_off_grid.erase (hills_off_grid.begin(), old_hills_off_grid_end);
01366   }
01367   
01368   has_data = true;
01369 
01370   if (comm != single_replica) {
01371     read_replica_files();
01372   }
01373 
01374   return is;
01375 }  
01376 
01377 
01378 std::istream & colvarbias_meta::read_hill (std::istream &is)
01379 {
01380   if (!is) return is; // do nothing if failbit is set
01381 
01382   size_t const start_pos = is.tellg();
01383 
01384   std::string data;
01385   if ( !(is >> read_block ("hill", data)) ) {
01386     is.clear();
01387     is.seekg (start_pos, std::ios::beg);
01388     is.setstate (std::ios::failbit);
01389     return is;
01390   }
01391 
01392   size_t h_it;
01393   get_keyval (data, "step", h_it, 0, parse_silent);
01394   if (h_it <= state_file_step) {
01395     if (cvm::debug())
01396       cvm::log ("Skipping a hill older than the state file for metadynamics bias \""+
01397                 this->name+"\""+
01398                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
01399     return is;
01400   }
01401 
01402   cvm::real h_weight;
01403   get_keyval (data, "weight", h_weight, hill_weight, parse_silent);
01404 
01405   std::vector<colvarvalue> h_centers (colvars.size());
01406   for (size_t i = 0; i < colvars.size(); i++) {
01407     h_centers[i].type ((colvars[i]->value()).type()); 
01408   }
01409   {
01410     // it is safer to read colvarvalue objects one at a time;
01411     // TODO: change this it later
01412     std::string centers_input;
01413     key_lookup (data, "centers", centers_input);
01414     std::istringstream centers_is (centers_input);
01415     for (size_t i = 0; i < colvars.size(); i++) {
01416       centers_is >> h_centers[i];
01417     }
01418   }
01419 
01420   std::vector<cvm::real> h_widths (colvars.size());
01421   get_keyval (data, "widths", h_widths,
01422               std::vector<cvm::real> (colvars.size(), (std::sqrt (2.0 * PI) / 2.0)),
01423               parse_silent);
01424   
01425   std::string h_replica = "";
01426   if (comm != single_replica) {
01427     get_keyval (data, "replicaID", h_replica, replica_id, parse_silent);
01428     if (h_replica != replica_id)
01429       cvm::fatal_error ("Error: trying to read a hill created by replica \""+h_replica+
01430                         "\" for replica \""+replica_id+
01431                         "\"; did you swap output files?\n");
01432   }
01433 
01434   hill_iter const hills_end = hills.end();
01435   hills.push_back (hill (h_it, h_weight, h_centers, h_widths, h_replica));
01436   if (new_hills_begin == hills_end) {
01437     // if new_hills_begin is unset, set it for the first time
01438     new_hills_begin = hills.end();
01439     new_hills_begin--;
01440   }
01441 
01442   if (use_grids) {
01443     // add this also to the list of hills that are off-grid, which will
01444     // be computed analytically
01445     cvm::real const min_dist =
01446       hills_energy->bin_distance_from_boundaries ((hills.back()).centers, true);
01447     if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
01448       hills_off_grid.push_back (hills.back());
01449     }
01450   }
01451 
01452   has_data = true;
01453   return is;
01454 }
01455 
01456 
01457 
01458 
01459 // **********************************************************************
01460 // output functions
01461 // **********************************************************************
01462 
01463 std::ostream & colvarbias_meta::write_restart (std::ostream& os)
01464 {
01465   os << "metadynamics {\n"
01466      << "  configuration {\n"
01467      << "    step " << cvm::step_absolute() << "\n"
01468      << "    name " << this->name << "\n";
01469   if (this->comm != single_replica)
01470     os << "    replicaID " << this->replica_id << "\n";
01471   os << "  }\n\n";
01472 
01473   if (use_grids) {
01474 
01475     // this is a very good time to project hills, if you haven't done
01476     // it already!
01477     project_hills (new_hills_begin, hills.end(),
01478                    hills_energy,    hills_energy_gradients);
01479     new_hills_begin = hills.end();
01480 
01481     // write down the grids to the restart file
01482     os << "  hills_energy\n";
01483     hills_energy->write_restart (os);
01484     os << "  hills_energy_gradients\n";
01485     hills_energy_gradients->write_restart (os);
01486   }
01487 
01488   if ( (!use_grids) || keep_hills ) {
01489     // write all hills currently in memory
01490     for (std::list<hill>::const_iterator h = this->hills.begin();
01491          h != this->hills.end();
01492          h++) {
01493       os << *h;
01494     }
01495   } else {
01496     // write just those that are near the grid boundaries
01497     for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
01498          h != this->hills_off_grid.end();
01499          h++) {
01500       os << *h;
01501     }
01502   }
01503 
01504   os << "}\n\n";
01505 
01506   if (comm != single_replica) {
01507     write_replica_state_file();
01508     // schedule to reread the state files of the other replicas (they
01509     // have also rewritten them)
01510     for (size_t ir = 0; ir < replicas.size(); ir++) {
01511       (replicas[ir])->replica_state_file_in_sync = false;
01512     }
01513   }
01514 
01515   if (dump_fes) {
01516     write_pmf();
01517   }
01518 
01519   return os;
01520 }  
01521 
01522 
01523 void colvarbias_meta::write_pmf()
01524 {
01525   // allocate a new grid to store the pmf
01526   colvar_grid_scalar *pmf = new colvar_grid_scalar (*hills_energy);
01527   pmf->create();
01528 
01529   std::string fes_file_name_prefix (cvm::output_prefix);
01530   
01531   if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) {
01532     // if this is not the only free energy integrator, append
01533     // this bias's name, to distinguish it from the output of the other
01534     // biases producing a .pmf file
01535     // TODO: fix for ABF with updateBias == no
01536     fes_file_name_prefix += ("."+this->name);
01537   }
01538 
01539   if ((comm == single_replica) || (dump_replica_fes)) {
01540     // output the PMF from this instance or replica
01541     pmf->reset();
01542     pmf->add_grid (*hills_energy);
01543     cvm::real const max = pmf->maximum_value();
01544     pmf->add_constant (-1.0 * max);
01545     pmf->multiply_constant (-1.0);
01546     if (well_tempered) {
01547       cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
01548       pmf->multiply_constant (well_temper_scale);
01549     }
01550     {
01551       std::string const fes_file_name (fes_file_name_prefix +
01552                                        ((comm != single_replica) ? ".partial" : "") +
01553                                        (dump_fes_save ?
01554                                         "."+cvm::to_str (cvm::step_absolute()) : "") +
01555                                        ".pmf");
01556       cvm::backup_file (fes_file_name.c_str());
01557       std::ofstream fes_os (fes_file_name.c_str());
01558       pmf->write_multicol (fes_os);
01559       fes_os.close();
01560     }
01561   }
01562   if (comm != single_replica) {
01563     // output the combined PMF from all replicas
01564     pmf->reset();
01565     pmf->add_grid (*hills_energy);
01566     for (size_t ir = 0; ir < replicas.size(); ir++) {
01567       pmf->add_grid (*(replicas[ir]->hills_energy));
01568     }
01569     cvm::real const max = pmf->maximum_value();
01570     pmf->add_constant (-1.0 * max);
01571     pmf->multiply_constant (-1.0);
01572     if (well_tempered) {
01573       cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
01574       pmf->multiply_constant (well_temper_scale);
01575     }
01576     std::string const fes_file_name (fes_file_name_prefix +
01577                                      (dump_fes_save ?
01578                                       "."+cvm::to_str (cvm::step_absolute()) : "") +
01579                                      ".pmf");
01580     cvm::backup_file (fes_file_name.c_str());
01581     std::ofstream fes_os (fes_file_name.c_str());
01582     pmf->write_multicol (fes_os);
01583     fes_os.close();
01584   }
01585 
01586   delete pmf;
01587 }
01588 
01589 
01590 
01591 void colvarbias_meta::write_replica_state_file()
01592 {
01593   // write down also the restart for the other replicas: TODO: this
01594   // is duplicated code, that could be cleaned up later
01595   cvm::backup_file (replica_state_file.c_str());
01596   std::ofstream rep_state_os (replica_state_file.c_str());
01597   if (!rep_state_os.good())
01598     cvm::fatal_error ("Error: in opening file \""+
01599                       replica_state_file+"\" for writing.\n");
01600 
01601   rep_state_os.setf (std::ios::scientific, std::ios::floatfield); 
01602   rep_state_os << "\n"
01603                << "metadynamics {\n"
01604                << "  configuration {\n"
01605                << "    name " << this->name << "\n"
01606                << "    step " << cvm::step_absolute() << "\n";
01607   if (this->comm != single_replica) {
01608     rep_state_os << "    replicaID " << this->replica_id << "\n";
01609   }
01610   rep_state_os << "  }\n\n";
01611   rep_state_os << "  hills_energy\n";
01612   rep_state_os << std::setprecision (cvm::cv_prec)
01613                << std::setw (cvm::cv_width);
01614   hills_energy->write_restart (rep_state_os);
01615   rep_state_os << "  hills_energy_gradients\n";
01616   rep_state_os << std::setprecision (cvm::cv_prec)
01617                << std::setw (cvm::cv_width);
01618   hills_energy_gradients->write_restart (rep_state_os);
01619 
01620   if ( (!use_grids) || keep_hills ) {
01621     // write all hills currently in memory
01622     for (std::list<hill>::const_iterator h = this->hills.begin();
01623          h != this->hills.end();
01624          h++) {
01625       rep_state_os << *h;
01626     }
01627   } else {
01628     // write just those that are near the grid boundaries
01629     for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
01630          h != this->hills_off_grid.end();
01631          h++) {
01632       rep_state_os << *h;
01633     }
01634   }
01635 
01636   rep_state_os << "}\n\n";
01637   rep_state_os.close();
01638 
01639   // reopen the hills file
01640   replica_hills_os.close();
01641   replica_hills_os.open (replica_hills_file.c_str());
01642   if (!replica_hills_os.good())
01643     cvm::fatal_error ("Error: in opening file \""+
01644                       replica_hills_file+"\" for writing.\n");
01645   replica_hills_os.setf (std::ios::scientific, std::ios::floatfield);
01646 }
01647 
01648 std::string colvarbias_meta::hill::output_traj()
01649 {
01650   std::ostringstream os;
01651   os.setf (std::ios::fixed, std::ios::floatfield);
01652   os << std::setw (cvm::it_width) << it << " ";
01653 
01654   os.setf (std::ios::scientific, std::ios::floatfield);
01655 
01656   os << "  ";
01657   for (size_t i = 0; i < centers.size(); i++) {
01658     os << " ";
01659     os << std::setprecision (cvm::cv_prec)
01660        << std::setw (cvm::cv_width)  << centers[i];
01661   }
01662 
01663   os << "  ";
01664   for (size_t i = 0; i < widths.size(); i++) {
01665     os << " ";
01666     os << std::setprecision (cvm::cv_prec)
01667        << std::setw (cvm::cv_width) << widths[i];
01668   }
01669 
01670   os << "  ";
01671   os << std::setprecision (cvm::en_prec)
01672      << std::setw (cvm::en_width) << W << "\n";
01673 
01674   return os.str();
01675 }  
01676 
01677 
01678 std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
01679 {
01680   os.setf (std::ios::scientific, std::ios::floatfield);
01681 
01682   os << "hill {\n";
01683   os << "  step " << std::setw (cvm::it_width) << h.it << "\n";
01684   os << "  weight   "
01685      << std::setprecision (cvm::en_prec)
01686      << std::setw (cvm::en_width)
01687      << h.W << "\n";
01688 
01689   if (h.replica.size())
01690     os << "  replicaID  " << h.replica << "\n";
01691 
01692   os << "  centers ";
01693   for (size_t i = 0; i < (h.centers).size(); i++) {
01694     os << " "
01695        << std::setprecision (cvm::cv_prec)
01696        << std::setw (cvm::cv_width)
01697        << h.centers[i];
01698   }
01699   os << "\n";
01700 
01701   os << "  widths  ";
01702   for (size_t i = 0; i < (h.widths).size(); i++) {
01703     os << " "
01704        << std::setprecision (cvm::cv_prec)
01705        << std::setw (cvm::cv_width)
01706        << h.widths[i];
01707   }
01708   os << "\n";
01709 
01710   os << "}\n";
01711 
01712   return os;
01713 }

Generated on Sun May 19 04:07:44 2013 for NAMD by  doxygen 1.3.9.1