00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <iostream>
00011
00012 #include "colvarmodule.h"
00013 #include "colvar.h"
00014 #include "colvarbias_abf.h"
00015
00016
00017 colvarbias_abf::colvarbias_abf(char const *key)
00018 : colvarbias(key),
00019 b_UI_estimator(false),
00020 b_CZAR_estimator(false),
00021 pabf_freq(0),
00022 system_force(NULL),
00023 gradients(NULL),
00024 samples(NULL),
00025 pmf(NULL),
00026 z_gradients(NULL),
00027 z_samples(NULL),
00028 czar_gradients(NULL),
00029 czar_pmf(NULL),
00030 last_gradients(NULL),
00031 last_samples(NULL)
00032 {
00033 colvarproxy *proxy = cvm::main()->proxy;
00034 if (!proxy->total_forces_same_step()) {
00035
00036 feature_states[f_cvb_step_zero_data].available = false;
00037 }
00038 }
00039
00040
00041 int colvarbias_abf::init(std::string const &conf)
00042 {
00043 colvarproxy *proxy = cvm::main()->proxy;
00044
00045 colvarbias::init(conf);
00046 cvm::main()->cite_feature("ABF colvar bias implementation");
00047
00048 enable(f_cvb_scalar_variables);
00049 enable(f_cvb_calc_pmf);
00050
00051 if ((proxy->target_temperature() == 0.0) && proxy->simulation_running()) {
00052 cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");
00053 }
00054
00055
00056
00057 get_keyval_feature((colvarparse *)this, conf, "applyBias", f_cvb_apply_force, true);
00058 if (!is_enabled(f_cvb_apply_force)){
00059 cvm::log("WARNING: ABF biases will *not* be applied!\n");
00060 }
00061
00062 get_keyval(conf, "updateBias", update_bias, true);
00063 if (update_bias) {
00064 enable(f_cvb_history_dependent);
00065 } else {
00066 cvm::log("WARNING: ABF biases will *not* be updated!\n");
00067 }
00068
00069 get_keyval(conf, "hideJacobian", hide_Jacobian, false);
00070 if (hide_Jacobian) {
00071 cvm::log("Jacobian (geometric) forces will be handled internally.\n");
00072 } else {
00073 cvm::log("Jacobian (geometric) forces will be included in reported free energy gradients.\n");
00074 }
00075
00076 get_keyval(conf, "fullSamples", full_samples, 200);
00077 if ( full_samples <= 1 ) full_samples = 1;
00078 min_samples = full_samples / 2;
00079
00080
00081 get_keyval(conf, "inputPrefix", input_prefix, std::vector<std::string>());
00082
00083 get_keyval(conf, "historyFreq", history_freq, 0);
00084 if (history_freq != 0) {
00085 if (output_freq == 0) {
00086 cvm::error("Error: historyFreq must be a multiple of outputFreq.\n",
00087 COLVARS_INPUT_ERROR);
00088 } else {
00089 if ((history_freq % output_freq) != 0) {
00090 cvm::error("Error: historyFreq must be a multiple of outputFreq.\n",
00091 COLVARS_INPUT_ERROR);
00092 }
00093 }
00094 }
00095 b_history_files = (history_freq > 0);
00096
00097
00098 get_keyval(conf, "shared", shared_on, false);
00099 if (shared_on) {
00100 cvm::main()->cite_feature("Multiple-walker ABF implementation");
00101 if ((proxy->replica_enabled() != COLVARS_OK) ||
00102 (proxy->num_replicas() <= 1)) {
00103 return cvm::error("Error: shared ABF requires more than one replica.",
00104 COLVARS_INPUT_ERROR);
00105 }
00106 cvm::log("shared ABF will be applied among "+
00107 cvm::to_str(proxy->num_replicas()) + " replicas.\n");
00108 if (cvm::proxy->smp_enabled() == COLVARS_OK) {
00109 cvm::error("Error: shared ABF is currently not available with SMP parallelism; "
00110 "please set \"SMP off\" at the top of the Colvars configuration file.\n",
00111 COLVARS_NOT_IMPLEMENTED);
00112 return COLVARS_NOT_IMPLEMENTED;
00113 }
00114
00115
00116 get_keyval(conf, "sharedFreq", shared_freq, output_freq);
00117 }
00118
00119
00120
00121 if (num_variables() == 0) {
00122 cvm::error("Error: no collective variables specified for the ABF bias.\n");
00123 return COLVARS_ERROR;
00124 }
00125
00126 if (update_bias) {
00127
00128 if(enable(f_cvb_get_total_force)) return cvm::get_error();
00129 }
00130
00131 bool b_extended = false;
00132 size_t i;
00133 for (i = 0; i < num_variables(); i++) {
00134
00135 if (colvars[i]->value().type() != colvarvalue::type_scalar) {
00136 cvm::error("Error: ABF bias can only use scalar-type variables.\n");
00137 }
00138 colvars[i]->enable(f_cv_grid);
00139 if (hide_Jacobian) {
00140 colvars[i]->enable(f_cv_hide_Jacobian);
00141 }
00142
00143
00144 if (colvars[i]->is_enabled(f_cv_extended_Lagrangian)
00145 && !colvars[i]->is_enabled(f_cv_external)) {
00146 b_extended = true;
00147 }
00148
00149
00150
00151
00152 if (colvars[i]->get_time_step_factor() != time_step_factor) {
00153 cvm::error("Error: " + colvars[i]->description + " has a value of timeStepFactor ("
00154 + cvm::to_str(colvars[i]->get_time_step_factor()) + ") different from that of "
00155 + description + " (" + cvm::to_str(time_step_factor) + ").\n");
00156 return COLVARS_ERROR;
00157 }
00158
00159
00160
00161 }
00162
00163 if (b_extended) {
00164 cvm::main()->cite_feature("eABF implementation");
00165 } else {
00166 cvm::main()->cite_feature("Internal-forces free energy estimator");
00167 }
00168
00169 if (get_keyval(conf, "maxForce", max_force)) {
00170 if (max_force.size() != num_variables()) {
00171 cvm::error("Error: Number of parameters to maxForce does not match number of colvars.");
00172 }
00173 for (i = 0; i < num_variables(); i++) {
00174 if (max_force[i] < 0.0) {
00175 cvm::error("Error: maxForce should be non-negative.");
00176 return COLVARS_ERROR;
00177 }
00178 }
00179 cap_force = true;
00180 } else {
00181 cap_force = false;
00182 }
00183
00184 bin.assign(num_variables(), 0);
00185 force_bin.assign(num_variables(), 0);
00186 system_force = new cvm::real [num_variables()];
00187
00188
00189 if (cvm::debug()) {
00190 cvm::log("Allocating count and free energy gradient grids.\n");
00191 }
00192
00193 samples = new colvar_grid_count(colvars);
00194 gradients = new colvar_grid_gradient(colvars);
00195 gradients->samples = samples;
00196 samples->has_parent_data = true;
00197
00198
00199 if ( b_extended ) {
00200 get_keyval(conf, "CZARestimator", b_CZAR_estimator, true);
00201 if ( b_CZAR_estimator ) {
00202 cvm::main()->cite_feature("CZAR eABF estimator");
00203 }
00204
00205 get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
00206 colvarparse::parse_silent);
00207
00208 z_bin.assign(num_variables(), 0);
00209 z_samples = new colvar_grid_count(colvars);
00210 z_samples->request_actual_value();
00211 z_gradients = new colvar_grid_gradient(colvars);
00212 z_gradients->request_actual_value();
00213 z_gradients->samples = z_samples;
00214 z_samples->has_parent_data = true;
00215 czar_gradients = new colvar_grid_gradient(colvars);
00216 }
00217
00218 get_keyval(conf, "integrate", b_integrate, num_variables() <= 3);
00219 if (b_integrate) {
00220
00221 if ( num_variables() > 3 ) {
00222 cvm::error("Error: cannot integrate free energy in dimension > 3.\n");
00223 return COLVARS_ERROR;
00224 }
00225 pmf = new integrate_potential(colvars, gradients);
00226 if ( b_CZAR_estimator ) {
00227 czar_pmf = new integrate_potential(colvars, czar_gradients);
00228 }
00229
00230 get_keyval(conf, "integrateMaxIterations", integrate_iterations, 10000, colvarparse::parse_silent);
00231 get_keyval(conf, "integrateTol", integrate_tol, 1e-6, colvarparse::parse_silent);
00232
00233 get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0, colvarparse::parse_silent);
00234 get_keyval(conf, "pABFintegrateMaxIterations", pabf_integrate_iterations, 100, colvarparse::parse_silent);
00235 get_keyval(conf, "pABFintegrateTol", pabf_integrate_tol, 1e-4, colvarparse::parse_silent);
00236 }
00237
00238
00239
00240
00241 last_samples = new colvar_grid_count(colvars);
00242 last_gradients = new colvar_grid_gradient(colvars);
00243 last_gradients->samples = last_samples;
00244 last_samples->has_parent_data = true;
00245 shared_last_step = -1;
00246
00247
00248 if ( input_prefix.size() > 0 ) {
00249 read_gradients_samples();
00250
00251 pmf->set_div();
00252 }
00253
00254
00255 if (b_extended) {
00256 get_keyval(conf, "UIestimator", b_UI_estimator, false);
00257
00258 if (b_UI_estimator) {
00259
00260 cvm::main()->cite_feature("Umbrella-integration eABF estimator");
00261 std::vector<double> UI_lowerboundary;
00262 std::vector<double> UI_upperboundary;
00263 std::vector<double> UI_width;
00264 std::vector<double> UI_krestr;
00265
00266 bool UI_restart = (input_prefix.size() > 0);
00267
00268 for (i = 0; i < num_variables(); i++) {
00269 UI_lowerboundary.push_back(colvars[i]->lower_boundary);
00270 UI_upperboundary.push_back(colvars[i]->upper_boundary);
00271 UI_width.push_back(colvars[i]->width);
00272 UI_krestr.push_back(colvars[i]->force_constant());
00273 }
00274 eabf_UI = UIestimator::UIestimator(UI_lowerboundary,
00275 UI_upperboundary,
00276 UI_width,
00277 UI_krestr,
00278 output_prefix,
00279 cvm::restart_out_freq,
00280 UI_restart,
00281 input_prefix,
00282 proxy->target_temperature());
00283 }
00284 }
00285
00286 cvm::log("Finished ABF setup.\n");
00287 return COLVARS_OK;
00288 }
00289
00291 colvarbias_abf::~colvarbias_abf()
00292 {
00293 if (samples) {
00294 delete samples;
00295 samples = NULL;
00296 }
00297
00298 if (gradients) {
00299 delete gradients;
00300 gradients = NULL;
00301 }
00302
00303 if (pmf) {
00304 delete pmf;
00305 pmf = NULL;
00306 }
00307
00308 if (z_samples) {
00309 delete z_samples;
00310 z_samples = NULL;
00311 }
00312
00313 if (z_gradients) {
00314 delete z_gradients;
00315 z_gradients = NULL;
00316 }
00317
00318 if (czar_gradients) {
00319 delete czar_gradients;
00320 czar_gradients = NULL;
00321 }
00322
00323 if (czar_pmf) {
00324 delete czar_pmf;
00325 czar_pmf = NULL;
00326 }
00327
00328
00329
00330
00331 if (last_samples) {
00332 delete last_samples;
00333 last_samples = NULL;
00334 }
00335
00336 if (last_gradients) {
00337 delete last_gradients;
00338 last_gradients = NULL;
00339 }
00340
00341 if (system_force) {
00342 delete [] system_force;
00343 system_force = NULL;
00344 }
00345 }
00346
00347
00350
00351 int colvarbias_abf::update()
00352 {
00353 if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
00354
00355 size_t i;
00356 for (i = 0; i < num_variables(); i++) {
00357 bin[i] = samples->current_bin_scalar(i);
00358 }
00359 if (cvm::proxy->total_forces_same_step()) {
00360
00361 force_bin = bin;
00362 }
00363
00364 if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
00365
00366 if (update_bias) {
00367
00368
00369
00370
00371
00372 if (samples->index_ok(force_bin)) {
00373
00374
00375 for (i = 0; i < num_variables(); i++) {
00376
00377
00378 update_system_force(i);
00379 }
00380 gradients->acc_force(force_bin, system_force);
00381 if ( b_integrate ) {
00382 pmf->update_div_neighbors(force_bin);
00383 }
00384 }
00385 }
00386
00387 if ( z_gradients && update_bias ) {
00388 for (i = 0; i < num_variables(); i++) {
00389 z_bin[i] = z_samples->current_bin_scalar(i);
00390 }
00391 if ( z_samples->index_ok(z_bin) ) {
00392 for (i = 0; i < num_variables(); i++) {
00393
00394
00395 update_system_force(i);
00396 }
00397 z_gradients->acc_force(z_bin, system_force);
00398 }
00399 }
00400
00401 if ( b_integrate ) {
00402 if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
00403 cvm::real err;
00404 int iter = pmf->integrate(pabf_integrate_iterations, pabf_integrate_tol, err);
00405 if ( iter == pabf_integrate_iterations ) {
00406 cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(pabf_integrate_tol)
00407 + " in " + cvm::to_str(pabf_integrate_iterations)
00408 + " steps. Residual error: " + cvm::to_str(err));
00409 }
00410 pmf->set_zero_minimum();
00411 }
00412 }
00413 }
00414
00415 if (!cvm::proxy->total_forces_same_step()) {
00416
00417
00418 force_bin = bin;
00419 }
00420
00421
00422 for (i = 0; i < num_variables(); i++) {
00423 colvar_forces[i].reset();
00424 }
00425
00426
00427 if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) {
00428
00429 cvm::real count = cvm::real(samples->value(bin));
00430 cvm::real fact = 1.0;
00431
00432
00433 if ( count < full_samples ) {
00434 fact = (count < min_samples) ? 0.0 :
00435 (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
00436 }
00437
00438 std::vector<cvm::real> grad(num_variables());
00439
00440 if ( pabf_freq ) {
00441
00442 pmf->vector_gradient_finite_diff(bin, grad);
00443 } else {
00444
00445 gradients->vector_value(bin, grad);
00446 }
00447
00448
00449
00450
00451
00452 if ( fact != 0.0 ) {
00453 if ( (num_variables() == 1) && colvars[0]->periodic_boundaries() ) {
00454
00455
00456
00457 colvar_forces[0].real_value = fact * (grad[0] - gradients->average ());
00458 } else {
00459 for (i = 0; i < num_variables(); i++) {
00460
00461 colvar_forces[i].real_value = fact * grad[i];
00462 }
00463 }
00464 if (cap_force) {
00465 for (i = 0; i < num_variables(); i++) {
00466 if ( colvar_forces[i].real_value * colvar_forces[i].real_value > max_force[i] * max_force[i] ) {
00467 colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]);
00468 }
00469 }
00470 }
00471 }
00472 }
00473
00474
00475 if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
00476
00477 output_prefix = cvm::output_prefix();
00478 } else {
00479 output_prefix = cvm::output_prefix() + "." + this->name;
00480 }
00481
00482 if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) {
00483
00484 replica_share();
00485 }
00486
00487
00488 if (shared_last_step < 0) {
00489
00490 last_gradients->copy_grid(*gradients);
00491 last_samples->copy_grid(*samples);
00492 shared_last_step = cvm::step_absolute();
00493 cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".\n");
00494 }
00495
00496
00497 if (b_UI_estimator)
00498 {
00499 std::vector<double> x(num_variables(),0);
00500 std::vector<double> y(num_variables(),0);
00501 for (i = 0; i < num_variables(); i++)
00502 {
00503 x[i] = colvars[i]->actual_value();
00504 y[i] = colvars[i]->value();
00505 }
00506 eabf_UI.update_output_filename(output_prefix);
00507 eabf_UI.update(cvm::step_absolute(), x, y);
00508 }
00509
00511 int error_code = calc_energy(NULL);
00512
00513 return error_code;
00514 }
00515
00516
00517 int colvarbias_abf::replica_share() {
00518
00519 colvarproxy *proxy = cvm::main()->proxy;
00520
00521 if (proxy->replica_enabled() != COLVARS_OK) {
00522 cvm::error("Error: shared ABF: No replicas.\n");
00523 return COLVARS_ERROR;
00524 }
00525
00526 if (shared_last_step < 0 ) {
00527 cvm::error("Error: shared ABF: Tried to apply shared ABF before any sampling had occurred.\n");
00528 return COLVARS_ERROR;
00529 }
00530
00531
00532 cvm::log("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) );
00533
00534
00535 size_t data_n = gradients->raw_data_num();
00536 size_t samp_start = data_n*sizeof(cvm::real);
00537 size_t msg_total = data_n*sizeof(size_t) + samp_start;
00538 char* msg_data = new char[msg_total];
00539
00540 if (proxy->replica_index() == 0) {
00541 int p;
00542
00543 for (p = 1; p < proxy->num_replicas(); p++) {
00544
00545 proxy->replica_comm_recv(msg_data, msg_total, p);
00546
00547
00548 last_gradients->raw_data_in((cvm::real*)(&msg_data[0]));
00549 last_samples->raw_data_in((size_t*)(&msg_data[samp_start]));
00550
00551
00552
00553 gradients->add_grid( *last_gradients );
00554 samples->add_grid( *last_samples );
00555 }
00556
00557
00558 gradients->raw_data_out((cvm::real*)(&msg_data[0]));
00559 samples->raw_data_out((size_t*)(&msg_data[samp_start]));
00560 for (p = 1; p < proxy->num_replicas(); p++) {
00561 proxy->replica_comm_send(msg_data, msg_total, p);
00562 }
00563
00564 } else {
00565
00566
00567 last_gradients->delta_grid(*gradients);
00568 last_samples->delta_grid(*samples);
00569
00570
00571 last_gradients->raw_data_out((cvm::real*)(&msg_data[0]));
00572 last_samples->raw_data_out((size_t*)(&msg_data[samp_start]));
00573 proxy->replica_comm_send(msg_data, msg_total, 0);
00574
00575
00576 proxy->replica_comm_recv(msg_data, msg_total, 0);
00577
00578 gradients->raw_data_in((cvm::real*)(&msg_data[0]));
00579 samples->raw_data_in((size_t*)(&msg_data[samp_start]));
00580 }
00581
00582
00583
00584 proxy->replica_comm_barrier();
00585
00586 delete[] msg_data;
00587
00588
00589 last_gradients->copy_grid(*gradients);
00590 last_samples->copy_grid(*samples);
00591 shared_last_step = cvm::step_absolute();
00592
00593 if (b_integrate) {
00594
00595 pmf->set_div();
00596 }
00597 return COLVARS_OK;
00598 }
00599
00600
00601 template <class T> int colvarbias_abf::write_grid_to_file(T const *grid,
00602 std::string const &filename,
00603 bool close) {
00604 std::ostream &os = cvm::proxy->output_stream(filename);
00605 if (!os) {
00606 return cvm::error("Error opening file " + filename + " for writing.\n", COLVARS_ERROR | COLVARS_FILE_ERROR);
00607 }
00608 grid->write_multicol(os);
00609 if (close) {
00610 cvm::proxy->close_output_stream(filename);
00611 } else {
00612
00613 os << std::endl;
00614 cvm::proxy->flush_output_stream(filename);
00615 }
00616
00617
00618
00619
00620 if (num_variables() > 2 && close) {
00621 std::string dx = filename + ".dx";
00622 std::ostream &dx_os = cvm::proxy->output_stream(dx);
00623 if (!dx_os) {
00624 return cvm::error("Error opening file " + dx + " for writing.\n", COLVARS_ERROR | COLVARS_FILE_ERROR);
00625 }
00626 grid->write_opendx(dx_os);
00627
00628 cvm::proxy->close_output_stream(dx);
00629
00630
00631
00632
00633
00634
00635 }
00636 return COLVARS_OK;
00637 }
00638
00639
00640 void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool close)
00641 {
00642 colvarproxy *proxy = cvm::main()->proxy;
00643
00644 write_grid_to_file<colvar_grid_count>(samples, prefix + ".count", close);
00645 write_grid_to_file<colvar_grid_gradient>(gradients, prefix + ".grad", close);
00646
00647 if (b_integrate) {
00648
00649 cvm::real err;
00650 pmf->integrate(integrate_iterations, integrate_tol, err);
00651 pmf->set_zero_minimum();
00652 write_grid_to_file<colvar_grid_scalar>(pmf, prefix + ".pmf", close);
00653 }
00654
00655 if (b_CZAR_estimator) {
00656
00657 write_grid_to_file<colvar_grid_count>(z_samples, prefix + ".zcount", close);
00658 if (b_czar_window_file) {
00659 write_grid_to_file<colvar_grid_gradient>(z_gradients, prefix + ".zgrad", close);
00660 }
00661
00662
00663 for (std::vector<int> ix = czar_gradients->new_index();
00664 czar_gradients->index_ok(ix); czar_gradients->incr(ix)) {
00665 for (size_t n = 0; n < czar_gradients->multiplicity(); n++) {
00666 czar_gradients->set_value(ix, z_gradients->value_output(ix, n)
00667 - proxy->target_temperature() * proxy->boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n);
00668 }
00669 }
00670 write_grid_to_file<colvar_grid_gradient>(czar_gradients, prefix + ".czar.grad", close);
00671
00672 if (b_integrate) {
00673
00674 cvm::real err;
00675 czar_pmf->set_div();
00676 czar_pmf->integrate(integrate_iterations, integrate_tol, err);
00677 czar_pmf->set_zero_minimum();
00678 write_grid_to_file<colvar_grid_scalar>(czar_pmf, prefix + ".czar.pmf", close);
00679 }
00680 }
00681 return;
00682 }
00683
00684
00685
00687
00688 return samples->number_of_points(0);
00689 }
00691 int colvarbias_abf::current_bin() {
00692 return samples->current_bin_scalar(0);
00693 }
00695 int colvarbias_abf::bin_count(int bin_index) {
00696 if (bin_index < 0 || bin_index >= bin_num()) {
00697 cvm::error("Error: Tried to get bin count from invalid bin index "+cvm::to_str(bin_index));
00698 return -1;
00699 }
00700 std::vector<int> ix(1,(int)bin_index);
00701 return samples->value(ix);
00702 }
00703
00704
00705 int colvarbias_abf::read_gradients_samples()
00706 {
00707 int error_code = COLVARS_OK;
00708
00709 std::string samples_in_name, gradients_in_name, z_samples_in_name, z_gradients_in_name;
00710
00711 for ( size_t i = 0; i < input_prefix.size(); i++ ) {
00712 samples_in_name = input_prefix[i] + ".count";
00713 gradients_in_name = input_prefix[i] + ".grad";
00714 z_samples_in_name = input_prefix[i] + ".zcount";
00715 z_gradients_in_name = input_prefix[i] + ".zgrad";
00716
00717 cvm::log("Reading sample count from " + samples_in_name +
00718 " and gradient from " + gradients_in_name);
00719
00720 error_code |= samples->read_multicol(samples_in_name,
00721 "ABF samples file",
00722 true);
00723
00724 error_code |= gradients->read_multicol(gradients_in_name,
00725 "ABF gradient file",
00726 true);
00727
00728 if (b_CZAR_estimator) {
00729
00730 cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);
00731 error_code |= z_samples->read_multicol(z_samples_in_name,
00732 "eABF z-histogram file",
00733 true);
00734 error_code |= z_gradients->read_multicol(z_gradients_in_name,
00735 "eABF z-gradient file",
00736 true);
00737 }
00738 }
00739
00740 return error_code;
00741 }
00742
00743
00744 std::ostream & colvarbias_abf::write_state_data(std::ostream& os)
00745 {
00746 std::ios::fmtflags flags(os.flags());
00747
00748 os.setf(std::ios::fmtflags(0), std::ios::floatfield);
00749 os << "\nsamples\n";
00750 samples->write_raw(os, 8);
00751 os.flags(flags);
00752
00753 os << "\ngradient\n";
00754 gradients->write_raw(os, 8);
00755
00756 if (b_CZAR_estimator) {
00757 os.setf(std::ios::fmtflags(0), std::ios::floatfield);
00758 os << "\nz_samples\n";
00759 z_samples->write_raw(os, 8);
00760 os.flags(flags);
00761 os << "\nz_gradient\n";
00762 z_gradients->write_raw(os, 8);
00763 }
00764
00765 os.flags(flags);
00766 return os;
00767 }
00768
00769
00770 std::istream & colvarbias_abf::read_state_data(std::istream& is)
00771 {
00772 if ( input_prefix.size() > 0 ) {
00773 cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", COLVARS_INPUT_ERROR);
00774 }
00775
00776 if (! read_state_data_key(is, "samples")) {
00777 return is;
00778 }
00779 if (! samples->read_raw(is)) {
00780 return is;
00781 }
00782
00783 if (! read_state_data_key(is, "gradient")) {
00784 return is;
00785 }
00786 if (! gradients->read_raw(is)) {
00787 return is;
00788 }
00789 if (b_integrate) {
00790
00791 pmf->set_div();
00792 }
00793
00794 if (b_CZAR_estimator) {
00795
00796 if (! read_state_data_key(is, "z_samples")) {
00797 return is;
00798 }
00799 if (! z_samples->read_raw(is)) {
00800 return is;
00801 }
00802
00803 if (! read_state_data_key(is, "z_gradient")) {
00804 return is;
00805 }
00806 if (! z_gradients->read_raw(is)) {
00807 return is;
00808 }
00809 }
00810
00811 return is;
00812 }
00813
00814
00815 int colvarbias_abf::write_output_files()
00816 {
00817 if (cvm::debug()) {
00818 cvm::log("ABF bias trying to write gradients and samples to disk");
00819 }
00820
00821 if (shared_on && cvm::main()->proxy->replica_index() > 0
00822 && ! (b_CZAR_estimator || b_UI_estimator) ) {
00823
00824
00825 return COLVARS_OK;
00826 }
00827
00828 write_gradients_samples(output_prefix);
00829 if (b_history_files) {
00830 if ((cvm::step_absolute() % history_freq) == 0) {
00831 write_gradients_samples(output_prefix + ".hist", false);
00832 }
00833 }
00834
00835 if (b_UI_estimator) {
00836 eabf_UI.calc_pmf();
00837 eabf_UI.write_files();
00838 }
00839
00840 return COLVARS_OK;
00841 }
00842
00843
00844 int colvarbias_abf::calc_energy(std::vector<colvarvalue> const *values)
00845 {
00846 bias_energy = 0.0;
00847
00848 if (num_variables() > 1 || values != NULL) {
00849
00850
00851 if (pmf != NULL) {
00852 std::vector<int> const curr_bin = values ?
00853 pmf->get_colvars_index(*values) :
00854 pmf->get_colvars_index();
00855
00856 if (pmf->index_ok(curr_bin)) {
00857 bias_energy = pmf->value(curr_bin);
00858 }
00859 }
00860 return COLVARS_OK;
00861 }
00862
00863
00864 int home0 = gradients->current_bin_scalar(0);
00865 if (home0 < 0) return COLVARS_OK;
00866 int gradient_len = (int)(gradients->number_of_points(0));
00867 int home = (home0 < gradient_len) ? home0 : (gradient_len-1);
00868
00869
00870 cvm::real sum = 0.0;
00871 for (int i = 0; i < home; i++) {
00872 std::vector<int> ix(1,i);
00873
00874
00875 unsigned int count = samples->value(ix);
00876 cvm::real fact = 1.0;
00877 if ( count < full_samples ) {
00878 fact = (count < min_samples) ? 0.0 :
00879 (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
00880 }
00881 if (count > 0) sum += fact*gradients->value(ix)/count*gradients->widths[0];
00882 }
00883
00884
00885 std::vector<int> ix(1,home);
00886 cvm::real frac = gradients->current_bin_scalar_fraction(0);
00887 unsigned int count = samples->value(ix);
00888 cvm::real fact = 1.0;
00889 if ( count < full_samples ) {
00890 fact = (count < min_samples) ? 0.0 :
00891 (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
00892 }
00893 if (count > 0)
00894 sum += fact*gradients->value(ix)/count*gradients->widths[0]*frac;
00895
00896
00897 bias_energy = -sum;
00898 return COLVARS_OK;
00899 }