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

colvarbias.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 <cstring>
00012 
00013 #include "colvarmodule.h"
00014 #include "colvarproxy.h"
00015 #include "colvarvalue.h"
00016 #include "colvarbias.h"
00017 #include "colvargrid.h"
00018 
00019 
00020 colvarbias::colvarbias(char const *key)
00021 {
00022   bias_type = colvarparse::to_lower_cppstr(key);
00023   state_keyword = bias_type;
00024 
00025   rank = -1;
00026   description = "uninitialized " + bias_type + " bias";
00027 
00028   colvarbias::init_dependencies();
00029 
00030   time_step_factor = 1;
00031 
00032   has_data = false;
00033   b_output_energy = false;
00034   output_freq = cvm::restart_out_freq;
00035 
00036   colvarbias::reset();
00037   state_file_step = 0L;
00038   matching_state = false;
00039   biasing_force_scaling_factors = NULL;
00040 }
00041 
00042 
00043 int colvarbias::init(std::string const &conf)
00044 {
00045   int error_code = COLVARS_OK;
00046 
00047   name = bias_type + cvm::to_str(rank);
00048   colvarparse::set_string(conf);
00049 
00050   size_t i = 0;
00051 
00052   if (num_variables() == 0) {
00053     // First initialization
00054 
00055     cvm::log("Initializing a new \""+bias_type+"\" instance.\n");
00056 
00057     // Only allow setting a non-default name on first init
00058     get_keyval(conf, "name", name, name);
00059 
00060     colvarbias *bias_with_name = cvm::bias_by_name(this->name);
00061     if (bias_with_name != NULL) {
00062       if ((bias_with_name->rank != this->rank) ||
00063           (bias_with_name->bias_type != this->bias_type)) {
00064         error_code |= cvm::error("Error: this bias cannot have the same name, \""+
00065                                  this->name+"\", as another bias.\n",
00066                                  COLVARS_INPUT_ERROR);
00067       }
00068     }
00069     description = "bias " + name;
00070 
00071     {
00072       // lookup the associated colvars
00073       std::vector<std::string> colvar_names;
00074       if (get_keyval(conf, "colvars", colvar_names)) {
00075         if (num_variables()) {
00076           error_code |= cvm::error("Error: cannot redefine the colvars that "
00077                                    "a bias was already defined on.\n",
00078                                    COLVARS_INPUT_ERROR);
00079         }
00080         for (i = 0; i < colvar_names.size(); i++) {
00081           add_colvar(colvar_names[i]);
00082         }
00083       }
00084     }
00085 
00086     if (!num_variables()) {
00087       error_code |= cvm::error("Error: no collective variables specified.\n",
00088                                COLVARS_INPUT_ERROR);
00089     }
00090 
00091   } else {
00092     cvm::log("Reinitializing bias \""+name+"\".\n");
00093   }
00094 
00095   colvar_values.resize(num_variables());
00096   for (i = 0; i < num_variables(); i++) {
00097     colvar_values[i].type(colvars[i]->value().type());
00098     colvar_forces[i].type(colvar_values[i].type());
00099     previous_colvar_forces[i].type(colvar_values[i].type());
00100   }
00101 
00102   output_prefix = cvm::output_prefix();
00103 
00104   get_keyval_feature(this, conf, "stepZeroData", f_cvb_step_zero_data, is_enabled(f_cvb_step_zero_data));
00105 
00106   // Write energy to traj file?
00107   get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
00108 
00109   // How often to write full output files?
00110   get_keyval(conf, "outputFreq", output_freq, output_freq);
00111 
00112   // Disabled by default in base class; default value can be overridden by derived class constructor
00113   get_keyval_feature(this, conf, "bypassExtendedLagrangian", f_cvb_bypass_ext_lagrangian, is_enabled(f_cvb_bypass_ext_lagrangian), parse_echo);
00114 
00115   get_keyval(conf, "timeStepFactor", time_step_factor, time_step_factor);
00116   if (time_step_factor < 1) {
00117     error_code |= cvm::error("Error: timeStepFactor must be 1 or greater.\n",
00118                              COLVARS_INPUT_ERROR);
00119   }
00120 
00121   // Use the scaling factors from a grid?
00122   get_keyval_feature(this, conf, "scaledBiasingForce",
00123                      f_cvb_scale_biasing_force,
00124                      is_enabled(f_cvb_scale_biasing_force), parse_echo);
00125   if (is_enabled(f_cvb_scale_biasing_force)) {
00126     std::string biasing_force_scaling_factors_in_filename;
00127     get_keyval(conf, "scaledBiasingForceFactorsGrid",
00128                biasing_force_scaling_factors_in_filename, std::string());
00129     biasing_force_scaling_factors = new colvar_grid_scalar(colvars);
00130     error_code |= biasing_force_scaling_factors->read_multicol(biasing_force_scaling_factors_in_filename,
00131                                                                "grid file");
00132     biasing_force_scaling_factors_bin.assign(num_variables(), 0);
00133   }
00134 
00135   // Now that children are defined, we can solve dependencies
00136   enable(f_cvb_active);
00137   if (cvm::debug()) print_state();
00138 
00139   return error_code;
00140 }
00141 
00142 
00143 int colvarbias::init_dependencies() {
00144   int i;
00145   if (features().size() == 0) {
00146     for (i = 0; i < f_cvb_ntot; i++) {
00147       modify_features().push_back(new feature);
00148     }
00149 
00150     init_feature(f_cvb_active, "active", f_type_dynamic);
00151     require_feature_children(f_cvb_active, f_cv_active);
00152 
00153     init_feature(f_cvb_awake, "awake", f_type_static);
00154     require_feature_self(f_cvb_awake, f_cvb_active);
00155 
00156     init_feature(f_cvb_step_zero_data, "step_zero_data", f_type_user);
00157 
00158     init_feature(f_cvb_apply_force, "apply_force", f_type_user);
00159     require_feature_children(f_cvb_apply_force, f_cv_gradient);
00160 
00161     init_feature(f_cvb_bypass_ext_lagrangian, "bypass_extended_Lagrangian_coordinates", f_type_user);
00162 
00163     // The exclusion below prevents the inconsistency where biasing forces are applied onto
00164     // the actual colvar, while total forces are measured on the extended coordinate
00165     exclude_feature_self(f_cvb_bypass_ext_lagrangian, f_cvb_get_total_force);
00166 
00167     init_feature(f_cvb_get_total_force, "obtain_total_force", f_type_dynamic);
00168     require_feature_children(f_cvb_get_total_force, f_cv_total_force);
00169 
00170     init_feature(f_cvb_output_acc_work, "output_accumulated_work", f_type_user);
00171     require_feature_self(f_cvb_output_acc_work, f_cvb_apply_force);
00172 
00173     init_feature(f_cvb_history_dependent, "history_dependent", f_type_static);
00174 
00175     init_feature(f_cvb_time_dependent, "time_dependent", f_type_static);
00176 
00177     init_feature(f_cvb_scalar_variables, "require_scalar_variables", f_type_static);
00178     require_feature_children(f_cvb_scalar_variables, f_cv_scalar);
00179 
00180     init_feature(f_cvb_calc_pmf, "calculate_a_PMF", f_type_static);
00181 
00182     init_feature(f_cvb_calc_ti_samples, "calculate_TI_samples", f_type_dynamic);
00183     require_feature_self(f_cvb_calc_ti_samples, f_cvb_get_total_force);
00184     require_feature_children(f_cvb_calc_ti_samples, f_cv_grid);
00185 
00186     init_feature(f_cvb_write_ti_samples, "write_TI_samples_", f_type_user);
00187     require_feature_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples);
00188 
00189     init_feature(f_cvb_write_ti_pmf, "write_TI_PMF", f_type_user);
00190     require_feature_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples);
00191 
00192     init_feature(f_cvb_scale_biasing_force, "scale_biasing_force", f_type_user);
00193     require_feature_children(f_cvb_scale_biasing_force, f_cv_grid);
00194 
00195     // check that everything is initialized
00196     for (i = 0; i < colvardeps::f_cvb_ntot; i++) {
00197       if (is_not_set(i)) {
00198         cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
00199       }
00200     }
00201   }
00202 
00203   // Initialize feature_states for each instance
00204   feature_states.reserve(f_cvb_ntot);
00205   for (i = 0; i < f_cvb_ntot; i++) {
00206     feature_states.push_back(feature_state(true, false));
00207     // Most features are available, so we set them so
00208     // and list exceptions below
00209   }
00210 
00211   // only compute TI samples when deriving from colvarbias_ti
00212   feature_states[f_cvb_calc_ti_samples].available = false;
00213 
00214   // The feature f_cvb_bypass_ext_lagrangian is only implemented by some derived classes
00215   // (initially, harmonicWalls)
00216   feature_states[f_cvb_bypass_ext_lagrangian].available = false;
00217   // disabled by default; can be changed by derived classes that implement it
00218   feature_states[f_cvb_bypass_ext_lagrangian].enabled = false;
00219 
00220   return COLVARS_OK;
00221 }
00222 
00223 
00224 int colvarbias::reset()
00225 {
00226   bias_energy = 0.0;
00227   for (size_t i = 0; i < num_variables(); i++) {
00228     colvar_forces[i].reset();
00229   }
00230   return COLVARS_OK;
00231 }
00232 
00233 
00234 colvarbias::colvarbias()
00235   : colvarparse(), has_data(false)
00236 {}
00237 
00238 
00239 colvarbias::~colvarbias()
00240 {
00241   colvarbias::clear();
00242 }
00243 
00244 
00245 int colvarbias::clear()
00246 {
00247   free_children_deps();
00248 
00249   // Remove references to this bias from colvars
00250   for (std::vector<colvar *>::iterator cvi = colvars.begin();
00251        cvi != colvars.end();
00252        ++cvi) {
00253     for (std::vector<colvarbias *>::iterator bi = (*cvi)->biases.begin();
00254          bi != (*cvi)->biases.end();
00255          ++bi) {
00256       if ( *bi == this) {
00257         (*cvi)->biases.erase(bi);
00258         break;
00259       }
00260     }
00261   }
00262 
00263   colvarmodule *cv = cvm::main();
00264   // ...and from the colvars module
00265   for (std::vector<colvarbias *>::iterator bi = cv->biases.begin();
00266        bi != cv->biases.end();
00267        ++bi) {
00268     if ( *bi == this) {
00269       cv->biases.erase(bi);
00270       break;
00271     }
00272   }
00273 
00274   if (biasing_force_scaling_factors != NULL) {
00275     delete biasing_force_scaling_factors;
00276     biasing_force_scaling_factors = NULL;
00277     biasing_force_scaling_factors_bin.clear();
00278   }
00279 
00280   cv->config_changed();
00281 
00282   return COLVARS_OK;
00283 }
00284 
00285 
00286 int colvarbias::clear_state_data()
00287 {
00288   // no mutable content to delete for base class
00289   return COLVARS_OK;
00290 }
00291 
00292 
00293 int colvarbias::add_colvar(std::string const &cv_name)
00294 {
00295   if (colvar *cv = cvm::colvar_by_name(cv_name)) {
00296 
00297     if (cvm::debug()) {
00298       cvm::log("Applying this bias to collective variable \""+
00299                cv->name+"\".\n");
00300     }
00301 
00302     colvars.push_back(cv);
00303     cv->biases.push_back(this); // add back-reference to this bias to colvar
00304 
00305     // Add dependency link. All biases need at least the value of each colvar
00306     // although possibly not at all timesteps
00307     add_child(cv);
00308 
00309     colvar_forces.push_back(colvarvalue());
00310     colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero
00311     colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force
00312     colvar_forces.back().reset();
00313     previous_colvar_forces.push_back(colvar_forces.back());
00314 
00315   } else {
00316     cvm::error("Error: cannot find a colvar named \""+
00317                cv_name+"\".\n", COLVARS_INPUT_ERROR);
00318     return COLVARS_INPUT_ERROR;
00319   }
00320 
00321   return COLVARS_OK;
00322 }
00323 
00324 
00325 int colvarbias::update()
00326 {
00327   if (cvm::debug()) {
00328     cvm::log("Updating the "+bias_type+" bias \""+this->name+"\".\n");
00329   }
00330 
00331   int error_code = COLVARS_OK;
00332 
00333   has_data = true;
00334 
00335   // Update the cached colvar values
00336   for (size_t i = 0; i < num_variables(); i++) {
00337     colvar_values[i] = colvars[i]->value();
00338   }
00339 
00340   error_code |= calc_energy(NULL);
00341   error_code |= calc_forces(NULL);
00342 
00343   return error_code;
00344 }
00345 
00346 
00347 bool colvarbias::can_accumulate_data()
00348 {
00349   colvarproxy *proxy = cvm::main()->proxy;
00350   if (((cvm::step_relative() > 0) && !proxy->simulation_continuing()) ||
00351       is_enabled(f_cvb_step_zero_data)) {
00352     return true;
00353   }
00354   return false;
00355 }
00356 
00357 
00358 int colvarbias::calc_energy(std::vector<colvarvalue> const *)
00359 {
00360   bias_energy = 0.0;
00361   return COLVARS_OK;
00362 }
00363 
00364 
00365 int colvarbias::calc_forces(std::vector<colvarvalue> const *)
00366 {
00367   for (size_t ir = 0; ir < num_variables(); ir++) {
00368     colvar_forces[ir].reset();
00369   }
00370   return COLVARS_OK;
00371 }
00372 
00373 
00374 int colvarbias::communicate_forces()
00375 {
00376   int error_code = COLVARS_OK;
00377   if (! is_enabled(f_cvb_apply_force)) {
00378     return error_code;
00379   }
00380   cvm::real biasing_force_factor = 1.0;
00381   size_t i = 0;
00382   if (is_enabled(f_cvb_scale_biasing_force)) {
00383     for (i = 0; i < num_variables(); i++) {
00384       biasing_force_scaling_factors_bin[i] = biasing_force_scaling_factors->current_bin_scalar(i);
00385     }
00386     if (biasing_force_scaling_factors->index_ok(biasing_force_scaling_factors_bin)) {
00387       biasing_force_factor *= biasing_force_scaling_factors->value(biasing_force_scaling_factors_bin);
00388     }
00389   }
00390   for (i = 0; i < num_variables(); i++) {
00391     if (cvm::debug()) {
00392       cvm::log("Communicating a force to colvar \""+
00393                variables(i)->name+"\".\n");
00394     }
00395     // Impulse-style multiple timestep
00396     // Note that biases with different values of time_step_factor
00397     // may send forces to the same colvar
00398     // which is why rescaling has to happen now: the colvar is not
00399     // aware of this bias' time_step_factor
00400     if (is_enabled(f_cvb_bypass_ext_lagrangian)) {
00401       variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i] * biasing_force_factor);
00402     } else {
00403       variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i] * biasing_force_factor);
00404     }
00405     previous_colvar_forces[i] = colvar_forces[i];
00406   }
00407   return error_code;
00408 }
00409 
00410 
00411 int colvarbias::end_of_step()
00412 {
00413   return COLVARS_OK;
00414 }
00415 
00416 
00417 int colvarbias::change_configuration(std::string const & /* conf */)
00418 {
00419   cvm::error("Error: change_configuration() not implemented.\n",
00420              COLVARS_NOT_IMPLEMENTED);
00421   return COLVARS_NOT_IMPLEMENTED;
00422 }
00423 
00424 
00425 cvm::real colvarbias::energy_difference(std::string const & /* conf */)
00426 {
00427   cvm::error("Error: energy_difference() not implemented.\n",
00428              COLVARS_NOT_IMPLEMENTED);
00429   return 0.0;
00430 }
00431 
00432 
00433 // So far, these are only implemented in colvarbias_abf
00434 int colvarbias::bin_num()
00435 {
00436   cvm::error("Error: bin_num() not implemented.\n");
00437   return COLVARS_NOT_IMPLEMENTED;
00438 }
00439 int colvarbias::current_bin()
00440 {
00441   cvm::error("Error: current_bin() not implemented.\n");
00442   return COLVARS_NOT_IMPLEMENTED;
00443 }
00444 int colvarbias::bin_count(int /* bin_index */)
00445 {
00446   cvm::error("Error: bin_count() not implemented.\n");
00447   return COLVARS_NOT_IMPLEMENTED;
00448 }
00449 int colvarbias::replica_share()
00450 {
00451   cvm::error("Error: replica_share() not implemented.\n");
00452   return COLVARS_NOT_IMPLEMENTED;
00453 }
00454 
00455 
00456 std::string const colvarbias::get_state_params() const
00457 {
00458   std::ostringstream os;
00459   os << "step " << cvm::step_absolute() << "\n"
00460      << "name " << this->name << "\n";
00461   return os.str();
00462 }
00463 
00464 
00465 int colvarbias::set_state_params(std::string const &conf)
00466 {
00467   matching_state = false;
00468 
00469   std::string check_name = "";
00470   colvarparse::get_keyval(conf, "name", check_name,
00471                           std::string(""), colvarparse::parse_silent);
00472 
00473   if (check_name.size() == 0) {
00474     cvm::error("Error: \""+bias_type+"\" block within the restart file "
00475                "has no identifiers.\n", COLVARS_INPUT_ERROR);
00476   }
00477 
00478   if (check_name != this->name) {
00479     if (cvm::debug()) {
00480       cvm::log("Ignoring state of bias \""+check_name+
00481                "\": this bias is named \""+name+"\".\n");
00482     }
00483     return COLVARS_OK;
00484   }
00485 
00486   matching_state = true;
00487 
00488   colvarparse::get_keyval(conf, "step", state_file_step,
00489                           cvm::step_absolute(), colvarparse::parse_silent);
00490 
00491   return COLVARS_OK;
00492 }
00493 
00494 
00495 std::ostream & colvarbias::write_state(std::ostream &os)
00496 {
00497   if (cvm::debug()) {
00498     cvm::log("Writing state file for bias \""+name+"\"\n");
00499   }
00500   os.setf(std::ios::scientific, std::ios::floatfield);
00501   os.precision(cvm::cv_prec);
00502   os << state_keyword << " {\n"
00503      << "  configuration {\n";
00504   std::istringstream is(get_state_params());
00505   std::string line;
00506   while (std::getline(is, line)) {
00507     os << "    " << line << "\n";
00508   }
00509   os << "  }\n";
00510   write_state_data(os);
00511   os << "}\n\n";
00512   return os;
00513 }
00514 
00515 
00516 std::istream & colvarbias::read_state(std::istream &is)
00517 {
00518   std::streampos const start_pos = is.tellg();
00519 
00520   std::string key, brace, conf;
00521   if ( !(is >> key)   || !(key == state_keyword || key == bias_type) ||
00522        !(is >> brace) || !(brace == "{") ||
00523        !(is >> colvarparse::read_block("configuration", &conf)) ||
00524        (set_state_params(conf) != COLVARS_OK) ) {
00525     cvm::error("Error: in reading state configuration for \""+bias_type+
00526                "\" bias \""+
00527                this->name+"\" at position "+
00528                cvm::to_str(static_cast<size_t>(is.tellg()))+
00529                " in stream.\n", COLVARS_INPUT_ERROR);
00530     is.clear();
00531     is.seekg(start_pos, std::ios::beg);
00532     is.setstate(std::ios::failbit);
00533     return is;
00534   }
00535 
00536   if (matching_state == false) {
00537     // This state is not for this bias
00538     is.seekg(start_pos, std::ios::beg);
00539     return is;
00540   }
00541 
00542   cvm::log("Restarting "+bias_type+" bias \""+name+"\" from step number "+
00543            cvm::to_str(state_file_step)+".\n");
00544 
00545   if (!read_state_data(is)) {
00546     cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+
00547                this->name+"\" at position "+
00548                cvm::to_str(static_cast<size_t>(is.tellg()))+
00549                " in stream.\n", COLVARS_INPUT_ERROR);
00550     is.clear();
00551     is.seekg(start_pos, std::ios::beg);
00552     is.setstate(std::ios::failbit);
00553   }
00554 
00555   is >> brace;
00556   if (brace != "}") {
00557     cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+
00558                this->name+"\": no matching brace at position "+
00559                cvm::to_str(static_cast<size_t>(is.tellg()))+
00560                " in stream.\n");
00561     is.setstate(std::ios::failbit);
00562   }
00563 
00564   return is;
00565 }
00566 
00567 
00568 int colvarbias::write_state_prefix(std::string const &prefix)
00569 {
00570   std::string const filename =
00571     cvm::state_file_prefix(prefix.c_str())+".colvars.state";
00572   std::ostream &os = cvm::proxy->output_stream(filename.c_str(), "bias state file");
00573   int error_code = COLVARS_OK;
00574   if (os) {
00575     os.setf(std::ios::scientific, std::ios::floatfield);
00576     error_code = write_state(os) ? COLVARS_OK : COLVARS_FILE_ERROR;
00577   } else {
00578     error_code = COLVARS_FILE_ERROR;
00579   }
00580   cvm::proxy->close_output_stream(filename.c_str());
00581   return error_code;
00582 }
00583 
00584 
00585 int colvarbias::write_state_string(std::string &output)
00586 {
00587   std::ostringstream os;
00588   if (!write_state(os)) {
00589     return cvm::error("Error: in writing state of bias \""+name+
00590                       "\" to buffer.\n", COLVARS_FILE_ERROR);
00591   }
00592   output = os.str();
00593   return COLVARS_OK;
00594 }
00595 
00596 
00597 int colvarbias::read_state_prefix(std::string const &prefix)
00598 {
00599   std::string filename(prefix+std::string(".colvars.state"));
00600   std::istream *is = &(cvm::main()->proxy->input_stream(filename,
00601                                                         "bias state file",
00602                                                         false));
00603   if (!*is) {
00604     filename = prefix;
00605     is = &(cvm::main()->proxy->input_stream(filename, "bias state file"));
00606   }
00607 
00608   if (read_state(*is)) {
00609     return cvm::main()->proxy->close_input_stream(filename);
00610   }
00611   return COLVARS_FILE_ERROR;
00612 }
00613 
00614 
00615 int colvarbias::read_state_string(char const *buffer)
00616 {
00617   if (buffer != NULL) {
00618     size_t const buffer_size = strlen(buffer);
00619     if (cvm::debug()) {
00620       cvm::log("colvarbias::read_state_string() with argument:\n");
00621       cvm::log(buffer);
00622     }
00623 
00624     if (buffer_size > 0) {
00625       std::istringstream is;
00626       is.rdbuf()->pubsetbuf(const_cast<char *>(buffer), buffer_size);
00627       return read_state(is).good() ? COLVARS_OK :
00628         cvm::error("Error: in reading state for \""+name+"\" from buffer.\n",
00629                    COLVARS_FILE_ERROR);
00630     }
00631     return COLVARS_OK;
00632   }
00633   return cvm::error("Error: NULL pointer for colvarbias::read_state_string()",
00634                     COLVARS_BUG_ERROR);
00635 }
00636 
00637 
00638 std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key)
00639 {
00640   std::streampos const start_pos = is.tellg();
00641   std::string key_in;
00642   if ( !(is >> key_in) ||
00643        !(to_lower_cppstr(key_in) == to_lower_cppstr(std::string(key))) ) {
00644     cvm::error("Error: in reading restart configuration for "+
00645                bias_type+" bias \""+this->name+"\" at position "+
00646                cvm::to_str(static_cast<size_t>(is.tellg()))+
00647                " in stream.\n", COLVARS_INPUT_ERROR);
00648     is.clear();
00649     is.seekg(start_pos, std::ios::beg);
00650     is.setstate(std::ios::failbit);
00651     return is;
00652   }
00653   return is;
00654 }
00655 
00656 
00657 
00658 std::ostream & colvarbias::write_traj_label(std::ostream &os)
00659 {
00660   os << " ";
00661   if (b_output_energy)
00662     os << " E_"
00663        << cvm::wrap_string(this->name, cvm::en_width-2);
00664   return os;
00665 }
00666 
00667 
00668 std::ostream & colvarbias::write_traj(std::ostream &os)
00669 {
00670   os << " ";
00671   if (b_output_energy)
00672     os << " "
00673        << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
00674        << bias_energy;
00675   return os;
00676 }
00677 
00678 
00679 
00680 colvarbias_ti::colvarbias_ti(char const *key)
00681   : colvarbias(key)
00682 {
00683   colvarproxy *proxy = cvm::main()->proxy;
00684   provide(f_cvb_calc_ti_samples);
00685   if (!proxy->total_forces_same_step()) {
00686     // Samples at step zero can not be collected
00687     feature_states[f_cvb_step_zero_data].available = false;
00688   }
00689   ti_avg_forces = NULL;
00690   ti_count = NULL;
00691 }
00692 
00693 
00694 colvarbias_ti::~colvarbias_ti()
00695 {
00696   colvarbias_ti::clear_state_data();
00697 }
00698 
00699 
00700 int colvarbias_ti::clear_state_data()
00701 {
00702   if (ti_avg_forces != NULL) {
00703     delete ti_avg_forces;
00704     ti_avg_forces = NULL;
00705   }
00706   if (ti_count != NULL) {
00707     delete ti_count;
00708     ti_count = NULL;
00709   }
00710   return COLVARS_OK;
00711 }
00712 
00713 
00714 int colvarbias_ti::init(std::string const &conf)
00715 {
00716   int error_code = COLVARS_OK;
00717 
00718   get_keyval_feature(this, conf, "writeTISamples",
00719                      f_cvb_write_ti_samples,
00720                      is_enabled(f_cvb_write_ti_samples));
00721 
00722   get_keyval_feature(this, conf, "writeTIPMF",
00723                      f_cvb_write_ti_pmf,
00724                      is_enabled(f_cvb_write_ti_pmf));
00725 
00726   if ((num_variables() > 1) && is_enabled(f_cvb_write_ti_pmf)) {
00727     return cvm::error("Error: only 1-dimensional PMFs can be written "
00728                       "on the fly.\n"
00729                       "Consider using writeTISamples instead and "
00730                       "post-processing the sampled free-energy gradients.\n",
00731                       COLVARS_NOT_IMPLEMENTED);
00732   } else {
00733     error_code |= init_grids();
00734   }
00735 
00736   if (is_enabled(f_cvb_write_ti_pmf)) {
00737     enable(f_cvb_write_ti_samples);
00738   }
00739 
00740   if (is_enabled(f_cvb_calc_ti_samples)) {
00741     std::vector<std::string> const time_biases =
00742       cvm::main()->time_dependent_biases();
00743     if (time_biases.size() > 0) {
00744       if ((time_biases.size() > 1) || (time_biases[0] != this->name)) {
00745         for (size_t i = 0; i < num_variables(); i++) {
00746           if (! variables(i)->is_enabled(f_cv_subtract_applied_force)) {
00747             return cvm::error("Error: cannot collect TI samples while other "
00748                               "time-dependent biases are active and not all "
00749                               "variables have subtractAppliedForces on.\n",
00750                               COLVARS_INPUT_ERROR);
00751           }
00752         }
00753       }
00754     }
00755   }
00756 
00757   if (is_enabled(f_cvb_write_ti_pmf) || is_enabled(f_cvb_write_ti_samples)) {
00758     cvm::main()->cite_feature("Internal-forces free energy estimator");
00759   }
00760 
00761   return error_code;
00762 }
00763 
00764 
00765 int colvarbias_ti::init_grids()
00766 {
00767   if (is_enabled(f_cvb_calc_ti_samples)) {
00768     if (ti_avg_forces == NULL) {
00769       ti_bin.resize(num_variables());
00770       ti_system_forces.resize(num_variables());
00771       for (size_t icv = 0; icv < num_variables(); icv++) {
00772         ti_system_forces[icv].type(variables(icv)->value());
00773         ti_system_forces[icv].is_derivative();
00774         ti_system_forces[icv].reset();
00775       }
00776       ti_avg_forces = new colvar_grid_gradient(colvars);
00777       ti_count = new colvar_grid_count(colvars);
00778       ti_avg_forces->samples = ti_count;
00779       ti_count->has_parent_data = true;
00780     }
00781   }
00782 
00783   return COLVARS_OK;
00784 }
00785 
00786 
00787 int colvarbias_ti::update()
00788 {
00789   return update_system_forces(NULL);
00790 }
00791 
00792 
00793 int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
00794                                         *subtract_forces)
00795 {
00796   if (! is_enabled(f_cvb_calc_ti_samples)) {
00797     return COLVARS_OK;
00798   }
00799 
00800   has_data = true;
00801 
00802   if (cvm::debug()) {
00803     cvm::log("Updating system forces for bias "+this->name+"\n");
00804   }
00805 
00806   colvarproxy *proxy = cvm::main()->proxy;
00807 
00808   size_t i;
00809 
00810   if (proxy->total_forces_same_step()) {
00811     for (i = 0; i < num_variables(); i++) {
00812       ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
00813     }
00814   }
00815 
00816   // Collect total colvar forces
00817   if ((cvm::step_relative() > 0) || proxy->total_forces_same_step()) {
00818     if (ti_avg_forces->index_ok(ti_bin)) {
00819       for (i = 0; i < num_variables(); i++) {
00820         if (variables(i)->is_enabled(f_cv_subtract_applied_force)) {
00821           // this colvar is already subtracting all applied forces
00822           ti_system_forces[i] = variables(i)->total_force();
00823         } else {
00824           ti_system_forces[i] = variables(i)->total_force() -
00825             ((subtract_forces != NULL) ?
00826              (*subtract_forces)[i] : previous_colvar_forces[i]);
00827         }
00828       }
00829       if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
00830         ti_avg_forces->acc_value(ti_bin, ti_system_forces);
00831       }
00832     }
00833   }
00834 
00835   if (!proxy->total_forces_same_step()) {
00836     // Set the index for use in the next iteration, when total forces come in
00837     for (i = 0; i < num_variables(); i++) {
00838       ti_bin[i] = ti_avg_forces->current_bin_scalar(i);
00839     }
00840   }
00841 
00842   return COLVARS_OK;
00843 }
00844 
00845 
00846 std::string const colvarbias_ti::get_state_params() const
00847 {
00848   return std::string("");
00849 }
00850 
00851 
00852 int colvarbias_ti::set_state_params(std::string const & /* state_conf */)
00853 {
00854   return COLVARS_OK;
00855 }
00856 
00857 
00858 std::ostream & colvarbias_ti::write_state_data(std::ostream &os)
00859 {
00860   if (! is_enabled(f_cvb_calc_ti_samples)) {
00861     return os;
00862   }
00863   os << "\nhistogram\n";
00864   ti_count->write_raw(os);
00865   os << "\nsystem_forces\n";
00866   ti_avg_forces->write_raw(os);
00867   return os;
00868 }
00869 
00870 
00871 std::istream & colvarbias_ti::read_state_data(std::istream &is)
00872 {
00873   if (! is_enabled(f_cvb_calc_ti_samples)) {
00874     return is;
00875   }
00876   if (cvm::debug()) {
00877     cvm::log("Reading state data for the TI estimator.\n");
00878   }
00879   if (! read_state_data_key(is, "histogram")) {
00880     return is;
00881   }
00882   if (! ti_count->read_raw(is)) {
00883     return is;
00884   }
00885   if (! read_state_data_key(is, "system_forces")) {
00886     return is;
00887   }
00888   if (! ti_avg_forces->read_raw(is)) {
00889     return is;
00890   }
00891   if (cvm::debug()) {
00892     cvm::log("Done reading state data for the TI estimator.\n");
00893   }
00894   return is;
00895 }
00896 
00897 
00898 int colvarbias_ti::write_output_files()
00899 {
00900   int error_code = COLVARS_OK;
00901 
00902   if (!has_data) {
00903     // nothing to write
00904     return COLVARS_OK;
00905   }
00906 
00907   std::string const ti_output_prefix = cvm::output_prefix()+"."+this->name;
00908 
00909   if (is_enabled(f_cvb_write_ti_samples)) {
00910     std::string const ti_count_file_name(ti_output_prefix+".ti.count");
00911     error_code |= ti_count->write_multicol(ti_count_file_name, "TI count file");
00912 
00913     std::string const ti_grad_file_name(ti_output_prefix+".ti.force");
00914     error_code |= ti_avg_forces->write_multicol(ti_grad_file_name, "TI gradient file");
00915   }
00916 
00917   if (is_enabled(f_cvb_write_ti_pmf)) {
00918     std::string const pmf_file_name(ti_output_prefix+".ti.pmf");
00919     cvm::log("Writing TI PMF to file \""+pmf_file_name+"\".\n");
00920     std::ostream &os = cvm::proxy->output_stream(pmf_file_name, "TI PMF");
00921     if (os) {
00922       // get the FE gradient
00923       ti_avg_forces->multiply_constant(-1.0);
00924       ti_avg_forces->write_1D_integral(os);
00925       ti_avg_forces->multiply_constant(-1.0);
00926       cvm::proxy->close_output_stream(pmf_file_name);
00927     } else {
00928       error_code |= COLVARS_FILE_ERROR;
00929     }
00930   }
00931 
00932   return error_code;
00933 }
00934 
00935 
00936 // Static members
00937 
00938 std::vector<colvardeps::feature *> colvarbias::cvb_features;

Generated on Fri Mar 29 02:44:46 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002