colvarbias_histogram_reweight_amd.C

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 }