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 #ifndef COLVARGRID_H 00011 #define COLVARGRID_H 00012 00013 #include <iostream> 00014 #include <iomanip> 00015 00016 #include "colvar.h" 00017 #include "colvarmodule.h" 00018 #include "colvarvalue.h" 00019 #include "colvarparse.h" 00020 00025 template <class T> class colvar_grid : public colvarparse { 00026 00027 protected: 00028 00030 size_t nd; 00031 00033 std::vector<int> nx; 00034 00036 std::vector<int> nxc; 00037 00040 size_t mult; 00041 00043 size_t nt; 00044 00046 std::vector<T> data; 00047 00049 std::vector<size_t> new_data; 00050 00052 std::vector<colvar *> cv; 00053 00055 std::vector<bool> use_actual_value; 00056 00058 inline size_t address(std::vector<int> const &ix) const 00059 { 00060 size_t addr = 0; 00061 for (size_t i = 0; i < nd; i++) { 00062 addr += ix[i]*static_cast<size_t>(nxc[i]); 00063 if (cvm::debug()) { 00064 if (ix[i] >= nx[i]) { 00065 cvm::error("Error: exceeding bounds in colvar_grid.\n", COLVARS_BUG_ERROR); 00066 return 0; 00067 } 00068 } 00069 } 00070 return addr; 00071 } 00072 00073 public: 00074 00076 std::vector<colvarvalue> lower_boundaries; 00077 00079 std::vector<colvarvalue> upper_boundaries; 00080 00082 std::vector<bool> periodic; 00083 00085 std::vector<bool> hard_lower_boundaries; 00086 00088 std::vector<bool> hard_upper_boundaries; 00089 00091 std::vector<cvm::real> widths; 00092 00094 bool has_parent_data; 00095 00097 bool has_data; 00098 00100 inline size_t num_variables() const 00101 { 00102 return nd; 00103 } 00104 00106 inline std::vector<int> const &number_of_points_vec() const 00107 { 00108 return nx; 00109 } 00110 00113 inline size_t number_of_points(int const icv = -1) const 00114 { 00115 if (icv < 0) { 00116 return nt; 00117 } else { 00118 return nx[icv]; 00119 } 00120 } 00121 00123 inline std::vector<int> const & sizes() const 00124 { 00125 return nx; 00126 } 00127 00129 inline void set_sizes(std::vector<int> const &new_sizes) 00130 { 00131 nx = new_sizes; 00132 } 00133 00135 inline size_t multiplicity() const 00136 { 00137 return mult; 00138 } 00139 00141 inline void request_actual_value(bool b = true) 00142 { 00143 size_t i; 00144 for (i = 0; i < use_actual_value.size(); i++) { 00145 use_actual_value[i] = b; 00146 } 00147 } 00148 00150 int setup(std::vector<int> const &nx_i, 00151 T const &t = T(), 00152 size_t const &mult_i = 1) 00153 { 00154 if (cvm::debug()) { 00155 cvm::log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+ 00156 ", dimensions = "+cvm::to_str(nx_i)+".\n"); 00157 } 00158 00159 mult = mult_i; 00160 00161 data.clear(); 00162 00163 nx = nx_i; 00164 nd = nx.size(); 00165 00166 nxc.resize(nd); 00167 00168 // setup dimensions 00169 nt = mult; 00170 for (int i = nd-1; i >= 0; i--) { 00171 if (nx[i] <= 0) { 00172 cvm::error("Error: providing an invalid number of grid points, "+ 00173 cvm::to_str(nx[i])+".\n", COLVARS_BUG_ERROR); 00174 return COLVARS_ERROR; 00175 } 00176 nxc[i] = nt; 00177 nt *= nx[i]; 00178 } 00179 00180 if (cvm::debug()) { 00181 cvm::log("Total number of grid elements = "+cvm::to_str(nt)+".\n"); 00182 } 00183 00184 data.reserve(nt); 00185 data.assign(nt, t); 00186 00187 return COLVARS_OK; 00188 } 00189 00191 int setup() 00192 { 00193 return setup(this->nx, T(), this->mult); 00194 } 00195 00197 void reset(T const &t = T()) 00198 { 00199 data.assign(nt, t); 00200 } 00201 00202 00204 colvar_grid() : has_data(false) 00205 { 00206 nd = nt = 0; 00207 mult = 1; 00208 has_parent_data = false; 00209 this->setup(); 00210 } 00211 00213 virtual ~colvar_grid() 00214 {} 00215 00219 colvar_grid(colvar_grid<T> const &g) : colvarparse(), 00220 nd(g.nd), 00221 nx(g.nx), 00222 mult(g.mult), 00223 data(), 00224 cv(g.cv), 00225 use_actual_value(g.use_actual_value), 00226 lower_boundaries(g.lower_boundaries), 00227 upper_boundaries(g.upper_boundaries), 00228 periodic(g.periodic), 00229 hard_lower_boundaries(g.hard_lower_boundaries), 00230 hard_upper_boundaries(g.hard_upper_boundaries), 00231 widths(g.widths), 00232 has_parent_data(false), 00233 has_data(false) 00234 {} 00235 00240 colvar_grid(std::vector<int> const &nx_i, 00241 T const &t = T(), 00242 size_t mult_i = 1) 00243 : has_parent_data(false), has_data(false) 00244 { 00245 this->setup(nx_i, t, mult_i); 00246 } 00247 00251 colvar_grid(std::vector<colvar *> const &colvars, 00252 T const &t = T(), 00253 size_t mult_i = 1, 00254 bool add_extra_bin = false) 00255 : has_parent_data(false), has_data(false) 00256 { 00257 (void) t; 00258 this->init_from_colvars(colvars, mult_i, add_extra_bin); 00259 } 00260 00261 int init_from_colvars(std::vector<colvar *> const &colvars, 00262 size_t mult_i = 1, 00263 bool add_extra_bin = false) 00264 { 00265 if (cvm::debug()) { 00266 cvm::log("Reading grid configuration from collective variables.\n"); 00267 } 00268 00269 cv = colvars; 00270 nd = colvars.size(); 00271 mult = mult_i; 00272 00273 size_t i; 00274 00275 if (cvm::debug()) { 00276 cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+ 00277 " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n"); 00278 } 00279 00280 for (i = 0; i < cv.size(); i++) { 00281 00282 if (cv[i]->value().type() != colvarvalue::type_scalar) { 00283 cvm::error("Colvar grids can only be automatically " 00284 "constructed for scalar variables. " 00285 "ABF and histogram can not be used; metadynamics " 00286 "can be used with useGrids disabled.\n", COLVARS_INPUT_ERROR); 00287 return COLVARS_ERROR; 00288 } 00289 00290 if (cv[i]->width <= 0.0) { 00291 cvm::error("Tried to initialize a grid on a " 00292 "variable with negative or zero width.\n", COLVARS_INPUT_ERROR); 00293 return COLVARS_ERROR; 00294 } 00295 00296 widths.push_back(cv[i]->width); 00297 hard_lower_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_lower_boundary)); 00298 hard_upper_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_upper_boundary)); 00299 periodic.push_back(cv[i]->periodic_boundaries()); 00300 00301 // By default, get reported colvar value (for extended Lagrangian colvars) 00302 use_actual_value.push_back(false); 00303 00304 // except if a colvar is specified twice in a row 00305 // then the first instance is the actual value 00306 // For histograms of extended-system coordinates 00307 if (i > 0 && cv[i-1] == cv[i]) { 00308 use_actual_value[i-1] = true; 00309 } 00310 00311 if (add_extra_bin) { 00312 if (periodic[i]) { 00313 // Shift the grid by half the bin width (values at edges instead of center of bins) 00314 lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]); 00315 upper_boundaries.push_back(cv[i]->upper_boundary.real_value - 0.5 * widths[i]); 00316 } else { 00317 // Make this grid larger by one bin width 00318 lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]); 00319 upper_boundaries.push_back(cv[i]->upper_boundary.real_value + 0.5 * widths[i]); 00320 } 00321 } else { 00322 lower_boundaries.push_back(cv[i]->lower_boundary); 00323 upper_boundaries.push_back(cv[i]->upper_boundary); 00324 } 00325 } 00326 00327 this->init_from_boundaries(); 00328 return this->setup(); 00329 } 00330 00331 int init_from_boundaries() 00332 { 00333 if (cvm::debug()) { 00334 cvm::log("Configuring grid dimensions from colvars boundaries.\n"); 00335 } 00336 00337 // these will have to be updated 00338 nx.clear(); 00339 nxc.clear(); 00340 nt = 0; 00341 00342 for (size_t i = 0; i < lower_boundaries.size(); i++) { 00343 // Re-compute periodicity using current grid boundaries 00344 periodic[i] = cv[i]->periodic_boundaries(lower_boundaries[i].real_value, 00345 upper_boundaries[i].real_value); 00346 00347 cvm::real nbins = ( upper_boundaries[i].real_value - 00348 lower_boundaries[i].real_value ) / widths[i]; 00349 int nbins_round = (int)(nbins+0.5); 00350 00351 if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) { 00352 cvm::log("Warning: grid interval("+ 00353 cvm::to_str(lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+ 00354 cvm::to_str(upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+ 00355 ") is not commensurate to its bin width("+ 00356 cvm::to_str(widths[i], cvm::cv_width, cvm::cv_prec)+").\n"); 00357 upper_boundaries[i].real_value = lower_boundaries[i].real_value + 00358 (nbins_round * widths[i]); 00359 } 00360 00361 if (cvm::debug()) 00362 cvm::log("Number of points is "+cvm::to_str((int) nbins_round)+ 00363 " for the colvar no. "+cvm::to_str(i+1)+".\n"); 00364 00365 nx.push_back(nbins_round); 00366 } 00367 00368 return COLVARS_OK; 00369 } 00370 00373 inline void wrap(std::vector<int> & ix) const 00374 { 00375 for (size_t i = 0; i < nd; i++) { 00376 if (periodic[i]) { 00377 ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result 00378 } else { 00379 if (ix[i] < 0 || ix[i] >= nx[i]) { 00380 cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: " 00381 + cvm::to_str(ix), COLVARS_BUG_ERROR); 00382 return; 00383 } 00384 } 00385 } 00386 } 00387 00390 inline bool wrap_edge(std::vector<int> & ix) const 00391 { 00392 bool edge = false; 00393 for (size_t i = 0; i < nd; i++) { 00394 if (periodic[i]) { 00395 ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result 00396 } else if (ix[i] == -1 || ix[i] == nx[i]) { 00397 edge = true; 00398 } 00399 } 00400 return edge; 00401 } 00402 00404 inline int current_bin_scalar(int const i) const 00405 { 00406 return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); 00407 } 00408 00411 inline int current_bin_scalar_bound(int const i) const 00412 { 00413 return value_to_bin_scalar_bound(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); 00414 } 00415 00417 inline int current_bin_scalar(int const i, int const iv) const 00418 { 00419 return value_to_bin_scalar(use_actual_value[i] ? 00420 cv[i]->actual_value().vector1d_value[iv] : 00421 cv[i]->value().vector1d_value[iv], i); 00422 } 00423 00426 inline int value_to_bin_scalar(colvarvalue const &value, const int i) const 00427 { 00428 return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); 00429 } 00430 00432 inline cvm::real current_bin_scalar_fraction(int const i) const 00433 { 00434 return value_to_bin_scalar_fraction(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); 00435 } 00436 00439 inline cvm::real value_to_bin_scalar_fraction(colvarvalue const &value, const int i) const 00440 { 00441 cvm::real x = (value.real_value - lower_boundaries[i].real_value) / widths[i]; 00442 return x - cvm::floor(x); 00443 } 00444 00447 inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const 00448 { 00449 int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); 00450 if (bin_index < 0) bin_index=0; 00451 if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1; 00452 return (int) bin_index; 00453 } 00454 00456 inline int value_to_bin_scalar(colvarvalue const &value, 00457 colvarvalue const &new_offset, 00458 cvm::real const &new_width) const 00459 { 00460 return (int) cvm::floor( (value.real_value - new_offset.real_value) / new_width ); 00461 } 00462 00465 inline colvarvalue bin_to_value_scalar(int const &i_bin, int const i) const 00466 { 00467 return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin); 00468 } 00469 00471 inline colvarvalue bin_to_value_scalar(int const &i_bin, 00472 colvarvalue const &new_offset, 00473 cvm::real const &new_width) const 00474 { 00475 return new_offset.real_value + new_width * (0.5 + i_bin); 00476 } 00477 00479 inline void set_value(std::vector<int> const &ix, 00480 T const &t, 00481 size_t const &imult = 0) 00482 { 00483 data[this->address(ix)+imult] = t; 00484 has_data = true; 00485 } 00486 00488 inline void set_value(size_t i, T const &t) 00489 { 00490 data[i] = t; 00491 } 00492 00497 void delta_grid(colvar_grid<T> const &other_grid) 00498 { 00499 00500 if (other_grid.multiplicity() != this->multiplicity()) { 00501 cvm::error("Error: trying to subtract two grids with " 00502 "different multiplicity.\n"); 00503 return; 00504 } 00505 00506 if (other_grid.data.size() != this->data.size()) { 00507 cvm::error("Error: trying to subtract two grids with " 00508 "different size.\n"); 00509 return; 00510 } 00511 00512 for (size_t i = 0; i < data.size(); i++) { 00513 data[i] = other_grid.data[i] - data[i]; 00514 } 00515 has_data = true; 00516 } 00517 00521 void copy_grid(colvar_grid<T> const &other_grid) 00522 { 00523 if (other_grid.multiplicity() != this->multiplicity()) { 00524 cvm::error("Error: trying to copy two grids with " 00525 "different multiplicity.\n"); 00526 return; 00527 } 00528 00529 if (other_grid.data.size() != this->data.size()) { 00530 cvm::error("Error: trying to copy two grids with " 00531 "different size.\n"); 00532 return; 00533 } 00534 00535 00536 for (size_t i = 0; i < data.size(); i++) { 00537 data[i] = other_grid.data[i]; 00538 } 00539 has_data = true; 00540 } 00541 00544 void raw_data_out(T* out_data) const 00545 { 00546 for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i]; 00547 } 00548 void raw_data_out(std::vector<T>& out_data) const 00549 { 00550 out_data = data; 00551 } 00553 void raw_data_in(const T* in_data) 00554 { 00555 for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i]; 00556 has_data = true; 00557 } 00558 void raw_data_in(const std::vector<T>& in_data) 00559 { 00560 data = in_data; 00561 has_data = true; 00562 } 00564 size_t raw_data_num() const { return data.size(); } 00565 00566 00569 inline T const & value(std::vector<int> const &ix, 00570 size_t const &imult = 0) const 00571 { 00572 return data[this->address(ix) + imult]; 00573 } 00574 00576 inline T const & value(size_t i) const 00577 { 00578 return data[i]; 00579 } 00580 00582 inline void add_constant(T const &t) 00583 { 00584 for (size_t i = 0; i < nt; i++) 00585 data[i] += t; 00586 has_data = true; 00587 } 00588 00590 inline void multiply_constant(cvm::real const &a) 00591 { 00592 for (size_t i = 0; i < nt; i++) 00593 data[i] *= a; 00594 } 00595 00597 inline void remove_small_values(cvm::real const &a) 00598 { 00599 for (size_t i = 0; i < nt; i++) 00600 if(data[i]<a) data[i] = a; 00601 } 00602 00603 00606 inline std::vector<int> const get_colvars_index(std::vector<colvarvalue> const &values) const 00607 { 00608 std::vector<int> index = new_index(); 00609 for (size_t i = 0; i < nd; i++) { 00610 index[i] = value_to_bin_scalar(values[i], i); 00611 } 00612 return index; 00613 } 00614 00617 inline std::vector<int> const get_colvars_index() const 00618 { 00619 std::vector<int> index = new_index(); 00620 for (size_t i = 0; i < nd; i++) { 00621 index[i] = current_bin_scalar(i); 00622 } 00623 return index; 00624 } 00625 00628 inline std::vector<int> const get_colvars_index_bound() const 00629 { 00630 std::vector<int> index = new_index(); 00631 for (size_t i = 0; i < nd; i++) { 00632 index[i] = current_bin_scalar_bound(i); 00633 } 00634 return index; 00635 } 00636 00640 inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values, 00641 bool skip_hard_boundaries = false) 00642 { 00643 cvm::real minimum = 1.0E+16; 00644 for (size_t i = 0; i < nd; i++) { 00645 00646 if (periodic[i]) continue; 00647 00648 cvm::real dl = cvm::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i]; 00649 cvm::real du = cvm::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i]; 00650 00651 if (values[i].real_value < lower_boundaries[i]) 00652 dl *= -1.0; 00653 if (values[i].real_value > upper_boundaries[i]) 00654 du *= -1.0; 00655 00656 if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum)) 00657 minimum = dl; 00658 if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum)) 00659 minimum = du; 00660 } 00661 00662 return minimum; 00663 } 00664 00665 00670 void map_grid(colvar_grid<T> const &other_grid) 00671 { 00672 if (other_grid.multiplicity() != this->multiplicity()) { 00673 cvm::error("Error: trying to merge two grids with values of " 00674 "different multiplicity.\n"); 00675 return; 00676 } 00677 00678 std::vector<colvarvalue> const &gb = this->lower_boundaries; 00679 std::vector<cvm::real> const &gw = this->widths; 00680 std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries; 00681 std::vector<cvm::real> const &ogw = other_grid.widths; 00682 00683 std::vector<int> ix = this->new_index(); 00684 std::vector<int> oix = other_grid.new_index(); 00685 00686 if (cvm::debug()) 00687 cvm::log("Remapping grid...\n"); 00688 for ( ; this->index_ok(ix); this->incr(ix)) { 00689 00690 for (size_t i = 0; i < nd; i++) { 00691 oix[i] = 00692 value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]), 00693 ogb[i], 00694 ogw[i]); 00695 } 00696 00697 if (! other_grid.index_ok(oix)) { 00698 continue; 00699 } 00700 00701 for (size_t im = 0; im < mult; im++) { 00702 this->set_value(ix, other_grid.value(oix, im), im); 00703 } 00704 } 00705 00706 has_data = true; 00707 if (cvm::debug()) 00708 cvm::log("Remapping done.\n"); 00709 } 00710 00713 void add_grid(colvar_grid<T> const &other_grid, 00714 cvm::real scale_factor = 1.0) 00715 { 00716 if (other_grid.multiplicity() != this->multiplicity()) { 00717 cvm::error("Error: trying to sum togetehr two grids with values of " 00718 "different multiplicity.\n"); 00719 return; 00720 } 00721 if (scale_factor != 1.0) 00722 for (size_t i = 0; i < data.size(); i++) { 00723 data[i] += static_cast<T>(scale_factor * other_grid.data[i]); 00724 } 00725 else 00726 // skip multiplication if possible 00727 for (size_t i = 0; i < data.size(); i++) { 00728 data[i] += other_grid.data[i]; 00729 } 00730 has_data = true; 00731 } 00732 00735 virtual T value_output(std::vector<int> const &ix, 00736 size_t const &imult = 0) const 00737 { 00738 return value(ix, imult); 00739 } 00740 00744 virtual void value_input(std::vector<int> const &ix, 00745 T const &t, 00746 size_t const &imult = 0, 00747 bool add = false) 00748 { 00749 if ( add ) 00750 data[address(ix) + imult] += t; 00751 else 00752 data[address(ix) + imult] = t; 00753 has_data = true; 00754 } 00755 00756 // /// Get the pointer to the binned value indexed by ix 00757 // inline T const *value_p (std::vector<int> const &ix) 00758 // { 00759 // return &(data[address (ix)]); 00760 // } 00761 00764 inline std::vector<int> const new_index() const 00765 { 00766 return std::vector<int> (nd, 0); 00767 } 00768 00771 inline bool index_ok(std::vector<int> const &ix) const 00772 { 00773 for (size_t i = 0; i < nd; i++) { 00774 if ( (ix[i] < 0) || (ix[i] >= int(nx[i])) ) 00775 return false; 00776 } 00777 return true; 00778 } 00779 00782 inline void incr(std::vector<int> &ix) const 00783 { 00784 for (int i = ix.size()-1; i >= 0; i--) { 00785 00786 ix[i]++; 00787 00788 if (ix[i] >= nx[i]) { 00789 00790 if (i > 0) { 00791 ix[i] = 0; 00792 continue; 00793 } else { 00794 // this is the last iteration, a non-valid index is being 00795 // set for the outer index, which will be caught by 00796 // index_ok() 00797 ix[0] = nx[0]; 00798 return; 00799 } 00800 } else { 00801 return; 00802 } 00803 } 00804 } 00805 00807 std::ostream & write_params(std::ostream &os) 00808 { 00809 size_t i; 00810 os << "grid_parameters {\n n_colvars " << nd << "\n"; 00811 00812 os << " lower_boundaries "; 00813 for (i = 0; i < nd; i++) 00814 os << " " << lower_boundaries[i]; 00815 os << "\n"; 00816 00817 os << " upper_boundaries "; 00818 for (i = 0; i < nd; i++) 00819 os << " " << upper_boundaries[i]; 00820 os << "\n"; 00821 00822 os << " widths "; 00823 for (i = 0; i < nd; i++) 00824 os << " " << widths[i]; 00825 os << "\n"; 00826 00827 os << " sizes "; 00828 for (i = 0; i < nd; i++) 00829 os << " " << nx[i]; 00830 os << "\n"; 00831 00832 os << "}\n"; 00833 return os; 00834 } 00835 00837 int parse_params(std::string const &conf, 00838 colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal) 00839 { 00840 if (cvm::debug()) cvm::log("Reading grid configuration from string.\n"); 00841 00842 std::vector<int> old_nx = nx; 00843 std::vector<colvarvalue> old_lb = lower_boundaries; 00844 std::vector<colvarvalue> old_ub = upper_boundaries; 00845 std::vector<cvm::real> old_w = widths; 00846 00847 { 00848 size_t nd_in = 0; 00849 // this is only used in state files 00850 colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent); 00851 if (nd_in != nd) { 00852 cvm::error("Error: trying to read data for a grid " 00853 "that contains a different number of colvars ("+ 00854 cvm::to_str(nd_in)+") than the grid defined " 00855 "in the configuration file("+cvm::to_str(nd)+ 00856 ").\n"); 00857 return COLVARS_ERROR; 00858 } 00859 } 00860 00861 // underscore keywords are used in state file 00862 colvarparse::get_keyval(conf, "lower_boundaries", 00863 lower_boundaries, lower_boundaries, colvarparse::parse_silent); 00864 colvarparse::get_keyval(conf, "upper_boundaries", 00865 upper_boundaries, upper_boundaries, colvarparse::parse_silent); 00866 00867 // camel case keywords are used in config file 00868 colvarparse::get_keyval(conf, "lowerBoundaries", 00869 lower_boundaries, lower_boundaries, parse_mode); 00870 colvarparse::get_keyval(conf, "upperBoundaries", 00871 upper_boundaries, upper_boundaries, parse_mode); 00872 00873 colvarparse::get_keyval(conf, "widths", widths, widths, parse_mode); 00874 00875 // only used in state file 00876 colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent); 00877 00878 if (nd < lower_boundaries.size()) nd = lower_boundaries.size(); 00879 00880 if (! use_actual_value.size()) use_actual_value.assign(nd, false); 00881 if (! periodic.size()) periodic.assign(nd, false); 00882 if (! widths.size()) widths.assign(nd, 1.0); 00883 00884 cvm::real eps = 1.e-10; 00885 00886 bool new_params = false; 00887 if (old_nx.size()) { 00888 for (size_t i = 0; i < nd; i++) { 00889 if (old_nx[i] != nx[i] || 00890 cvm::sqrt(cv[i]->dist2(old_lb[i], lower_boundaries[i])) > eps || 00891 cvm::sqrt(cv[i]->dist2(old_ub[i], upper_boundaries[i])) > eps || 00892 cvm::fabs(old_w[i] - widths[i]) > eps) { 00893 new_params = true; 00894 } 00895 } 00896 } else { 00897 new_params = true; 00898 } 00899 00900 // reallocate the array in case the grid params have just changed 00901 if (new_params) { 00902 init_from_boundaries(); 00903 // data.clear(); // no longer needed: setup calls clear() 00904 return this->setup(nx, T(), mult); 00905 } 00906 00907 return COLVARS_OK; 00908 } 00909 00913 void check_consistency() 00914 { 00915 for (size_t i = 0; i < nd; i++) { 00916 if ( (cvm::sqrt(cv[i]->dist2(cv[i]->lower_boundary, 00917 lower_boundaries[i])) > 1.0E-10) || 00918 (cvm::sqrt(cv[i]->dist2(cv[i]->upper_boundary, 00919 upper_boundaries[i])) > 1.0E-10) || 00920 (cvm::sqrt(cv[i]->dist2(cv[i]->width, 00921 widths[i])) > 1.0E-10) ) { 00922 cvm::error("Error: restart information for a grid is " 00923 "inconsistent with that of its colvars.\n"); 00924 return; 00925 } 00926 } 00927 } 00928 00929 00932 void check_consistency(colvar_grid<T> const &other_grid) 00933 { 00934 for (size_t i = 0; i < nd; i++) { 00935 // we skip dist2(), because periodicities and the like should 00936 // matter: boundaries should be EXACTLY the same (otherwise, 00937 // map_grid() should be used) 00938 if ( (cvm::fabs(other_grid.lower_boundaries[i] - 00939 lower_boundaries[i]) > 1.0E-10) || 00940 (cvm::fabs(other_grid.upper_boundaries[i] - 00941 upper_boundaries[i]) > 1.0E-10) || 00942 (cvm::fabs(other_grid.widths[i] - 00943 widths[i]) > 1.0E-10) || 00944 (data.size() != other_grid.data.size()) ) { 00945 cvm::error("Error: inconsistency between " 00946 "two grids that are supposed to be equal, " 00947 "aside from the data stored.\n"); 00948 return; 00949 } 00950 } 00951 } 00952 00953 00955 std::istream & read_restart(std::istream &is) 00956 { 00957 std::streampos const start_pos = is.tellg(); 00958 std::string key, conf; 00959 if ((is >> key) && (key == std::string("grid_parameters"))) { 00960 is.seekg(start_pos, std::ios::beg); 00961 is >> colvarparse::read_block("grid_parameters", &conf); 00962 parse_params(conf, colvarparse::parse_silent); 00963 } else { 00964 cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n"); 00965 is.seekg(start_pos, std::ios::beg); 00966 } 00967 read_raw(is); 00968 return is; 00969 } 00970 00972 std::ostream & write_restart(std::ostream &os) 00973 { 00974 write_params(os); 00975 write_raw(os); 00976 return os; 00977 } 00978 00979 00983 std::ostream & write_raw(std::ostream &os, 00984 size_t const buf_size = 3) const 00985 { 00986 std::streamsize const w = os.width(); 00987 std::streamsize const p = os.precision(); 00988 00989 std::vector<int> ix = new_index(); 00990 size_t count = 0; 00991 for ( ; index_ok(ix); incr(ix)) { 00992 for (size_t imult = 0; imult < mult; imult++) { 00993 os << " " 00994 << std::setw(w) << std::setprecision(p) 00995 << value_output(ix, imult); 00996 if (((++count) % buf_size) == 0) 00997 os << "\n"; 00998 } 00999 } 01000 // write a final newline only if buffer is not empty 01001 if ((count % buf_size) != 0) 01002 os << "\n"; 01003 01004 return os; 01005 } 01006 01008 std::istream & read_raw(std::istream &is) 01009 { 01010 std::streampos const start_pos = is.tellg(); 01011 01012 for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) { 01013 for (size_t imult = 0; imult < mult; imult++) { 01014 T new_value; 01015 if (is >> new_value) { 01016 value_input(ix, new_value, imult); 01017 } else { 01018 is.clear(); 01019 is.seekg(start_pos, std::ios::beg); 01020 is.setstate(std::ios::failbit); 01021 cvm::error("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n"); 01022 return is; 01023 } 01024 } 01025 } 01026 01027 has_data = true; 01028 return is; 01029 } 01030 01032 std::istream & read_multicol(std::istream &is, bool add = false); 01033 01035 int read_multicol(std::string const &filename, 01036 std::string description = "grid file", 01037 bool add = false); 01038 01040 std::ostream & write_multicol(std::ostream &os) const; 01041 01043 int write_multicol(std::string const &filename, 01044 std::string description = "grid file") const; 01045 01047 std::ostream & write_opendx(std::ostream &os) const; 01048 01050 int write_opendx(std::string const &filename, 01051 std::string description = "grid file") const; 01052 }; 01053 01054 01055 01058 class colvar_grid_count : public colvar_grid<size_t> 01059 { 01060 public: 01061 01063 colvar_grid_count(); 01064 01066 virtual ~colvar_grid_count() 01067 {} 01068 01070 colvar_grid_count(std::vector<int> const &nx_i, 01071 size_t const &def_count = 0); 01072 01074 colvar_grid_count(std::vector<colvar *> &colvars, 01075 size_t const &def_count = 0, 01076 bool add_extra_bin = false); 01077 01079 inline void incr_count(std::vector<int> const &ix) 01080 { 01081 ++(data[this->address(ix)]); 01082 } 01083 01085 inline size_t const & new_count(std::vector<int> const &ix, 01086 size_t const &imult = 0) 01087 { 01088 return new_data[address(ix) + imult]; 01089 } 01090 01092 std::istream & read_multicol(std::istream &is, bool add = false); 01093 01095 int read_multicol(std::string const &filename, 01096 std::string description = "grid file", 01097 bool add = false); 01098 01100 std::ostream & write_multicol(std::ostream &os) const; 01101 01103 int write_multicol(std::string const &filename, 01104 std::string description = "grid file") const; 01105 01107 std::ostream & write_opendx(std::ostream &os) const; 01108 01110 int write_opendx(std::string const &filename, 01111 std::string description = "grid file") const; 01112 01114 virtual void value_input(std::vector<int> const &ix, 01115 size_t const &t, 01116 size_t const &imult = 0, 01117 bool add = false) 01118 { 01119 (void) imult; 01120 if (add) { 01121 data[address(ix)] += t; 01122 if (this->has_parent_data) { 01123 // save newly read data for inputting parent grid 01124 new_data[address(ix)] = t; 01125 } 01126 } else { 01127 data[address(ix)] = t; 01128 } 01129 has_data = true; 01130 } 01131 01134 inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0, 01135 int n = 0) 01136 { 01137 cvm::real A0, A1, A2; 01138 std::vector<int> ix = ix0; 01139 01140 // TODO this can be rewritten more concisely with wrap_edge() 01141 if (periodic[n]) { 01142 ix[n]--; wrap(ix); 01143 A0 = value(ix); 01144 ix = ix0; 01145 ix[n]++; wrap(ix); 01146 A1 = value(ix); 01147 if (A0 * A1 == 0) { 01148 return 0.; // can't handle empty bins 01149 } else { 01150 return (cvm::logn(A1) - cvm::logn(A0)) 01151 / (widths[n] * 2.); 01152 } 01153 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge 01154 ix[n]--; 01155 A0 = value(ix); 01156 ix = ix0; 01157 ix[n]++; 01158 A1 = value(ix); 01159 if (A0 * A1 == 0) { 01160 return 0.; // can't handle empty bins 01161 } else { 01162 return (cvm::logn(A1) - cvm::logn(A0)) 01163 / (widths[n] * 2.); 01164 } 01165 } else { 01166 // edge: use 2nd order derivative 01167 int increment = (ix[n] == 0 ? 1 : -1); 01168 // move right from left edge, or the other way around 01169 A0 = value(ix); 01170 ix[n] += increment; A1 = value(ix); 01171 ix[n] += increment; A2 = value(ix); 01172 if (A0 * A1 * A2 == 0) { 01173 return 0.; // can't handle empty bins 01174 } else { 01175 return (-1.5 * cvm::logn(A0) + 2. * cvm::logn(A1) 01176 - 0.5 * cvm::logn(A2)) * increment / widths[n]; 01177 } 01178 } 01179 } 01180 01183 inline cvm::real gradient_finite_diff(const std::vector<int> &ix0, 01184 int n = 0) 01185 { 01186 cvm::real A0, A1, A2; 01187 std::vector<int> ix = ix0; 01188 01189 // FIXME this can be rewritten more concisely with wrap_edge() 01190 if (periodic[n]) { 01191 ix[n]--; wrap(ix); 01192 A0 = value(ix); 01193 ix = ix0; 01194 ix[n]++; wrap(ix); 01195 A1 = value(ix); 01196 if (A0 * A1 == 0) { 01197 return 0.; // can't handle empty bins 01198 } else { 01199 return (A1 - A0) / (widths[n] * 2.); 01200 } 01201 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge 01202 ix[n]--; 01203 A0 = value(ix); 01204 ix = ix0; 01205 ix[n]++; 01206 A1 = value(ix); 01207 if (A0 * A1 == 0) { 01208 return 0.; // can't handle empty bins 01209 } else { 01210 return (A1 - A0) / (widths[n] * 2.); 01211 } 01212 } else { 01213 // edge: use 2nd order derivative 01214 int increment = (ix[n] == 0 ? 1 : -1); 01215 // move right from left edge, or the other way around 01216 A0 = value(ix); 01217 ix[n] += increment; A1 = value(ix); 01218 ix[n] += increment; A2 = value(ix); 01219 return (-1.5 * A0 + 2. * A1 01220 - 0.5 * A2) * increment / widths[n]; 01221 } 01222 } 01223 }; 01224 01225 01227 class colvar_grid_scalar : public colvar_grid<cvm::real> 01228 { 01229 public: 01230 01233 colvar_grid_count *samples; 01234 01236 colvar_grid_scalar(); 01237 01239 colvar_grid_scalar(colvar_grid_scalar const &g); 01240 01242 virtual ~colvar_grid_scalar(); 01243 01245 colvar_grid_scalar(std::vector<int> const &nx_i); 01246 01248 colvar_grid_scalar(std::vector<colvar *> &colvars, 01249 bool add_extra_bin = false); 01250 01252 inline void acc_value(std::vector<int> const &ix, 01253 cvm::real const &new_value, 01254 size_t const &imult = 0) 01255 { 01256 (void) imult; 01257 // only legal value of imult here is 0 01258 data[address(ix)] += new_value; 01259 if (samples) 01260 samples->incr_count(ix); 01261 has_data = true; 01262 } 01263 01265 std::istream & read_multicol(std::istream &is, bool add = false); 01266 01268 int read_multicol(std::string const &filename, 01269 std::string description = "grid file", 01270 bool add = false); 01271 01273 std::ostream & write_multicol(std::ostream &os) const; 01274 01276 int write_multicol(std::string const &filename, 01277 std::string description = "grid file") const; 01278 01280 std::ostream & write_opendx(std::ostream &os) const; 01281 01283 int write_opendx(std::string const &filename, 01284 std::string description = "grid file") const; 01285 01290 inline void vector_gradient_finite_diff( const std::vector<int> &ix0, std::vector<cvm::real> &grad) 01291 { 01292 cvm::real A0, A1; 01293 std::vector<int> ix; 01294 size_t i, j, k, n; 01295 01296 if (nd == 2) { 01297 for (n = 0; n < 2; n++) { 01298 ix = ix0; 01299 A0 = value(ix); 01300 ix[n]++; wrap(ix); 01301 A1 = value(ix); 01302 ix[1-n]++; wrap(ix); 01303 A1 += value(ix); 01304 ix[n]--; wrap(ix); 01305 A0 += value(ix); 01306 grad[n] = 0.5 * (A1 - A0) / widths[n]; 01307 } 01308 } else if (nd == 3) { 01309 01310 cvm::real p[8]; // potential values within cube, indexed in binary (4 i + 2 j + k) 01311 ix = ix0; 01312 int index = 0; 01313 for (i = 0; i<2; i++) { 01314 ix[1] = ix0[1]; 01315 for (j = 0; j<2; j++) { 01316 ix[2] = ix0[2]; 01317 for (k = 0; k<2; k++) { 01318 wrap(ix); 01319 p[index++] = value(ix); 01320 ix[2]++; 01321 } 01322 ix[1]++; 01323 } 01324 ix[0]++; 01325 } 01326 01327 // The following would be easier to read using binary literals 01328 // 100 101 110 111 000 001 010 011 01329 grad[0] = 0.25 * ((p[4] + p[5] + p[6] + p[7]) - (p[0] + p[1] + p[2] + p[3])) / widths[0]; 01330 // 010 011 110 111 000 001 100 101 01331 grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[1]; 01332 // 001 011 101 111 000 010 100 110 01333 grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[2]; 01334 } else { 01335 cvm::error("Finite differences available in dimension 2 and 3 only."); 01336 } 01337 } 01338 01341 inline cvm::real gradient_finite_diff(const std::vector<int> &ix0, 01342 int n = 0) 01343 { 01344 cvm::real A0, A1, A2; 01345 std::vector<int> ix = ix0; 01346 01347 // FIXME this can be rewritten more concisely with wrap_edge() 01348 if (periodic[n]) { 01349 ix[n]--; wrap(ix); 01350 A0 = value(ix); 01351 ix = ix0; 01352 ix[n]++; wrap(ix); 01353 A1 = value(ix); 01354 if (A0 * A1 == 0) { 01355 return 0.; // can't handle empty bins 01356 } else { 01357 return (A1 - A0) / (widths[n] * 2.); 01358 } 01359 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge 01360 ix[n]--; 01361 A0 = value(ix); 01362 ix = ix0; 01363 ix[n]++; 01364 A1 = value(ix); 01365 if (A0 * A1 == 0) { 01366 return 0.; // can't handle empty bins 01367 } else { 01368 return (A1 - A0) / (widths[n] * 2.); 01369 } 01370 } else { 01371 // edge: use 2nd order derivative 01372 int increment = (ix[n] == 0 ? 1 : -1); 01373 // move right from left edge, or the other way around 01374 A0 = value(ix); 01375 ix[n] += increment; A1 = value(ix); 01376 ix[n] += increment; A2 = value(ix); 01377 return (-1.5 * A0 + 2. * A1 01378 - 0.5 * A2) * increment / widths[n]; 01379 } 01380 } 01381 01384 virtual cvm::real value_output(std::vector<int> const &ix, 01385 size_t const &imult = 0) const 01386 { 01387 if (imult > 0) { 01388 cvm::error("Error: trying to access a component " 01389 "larger than 1 in a scalar data grid.\n"); 01390 return 0.; 01391 } 01392 if (samples) { 01393 return (samples->value(ix) > 0) ? 01394 (data[address(ix)] / cvm::real(samples->value(ix))) : 01395 0.0; 01396 } else { 01397 return data[address(ix)]; 01398 } 01399 } 01400 01402 virtual void value_input(std::vector<int> const &ix, 01403 cvm::real const &new_value, 01404 size_t const &imult = 0, 01405 bool add = false) 01406 { 01407 if (imult > 0) { 01408 cvm::error("Error: trying to access a component " 01409 "larger than 1 in a scalar data grid.\n"); 01410 return; 01411 } 01412 if (add) { 01413 if (samples) 01414 data[address(ix)] += new_value * samples->new_count(ix); 01415 else 01416 data[address(ix)] += new_value; 01417 } else { 01418 if (samples) 01419 data[address(ix)] = new_value * samples->value(ix); 01420 else 01421 data[address(ix)] = new_value; 01422 } 01423 has_data = true; 01424 } 01425 01427 cvm::real maximum_value() const; 01428 01430 cvm::real minimum_value() const; 01431 01433 cvm::real minimum_pos_value() const; 01434 01436 cvm::real integral() const; 01437 01440 cvm::real entropy() const; 01441 }; 01442 01443 01444 01446 class colvar_grid_gradient : public colvar_grid<cvm::real> 01447 { 01448 public: 01449 01452 colvar_grid_count *samples; 01453 01456 colvar_grid_scalar *weights; 01457 01459 colvar_grid_gradient(); 01460 01462 virtual ~colvar_grid_gradient() 01463 {} 01464 01466 colvar_grid_gradient(std::vector<int> const &nx_i); 01467 01469 colvar_grid_gradient(std::vector<colvar *> &colvars); 01470 01472 colvar_grid_gradient(std::string &filename); 01473 01475 virtual std::istream & read_multicol(std::istream &is, bool add = false); 01476 01478 virtual int read_multicol(std::string const &filename, 01479 std::string description = "grid file", 01480 bool add = false); 01481 01483 virtual std::ostream & write_multicol(std::ostream &os) const; 01484 01486 virtual int write_multicol(std::string const &filename, 01487 std::string description = "grid file") const; 01488 01490 virtual std::ostream & write_opendx(std::ostream &os) const; 01491 01493 virtual int write_opendx(std::string const &filename, 01494 std::string description = "grid file") const; 01495 01497 inline void vector_value(std::vector<int> const &ix, std::vector<cvm::real> &v) const 01498 { 01499 cvm::real const * p = &value(ix); 01500 if (samples) { 01501 int count = samples->value(ix); 01502 if (count) { 01503 cvm::real invcount = 1.0 / count; 01504 for (size_t i = 0; i < mult; i++) { 01505 v[i] = invcount * p[i]; 01506 } 01507 } else { 01508 for (size_t i = 0; i < mult; i++) { 01509 v[i] = 0.0; 01510 } 01511 } 01512 } else { 01513 for (size_t i = 0; i < mult; i++) { 01514 v[i] = p[i]; 01515 } 01516 } 01517 } 01518 01520 inline void acc_value(std::vector<int> const &ix, std::vector<colvarvalue> const &values) { 01521 for (size_t imult = 0; imult < mult; imult++) { 01522 data[address(ix) + imult] += values[imult].real_value; 01523 } 01524 if (samples) 01525 samples->incr_count(ix); 01526 } 01527 01530 inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) { 01531 for (size_t imult = 0; imult < mult; imult++) { 01532 data[address(ix) + imult] -= forces[imult]; 01533 } 01534 if (samples) 01535 samples->incr_count(ix); 01536 } 01537 01540 inline void acc_force_weighted(std::vector<int> const &ix, 01541 cvm::real const *forces, 01542 cvm::real weight) { 01543 for (size_t imult = 0; imult < mult; imult++) { 01544 data[address(ix) + imult] -= forces[imult] * weight; 01545 } 01546 weights->acc_value(ix, weight); 01547 } 01548 01551 virtual cvm::real value_output(std::vector<int> const &ix, 01552 size_t const &imult = 0) const 01553 { 01554 if (samples) 01555 return (samples->value(ix) > 0) ? 01556 (data[address(ix) + imult] / cvm::real(samples->value(ix))) : 01557 0.0; 01558 else 01559 return data[address(ix) + imult]; 01560 } 01561 01565 virtual void value_input(std::vector<int> const &ix, 01566 cvm::real const &new_value, 01567 size_t const &imult = 0, 01568 bool add = false) 01569 { 01570 if (add) { 01571 if (samples) 01572 data[address(ix) + imult] += new_value * samples->new_count(ix); 01573 else 01574 data[address(ix) + imult] += new_value; 01575 } else { 01576 if (samples) 01577 data[address(ix) + imult] = new_value * samples->value(ix); 01578 else 01579 data[address(ix) + imult] = new_value; 01580 } 01581 has_data = true; 01582 } 01583 01584 01586 inline cvm::real average() 01587 { 01588 size_t n = 0; 01589 01590 if (nd != 1 || nx[0] == 0) { 01591 return 0.0; 01592 } 01593 01594 cvm::real sum = 0.0; 01595 std::vector<int> ix = new_index(); 01596 if (samples) { 01597 for ( ; index_ok(ix); incr(ix)) { 01598 if ( (n = samples->value(ix)) ) 01599 sum += value(ix) / n; 01600 } 01601 } else { 01602 for ( ; index_ok(ix); incr(ix)) { 01603 sum += value(ix); 01604 } 01605 } 01606 return (sum / cvm::real(nx[0])); 01607 } 01608 01611 void write_1D_integral(std::ostream &os); 01612 01613 }; 01614 01615 01616 01618 01619 class integrate_potential : public colvar_grid_scalar 01620 { 01621 public: 01622 01623 integrate_potential(); 01624 01625 virtual ~integrate_potential() 01626 {} 01627 01629 integrate_potential(std::vector<colvar *> &colvars, colvar_grid_gradient * gradients); 01630 01632 integrate_potential(colvar_grid_gradient * gradients); 01633 01635 int integrate(const int itmax, const cvm::real & tol, cvm::real & err); 01636 01639 void update_div_neighbors(const std::vector<int> &ix); 01640 01643 void set_div(); 01644 01647 inline void set_zero_minimum() { 01648 add_constant(-1.0 * minimum_value()); 01649 } 01650 01651 protected: 01652 01653 // Reference to gradient grid 01654 colvar_grid_gradient *gradients; 01655 01657 std::vector<cvm::real> divergence; 01658 01659 // std::vector<cvm::real> inv_lap_diag; // Inverse of the diagonal of the Laplacian; for conditioning 01660 01663 void update_div_local(const std::vector<int> &ix); 01664 01668 void get_grad(cvm::real * g, std::vector<int> &ix); 01669 01671 void nr_linbcg_sym(const std::vector<cvm::real> &b, std::vector<cvm::real> &x, 01672 const cvm::real &tol, const int itmax, int &iter, cvm::real &err); 01673 01675 cvm::real l2norm(const std::vector<cvm::real> &x); 01676 01678 void atimes(const std::vector<cvm::real> &x, std::vector<cvm::real> &r); 01679 01680 // /// Inversion of preconditioner matrix 01681 // void asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x); 01682 }; 01683 01684 #endif 01685