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

colvarbias_abf.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 
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     // Samples at step zero can not be collected
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   // ************* parsing general ABF options ***********************
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   // full_samples - min_samples >= 1 is guaranteed
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   // shared ABF
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     // If shared_freq is not set, we default to output_freq
00116     get_keyval(conf, "sharedFreq", shared_freq, output_freq);
00117   }
00118 
00119   // ************* checking the associated colvars *******************
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     // Request calculation of total force
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); // Could be a child dependency of a f_cvb_use_grids feature
00139     if (hide_Jacobian) {
00140       colvars[i]->enable(f_cv_hide_Jacobian);
00141     }
00142 
00143     // If any colvar is extended-system (restrained style, not external with constraint), we are running eABF
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     // Cannot mix and match coarse time steps with ABF because it gives
00150     // wrong total force averages - total force needs to be averaged over
00151     // every time step
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     // Here we could check for orthogonality of the Cartesian coordinates
00160     // and make it just a warning if some parameter is set?
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   // Construct empty grids based on the colvars
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   // Data for eAB F z-based estimator
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     // CZAR output files for stratified eABF
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); // Integrate for output if d<=3
00219   if (b_integrate) {
00220     // For now, we integrate on-the-fly iff the grid is < 3D
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     // Parameters for integrating initial (and final) gradient data
00230     get_keyval(conf, "integrateMaxIterations", integrate_iterations, 10000, colvarparse::parse_silent);
00231     get_keyval(conf, "integrateTol", integrate_tol, 1e-6, colvarparse::parse_silent);
00232     // Projected ABF, updating the integrated PMF on the fly
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   // For shared ABF, we store a second set of grids.
00239   // This used to be only if "shared" was defined,
00240   // but now we allow calling share externally (e.g. from Tcl).
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   // If custom grids are provided, read them
00248   if ( input_prefix.size() > 0 ) {
00249     read_gradients_samples();
00250     // Update divergence to account for input data
00251     pmf->set_div();
00252   }
00253 
00254   // if extendedLangrangian is on, then call UI estimator
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,                // force constant in eABF
00278                                          output_prefix,              // the prefix of output files
00279                                          cvm::restart_out_freq,
00280                                          UI_restart,                    // whether restart from a .count and a .grad file
00281                                          input_prefix,   // the prefixes of input files
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   // shared ABF
00329   // We used to only do this if "shared" was defined,
00330   // but now we can call shared externally
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     // e.g. in LAMMPS, total forces are current
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 //       if (b_adiabatic_reweighting) {
00368 //         // Update gradients non-locally based on conditional distribution of
00369 //         // fictitious variable TODO
00370 //
00371 //       } else
00372       if (samples->index_ok(force_bin)) {
00373         // Only if requested and within bounds of the grid...
00374 
00375         for (i = 0; i < num_variables(); i++) {
00376           // get total forces (lagging by 1 timestep) from colvars
00377           // and subtract previous ABF force if necessary
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           // If we are outside the range of xi, the force has not been obtained above
00394           // the function is just an accessor, so cheap to call again anyway
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(); // TODO: do this only when necessary
00411       }
00412     }
00413   }
00414 
00415   if (!cvm::proxy->total_forces_same_step()) {
00416     // e.g. in NAMD, total forces will be available for next timestep
00417     // hence we store the current colvar bin
00418     force_bin = bin;
00419   }
00420 
00421   // Reset biasing forces from previous timestep
00422   for (i = 0; i < num_variables(); i++) {
00423     colvar_forces[i].reset();
00424   }
00425 
00426   // Compute and apply the new bias, if applicable
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     // Factor that ensures smooth introduction of the force
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       // In projected ABF, the force is the PMF gradient estimate
00442       pmf->vector_gradient_finite_diff(bin, grad);
00443     } else {
00444       // Normal ABF
00445       gradients->vector_value(bin, grad);
00446     }
00447 
00448 //     if ( b_adiabatic_reweighting) {
00449 //       // Average of force according to conditional distribution of fictitious variable
00450 //       // need freshly integrated PMF, gradient TODO
00451 //     } else
00452     if ( fact != 0.0 ) {
00453       if ( (num_variables() == 1) && colvars[0]->periodic_boundaries() ) {
00454         // Enforce a zero-mean bias on periodic, 1D coordinates
00455         // in other words: boundary condition is that the biasing potential is periodic
00456         // This is enforced naturally if using integrated PMF
00457         colvar_forces[0].real_value = fact * (grad[0] - gradients->average ());
00458       } else {
00459         for (i = 0; i < num_variables(); i++) {
00460           // subtracting the mean force (opposite of the FE gradient) means adding the gradient
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   // update the output prefix; TODO: move later to setup_output() function
00475   if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
00476     // This is the only bias computing PMFs
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     // Share gradients and samples for shared ABF.
00484     replica_share();
00485   }
00486 
00487   // Prepare for the first sharing.
00488   if (shared_last_step < 0) {
00489     // Copy the current gradient and count values into last.
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   // update UI estimator every step
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   // We must have stored the last_gradients and last_samples.
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   // Share gradients for shared ABF.
00532   cvm::log("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) );
00533 
00534   // Count of data items.
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     // Replica 0 collects the delta gradient and count from the others.
00543     for (p = 1; p < proxy->num_replicas(); p++) {
00544       // Receive the deltas.
00545       proxy->replica_comm_recv(msg_data, msg_total, p);
00546 
00547       // Map the deltas from the others into the grids.
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       // Combine the delta gradient and count of the other replicas
00552       // with Replica 0's current state (including its delta).
00553       gradients->add_grid( *last_gradients );
00554       samples->add_grid( *last_samples );
00555     }
00556 
00557     // Now we must send the combined gradient to the other replicas.
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     // All other replicas send their delta gradient and count.
00566     // Calculate the delta gradient and count.
00567     last_gradients->delta_grid(*gradients);
00568     last_samples->delta_grid(*samples);
00569 
00570     // Cast the raw char data to the gradient and samples.
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     // We now receive the combined gradient from Replica 0.
00576     proxy->replica_comm_recv(msg_data, msg_total, 0);
00577     // We sync to the combined gradient computed by Replica 0.
00578     gradients->raw_data_in((cvm::real*)(&msg_data[0]));
00579     samples->raw_data_in((size_t*)(&msg_data[samp_start]));
00580   }
00581 
00582   // Without a barrier it's possible that one replica starts
00583   // share 2 when other replicas haven't finished share 1.
00584   proxy->replica_comm_barrier();
00585   // Done syncing the replicas.
00586   delete[] msg_data;
00587 
00588   // Copy the current gradient and count values into last.
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     // Update divergence to account for newly shared gradients
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     // Insert empty line between frames in history files
00613     os << std::endl;
00614     cvm::proxy->flush_output_stream(filename);
00615   }
00616 
00617   // In dimension higher than 2, dx is easier to handle and visualize
00618   // but we cannot write multiple frames in a dx file now
00619   // (could be implemented as multiple dx files)
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     // if (close) {
00628       cvm::proxy->close_output_stream(dx);
00629     // }
00630     // else {
00631     //   // TODO, decide convention for multiple datasets in dx file
00632     //   *dx_os << std::endl;
00633     //   dx_os->flush();
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     // Do numerical integration (to high precision) and output a PMF
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     // Write eABF CZAR-related quantities
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     // Calculate CZAR estimator of gradients
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       // Do numerical integration (to high precision) and output a PMF
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 // For Tcl implementation of selection rules.
00687 int colvarbias_abf::bin_num() {
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     // For user-provided files, the per-bias naming scheme may not apply
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       // Read eABF z-averaged data for CZAR
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); // default floating-point format
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); // default floating-point format
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     // Update divergence to account for restart data
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     // No need to report the same data as replica 0, let it do the I/O job
00824     // except if using an eABF FE estimator
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; // default value, overridden if a value can be calculated
00847 
00848   if (num_variables() > 1 || values != NULL) {
00849     // Use simple estimate: neglect effect of fullSamples,
00850     // return value at center of bin
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   // Get the home bin.
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   // Integrate the gradient up to the home bin.
00870   cvm::real sum = 0.0;
00871   for (int i = 0; i < home; i++) {
00872     std::vector<int> ix(1,i);
00873 
00874     // Include the full_samples factor if necessary.
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   // Integrate the gradient up to the current position in the home interval, a fractional portion of a bin.
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   // The applied potential is the negative integral of force samples.
00897   bias_energy = -sum;
00898   return COLVARS_OK;
00899 }

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