Update Colvars to version 2018-12-14
[namd.git] / colvars / src / colvarbias_abf.cpp
1 // -*- c++ -*-
2
3 // This file is part of the Collective Variables module (Colvars).
4 // The original version of Colvars and its updates are located at:
5 // https://github.com/colvars/colvars
6 // Please update all Colvars source files before making any changes.
7 // If you wish to distribute your changes, please submit them to the
8 // Colvars repository at GitHub.
9
10 #include "colvarmodule.h"
11 #include "colvarproxy.h"
12 #include "colvar.h"
13 #include "colvarbias_abf.h"
14
15
16 colvarbias_abf::colvarbias_abf(char const *key)
17   : colvarbias(key),
18     b_UI_estimator(false),
19     b_CZAR_estimator(false),
20     pabf_freq(0),
21     system_force(NULL),
22     gradients(NULL),
23     samples(NULL),
24     pmf(NULL),
25     z_gradients(NULL),
26     z_samples(NULL),
27     czar_gradients(NULL),
28     czar_pmf(NULL),
29     last_gradients(NULL),
30     last_samples(NULL)
31 {
32 }
33
34 int colvarbias_abf::init(std::string const &conf)
35 {
36   colvarbias::init(conf);
37
38   enable(f_cvb_scalar_variables);
39   enable(f_cvb_calc_pmf);
40
41   // TODO relax this in case of VMD plugin
42   if (cvm::temperature() == 0.0)
43     cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");
44
45   // ************* parsing general ABF options ***********************
46
47   get_keyval_feature((colvarparse *)this, conf, "applyBias",  f_cvb_apply_force, true);
48   if (!is_enabled(f_cvb_apply_force)){
49     cvm::log("WARNING: ABF biases will *not* be applied!\n");
50   }
51
52   get_keyval(conf, "updateBias",  update_bias, true);
53   if (update_bias) {
54     enable(f_cvb_history_dependent);
55   } else {
56     cvm::log("WARNING: ABF biases will *not* be updated!\n");
57   }
58
59   get_keyval(conf, "hideJacobian", hide_Jacobian, false);
60   if (hide_Jacobian) {
61     cvm::log("Jacobian (geometric) forces will be handled internally.\n");
62   } else {
63     cvm::log("Jacobian (geometric) forces will be included in reported free energy gradients.\n");
64   }
65
66   get_keyval(conf, "fullSamples", full_samples, 200);
67   if ( full_samples <= 1 ) full_samples = 1;
68   min_samples = full_samples / 2;
69   // full_samples - min_samples >= 1 is guaranteed
70
71   get_keyval(conf, "inputPrefix",  input_prefix, std::vector<std::string>());
72   get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
73   get_keyval(conf, "historyFreq", history_freq, 0);
74   b_history_files = (history_freq > 0);
75
76   // shared ABF
77   get_keyval(conf, "shared", shared_on, false);
78   if (shared_on) {
79     if (!cvm::replica_enabled() || cvm::replica_num() <= 1) {
80       cvm::error("Error: shared ABF requires more than one replica.");
81       return COLVARS_ERROR;
82     }
83     cvm::log("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n");
84     if (cvm::proxy->smp_enabled() == COLVARS_OK) {
85       cvm::error("Error: shared ABF is currently not available with SMP parallelism; "
86                  "please set \"SMP off\" at the top of the Colvars configuration file.\n",
87                  COLVARS_NOT_IMPLEMENTED);
88       return COLVARS_NOT_IMPLEMENTED;
89     }
90
91     // If shared_freq is not set, we default to output_freq
92     get_keyval(conf, "sharedFreq", shared_freq, output_freq);
93   }
94
95   // ************* checking the associated colvars *******************
96
97   if (num_variables() == 0) {
98     cvm::error("Error: no collective variables specified for the ABF bias.\n");
99     return COLVARS_ERROR;
100   }
101
102   if (update_bias) {
103     // Request calculation of total force
104     if(enable(f_cvb_get_total_force)) return cvm::get_error();
105   }
106
107   bool b_extended = false;
108   size_t i;
109   for (i = 0; i < num_variables(); i++) {
110
111     if (colvars[i]->value().type() != colvarvalue::type_scalar) {
112       cvm::error("Error: ABF bias can only use scalar-type variables.\n");
113     }
114     colvars[i]->enable(f_cv_grid);
115     if (hide_Jacobian) {
116       colvars[i]->enable(f_cv_hide_Jacobian);
117     }
118
119     // If any colvar is extended-system, we need to collect the extended
120     // system gradient
121     if (colvars[i]->is_enabled(f_cv_extended_Lagrangian))
122       b_extended = true;
123
124     // Cannot mix and match coarse time steps with ABF because it gives
125     // wrong total force averages - total force needs to be averaged over
126     // every time step
127     if (colvars[i]->get_time_step_factor() != time_step_factor) {
128       cvm::error("Error: " + colvars[i]->description + " has a value of timeStepFactor ("
129         + cvm::to_str(colvars[i]->get_time_step_factor()) + ") different from that of "
130         + description + " (" + cvm::to_str(time_step_factor) + ").\n");
131       return COLVARS_ERROR;
132     }
133
134     // Here we could check for orthogonality of the Cartesian coordinates
135     // and make it just a warning if some parameter is set?
136   }
137
138   if (get_keyval(conf, "maxForce", max_force)) {
139     if (max_force.size() != num_variables()) {
140       cvm::error("Error: Number of parameters to maxForce does not match number of colvars.");
141     }
142     for (i = 0; i < num_variables(); i++) {
143       if (max_force[i] < 0.0) {
144         cvm::error("Error: maxForce should be non-negative.");
145       }
146     }
147     cap_force = true;
148   } else {
149     cap_force = false;
150   }
151
152   bin.assign(num_variables(), 0);
153   force_bin.assign(num_variables(), 0);
154   system_force = new cvm::real [num_variables()];
155
156   // Construct empty grids based on the colvars
157   if (cvm::debug()) {
158     cvm::log("Allocating count and free energy gradient grids.\n");
159   }
160
161   samples   = new colvar_grid_count(colvars);
162   gradients = new colvar_grid_gradient(colvars);
163   gradients->samples = samples;
164   samples->has_parent_data = true;
165
166   // Data for eAB F z-based estimator
167   if ( b_extended ) {
168     get_keyval(conf, "CZARestimator", b_CZAR_estimator, true);
169     // CZAR output files for stratified eABF
170     get_keyval(conf, "writeCZARwindowFile", b_czar_window_file, false,
171                colvarparse::parse_silent);
172
173     z_bin.assign(num_variables(), 0);
174     z_samples   = new colvar_grid_count(colvars);
175     z_samples->request_actual_value();
176     z_gradients = new colvar_grid_gradient(colvars);
177     z_gradients->request_actual_value();
178     z_gradients->samples = z_samples;
179     z_samples->has_parent_data = true;
180     czar_gradients = new colvar_grid_gradient(colvars);
181   }
182
183   // For now, we integrate on-the-fly iff the grid is < 3D
184   if ( num_variables() <= 3 ) {
185     pmf = new integrate_potential(colvars, gradients);
186     if ( b_CZAR_estimator ) {
187       czar_pmf = new integrate_potential(colvars, czar_gradients);
188     }
189     get_keyval(conf, "integrate", b_integrate, true); // Integrate for output
190     if ( num_variables() > 1 ) {
191       // Projected ABF
192       get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0);
193       // Parameters for integrating initial (and final) gradient data
194       get_keyval(conf, "integrateInitSteps", integrate_initial_steps, 1e4);
195       get_keyval(conf, "integrateInitTol", integrate_initial_tol, 1e-6);
196       // for updating the integrated PMF on the fly
197       get_keyval(conf, "integrateSteps", integrate_steps, 100);
198       get_keyval(conf, "integrateTol", integrate_tol, 1e-4);
199     }
200   } else {
201     b_integrate = false;
202   }
203
204   // For shared ABF, we store a second set of grids.
205   // This used to be only if "shared" was defined,
206   // but now we allow calling share externally (e.g. from Tcl).
207   last_samples   = new colvar_grid_count(colvars);
208   last_gradients = new colvar_grid_gradient(colvars);
209   last_gradients->samples = last_samples;
210   last_samples->has_parent_data = true;
211   shared_last_step = -1;
212
213   // If custom grids are provided, read them
214   if ( input_prefix.size() > 0 ) {
215     read_gradients_samples();
216     // Update divergence to account for input data
217     pmf->set_div();
218   }
219
220   // if extendedLangrangian is on, then call UI estimator
221   if (b_extended) {
222     get_keyval(conf, "UIestimator", b_UI_estimator, false);
223
224     if (b_UI_estimator) {
225     std::vector<double> UI_lowerboundary;
226     std::vector<double> UI_upperboundary;
227     std::vector<double> UI_width;
228     std::vector<double> UI_krestr;
229
230     bool UI_restart = (input_prefix.size() > 0);
231
232     for (i = 0; i < num_variables(); i++)
233     {
234       UI_lowerboundary.push_back(colvars[i]->lower_boundary);
235       UI_upperboundary.push_back(colvars[i]->upper_boundary);
236       UI_width.push_back(colvars[i]->width);
237       UI_krestr.push_back(colvars[i]->force_constant());
238     }
239       eabf_UI = UIestimator::UIestimator(UI_lowerboundary,
240                                          UI_upperboundary,
241                                          UI_width,
242                                          UI_krestr,                // force constant in eABF
243                                          output_prefix,              // the prefix of output files
244                                          cvm::restart_out_freq,
245                                          UI_restart,                    // whether restart from a .count and a .grad file
246                                          input_prefix,   // the prefixes of input files
247                                          cvm::temperature());
248     }
249   }
250
251   cvm::log("Finished ABF setup.\n");
252   return COLVARS_OK;
253 }
254
255 /// Destructor
256 colvarbias_abf::~colvarbias_abf()
257 {
258   if (samples) {
259     delete samples;
260     samples = NULL;
261   }
262
263   if (gradients) {
264     delete gradients;
265     gradients = NULL;
266   }
267
268   if (pmf) {
269     delete pmf;
270     pmf = NULL;
271   }
272
273   if (z_samples) {
274     delete z_samples;
275     z_samples = NULL;
276   }
277
278   if (z_gradients) {
279     delete z_gradients;
280     z_gradients = NULL;
281   }
282
283   if (czar_gradients) {
284     delete czar_gradients;
285     czar_gradients = NULL;
286   }
287
288   if (czar_pmf) {
289     delete czar_pmf;
290     czar_pmf = NULL;
291   }
292
293   // shared ABF
294   // We used to only do this if "shared" was defined,
295   // but now we can call shared externally
296   if (last_samples) {
297     delete last_samples;
298     last_samples = NULL;
299   }
300
301   if (last_gradients) {
302     delete last_gradients;
303     last_gradients = NULL;
304   }
305
306   if (system_force) {
307     delete [] system_force;
308     system_force = NULL;
309   }
310 }
311
312
313 /// Update the FE gradient, compute and apply biasing force
314 /// also output data to disk if needed
315
316 int colvarbias_abf::update()
317 {
318   if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
319
320   size_t i;
321   for (i = 0; i < num_variables(); i++) {
322     bin[i] = samples->current_bin_scalar(i);
323   }
324   if (cvm::proxy->total_forces_same_step()) {
325     // e.g. in LAMMPS, total forces are current
326     force_bin = bin;
327   }
328
329   if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) {
330
331     if (update_bias) {
332 //       if (b_adiabatic_reweighting) {
333 //         // Update gradients non-locally based on conditional distribution of
334 //         // fictitious variable TODO
335 //
336 //       } else
337       if (samples->index_ok(force_bin)) {
338         // Only if requested and within bounds of the grid...
339
340         for (i = 0; i < num_variables(); i++) {
341           // get total forces (lagging by 1 timestep) from colvars
342           // and subtract previous ABF force if necessary
343           update_system_force(i);
344         }
345         gradients->acc_force(force_bin, system_force);
346         if ( b_integrate ) {
347           pmf->update_div_neighbors(force_bin);
348         }
349       }
350     }
351
352     if ( z_gradients && update_bias ) {
353       for (i = 0; i < num_variables(); i++) {
354         z_bin[i] = z_samples->current_bin_scalar(i);
355       }
356       if ( z_samples->index_ok(z_bin) ) {
357         for (i = 0; i < num_variables(); i++) {
358           // If we are outside the range of xi, the force has not been obtained above
359           // the function is just an accessor, so cheap to call again anyway
360           update_system_force(i);
361         }
362         z_gradients->acc_force(z_bin, system_force);
363       }
364     }
365
366     if ( b_integrate ) {
367       if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
368         cvm::real err;
369         int iter = pmf->integrate(integrate_steps, integrate_tol, err);
370         if ( iter == integrate_steps ) {
371           cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol)
372             + " in " + cvm::to_str(integrate_steps)
373             + " steps. Residual error: " +  cvm::to_str(err));
374         }
375         pmf->set_zero_minimum(); // TODO: do this only when necessary
376       }
377     }
378   }
379
380   if (!cvm::proxy->total_forces_same_step()) {
381     // e.g. in NAMD, total forces will be available for next timestep
382     // hence we store the current colvar bin
383     force_bin = bin;
384   }
385
386   // Reset biasing forces from previous timestep
387   for (i = 0; i < num_variables(); i++) {
388     colvar_forces[i].reset();
389   }
390
391   // Compute and apply the new bias, if applicable
392   if (is_enabled(f_cvb_apply_force) && samples->index_ok(bin)) {
393
394     cvm::real count = samples->value(bin);
395     cvm::real fact = 1.0;
396
397     // Factor that ensures smooth introduction of the force
398     if ( count < full_samples ) {
399       fact = (count < min_samples) ? 0.0 :
400         (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples));
401     }
402
403     std::vector<cvm::real>  grad(num_variables());
404
405     if ( pabf_freq ) {
406       // In projected ABF, the force is the PMF gradient estimate
407       pmf->vector_gradient_finite_diff(bin, grad);
408     } else {
409       // Normal ABF
410       gradients->vector_value(bin, grad);
411     }
412
413 //     if ( b_adiabatic_reweighting) {
414 //       // Average of force according to conditional distribution of fictitious variable
415 //       // need freshly integrated PMF, gradient TODO
416 //     } else
417     if ( fact != 0.0 ) {
418       if ( (num_variables() == 1) && colvars[0]->periodic_boundaries() ) {
419         // Enforce a zero-mean bias on periodic, 1D coordinates
420         // in other words: boundary condition is that the biasing potential is periodic
421         // This is enforced naturally if using integrated PMF
422         colvar_forces[0].real_value = fact * (grad[0] - gradients->average ());
423       } else {
424         for (size_t i = 0; i < num_variables(); i++) {
425           // subtracting the mean force (opposite of the FE gradient) means adding the gradient
426           colvar_forces[i].real_value = fact * grad[i];
427         }
428       }
429       if (cap_force) {
430         for (size_t i = 0; i < num_variables(); i++) {
431           if ( colvar_forces[i].real_value * colvar_forces[i].real_value > max_force[i] * max_force[i] ) {
432             colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]);
433           }
434         }
435       }
436     }
437   }
438
439   // update the output prefix; TODO: move later to setup_output() function
440   if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) == 1) {
441     // This is the only bias computing PMFs
442     output_prefix = cvm::output_prefix();
443   } else {
444     output_prefix = cvm::output_prefix() + "." + this->name;
445   }
446
447   if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
448     if (cvm::debug()) cvm::log("ABF bias trying to write gradients and samples to disk");
449     write_gradients_samples(output_prefix);
450   }
451
452   if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
453     // file already exists iff cvm::step_relative() > 0
454     // otherwise, backup and replace
455     write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0));
456   }
457
458   if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) {
459     // Share gradients and samples for shared ABF.
460     replica_share();
461   }
462
463   // Prepare for the first sharing.
464   if (shared_last_step < 0) {
465     // Copy the current gradient and count values into last.
466     last_gradients->copy_grid(*gradients);
467     last_samples->copy_grid(*samples);
468     shared_last_step = cvm::step_absolute();
469     cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
470   }
471
472   // update UI estimator every step
473   if (b_UI_estimator)
474   {
475     std::vector<double> x(num_variables(),0);
476     std::vector<double> y(num_variables(),0);
477     for (size_t i = 0; i < num_variables(); i++)
478     {
479       x[i] = colvars[i]->actual_value();
480       y[i] = colvars[i]->value();
481     }
482     eabf_UI.update_output_filename(output_prefix);
483     eabf_UI.update(cvm::step_absolute(), x, y);
484   }
485
486   return COLVARS_OK;
487 }
488
489
490 int colvarbias_abf::replica_share() {
491
492   if ( !cvm::replica_enabled() ) {
493     cvm::error("Error: shared ABF: No replicas.\n");
494     return COLVARS_ERROR;
495   }
496   // We must have stored the last_gradients and last_samples.
497   if (shared_last_step < 0 ) {
498     cvm::error("Error: shared ABF: Tried to apply shared ABF before any sampling had occurred.\n");
499     return COLVARS_ERROR;
500   }
501
502   // Share gradients for shared ABF.
503   cvm::log("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) );
504
505   // Count of data items.
506   size_t data_n = gradients->raw_data_num();
507   size_t samp_start = data_n*sizeof(cvm::real);
508   size_t msg_total = data_n*sizeof(size_t) + samp_start;
509   char* msg_data = new char[msg_total];
510
511   if (cvm::replica_index() == 0) {
512     int p;
513     // Replica 0 collects the delta gradient and count from the others.
514     for (p = 1; p < cvm::replica_num(); p++) {
515       // Receive the deltas.
516       cvm::replica_comm_recv(msg_data, msg_total, p);
517
518       // Map the deltas from the others into the grids.
519       last_gradients->raw_data_in((cvm::real*)(&msg_data[0]));
520       last_samples->raw_data_in((size_t*)(&msg_data[samp_start]));
521
522       // Combine the delta gradient and count of the other replicas
523       // with Replica 0's current state (including its delta).
524       gradients->add_grid( *last_gradients );
525       samples->add_grid( *last_samples );
526     }
527
528     // Now we must send the combined gradient to the other replicas.
529     gradients->raw_data_out((cvm::real*)(&msg_data[0]));
530     samples->raw_data_out((size_t*)(&msg_data[samp_start]));
531     for (p = 1; p < cvm::replica_num(); p++) {
532       cvm::replica_comm_send(msg_data, msg_total, p);
533     }
534
535   } else {
536     // All other replicas send their delta gradient and count.
537     // Calculate the delta gradient and count.
538     last_gradients->delta_grid(*gradients);
539     last_samples->delta_grid(*samples);
540
541     // Cast the raw char data to the gradient and samples.
542     last_gradients->raw_data_out((cvm::real*)(&msg_data[0]));
543     last_samples->raw_data_out((size_t*)(&msg_data[samp_start]));
544     cvm::replica_comm_send(msg_data, msg_total, 0);
545
546     // We now receive the combined gradient from Replica 0.
547     cvm::replica_comm_recv(msg_data, msg_total, 0);
548     // We sync to the combined gradient computed by Replica 0.
549     gradients->raw_data_in((cvm::real*)(&msg_data[0]));
550     samples->raw_data_in((size_t*)(&msg_data[samp_start]));
551   }
552
553   // Without a barrier it's possible that one replica starts
554   // share 2 when other replicas haven't finished share 1.
555   cvm::replica_comm_barrier();
556   // Done syncing the replicas.
557   delete[] msg_data;
558
559   // Copy the current gradient and count values into last.
560   last_gradients->copy_grid(*gradients);
561   last_samples->copy_grid(*samples);
562   shared_last_step = cvm::step_absolute();
563
564   return COLVARS_OK;
565 }
566
567
568
569 void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool append)
570 {
571   std::string  samples_out_name = prefix + ".count";
572   std::string  gradients_out_name = prefix + ".grad";
573   std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
574
575   std::ostream *samples_os =
576     cvm::proxy->output_stream(samples_out_name, mode);
577   if (!samples_os) return;
578   samples->write_multicol(*samples_os);
579   cvm::proxy->close_output_stream(samples_out_name);
580
581   // In dimension higher than 2, dx is easier to handle and visualize
582   if (num_variables() > 2) {
583     std::string  samples_dx_out_name = prefix + ".count.dx";
584     std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
585     if (!samples_os) return;
586     samples->write_opendx(*samples_dx_os);
587     *samples_dx_os << std::endl;
588     cvm::proxy->close_output_stream(samples_dx_out_name);
589   }
590
591   std::ostream *gradients_os =
592     cvm::proxy->output_stream(gradients_out_name, mode);
593   if (!gradients_os) return;
594   gradients->write_multicol(*gradients_os);
595   cvm::proxy->close_output_stream(gradients_out_name);
596
597   if (b_integrate) {
598     // Do numerical integration (to high precision) and output a PMF
599     cvm::real err;
600     pmf->integrate(integrate_initial_steps, integrate_initial_tol, err);
601     pmf->set_zero_minimum();
602
603     std::string  pmf_out_name = prefix + ".pmf";
604     std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
605     if (!pmf_os) return;
606     pmf->write_multicol(*pmf_os);
607
608     // In dimension higher than 2, dx is easier to handle and visualize
609     if (num_variables() > 2) {
610       std::string  pmf_dx_out_name = prefix + ".pmf.dx";
611       std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
612       if (!pmf_dx_os) return;
613       pmf->write_opendx(*pmf_dx_os);
614       *pmf_dx_os << std::endl;
615       cvm::proxy->close_output_stream(pmf_dx_out_name);
616     }
617
618     *pmf_os << std::endl;
619     cvm::proxy->close_output_stream(pmf_out_name);
620   }
621
622   if (b_CZAR_estimator) {
623     // Write eABF CZAR-related quantities
624
625     std::string  z_samples_out_name = prefix + ".zcount";
626
627     std::ostream *z_samples_os =
628       cvm::proxy->output_stream(z_samples_out_name, mode);
629     if (!z_samples_os) return;
630     z_samples->write_multicol(*z_samples_os);
631     cvm::proxy->close_output_stream(z_samples_out_name);
632
633     if (b_czar_window_file) {
634       std::string  z_gradients_out_name = prefix + ".zgrad";
635
636       std::ostream *z_gradients_os =
637         cvm::proxy->output_stream(z_gradients_out_name, mode);
638       if (!z_gradients_os) return;
639       z_gradients->write_multicol(*z_gradients_os);
640       cvm::proxy->close_output_stream(z_gradients_out_name);
641     }
642
643     // Calculate CZAR estimator of gradients
644     for (std::vector<int> ix = czar_gradients->new_index();
645           czar_gradients->index_ok(ix); czar_gradients->incr(ix)) {
646       for (size_t n = 0; n < czar_gradients->multiplicity(); n++) {
647         czar_gradients->set_value(ix, z_gradients->value_output(ix, n)
648           - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n);
649       }
650     }
651
652     std::string  czar_gradients_out_name = prefix + ".czar.grad";
653
654     std::ostream *czar_gradients_os =
655       cvm::proxy->output_stream(czar_gradients_out_name, mode);
656     if (!czar_gradients_os) return;
657     czar_gradients->write_multicol(*czar_gradients_os);
658     cvm::proxy->close_output_stream(czar_gradients_out_name);
659
660     if (b_integrate) {
661       // Do numerical integration (to high precision) and output a PMF
662       cvm::real err;
663       czar_pmf->set_div();
664       czar_pmf->integrate(integrate_initial_steps, integrate_initial_tol, err);
665       czar_pmf->set_zero_minimum();
666
667       std::string  czar_pmf_out_name = prefix + ".czar.pmf";
668       std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
669       if (!czar_pmf_os) return;
670       czar_pmf->write_multicol(*czar_pmf_os);
671
672       // In dimension higher than 2, dx is easier to handle and visualize
673       if (num_variables() > 2) {
674         std::string  czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
675         std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
676         if (!czar_pmf_dx_os) return;
677         czar_pmf->write_opendx(*czar_pmf_dx_os);
678         *czar_pmf_dx_os << std::endl;
679         cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
680       }
681
682       *czar_pmf_os << std::endl;
683       cvm::proxy->close_output_stream(czar_pmf_out_name);
684     }
685   }
686   return;
687 }
688
689
690 // For Tcl implementation of selection rules.
691 /// Give the total number of bins for a given bias.
692 int colvarbias_abf::bin_num() {
693   return samples->number_of_points(0);
694 }
695 /// Calculate the bin index for a given bias.
696 int colvarbias_abf::current_bin() {
697   return samples->current_bin_scalar(0);
698 }
699 /// Give the count at a given bin index.
700 int colvarbias_abf::bin_count(int bin_index) {
701   if (bin_index < 0 || bin_index >= bin_num()) {
702     cvm::error("Error: Tried to get bin count from invalid bin index "+cvm::to_str(bin_index));
703     return -1;
704   }
705   std::vector<int> ix(1,(int)bin_index);
706   return samples->value(ix);
707 }
708
709
710 void colvarbias_abf::read_gradients_samples()
711 {
712   std::string samples_in_name, gradients_in_name, z_samples_in_name, z_gradients_in_name;
713
714   for ( size_t i = 0; i < input_prefix.size(); i++ ) {
715     samples_in_name = input_prefix[i] + ".count";
716     gradients_in_name = input_prefix[i] + ".grad";
717     z_samples_in_name = input_prefix[i] + ".zcount";
718     z_gradients_in_name = input_prefix[i] + ".zgrad";
719     // For user-provided files, the per-bias naming scheme may not apply
720
721     std::ifstream is;
722
723     cvm::log("Reading sample count from " + samples_in_name + " and gradient from " + gradients_in_name);
724     is.open(samples_in_name.c_str());
725     if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
726     samples->read_multicol(is, true);
727     is.close();
728     is.clear();
729
730     is.open(gradients_in_name.c_str());
731     if (!is.is_open()) {
732       cvm::error("Error opening ABF gradient file " +
733                  gradients_in_name + " for reading", INPUT_ERROR);
734     } else {
735       gradients->read_multicol(is, true);
736       is.close();
737     }
738
739     if (b_CZAR_estimator) {
740       // Read eABF z-averaged data for CZAR
741       cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name);
742
743       is.clear();
744       is.open(z_samples_in_name.c_str());
745       if (!is.is_open())  cvm::error("Error opening eABF z-histogram file " + z_samples_in_name + " for reading");
746       z_samples->read_multicol(is, true);
747       is.close();
748       is.clear();
749
750       is.open(z_gradients_in_name.c_str());
751       if (!is.is_open())  cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading");
752       z_gradients->read_multicol(is, true);
753       is.close();
754     }
755   }
756   return;
757 }
758
759
760 std::ostream & colvarbias_abf::write_state_data(std::ostream& os)
761 {
762   std::ios::fmtflags flags(os.flags());
763
764   os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
765   os << "\nsamples\n";
766   samples->write_raw(os, 8);
767   os.flags(flags);
768
769   os << "\ngradient\n";
770   gradients->write_raw(os, 8);
771
772   if (b_CZAR_estimator) {
773     os.setf(std::ios::fmtflags(0), std::ios::floatfield); // default floating-point format
774     os << "\nz_samples\n";
775     z_samples->write_raw(os, 8);
776     os.flags(flags);
777     os << "\nz_gradient\n";
778     z_gradients->write_raw(os, 8);
779   }
780
781   os.flags(flags);
782   return os;
783 }
784
785
786 std::istream & colvarbias_abf::read_state_data(std::istream& is)
787 {
788   if ( input_prefix.size() > 0 ) {
789     cvm::error("ERROR: cannot provide both inputPrefix and a colvars state file.\n", INPUT_ERROR);
790   }
791
792   if (! read_state_data_key(is, "samples")) {
793     return is;
794   }
795   if (! samples->read_raw(is)) {
796     return is;
797   }
798
799   if (! read_state_data_key(is, "gradient")) {
800     return is;
801   }
802   if (! gradients->read_raw(is)) {
803     return is;
804   }
805   if (b_integrate) {
806     // Update divergence to account for restart data
807     pmf->set_div();
808   }
809
810   if (b_CZAR_estimator) {
811
812     if (! read_state_data_key(is, "z_samples")) {
813       return is;
814     }
815     if (! z_samples->read_raw(is)) {
816       return is;
817     }
818
819     if (! read_state_data_key(is, "z_gradient")) {
820       return is;
821     }
822     if (! z_gradients->read_raw(is)) {
823       return is;
824     }
825   }
826
827   return is;
828 }
829
830 int colvarbias_abf::write_output_files()
831 {
832   write_gradients_samples(output_prefix);
833   return COLVARS_OK;
834 }