Difference for src/colvarbias_abf.C from version 1.19 to 1.20

version 1.19version 1.20
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include "colvarmodule.h" #include "colvarmodule.h"
 #include "colvar.h" #include "colvar.h"
 #include "colvarbias_abf.h" #include "colvarbias_abf.h"
  
 /// ABF bias constructor; parses the config file 
  
 colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) colvarbias_abf::colvarbias_abf(char const *key)
   : colvarbias(conf, key),   : colvarbias(key),
      system_force(NULL),
     gradients(NULL),     gradients(NULL),
     samples(NULL)     samples(NULL),
      z_gradients(NULL),
      z_samples(NULL),
      czar_gradients(NULL),
      last_gradients(NULL),
      last_samples(NULL)
 { {
  }
  
  
  int colvarbias_abf::init(std::string const &conf)
  {
    colvarbias::init(conf);
  
    provide(f_cvb_history_dependent);
  
   // TODO relax this in case of VMD plugin   // TODO relax this in case of VMD plugin
   if (cvm::temperature() == 0.0)   if (cvm::temperature() == 0.0)
     cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");     cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");
  
   // ************* parsing general ABF options ***********************   // ************* parsing general ABF options ***********************
  
   get_keyval(conf, "applyBias",  apply_bias, true);   get_keyval_feature((colvarparse *)this, conf, "applyBias",  f_cvb_apply_force, true);
   if (!apply_bias) cvm::log("WARNING: ABF biases will *not* be applied!\n");   if (!is_enabled(f_cvb_apply_force)){
      cvm::log("WARNING: ABF biases will *not* be applied!\n");
    }
  
   get_keyval(conf, "updateBias",  update_bias, true);   get_keyval(conf, "updateBias",  update_bias, true);
   if (!update_bias) cvm::log("WARNING: ABF biases will *not* be updated!\n");   if (update_bias) {
      enable(f_cvb_history_dependent);
    } else {
      cvm::log("WARNING: ABF biases will *not* be updated!\n");
    }
  
   get_keyval(conf, "hideJacobian", hide_Jacobian, false);   get_keyval(conf, "hideJacobian", hide_Jacobian, false);
   if (hide_Jacobian) {   if (hide_Jacobian) {
Line 58
Line 78
     cvm::error("Error: no collective variables specified for the ABF bias.\n");     cvm::error("Error: no collective variables specified for the ABF bias.\n");
   }   }
  
    if (update_bias) {
    // Request calculation of total force (which also checks for availability)
      if(enable(f_cvb_get_total_force)) return cvm::get_error();
    }
  
    bool b_extended = false;
   for (size_t i = 0; i < colvars.size(); i++) {   for (size_t i = 0; i < colvars.size(); i++) {
  
     if (colvars[i]->value().type() != colvarvalue::type_scalar) {     if (colvars[i]->value().type() != colvarvalue::type_scalar) {
       cvm::error("Error: ABF bias can only use scalar-type variables.\n");       cvm::error("Error: ABF bias can only use scalar-type variables.\n");
     }     }
      colvars[i]->enable(f_cv_grid);
     colvars[i]->enable(colvar::task_gradients); 
  
     if (update_bias) { 
       // Request calculation of system force (which also checks for availability) 
       colvars[i]->enable(colvar::task_system_force); 
  
       if (!colvars[i]->tasks[colvar::task_extended_lagrangian]) { 
         // request computation of Jacobian force 
         // ultimately, will be done regardless of extended Lagrangian 
         // and colvar should then just report zero Jacobian force 
         colvars[i]->enable(colvar::task_Jacobian_force); 
  
         // request Jacobian force as part as system force 
         // except if the user explicitly requires the "silent" Jacobian 
         // correction AND the colvar has a single component 
         if (hide_Jacobian) {         if (hide_Jacobian) {
           if (colvars[i]->n_components() > 1) {       colvars[i]->enable(f_cv_hide_Jacobian);
             cvm::log("WARNING: colvar \"" + colvars[i]->name 
             + "\" has multiple components; reporting its Jacobian forces\n"); 
             colvars[i]->enable(colvar::task_report_Jacobian_force); 
           } 
         } else { 
           colvars[i]->enable(colvar::task_report_Jacobian_force); 
         } 
       } 
     }     }
  
      // If any colvar is extended-system, we need to collect the extended
      // system gradient
      if (colvars[i]->is_enabled(f_cv_extended_Lagrangian))
        b_extended = true;
  
     // Here we could check for orthogonality of the Cartesian coordinates     // Here we could check for orthogonality of the Cartesian coordinates
     // and make it just a warning if some parameter is set?     // and make it just a warning if some parameter is set?
   }   }
Line 111
Line 119
  
   bin.assign(colvars.size(), 0);   bin.assign(colvars.size(), 0);
   force_bin.assign(colvars.size(), 0);   force_bin.assign(colvars.size(), 0);
   force = new cvm::real [colvars.size()];   system_force = new cvm::real [colvars.size()];
  
   // Construct empty grids based on the colvars   // Construct empty grids based on the colvars
   if (cvm::debug()) {   if (cvm::debug()) {
Line 123
Line 131
   gradients->samples = samples;   gradients->samples = samples;
   samples->has_parent_data = true;   samples->has_parent_data = true;
  
    // Data for eABF z-based estimator
    if (b_extended) {
      // CZAR output files for stratified eABF
      get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
                 colvarparse::parse_silent);
  
      z_bin.assign(colvars.size(), 0);
      z_samples   = new colvar_grid_count(colvars);
      z_samples->request_actual_value();
      z_gradients = new colvar_grid_gradient(colvars);
      z_gradients->request_actual_value();
      z_gradients->samples = z_samples;
      z_samples->has_parent_data = true;
      czar_gradients = new colvar_grid_gradient(colvars);
    }
  
   // For shared ABF, we store a second set of grids.   // For shared ABF, we store a second set of grids.
   // This used to be only if "shared" was defined,   // This used to be only if "shared" was defined,
   // but now we allow calling share externally (e.g. from Tcl).   // but now we allow calling share externally (e.g. from Tcl).
Line 138
Line 162
   }   }
  
   cvm::log("Finished ABF setup.\n");   cvm::log("Finished ABF setup.\n");
  
    return COLVARS_OK;
 } }
  
 /// Destructor /// Destructor
Line 153
Line 179
     gradients = NULL;     gradients = NULL;
   }   }
  
    if (z_samples) {
      delete z_samples;
      z_samples = NULL;
    }
  
    if (z_gradients) {
      delete z_gradients;
      z_gradients = NULL;
    }
  
    if (czar_gradients) {
      delete czar_gradients;
      czar_gradients = NULL;
    }
  
   // shared ABF   // shared ABF
   // We used to only do this if "shared" was defined,   // We used to only do this if "shared" was defined,
   // but now we can call shared externally   // but now we can call shared externally
Line 166
Line 207
     last_gradients = NULL;     last_gradients = NULL;
   }   }
  
   delete [] force;   if (system_force) {
      delete [] system_force;
      system_force = NULL;
    }
  
   if (cvm::n_abf_biases > 0)   if (cvm::n_abf_biases > 0)
     cvm::n_abf_biases -= 1;     cvm::n_abf_biases -= 1;
Line 176
Line 220
 /// Update the FE gradient, compute and apply biasing force /// Update the FE gradient, compute and apply biasing force
 /// also output data to disk if needed /// also output data to disk if needed
  
 cvm::real colvarbias_abf::update() int colvarbias_abf::update()
 { {
   if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);   if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
  
Line 199
Line 243
     if ( update_bias && samples->index_ok(force_bin) ) {     if ( update_bias && samples->index_ok(force_bin) ) {
       // Only if requested and within bounds of the grid...       // Only if requested and within bounds of the grid...
  
       for (size_t i=0; i<colvars.size(); i++) {   // get forces(lagging by 1 timestep) from colvars       for (size_t i = 0; i < colvars.size(); i++) {
         force[i] = colvars[i]->system_force();         // get total forces (lagging by 1 timestep) from colvars
          // and subtract previous ABF force
          update_system_force(i);
        }
        gradients->acc_force(force_bin, system_force);
      }
      if ( z_gradients && update_bias ) {
        for (size_t i = 0; i < colvars.size(); i++) {
          z_bin[i] = z_samples->current_bin_scalar(i);
        }
        if ( z_samples->index_ok(z_bin) ) {
          for (size_t i = 0; i < colvars.size(); i++) {
            // If we are outside the range of xi, the force has not been obtained above
            // the function is just an accessor, so cheap to call again anyway
            update_system_force(i);
          }
          z_gradients->acc_force(z_bin, system_force);
       }       }
       gradients->acc_force(force_bin, force); 
     }     }
   }   }
  
Line 215
Line 274
   }   }
  
   // Compute and apply the new bias, if applicable   // Compute and apply the new bias, if applicable
   if ( apply_bias && samples->index_ok(bin) ) {   if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) {
  
     size_t  count = samples->value(bin);     size_t  count = samples->value(bin);
     cvm::real fact = 1.0;     cvm::real fact = 1.0;
Line 263
Line 322
   }   }
  
   if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {   if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
     cvm::log("ABFHISTORYFILE "+cvm::to_str(cvm::step_absolute())); 
     // file already exists iff cvm::step_relative() > 0     // file already exists iff cvm::step_relative() > 0
     // otherwise, backup and replace     // otherwise, backup and replace
     write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0));     write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0));
Line 283
Line 341
     cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");     cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
   }   }
  
   return 0.0;   return COLVARS_OK;
 } }
  
  
 int colvarbias_abf::replica_share() { int colvarbias_abf::replica_share() {
   int p;   int p;
  
Line 399
Line 458
     pmf_os << std::endl;     pmf_os << std::endl;
     pmf_os.close();     pmf_os.close();
   }   }
  
    if (z_gradients) {
      // Write eABF-related quantities
  
      std::string  z_samples_out_name = prefix + ".zcount";
      cvm::ofstream z_samples_os;
  
      if (!append) cvm::backup_file(z_samples_out_name.c_str());
      z_samples_os.open(z_samples_out_name.c_str(), mode);
      if (!z_samples_os.is_open()) {
        cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
      }
      z_samples->write_multicol(z_samples_os);
      z_samples_os.close();
  
      if (b_czar_window_file) {
        std::string  z_gradients_out_name = prefix + ".zgrad";
        cvm::ofstream z_gradients_os;
  
        if (!append) cvm::backup_file(z_gradients_out_name.c_str());
        z_gradients_os.open(z_gradients_out_name.c_str(), mode);
        if (!z_gradients_os.is_open()) {
          cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
        }
        z_gradients->write_multicol(z_gradients_os);
        z_gradients_os.close();
      }
  
      // Calculate CZAR estimator of gradients
      for (std::vector<int> ix = czar_gradients->new_index();
            czar_gradients->index_ok(ix); czar_gradients->incr(ix)) {
        for (size_t n = 0; n < czar_gradients->multiplicity(); n++) {
          czar_gradients->set_value(ix, z_gradients->value_output(ix, n)
              - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n),
              n);
        }
      }
  
      std::string  czar_gradients_out_name = prefix + ".czar.grad";
      cvm::ofstream czar_gradients_os;
  
      if (!append) cvm::backup_file(czar_gradients_out_name.c_str());
      czar_gradients_os.open(czar_gradients_out_name.c_str(), mode);
      if (!czar_gradients_os.is_open()) {
        cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
      }
      czar_gradients->write_multicol(czar_gradients_os);
      czar_gradients_os.close();
  
      if (colvars.size() == 1) {
        std::string  czar_pmf_out_name = prefix + ".czar.pmf";
        if (!append) cvm::backup_file(czar_pmf_out_name.c_str());
        cvm::ofstream czar_pmf_os;
        // Do numerical integration and output a PMF
        czar_pmf_os.open(czar_pmf_out_name.c_str(), mode);
        if (!czar_pmf_os.is_open())  cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
        czar_gradients->write_1D_integral(czar_pmf_os);
        czar_pmf_os << std::endl;
        czar_pmf_os.close();
      }
    }
   return;   return;
 } }
  
Line 425
Line 545
  
 void colvarbias_abf::read_gradients_samples() void colvarbias_abf::read_gradients_samples()
 { {
   std::string samples_in_name, gradients_in_name;   std::string samples_in_name, gradients_in_name, z_samples_in_name, z_gradients_in_name;
  
   for ( size_t i = 0; i < input_prefix.size(); i++ ) {   for ( size_t i = 0; i < input_prefix.size(); i++ ) {
     samples_in_name = input_prefix[i] + ".count";     samples_in_name = input_prefix[i] + ".count";
     gradients_in_name = input_prefix[i] + ".grad";     gradients_in_name = input_prefix[i] + ".grad";
      z_samples_in_name = input_prefix[i] + ".zcount";
      z_gradients_in_name = input_prefix[i] + ".zgrad";
     // For user-provided files, the per-bias naming scheme may not apply     // For user-provided files, the per-bias naming scheme may not apply
  
     std::ifstream is;     std::ifstream is;
  
     cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);     cvm::log("Reading sample count from " + samples_in_name + " and gradient from " + gradients_in_name);
     is.open(samples_in_name.c_str());     is.open(samples_in_name.c_str());
     if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");     if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
     samples->read_multicol(is, true);     samples->read_multicol(is, true);
Line 445
Line 567
     if (!is.is_open()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading");     if (!is.is_open()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading");
     gradients->read_multicol(is, true);     gradients->read_multicol(is, true);
     is.close();     is.close();
  
      if (z_gradients) {
        // Read eABF z-averaged data for CZAR
        cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);
  
        is.clear();
        is.open(z_samples_in_name.c_str());
        if (!is.is_open())  cvm::error("Error opening eABF z-histogram file " + z_samples_in_name + " for reading");
        z_samples->read_multicol(is, true);
        is.close();
        is.clear();
  
        is.open(z_gradients_in_name.c_str());
        if (!is.is_open())  cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading");
        z_gradients->read_multicol(is, true);
        is.close();
      }
   }   }
   return;   return;
 } }
Line 454
Line 593
 { {
  
   std::ios::fmtflags flags(os.flags());   std::ios::fmtflags flags(os.flags());
   os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format 
  
   os << "abf {\n"   os << "abf {\n"
      << "  configuration {\n"      << "  configuration {\n"
      << "    name " << this->name << "\n";      << "    name " << this->name << "\n";
   os << "  }\n";   os << "  }\n";
  
   os << "samples\n";   os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
    os << "\nsamples\n";
   samples->write_raw(os, 8);   samples->write_raw(os, 8);
    os.flags(flags);
  
   os << "\ngradient\n";   os << "\ngradient\n";
   gradients->write_raw(os);   gradients->write_raw(os, 8);
  
    if (z_gradients) {
      os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
      os << "\nz_samples\n";
      z_samples->write_raw(os, 8);
      os.flags(flags);
      os << "\nz_gradient\n";
      z_gradients->write_raw(os, 8);
    }
  
   os << "}\n\n";   os << "}\n\n";
  
Line 477
Line 626
 std::istream & colvarbias_abf::read_restart(std::istream& is) std::istream & colvarbias_abf::read_restart(std::istream& is)
 { {
   if ( input_prefix.size() > 0 ) {   if ( input_prefix.size() > 0 ) {
     cvm::error("ERROR: cannot provide both inputPrefix and restart information(colvarsInput)");     cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", INPUT_ERROR);
   }   }
  
   size_t const start_pos = is.tellg();   size_t const start_pos = is.tellg();
Line 539
Line 688
     return is;     return is;
   }   }
  
    if (z_gradients) {
      if ( !(is >> key) || !(key == "z_samples")) {
        cvm::log("Error: in reading restart configuration for ABF bias \""+
                  this->name+"\" at position "+
                  cvm::to_str(is.tellg())+" in stream.\n");
        is.clear();
        is.seekg(start_pos, std::ios::beg);
        is.setstate(std::ios::failbit);
        return is;
      }
      if (! z_samples->read_raw(is)) {
        is.clear();
        is.seekg(start_pos, std::ios::beg);
        is.setstate(std::ios::failbit);
        return is;
      }
  
      if ( !(is >> key)   || !(key == "z_gradient")) {
        cvm::log("Error: in reading restart configuration for ABF bias \""+
                  this->name+"\" at position "+
                  cvm::to_str(is.tellg())+" in stream.\n");
        is.clear();
        is.seekg(start_pos, std::ios::beg);
        is.setstate(std::ios::failbit);
        return is;
      }
      if (! z_gradients->read_raw(is)) {
        is.clear();
        is.seekg(start_pos, std::ios::beg);
        is.setstate(std::ios::failbit);
        return is;
      }
    }
   is >> brace;   is >> brace;
   if (brace != "}") {   if (brace != "}") {
     cvm::error("Error: corrupt restart information for ABF bias \""+     cvm::error("Error: corrupt restart information for ABF bias \""+

Legend:
Removed in v.1.19 
changed lines
 Added in v.1.20


Made by using version 1.53 of cvs2html