00001
00002
00003 #ifndef COLVARGRID_H
00004 #define COLVARGRID_H
00005
00006 #include <iostream>
00007 #include <iomanip>
00008 #include <cmath>
00009
00010 #include "colvar.h"
00011 #include "colvarmodule.h"
00012 #include "colvarvalue.h"
00013 #include "colvarparse.h"
00014
00019 template <class T> class colvar_grid : public colvarparse {
00020
00021 protected:
00022
00024 size_t nd;
00025
00027 std::vector<int> nx;
00028
00030 std::vector<int> nxc;
00031
00034 size_t mult;
00035
00037 size_t nt;
00038
00040 std::vector<T> data;
00041
00043 std::vector<size_t> new_data;
00044
00046 std::vector<colvar *> cv;
00047
00049 inline size_t address (std::vector<int> const &ix) const
00050 {
00051 size_t addr = 0;
00052 for (size_t i = 0; i < nd; i++) {
00053 addr += ix[i]*nxc[i];
00054 if (cvm::debug()) {
00055 if (ix[i] >= nx[i])
00056 cvm::fatal_error ("Error: exceeding bounds in colvar_grid.\n");
00057 }
00058 }
00059 return addr;
00060 }
00061
00062 public:
00063
00065 std::vector<colvarvalue> lower_boundaries;
00066
00068 std::vector<colvarvalue> upper_boundaries;
00069
00071 std::vector<bool> periodic;
00072
00074 std::vector<cvm::real> widths;
00075
00077 bool has_parent_data;
00078
00080 inline size_t number_of_colvars() const
00081 {
00082 return nd;
00083 }
00084
00087 inline size_t number_of_points (int const icv = -1) const
00088 {
00089 if (icv < 0) {
00090 return nt;
00091 } else {
00092 return nx[icv];
00093 }
00094 }
00095
00097 inline std::vector<int> const & sizes() const
00098 {
00099 return nx;
00100 }
00101
00103 inline void set_sizes (std::vector<int> const &new_sizes)
00104 {
00105 nx = new_sizes;
00106 }
00107
00109 inline size_t multiplicity() const
00110 {
00111 return mult;
00112 }
00113
00115 void create (std::vector<int> const &nx_i,
00116 T const &t = T(),
00117 size_t const &mult_i = 1)
00118 {
00119 mult = mult_i;
00120 nd = nx_i.size();
00121 nxc.resize (nd);
00122 nx = nx_i;
00123
00124 nt = mult;
00125 for (int i = nd-1; i >= 0; i--) {
00126 if (nx_i[i] <= 0)
00127 cvm::fatal_error ("Error: providing an invalid number of points, "+
00128 cvm::to_str (nx_i[i])+".\n");
00129 nxc[i] = nt;
00130 nt *= nx[i];
00131 }
00132
00133 data.reserve (nt);
00134 data.assign (nt, t);
00135 }
00136
00138 void create()
00139 {
00140 create (this->nx, T(), this->mult);
00141 }
00142
00143
00145 colvar_grid()
00146 {
00147 save_delimiters = false;
00148 nd = nt = 0;
00149 }
00150
00152 virtual ~colvar_grid()
00153 {}
00154
00158 colvar_grid (colvar_grid<T> const &g)
00159 : nd (g.nd),
00160 mult (g.mult),
00161 cv (g.cv),
00162 lower_boundaries (g.lower_boundaries),
00163 upper_boundaries (g.upper_boundaries),
00164 periodic (g.periodic),
00165 widths (g.widths),
00166 data()
00167 {
00168 save_delimiters = false;
00169 }
00170
00175 colvar_grid (std::vector<int> const &nx_i,
00176 T const &t = T(),
00177 size_t const &mult_i = 1)
00178 {
00179 save_delimiters = false;
00180 this->create (nx_i, t, mult_i);
00181 }
00182
00184 colvar_grid (std::vector<colvar *> const &colvars,
00185 T const &t = T(),
00186 size_t const &mult_i = 1,
00187 bool margin = false)
00188 : cv (colvars)
00189 {
00190 save_delimiters = false;
00191
00192 std::vector<int> nx_i;
00193
00194 if (cvm::debug())
00195 cvm::log ("Allocating a grid for "+cvm::to_str (colvars.size())+
00196 " collective variables.\n");
00197
00198 for (size_t i = 0; i < cv.size(); i++) {
00199
00200 if (cv[i]->type() != colvarvalue::type_scalar) {
00201 cvm::fatal_error ("Colvar grids can only be automatically "
00202 "constructed for scalar variables. "
00203 "ABF and histogram can not be used; metadynamics "
00204 "can be used with useGrids disabled.\n");
00205 }
00206
00207 if (cv[i]->width <= 0.0) {
00208 cvm::fatal_error ("Tried to initialize a grid on a"
00209 "variable with negative or zero width.\n");
00210 }
00211
00212 if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) {
00213 cvm::fatal_error ("Tried to initialize a grid on a"
00214 "variable with undefined boundaries.\n");
00215 }
00216
00217 widths.push_back (cv[i]->width);
00218 periodic.push_back (cv[i]->periodic_boundaries());
00219
00220 if (margin) {
00221 if (periodic[i]) {
00222
00223 lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
00224 upper_boundaries.push_back (cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
00225 } else {
00226
00227 lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
00228 upper_boundaries.push_back (cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
00229 }
00230 } else {
00231 lower_boundaries.push_back (cv[i]->lower_boundary);
00232 upper_boundaries.push_back (cv[i]->upper_boundary);
00233 }
00234
00235
00236 {
00237 cvm::real nbins = ( upper_boundaries[i].real_value -
00238 lower_boundaries[i].real_value ) / widths[i];
00239
00240 int nbins_round = (int)(nbins+0.5);
00241
00242 if (::fabs (nbins - cvm::real (nbins_round)) > 1.0E-10) {
00243 cvm::log ("Warning: grid interval ("+
00244 cvm::to_str (lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
00245 cvm::to_str (upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
00246 ") is not commensurate to its bin width ("+
00247 cvm::to_str (widths[i], cvm::cv_width, cvm::cv_prec)+").\n");
00248 upper_boundaries[i].real_value = lower_boundaries[i].real_value +
00249 (nbins_round * widths[i]);
00250 }
00251
00252 if (cvm::debug())
00253 cvm::log ("Number of points is "+cvm::to_str ((int) nbins_round)+
00254 " for the colvar no. "+cvm::to_str (i+1)+".\n");
00255
00256 nx_i.push_back (nbins_round);
00257 }
00258
00259 }
00260
00261 create (nx_i, t, mult_i);
00262 }
00263
00264
00267 inline void wrap (std::vector<int> & ix) const
00268 {
00269 for (size_t i = 0; i < nd; i++) {
00270 if (periodic[i]) {
00271 ix[i] = (ix[i] + nx[i]) % nx[i];
00272 } else {
00273 if (ix[i] < 0 || ix[i] >= nx[i])
00274 cvm::fatal_error ("Trying to wrap illegal index vector (non-PBC): "
00275 + cvm::to_str (ix));
00276 }
00277 }
00278 }
00279
00280
00281
00284 inline int value_to_bin_scalar (colvarvalue const &value, const int i) const
00285 {
00286 return (int) ::floor ( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
00287 }
00288
00291 inline colvarvalue bin_to_value_scalar (int const &i_bin, int const i) const
00292 {
00293 return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
00294 }
00295
00297 inline void set_value (std::vector<int> const &ix,
00298 T const &t,
00299 size_t const &imult = 0)
00300 {
00301 data[this->address (ix)+imult] = t;
00302 }
00303
00304
00307 inline T const & value (std::vector<int> const &ix,
00308 size_t const &imult = 0) const
00309 {
00310 return data[this->address (ix) + imult];
00311 }
00312
00313
00315 inline void add_constant (T const &t)
00316 {
00317 for (size_t i = 0; i < nt; i++)
00318 data[i] += t;
00319 }
00320
00322 inline void multiply_constant (cvm::real const &a)
00323 {
00324 for (size_t i = 0; i < nt; i++)
00325 data[i] *= a;
00326 }
00327
00328
00331 inline std::vector<int> const get_colvars_index (std::vector<colvarvalue> const &values) const
00332 {
00333 std::vector<int> index = new_index();
00334 for (size_t i = 0; i < nd; i++) {
00335 index[i] = cv[i]->value_to_bin_scalar (values[i], lower_boundaries[i], widths[i]);
00336 }
00337 return index;
00338 }
00339
00342 inline std::vector<int> const get_colvars_index() const
00343 {
00344 std::vector<colvarvalue> current_values (nd);
00345 for (size_t i = 0; i < nd; i++) {
00346 current_values[i] = cv[i]->value();
00347 }
00348 return get_colvars_index (current_values);
00349 }
00350
00353 inline cvm::real bin_distance_from_boundaries (std::vector<colvarvalue> const &values)
00354 {
00355 cvm::real minimum = 1.0E+16;
00356 for (size_t i = 0; i < nd; i++) {
00357
00358 cvm::real dl = ::sqrt (cv[i]->dist2 (values[i], lower_boundaries[i])) / widths[i];
00359 cvm::real du = ::sqrt (cv[i]->dist2 (values[i], upper_boundaries[i])) / widths[i];
00360
00361 if (! periodic[i]) {
00362 if (values[i].real_value < lower_boundaries[i])
00363 dl *= -1.0;
00364 if (values[i].real_value > upper_boundaries[i])
00365 du *= -1.0;
00366 }
00367
00368 if (dl < minimum)
00369 minimum = dl;
00370 if (du < minimum)
00371 minimum = du;
00372 }
00373
00374 return minimum;
00375 }
00376
00377
00385 void map_grid (colvar_grid<T> const &other_grid)
00386 {
00387 if (other_grid.multiplicity() != this->multiplicity())
00388 cvm::fatal_error ("Error: trying to merge two grids with values of "
00389 "different multiplicity.\n");
00390
00391 std::vector<colvarvalue> const &gb = this->lower_boundaries;
00392 std::vector<cvm::real> const &gw = this->widths;
00393 std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
00394 std::vector<cvm::real> const &ogw = other_grid.widths;
00395
00396 std::vector<int> ix = this->new_index();
00397 std::vector<int> oix = other_grid.new_index();
00398
00399 if (cvm::debug())
00400 cvm::log ("Remapping grid...\n");
00401 for ( ; this->index_ok (ix); this->incr (ix)) {
00402
00403 for (size_t i = 0; i < nd; i++) {
00404 oix[i] =
00405 cv[i]->value_to_bin_scalar (cv[i]->bin_to_value_scalar (ix[i], gb[i], gw[i]),
00406 ogb[i],
00407 ogw[i]);
00408 }
00409
00410 if (! other_grid.index_ok (oix)) {
00411 continue;
00412 }
00413
00414 for (size_t im = 0; im < mult; im++) {
00415 this->set_value (ix, other_grid.value (oix, im), im);
00416 }
00417 }
00418 if (cvm::debug())
00419 cvm::log ("Remapping done.\n");
00420 }
00421
00424 virtual inline T value_output (std::vector<int> const &ix,
00425 size_t const &imult = 0)
00426 {
00427 return value (ix, imult);
00428 }
00429
00433 virtual inline void value_input (std::vector<int> const &ix,
00434 T const &t,
00435 size_t const &imult = 0,
00436 bool add = false)
00437 {
00438 if ( add )
00439 data[address (ix) + imult] += t;
00440 else
00441 data[address (ix) + imult] = t;
00442 }
00443
00444
00445
00446
00447
00448
00449
00452 inline std::vector<int> const new_index() const
00453 {
00454 return std::vector<int> (nd, 0);
00455 }
00456
00459 inline bool index_ok (std::vector<int> const &ix) const
00460 {
00461 for (size_t i = 0; i < nd; i++) {
00462 if ( (ix[i] < 0) || (ix[i] >= int (nx[i])) )
00463 return false;
00464 }
00465 return true;
00466 }
00467
00470 inline void incr (std::vector<int> &ix) const
00471 {
00472 for (int i = ix.size()-1; i >= 0; i--) {
00473
00474 ix[i]++;
00475
00476 if (ix[i] >= nx[i]) {
00477
00478 if (i > 0) {
00479 ix[i] = 0;
00480 continue;
00481 } else {
00482
00483
00484
00485 ix[0] = nx[0];
00486 return;
00487 }
00488 } else {
00489 return;
00490 }
00491 }
00492 }
00493
00495 std::ostream & write_params (std::ostream &os)
00496 {
00497 os << "grid_parameters {\n n_colvars " << nd << "\n";
00498
00499 os << " lower_boundaries ";
00500 for (size_t i = 0; i < nd; i++)
00501 os << " " << lower_boundaries[i];
00502 os << "\n";
00503
00504 os << " upper_boundaries ";
00505 for (size_t i = 0; i < nd; i++)
00506 os << " " << upper_boundaries[i];
00507 os << "\n";
00508
00509 os << " widths ";
00510 for (size_t i = 0; i < nd; i++)
00511 os << " " << widths[i];
00512 os << "\n";
00513
00514 os << " sizes ";
00515 for (size_t i = 0; i < nd; i++)
00516 os << " " << nx[i];
00517 os << "\n";
00518
00519 os << "}\n";
00520 return os;
00521 }
00522
00523
00524 bool parse_params (std::string const &conf)
00525 {
00526 std::vector<int> old_nx = nx;
00527 std::vector<colvarvalue> old_lb = lower_boundaries;
00528
00529 {
00530 size_t nd_in = 0;
00531 colvarparse::get_keyval (conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
00532 if (nd_in != nd)
00533 cvm::fatal_error ("Error: trying to read data for a grid "
00534 "that contains a different number of colvars ("+
00535 cvm::to_str (nd_in)+") than the grid defined "
00536 "in the configuration file ("+cvm::to_str (nd)+
00537 ").\n");
00538 }
00539
00540 colvarparse::get_keyval (conf, "lower_boundaries",
00541 lower_boundaries, lower_boundaries, colvarparse::parse_silent);
00542
00543 colvarparse::get_keyval (conf, "upper_boundaries",
00544 upper_boundaries, upper_boundaries, colvarparse::parse_silent);
00545
00546 colvarparse::get_keyval (conf, "widths", widths, widths, colvarparse::parse_silent);
00547
00548 colvarparse::get_keyval (conf, "sizes", nx, nx, colvarparse::parse_silent);
00549
00550 bool new_params = false;
00551 for (size_t i = 0; i < nd; i++) {
00552 if ( (old_nx[i] != nx[i]) ||
00553 (::sqrt (cv[i]->dist2 (old_lb[i],
00554 lower_boundaries[i])) > 1.0E-10) ) {
00555 new_params = true;
00556 }
00557 }
00558
00559
00560 if (new_params) {
00561 data.resize (0);
00562 this->create (nx, T(), mult);
00563 }
00564
00565 return true;
00566 }
00567
00571 void check_consistency()
00572 {
00573 for (size_t i = 0; i < nd; i++) {
00574 if ( (::sqrt (cv[i]->dist2 (cv[i]->lower_boundary,
00575 lower_boundaries[i])) > 1.0E-10) ||
00576 (::sqrt (cv[i]->dist2 (cv[i]->upper_boundary,
00577 upper_boundaries[i])) > 1.0E-10) ||
00578 (::sqrt (cv[i]->dist2 (cv[i]->width,
00579 widths[i])) > 1.0E-10) ) {
00580 cvm::fatal_error ("Error: restart information for a grid is "
00581 "inconsistent with that of its colvars.\n");
00582 }
00583 }
00584 }
00585
00586
00590 std::ostream & write_raw (std::ostream &os,
00591 size_t const buf_size = 3)
00592 {
00593 std::streamsize const w = os.width();
00594 std::streamsize const p = os.precision();
00595
00596 std::vector<int> ix = new_index();
00597 size_t count = 0;
00598 for ( ; index_ok (ix); incr (ix)) {
00599 for (size_t imult = 0; imult < mult; imult++) {
00600 os << " "
00601 << std::setw (w) << std::setprecision (p)
00602 << value_output (ix, imult);
00603 if (((++count) % buf_size) == 0)
00604 os << "\n";
00605 }
00606 }
00607
00608 if ((count % buf_size) != 0)
00609 os << "\n";
00610
00611 return os;
00612 }
00613
00615 std::istream & read_raw (std::istream &is)
00616 {
00617 std::vector<int> ix = new_index();
00618 size_t count = 0;
00619 for ( ; index_ok (ix); incr (ix)) {
00620 for (size_t imult = 0; imult < mult; imult++) {
00621 T new_value;
00622 if (is >> new_value) {
00623 value_input (ix, new_value, imult);
00624 count++;
00625 } else {
00626 cvm::fatal_error ("Error: malformed grid data in input.\n");
00627 }
00628 }
00629 }
00630
00631 if (count != number_of_points()) {
00632 cvm::log ("Error: number of grid points read from file ("+
00633 cvm::to_str (count)+") is wrong (should be "+
00634 cvm::to_str (number_of_points())+").\n");
00635 cvm::fatal_error ("Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt.\n");
00636 }
00637
00638 return is;
00639 }
00640
00643 void write_multicol (std::ostream &os)
00644 {
00645 std::streamsize const w = os.width();
00646 std::streamsize const p = os.precision();
00647
00648
00649
00650
00651 os << std::setw (2) << "# " << nd << "\n";
00652 for (size_t i = 0; i < nd; i++) {
00653 os << "# "
00654 << std::setw (10) << lower_boundaries[i]
00655 << std::setw (10) << widths[i]
00656 << std::setw (10) << nx[i] << " "
00657 << periodic[i] << "\n";
00658 }
00659
00660 for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
00661
00662 if (ix.back() == 0) {
00663
00664 os << "\n";
00665 }
00666
00667 for (size_t i = 0; i < nd; i++) {
00668 os << " "
00669 << std::setw (w) << std::setprecision (p)
00670 << bin_to_value_scalar (ix[i], i);
00671 }
00672 os << " ";
00673 for (size_t imult = 0; imult < mult; imult++) {
00674 os << " "
00675 << std::setw (w) << std::setprecision (p)
00676 << value_output (ix, imult);
00677 }
00678 os << "\n";
00679 }
00680 }
00681
00684 void read_multicol (std::istream &is, bool add = false)
00685 {
00686
00687
00688
00689 std::string hash;
00690 cvm::real lower, width, x;
00691 size_t n, periodic;
00692 bool remap;
00693 std::vector<T> new_value;
00694 std::vector<int> nx_read;
00695 std::vector<int> bin;
00696
00697 if ( cv.size() != nd ) {
00698 cvm::fatal_error ("Cannot read grid file: missing reference to colvars.");
00699 }
00700
00701 if ( !(is >> hash) || (hash != "#") ) {
00702 cvm::fatal_error ("Error reading grid at position "+
00703 cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
00704 }
00705
00706 is >> n;
00707 if ( n != nd ) {
00708 cvm::fatal_error ("Error reading grid: wrong number of collective variables.\n");
00709 }
00710
00711 nx_read.resize (n);
00712 bin.resize (n);
00713 new_value.resize (mult);
00714
00715 if (this->has_parent_data && add) {
00716 new_data.resize (data.size());
00717 }
00718
00719 remap = false;
00720 for (size_t i = 0; i < nd; i++ ) {
00721 if ( !(is >> hash) || (hash != "#") ) {
00722 cvm::fatal_error ("Error reading grid at position "+
00723 cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
00724 }
00725
00726 is >> lower >> width >> nx_read[i] >> periodic;
00727
00728
00729 if ( (::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) ||
00730 (::fabs (width - widths[i] ) > 1.0e-10) ||
00731 (nx_read[i] != nx[i]) ) {
00732 cvm::log ("Warning: reading from different grid definition (colvar "
00733 + cvm::to_str (i+1) + "); remapping data on new grid.\n");
00734 remap = true;
00735 }
00736 }
00737
00738 if ( remap ) {
00739
00740 while (is.good()) {
00741 bool end_of_file = false;
00742
00743 for (size_t i = 0; i < nd; i++ ) {
00744 if ( !(is >> x) ) end_of_file = true;
00745 bin[i] = value_to_bin_scalar (x, i);
00746 }
00747 if (end_of_file) break;
00748
00749 for (size_t imult = 0; imult < mult; imult++) {
00750 is >> new_value[imult];
00751 }
00752
00753 if ( index_ok(bin) ) {
00754 for (size_t imult = 0; imult < mult; imult++) {
00755 value_input (bin, new_value[imult], imult, add);
00756 }
00757 }
00758 }
00759 } else {
00760
00761 for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
00762 for (size_t i = 0; i < nd; i++ ) {
00763 is >> x;
00764 }
00765 for (size_t imult = 0; imult < mult; imult++) {
00766 is >> new_value[imult];
00767 value_input (ix, new_value[imult], imult, add);
00768 }
00769 }
00770 }
00771 return;
00772 }
00773
00774 };
00775
00776
00777
00780 class colvar_grid_count : public colvar_grid<size_t>
00781 {
00782 public:
00783
00785 colvar_grid_count();
00786
00788 virtual inline ~colvar_grid_count()
00789 {}
00790
00792 colvar_grid_count (std::vector<int> const &nx_i,
00793 size_t const &def_count = 0);
00794
00796 colvar_grid_count (std::vector<colvar *> &colvars,
00797 size_t const &def_count = 0);
00798
00800 inline void incr_count (std::vector<int> const &ix)
00801 {
00802 ++(data[this->address (ix)]);
00803 }
00804
00806 inline size_t const & new_count (std::vector<int> const &ix,
00807 size_t const &imult = 0)
00808 {
00809 return new_data[address (ix) + imult];
00810 }
00811
00813 std::istream & read_restart (std::istream &is);
00814
00816 std::ostream & write_restart (std::ostream &os);
00817
00821 virtual inline void value_input (std::vector<int> const &ix,
00822 size_t const &t,
00823 size_t const &imult = 0,
00824 bool add = false)
00825 {
00826 if (add) {
00827 data[address (ix)] += t;
00828 if (this->has_parent_data) {
00829
00830 new_data[address (ix)] = t;
00831 }
00832 } else {
00833 data[address (ix)] = t;
00834 }
00835 }
00836 };
00837
00838
00840 class colvar_grid_scalar : public colvar_grid<cvm::real>
00841 {
00842 public:
00843
00846 colvar_grid_count *samples;
00847
00849 colvar_grid_scalar();
00850
00852 colvar_grid_scalar (colvar_grid_scalar const &g);
00853
00855 ~colvar_grid_scalar();
00856
00858 colvar_grid_scalar (std::vector<int> const &nx_i);
00859
00861 colvar_grid_scalar (std::vector<colvar *> &colvars,
00862 bool margin = 0);
00863
00865 inline void acc_value (std::vector<int> const &ix,
00866 cvm::real const &new_value,
00867 size_t const &imult = 0)
00868 {
00869
00870 data[address (ix)] += new_value;
00871 if (samples)
00872 samples->incr_count (ix);
00873 }
00874
00876 inline const cvm::real * gradient_finite_diff ( const std::vector<int> &ix0 )
00877 {
00878 cvm::real A0, A1;
00879 std::vector<int> ix;
00880 if (nd != 2) cvm::fatal_error ("Finite differences available in dimension 2 only.");
00881 for (int n = 0; n < nd; n++) {
00882 ix = ix0;
00883 A0 = data[address (ix)];
00884 ix[n]++; wrap (ix);
00885 A1 = data[address (ix)];
00886 ix[1-n]++; wrap (ix);
00887 A1 += data[address (ix)];
00888 ix[n]--; wrap (ix);
00889 A0 += data[address (ix)];
00890 grad[n] = 0.5 * (A1 - A0) / widths[n];
00891 }
00892 return grad;
00893 }
00894
00897 virtual cvm::real value_output (std::vector<int> const &ix,
00898 size_t const &imult = 0)
00899 {
00900 if (imult > 0)
00901 cvm::fatal_error ("Error: trying to access a component "
00902 "larger than 1 in a scalar data grid.\n");
00903 if (samples)
00904 return (samples->value (ix) > 0) ?
00905 (data[address (ix)] / cvm::real (samples->value (ix))) :
00906 0.0;
00907 else
00908 return data[address (ix)];
00909 }
00910
00914 virtual void value_input (std::vector<int> const &ix,
00915 cvm::real const &new_value,
00916 size_t const &imult = 0,
00917 bool add = false)
00918 {
00919 if (imult > 0)
00920 cvm::fatal_error ("Error: trying to access a component "
00921 "larger than 1 in a scalar data grid.\n");
00922 if (add) {
00923 if (samples)
00924 data[address (ix)] += new_value * samples->new_count (ix);
00925 else
00926 data[address (ix)] += new_value;
00927 } else {
00928 if (samples)
00929 data[address (ix)] = new_value * samples->value (ix);
00930 else
00931 data[address (ix)] = new_value;
00932 }
00933 }
00934
00936 std::istream & read_restart (std::istream &is);
00937
00939 std::ostream & write_restart (std::ostream &os);
00940
00942 inline cvm::real maximum_value()
00943 {
00944 cvm::real max = data[0];
00945 for (size_t i = 0; i < nt; i++) {
00946 if (data[i] > max) max = data[i];
00947 }
00948 return max;
00949 }
00950
00952 inline cvm::real minimum_value()
00953 {
00954 cvm::real min = data[0];
00955 for (size_t i = 0; i < nt; i++) {
00956 if (data[i] < min) min = data[i];
00957 }
00958 return min;
00959 }
00960
00961 private:
00962
00963 cvm::real * grad;
00964 };
00965
00966
00967
00969 class colvar_grid_gradient : public colvar_grid<cvm::real>
00970 {
00971 public:
00972
00975 colvar_grid_count *samples;
00976
00978 colvar_grid_gradient();
00979
00981 virtual inline ~colvar_grid_gradient()
00982 {}
00983
00985 colvar_grid_gradient (std::vector<int> const &nx_i);
00986
00988 colvar_grid_gradient (std::vector<colvar *> &colvars);
00989
00991 inline void acc_grad (std::vector<int> const &ix, cvm::real const *grads) {
00992 for (size_t imult = 0; imult < mult; imult++) {
00993 data[address (ix) + imult] += grads[imult];
00994 }
00995 if (samples)
00996 samples->incr_count (ix);
00997 }
00998
01001 inline void acc_force (std::vector<int> const &ix, cvm::real const *forces) {
01002 for (size_t imult = 0; imult < mult; imult++) {
01003 data[address (ix) + imult] -= forces[imult];
01004 }
01005 if (samples)
01006 samples->incr_count (ix);
01007 }
01008
01011 virtual inline cvm::real value_output (std::vector<int> const &ix,
01012 size_t const &imult = 0)
01013 {
01014 if (samples)
01015 return (samples->value (ix) > 0) ?
01016 (data[address (ix) + imult] / cvm::real (samples->value (ix))) :
01017 0.0;
01018 else
01019 return data[address (ix) + imult];
01020 }
01021
01025 virtual inline void value_input (std::vector<int> const &ix,
01026 cvm::real const &new_value,
01027 size_t const &imult = 0,
01028 bool add = false)
01029 {
01030 if (add) {
01031 if (samples)
01032 data[address (ix) + imult] += new_value * samples->new_count (ix);
01033 else
01034 data[address (ix) + imult] += new_value;
01035 } else {
01036 if (samples)
01037 data[address (ix) + imult] = new_value * samples->value (ix);
01038 else
01039 data[address (ix) + imult] = new_value;
01040 }
01041 }
01042
01043
01045 std::istream & read_restart (std::istream &is);
01046
01048 std::ostream & write_restart (std::ostream &os);
01049
01051 inline cvm::real average()
01052 {
01053 size_t n = 0;
01054
01055 if (nd != 1 || nx[0] == 0) {
01056 return 0.0;
01057 }
01058
01059 cvm::real sum = 0.0;
01060 std::vector<int> ix = new_index();
01061 if (samples) {
01062 for ( ; index_ok (ix); incr (ix)) {
01063 if ( (n = samples->value (ix)) )
01064 sum += value (ix) / n;
01065 }
01066 } else {
01067 for ( ; index_ok (ix); incr (ix)) {
01068 sum += value (ix);
01069 }
01070 }
01071 return (sum / cvm::real (nx[0]));
01072 }
01073
01076 void write_1D_integral (std::ostream &os);
01077
01078 };
01079
01080
01081 #endif
01082