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

colvarbias_histogram_reweight_amd.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 "colvarbias_histogram_reweight_amd.h"
00011 #include "colvarproxy.h"
00012 
00013 colvarbias_reweightaMD::colvarbias_reweightaMD(char const *key)
00014   : colvarbias_histogram(key), grid_count(NULL), grid_dV(NULL),
00015     grid_dV_square(NULL), pmf_grid_exp_avg(NULL), pmf_grid_cumulant(NULL),
00016     grad_grid_exp_avg(NULL), grad_grid_cumulant(NULL)
00017 {
00018 }
00019 
00020 colvarbias_reweightaMD::~colvarbias_reweightaMD() {
00021   if (grid_dV) {
00022     delete grid_dV;
00023     grid_dV = NULL;
00024   }
00025   if (grid_dV_square) {
00026     delete grid_dV_square;
00027     grid_dV_square = NULL;
00028   }
00029   if (grid_count) {
00030     delete grid_count;
00031     grid_count = NULL;
00032   }
00033   if (pmf_grid_exp_avg) {
00034     delete pmf_grid_exp_avg;
00035     pmf_grid_exp_avg = NULL;
00036   }
00037   if (pmf_grid_cumulant) {
00038     delete pmf_grid_cumulant;
00039     pmf_grid_cumulant = NULL;
00040   }
00041   if (grad_grid_exp_avg) {
00042     delete grad_grid_exp_avg;
00043     grad_grid_exp_avg = NULL;
00044   }
00045   if (grad_grid_cumulant) {
00046     delete grad_grid_cumulant;
00047     grad_grid_cumulant = NULL;
00048   }
00049 }
00050 
00051 int colvarbias_reweightaMD::init(std::string const &conf) {
00052   if (cvm::proxy->accelMD_enabled() == false) {
00053     cvm::error("Error: accelerated MD in your MD engine is not enabled.\n", COLVARS_INPUT_ERROR);
00054   }
00055   cvm::main()->cite_feature("reweightaMD colvar bias implementation (NAMD)");
00056   int baseclass_init_code = colvarbias_histogram::init(conf);
00057   get_keyval(conf, "CollectAfterSteps", start_after_steps, 0);
00058   get_keyval(conf, "CumulantExpansion", b_use_cumulant_expansion, true);
00059   get_keyval(conf, "WritePMFGradients", b_write_gradients, true);
00060   get_keyval(conf, "historyFreq", history_freq, 0);
00061   b_history_files = (history_freq > 0);
00062   grid_count = new colvar_grid_scalar(colvars);
00063   grid_count->request_actual_value();
00064   grid->request_actual_value();
00065   pmf_grid_exp_avg = new colvar_grid_scalar(colvars);
00066   if (b_write_gradients) {
00067     grad_grid_exp_avg = new colvar_grid_gradient(colvars);
00068   }
00069   if (b_use_cumulant_expansion) {
00070     grid_dV = new colvar_grid_scalar(colvars);
00071     grid_dV_square = new colvar_grid_scalar(colvars);
00072     pmf_grid_cumulant = new colvar_grid_scalar(colvars);
00073     grid_dV->request_actual_value();
00074     grid_dV_square->request_actual_value();
00075     if (b_write_gradients) {
00076       grad_grid_cumulant = new colvar_grid_gradient(colvars);
00077     }
00078   }
00079   previous_bin.assign(num_variables(), -1);
00080   return baseclass_init_code;
00081 }
00082 
00083 int colvarbias_reweightaMD::update() {
00084   colvarproxy *proxy = cvm::main()->proxy;
00085   int error_code = COLVARS_OK;
00086   if (cvm::step_relative() >= start_after_steps) {
00087     // update base class
00088     error_code |= colvarbias::update();
00089 
00090     if (cvm::debug()) {
00091       cvm::log("Updating histogram bias " + this->name);
00092     }
00093 
00094     if (cvm::step_relative() > 0) {
00095       previous_bin = bin;
00096     }
00097 
00098     // assign a valid bin size
00099     bin.assign(num_variables(), 0);
00100 
00101     if (colvar_array_size == 0) {
00102       // update indices for scalar values
00103       size_t i;
00104       for (i = 0; i < num_variables(); i++) {
00105         bin[i] = grid->current_bin_scalar(i);
00106       }
00107 
00108       if (grid->index_ok(previous_bin) && cvm::step_relative() > 0) {
00109         const cvm::real reweighting_factor = cvm::proxy->get_accelMD_factor();
00110         grid_count->acc_value(previous_bin, 1.0);
00111         grid->acc_value(previous_bin, reweighting_factor);
00112         if (b_use_cumulant_expansion) {
00113           const cvm::real dV = cvm::logn(reweighting_factor) *
00114             proxy->target_temperature() * proxy->boltzmann();
00115           grid_dV->acc_value(previous_bin, dV);
00116           grid_dV_square->acc_value(previous_bin, dV * dV);
00117         }
00118       }
00119     } else {
00120       // update indices for vector/array values
00121       size_t iv, i;
00122       for (iv = 0; iv < colvar_array_size; iv++) {
00123         for (i = 0; i < num_variables(); i++) {
00124           bin[i] = grid->current_bin_scalar(i, iv);
00125         }
00126 
00127       if (grid->index_ok(previous_bin) && cvm::step_relative() > 0) {
00128           const cvm::real reweighting_factor = cvm::proxy->get_accelMD_factor();
00129           grid_count->acc_value(previous_bin, 1.0);
00130           grid->acc_value(previous_bin, reweighting_factor);
00131           if (b_use_cumulant_expansion) {
00132             const cvm::real dV = cvm::logn(reweighting_factor) *
00133               proxy->target_temperature() * proxy->boltzmann();
00134             grid_dV->acc_value(previous_bin, dV);
00135             grid_dV_square->acc_value(previous_bin, dV * dV);
00136           }
00137         }
00138       }
00139     }
00140     previous_bin.assign(num_variables(), 0);
00141 
00142     error_code |= cvm::get_error();
00143   }
00144   return error_code;
00145 }
00146 
00147 int colvarbias_reweightaMD::write_output_files() {
00148   int error_code = COLVARS_OK;
00149   // error_code |= colvarbias_histogram::write_output_files();
00150   const std::string out_name_pmf = cvm::output_prefix() + "." +
00151                                    this->name + ".reweight";
00152   error_code |= write_exponential_reweighted_pmf(out_name_pmf);
00153   const std::string out_count_prefix = cvm::output_prefix() + "." +
00154                                        this->name;
00155   error_code |= write_count(out_count_prefix);
00156   const bool write_history = b_history_files &&
00157                              (cvm::step_absolute() % history_freq) == 0;
00158   if (write_history) {
00159     error_code |= write_exponential_reweighted_pmf(
00160       out_name_pmf + ".hist", true);
00161     error_code |= write_count(out_count_prefix + ".hist",
00162                               (cvm::step_relative() > 0));
00163   }
00164   if (b_use_cumulant_expansion) {
00165     const std::string out_name_cumulant_pmf = cvm::output_prefix() + "." +
00166                                               this->name + ".cumulant";
00167     error_code |= write_cumulant_expansion_pmf(out_name_cumulant_pmf);
00168     if (write_history) {
00169       error_code |= write_cumulant_expansion_pmf(
00170         out_name_cumulant_pmf + ".hist", true);
00171     }
00172   }
00173   error_code |= cvm::get_error();
00174   return error_code;
00175 }
00176 
00177 int colvarbias_reweightaMD::write_exponential_reweighted_pmf(
00178   const std::string& p_output_prefix, bool keep_open) {
00179   const std::string output_pmf = p_output_prefix + ".pmf";
00180 
00181   cvm::log("Writing the accelerated MD PMF file \"" + output_pmf + "\".\n");
00182   std::ostream &pmf_grid_os = cvm::proxy->output_stream(output_pmf, "PMF file");
00183   if (!pmf_grid_os) {
00184     return COLVARS_FILE_ERROR;
00185   }
00186   pmf_grid_exp_avg->copy_grid(*grid);
00187   // compute the average
00188   for (size_t i = 0; i < pmf_grid_exp_avg->raw_data_num(); ++i) {
00189     const double count = grid_count->value(i);
00190     if (count > 0) {
00191       const double tmp = pmf_grid_exp_avg->value(i);
00192       pmf_grid_exp_avg->set_value(i, tmp / count);
00193     }
00194   }
00195   hist_to_pmf(pmf_grid_exp_avg, grid_count);
00196   pmf_grid_exp_avg->write_multicol(pmf_grid_os);
00197   if (!keep_open) {
00198     cvm::proxy->close_output_stream(output_pmf);
00199   }
00200 
00201   if (b_write_gradients) {
00202     const std::string output_grad = p_output_prefix + ".grad";
00203     cvm::log("Writing the accelerated MD gradients file \"" + output_grad +
00204              "\".\n");
00205     std::ostream &grad_grid_os = cvm::proxy->output_stream(output_grad, "gradient file");
00206     if (!grad_grid_os) {
00207       return COLVARS_FILE_ERROR;
00208     }
00209     for (std::vector<int> ix = grad_grid_exp_avg->new_index();
00210           grad_grid_exp_avg->index_ok(ix); grad_grid_exp_avg->incr(ix)) {
00211       for (size_t n = 0; n < grad_grid_exp_avg->multiplicity(); n++) {
00212         grad_grid_exp_avg->set_value(
00213           ix, pmf_grid_exp_avg->gradient_finite_diff(ix, n), n);
00214       }
00215     }
00216     grad_grid_exp_avg->write_multicol(grad_grid_os);
00217     if (!keep_open) {
00218       cvm::proxy->close_output_stream(output_grad);
00219     }
00220   }
00221 
00222   return COLVARS_OK;
00223 }
00224 
00225 int colvarbias_reweightaMD::write_cumulant_expansion_pmf(
00226   const std::string& p_output_prefix, bool keep_open) {
00227   const std::string output_pmf = p_output_prefix + ".pmf";
00228   cvm::log("Writing the accelerated MD PMF file using cumulant expansion: \"" + output_pmf + "\".\n");
00229   std::ostream &pmf_grid_cumulant_os = cvm::proxy->output_stream(output_pmf, "PMF file");
00230   if (!pmf_grid_cumulant_os) {
00231     return COLVARS_FILE_ERROR;
00232   }
00233   compute_cumulant_expansion_factor(grid_dV, grid_dV_square,
00234                                     grid_count, pmf_grid_cumulant);
00235   hist_to_pmf(pmf_grid_cumulant, grid_count);
00236   pmf_grid_cumulant->write_multicol(pmf_grid_cumulant_os);
00237   if (!keep_open) {
00238     cvm::proxy->close_output_stream(output_pmf);
00239   }
00240 
00241   if (b_write_gradients) {
00242     const std::string output_grad = p_output_prefix + ".grad";
00243     cvm::log("Writing the accelerated MD gradients file \"" + output_grad + "\".\n");
00244     std::ostream &grad_grid_os = cvm::proxy->output_stream(output_grad, "grad file");
00245     if (!grad_grid_os) {
00246       return COLVARS_FILE_ERROR;
00247     }
00248     for (std::vector<int> ix = grad_grid_cumulant->new_index();
00249           grad_grid_cumulant->index_ok(ix); grad_grid_cumulant->incr(ix)) {
00250       for (size_t n = 0; n < grad_grid_cumulant->multiplicity(); n++) {
00251         grad_grid_cumulant->set_value(
00252           ix, pmf_grid_cumulant->gradient_finite_diff(ix, n), n);
00253       }
00254     }
00255     grad_grid_cumulant->write_multicol(grad_grid_os);
00256     cvm::proxy->close_output_stream(output_grad);
00257   }
00258   return COLVARS_OK;
00259 }
00260 
00261 int colvarbias_reweightaMD::write_count(const std::string& p_output_prefix, bool keep_open) {
00262   const std::string output_name = p_output_prefix + ".count";
00263   cvm::log("Writing the accelerated MD count file \""+output_name+"\".\n");
00264   std::ostream &grid_count_os = cvm::proxy->output_stream(output_name, "count file");
00265   if (!grid_count_os) {
00266     return COLVARS_FILE_ERROR;
00267   }
00268   grid_count->write_multicol(grid_count_os);
00269   if (!keep_open) {
00270     cvm::proxy->close_output_stream(output_name);
00271   }
00272   return COLVARS_OK;
00273 }
00274 
00275 void colvarbias_reweightaMD::hist_to_pmf(
00276   colvar_grid_scalar* hist,
00277   const colvar_grid_scalar* hist_count) const
00278 {
00279   colvarproxy *proxy = cvm::main()->proxy;
00280   if (hist->raw_data_num() == 0) return;
00281   const cvm::real kbt = proxy->boltzmann() * proxy->target_temperature();
00282   bool first_min_element = true;
00283   bool first_max_element = true;
00284   cvm::real min_element = 0;
00285   cvm::real max_element = 0;
00286   size_t i = 0;
00287   // the first loop: using logarithm to compute PMF
00288   for (i = 0; i < hist->raw_data_num(); ++i) {
00289     const cvm::real count = hist_count->value(i);
00290     if (count > 0) {
00291       const cvm::real x = hist->value(i);
00292       const cvm::real pmf_value = -1.0 * kbt * cvm::logn(x);
00293       hist->set_value(i, pmf_value);
00294       // find the minimum PMF value
00295       if (first_min_element) {
00296         // assign current PMF value to min_element at the first time
00297         min_element = pmf_value;
00298         first_min_element = false;
00299       } else {
00300         // if this is not the first time, then
00301         min_element = (pmf_value < min_element) ? pmf_value : min_element;
00302       }
00303       // do the same to the maximum
00304       if (first_max_element) {
00305         max_element = pmf_value;
00306         first_max_element = false;
00307       } else {
00308         max_element = (pmf_value > max_element) ? pmf_value : max_element;
00309       }
00310     }
00311   }
00312   // the second loop: bringing the minimum PMF value to zero
00313   for (i = 0; i < hist->raw_data_num(); ++i) {
00314     const cvm::real count = hist_count->value(i);
00315     if (count > 0) {
00316       // bins that have samples
00317       const cvm::real x = hist->value(i);
00318       hist->set_value(i, x - min_element);
00319     } else {
00320       hist->set_value(i, max_element - min_element);
00321     }
00322   }
00323 }
00324 
00325 
00326 void colvarbias_reweightaMD::compute_cumulant_expansion_factor(
00327   const colvar_grid_scalar* hist_dV,
00328   const colvar_grid_scalar* hist_dV_square,
00329   const colvar_grid_scalar* hist_count,
00330   colvar_grid_scalar* cumulant_expansion_factor) const
00331 {
00332   colvarproxy *proxy = cvm::main()->proxy;
00333   const cvm::real beta = 1.0 / (proxy->boltzmann() * proxy->target_temperature());
00334   size_t i = 0;
00335   for (i = 0; i < hist_dV->raw_data_num(); ++i) {
00336     const cvm::real count = hist_count->value(i);
00337     if (count > 0) {
00338       const cvm::real dV_avg = hist_dV->value(i) / count;
00339       const cvm::real dV_square_avg = hist_dV_square->value(i) / count;
00340       const cvm::real factor = cvm::exp(beta * dV_avg + 0.5 * beta * beta * (dV_square_avg - dV_avg * dV_avg));
00341       cumulant_expansion_factor->set_value(i, factor);
00342     }
00343   }
00344 }
00345 
00346 std::ostream & colvarbias_reweightaMD::write_state_data(std::ostream& os)
00347 {
00348   std::ios::fmtflags flags(os.flags());
00349   os.setf(std::ios::fmtflags(0), std::ios::floatfield);
00350   os << "grid\n";
00351   grid->write_raw(os, 8);
00352   os << "grid_count\n";
00353   grid_count->write_raw(os, 8);
00354   os << "grid_dV\n";
00355   grid_dV->write_raw(os, 8);
00356   os << "grid_dV_square\n";
00357   grid_dV_square->write_raw(os, 8);
00358   os.flags(flags);
00359   return os;
00360 }
00361 
00362 std::istream & colvarbias_reweightaMD::read_state_data(std::istream& is)
00363 {
00364   if (! read_state_data_key(is, "grid")) {
00365     return is;
00366   }
00367   if (! grid->read_raw(is)) {
00368     return is;
00369   }
00370   if (! read_state_data_key(is, "grid_count")) {
00371     return is;
00372   }
00373   if (! grid_count->read_raw(is)) {
00374     return is;
00375   }
00376   if (! read_state_data_key(is, "grid_dV")) {
00377     return is;
00378   }
00379   if (! grid_dV->read_raw(is)) {
00380     return is;
00381   }
00382   if (! read_state_data_key(is, "grid_dV_square")) {
00383     return is;
00384   }
00385   if (! grid_dV_square->read_raw(is)) {
00386     return is;
00387   }
00388   return is;
00389 }

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