00001 // -*- Mode:c++; c-basic-offset: 4; -*- 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 COLVAR_UIESTIMATOR_H 00011 #define COLVAR_UIESTIMATOR_H 00012 00013 #include <cmath> 00014 #include <vector> 00015 #include <iostream> 00016 #include <fstream> 00017 #include <string> 00018 00019 #include <typeinfo> 00020 00021 // only for colvar module! 00022 // when integrated into other code, just remove this line and "...cvm::backup_file(...)" 00023 #include "colvarmodule.h" 00024 #include "colvarproxy.h" 00025 00026 namespace UIestimator { 00027 const int Y_SIZE = 21; // defines the range of extended CV with respect to a given CV 00028 // For example, CV=10, width=1, Y_SIZE=21, then eCV=[0-20], having a size of 21 00029 const int HALF_Y_SIZE = 10; 00030 const int EXTENDED_X_SIZE = HALF_Y_SIZE; 00031 const double EPSILON = 0.000001; // for comparison of float numbers 00032 00033 class n_matrix { // Stores the distribution matrix of n(x,y) 00034 00035 public: 00036 n_matrix() {} 00037 n_matrix(const std::vector<double> & lowerboundary_input, // lowerboundary of x 00038 const std::vector<double> & upperboundary_input, // upperboundary of 00039 const std::vector<double> & width_input, // width of x 00040 const int y_size_input) { // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered 00041 00042 int i; 00043 00044 this->lowerboundary = lowerboundary_input; 00045 this->upperboundary = upperboundary_input; 00046 this->width = width_input; 00047 this->dimension = lowerboundary_input.size(); 00048 this->y_size = y_size_input; // keep in mind the internal (spare) matrix is stored in diagonal form 00049 this->y_total_size = int(cvm::pow(double(y_size_input), double(dimension)) + EPSILON); 00050 00051 // the range of the matrix is [lowerboundary, upperboundary] 00052 x_total_size = 1; 00053 for (i = 0; i < dimension; i++) { 00054 x_size.push_back(int((upperboundary_input[i] - lowerboundary_input[i]) / width_input[i] + EPSILON)); 00055 x_total_size *= x_size[i]; 00056 } 00057 00058 // initialize the internal matrix 00059 matrix.reserve(x_total_size); 00060 for (i = 0; i < x_total_size; i++) { 00061 matrix.push_back(std::vector<int>(y_total_size, 0)); 00062 } 00063 00064 temp.resize(dimension); 00065 } 00066 00067 int get_value(const std::vector<double> & x, const std::vector<double> & y) { 00068 return matrix[convert_x(x)][convert_y(x, y)]; 00069 } 00070 00071 void set_value(const std::vector<double> & x, const std::vector<double> & y, const int value) { 00072 matrix[convert_x(x)][convert_y(x,y)] = value; 00073 } 00074 00075 void increase_value(const std::vector<double> & x, const std::vector<double> & y, const int value) { 00076 matrix[convert_x(x)][convert_y(x,y)] += value; 00077 } 00078 00079 private: 00080 std::vector<double> lowerboundary; 00081 std::vector<double> upperboundary; 00082 std::vector<double> width; 00083 int dimension; 00084 std::vector<int> x_size; // the size of x in each dimension 00085 int x_total_size; // the size of x of the internal matrix 00086 int y_size; // the size of y in each dimension 00087 int y_total_size; // the size of y of the internal matrix 00088 00089 std::vector<std::vector<int> > matrix; // the internal matrix 00090 00091 std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource 00092 00093 int convert_x(const std::vector<double> & x) { // convert real x value to its interal index 00094 00095 int i, j; 00096 00097 for (i = 0; i < dimension; i++) { 00098 temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON); 00099 } 00100 00101 int index = 0; 00102 for (i = 0; i < dimension; i++) { 00103 if (i + 1 < dimension) { 00104 int x_temp = 1; 00105 for (j = i + 1; j < dimension; j++) 00106 x_temp *= x_size[j]; 00107 index += temp[i] * x_temp; 00108 } 00109 else 00110 index += temp[i]; 00111 } 00112 return index; 00113 } 00114 00115 int convert_y(const std::vector<double> & x, const std::vector<double> & y) { // convert real y value to its interal index 00116 00117 int i; 00118 00119 for (i = 0; i < dimension; i++) { 00120 temp[i] = int(round((round(y[i] / width[i] + EPSILON) - round(x[i] / width[i] + EPSILON)) + (y_size - 1) / 2 + EPSILON)); 00121 } 00122 00123 int index = 0; 00124 for (i = 0; i < dimension; i++) { 00125 if (i + 1 < dimension) 00126 index += temp[i] * int(cvm::pow(double(y_size), double(dimension - i - 1)) + EPSILON); 00127 else 00128 index += temp[i]; 00129 } 00130 return index; 00131 } 00132 00133 double round(double r) { 00134 return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5); 00135 } 00136 }; 00137 00138 // vector, store the sum_x, sum_x_square, count_y 00139 template <typename T> 00140 class n_vector { 00141 00142 public: 00143 n_vector() {} 00144 n_vector(const std::vector<double> & lowerboundary_input, // lowerboundary of x 00145 const std::vector<double> & upperboundary_input, // upperboundary of 00146 const std::vector<double> & width_input, // width of x 00147 const int y_size_input, // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered 00148 const T & default_value) { // the default value of T 00149 00150 this->width = width_input; 00151 this->dimension = lowerboundary_input.size(); 00152 00153 x_total_size = 1; 00154 for (int i = 0; i < dimension; i++) { 00155 this->lowerboundary.push_back(lowerboundary_input[i] - (y_size_input - 1) / 2 * width_input[i] - EPSILON); 00156 this->upperboundary.push_back(upperboundary_input[i] + (y_size_input - 1) / 2 * width_input[i] + EPSILON); 00157 00158 x_size.push_back(int((this->upperboundary[i] - this->lowerboundary[i]) / this->width[i] + EPSILON)); 00159 x_total_size *= x_size[i]; 00160 } 00161 00162 // initialize the internal vector 00163 vector.resize(x_total_size, default_value); 00164 00165 temp.resize(dimension); 00166 } 00167 00168 T & get_value(const std::vector<double> & x) { 00169 return vector[convert_x(x)]; 00170 } 00171 00172 void set_value(const std::vector<double> & x, const T value) { 00173 vector[convert_x(x)] = value; 00174 } 00175 00176 void increase_value(const std::vector<double> & x, const T value) { 00177 vector[convert_x(x)] += value; 00178 } 00179 00180 private: 00181 std::vector<double> lowerboundary; 00182 std::vector<double> upperboundary; 00183 std::vector<double> width; 00184 int dimension; 00185 std::vector<int> x_size; // the size of x in each dimension 00186 int x_total_size; // the size of x of the internal matrix 00187 00188 std::vector<T> vector; // the internal vector 00189 00190 std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource 00191 00192 int convert_x(const std::vector<double> & x) { // convert real x value to its interal index 00193 00194 int i, j; 00195 00196 for (i = 0; i < dimension; i++) { 00197 temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON); 00198 } 00199 00200 int index = 0; 00201 for (i = 0; i < dimension; i++) { 00202 if (i + 1 < dimension) { 00203 int x_temp = 1; 00204 for (j = i + 1; j < dimension; j++) 00205 x_temp *= x_size[j]; 00206 index += temp[i] * x_temp; 00207 } 00208 else 00209 index += temp[i]; 00210 } 00211 return index; 00212 } 00213 }; 00214 00215 class UIestimator { // the implemension of UI estimator 00216 00217 public: 00218 UIestimator() {} 00219 00220 //called when (re)start an eabf simulation 00221 UIestimator(const std::vector<double> & lowerboundary_input, 00222 const std::vector<double> & upperboundary_input, 00223 const std::vector<double> & width_input, 00224 const std::vector<double> & krestr_input, // force constant in eABF 00225 const std::string & output_filename_input, // the prefix of output files 00226 const int output_freq_input, 00227 const bool restart_input, // whether restart from a .count and a .grad file 00228 const std::vector<std::string> & input_filename_input, // the prefixes of input files 00229 const double temperature_input) { 00230 00231 // initialize variables 00232 this->lowerboundary = lowerboundary_input; 00233 this->upperboundary = upperboundary_input; 00234 this->width = width_input; 00235 this->krestr = krestr_input; 00236 this->output_filename = output_filename_input; 00237 this->output_freq = output_freq_input; 00238 this->restart = restart_input; 00239 this->input_filename = input_filename_input; 00240 this->temperature = temperature_input; 00241 00242 int i, j; 00243 00244 dimension = lowerboundary.size(); 00245 00246 for (i = 0; i < dimension; i++) { 00247 sum_x.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0)); 00248 sum_x_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0)); 00249 00250 x_av.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0)); 00251 sigma_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0)); 00252 } 00253 00254 count_y = n_vector<int>(lowerboundary, upperboundary, width, Y_SIZE, 0); 00255 distribution_x_y = n_matrix(lowerboundary, upperboundary, width, Y_SIZE); 00256 00257 grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0)); 00258 count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0); 00259 00260 written = false; 00261 written_1D = false; 00262 00263 if (dimension == 1) { 00264 std::vector<double> upperboundary_temp = upperboundary; 00265 upperboundary_temp[0] = upperboundary[0] + width[0]; 00266 oneD_pmf = n_vector<double>(lowerboundary, upperboundary_temp, width, 1, 0.0); 00267 } 00268 00269 if (restart == true) { 00270 input_grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0)); 00271 input_count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0); 00272 00273 // initialize input_Grad and input_count 00274 // the loop_flag is a n-dimensional vector, increae from lowerboundary to upperboundary when looping 00275 std::vector<double> loop_flag(dimension, 0); 00276 for (i = 0; i < dimension; i++) { 00277 loop_flag[i] = lowerboundary[i]; 00278 } 00279 00280 i = 0; 00281 while (i >= 0) { 00282 for (j = 0; j < dimension; j++) { 00283 input_grad.set_value(loop_flag, std::vector<double>(dimension,0)); 00284 } 00285 input_count.set_value(loop_flag, 0); 00286 00287 // iterate over any dimensions 00288 i = dimension - 1; 00289 while (i >= 0) { 00290 loop_flag[i] += width[i]; 00291 if (loop_flag[i] > upperboundary[i] - width[i] + EPSILON) { 00292 loop_flag[i] = lowerboundary[i]; 00293 i--; 00294 } 00295 else 00296 break; 00297 } 00298 } 00299 read_inputfiles(input_filename); 00300 } 00301 } 00302 00303 ~UIestimator() {} 00304 00305 // called from MD engine every step 00306 bool update(cvm::step_number /* step */, 00307 std::vector<double> x, std::vector<double> y) { 00308 00309 int i; 00310 00311 for (i = 0; i < dimension; i++) { 00312 // for dihedral RC, it is possible that x = 179 and y = -179, should correct it 00313 // may have problem, need to fix 00314 if (x[i] > 150 && y[i] < -150) { 00315 y[i] += 360; 00316 } 00317 if (x[i] < -150 && y[i] > 150) { 00318 y[i] -= 360; 00319 } 00320 00321 if (x[i] < lowerboundary[i] - EXTENDED_X_SIZE * width[i] + EPSILON || x[i] > upperboundary[i] + EXTENDED_X_SIZE * width[i] - EPSILON \ 00322 || y[i] - x[i] < -HALF_Y_SIZE * width[i] + EPSILON || y[i] - x[i] > HALF_Y_SIZE * width[i] - EPSILON \ 00323 || y[i] - lowerboundary[i] < -HALF_Y_SIZE * width[i] + EPSILON || y[i] - upperboundary[i] > HALF_Y_SIZE * width[i] - EPSILON) 00324 return false; 00325 } 00326 00327 for (i = 0; i < dimension; i++) { 00328 sum_x[i].increase_value(y, x[i]); 00329 sum_x_square[i].increase_value(y, x[i] * x[i]); 00330 } 00331 count_y.increase_value(y, 1); 00332 00333 for (i = 0; i < dimension; i++) { 00334 // adapt colvars precision 00335 if (x[i] < lowerboundary[i] + EPSILON || x[i] > upperboundary[i] - EPSILON) 00336 return false; 00337 } 00338 distribution_x_y.increase_value(x, y, 1); 00339 00340 return true; 00341 } 00342 00343 // update the output_filename 00344 void update_output_filename(const std::string& filename) { 00345 output_filename = filename; 00346 } 00347 00348 private: 00349 std::vector<n_vector<double> > sum_x; // the sum of x in each y bin 00350 std::vector<n_vector<double> > sum_x_square; // the sum of x in each y bin 00351 n_vector<int> count_y; // the distribution of y 00352 n_matrix distribution_x_y; // the distribution of <x, y> pair 00353 00354 int dimension; 00355 00356 std::vector<double> lowerboundary; 00357 std::vector<double> upperboundary; 00358 std::vector<double> width; 00359 std::vector<double> krestr; 00360 std::string output_filename; 00361 int output_freq; 00362 bool restart; 00363 std::vector<std::string> input_filename; 00364 double temperature; 00365 00366 n_vector<std::vector<double> > grad; 00367 n_vector<int> count; 00368 00369 n_vector<double> oneD_pmf; 00370 00371 n_vector<std::vector<double> > input_grad; 00372 n_vector<int> input_count; 00373 00374 // used in double integration 00375 std::vector<n_vector<double> > x_av; 00376 std::vector<n_vector<double> > sigma_square; 00377 00378 bool written; 00379 bool written_1D; 00380 00381 public: 00382 // calculate gradients from the internal variables 00383 void calc_pmf() { 00384 colvarproxy *proxy = cvm::main()->proxy; 00385 00386 int norm; 00387 int i, j, k; 00388 00389 std::vector<double> loop_flag(dimension, 0); 00390 for (i = 0; i < dimension; i++) { 00391 loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i]; 00392 } 00393 00394 i = 0; 00395 while (i >= 0) { 00396 norm = count_y.get_value(loop_flag) > 0 ? count_y.get_value(loop_flag) : 1; 00397 for (j = 0; j < dimension; j++) { 00398 x_av[j].set_value(loop_flag, sum_x[j].get_value(loop_flag) / norm); 00399 sigma_square[j].set_value(loop_flag, sum_x_square[j].get_value(loop_flag) / norm - x_av[j].get_value(loop_flag) * x_av[j].get_value(loop_flag)); 00400 } 00401 00402 // iterate over any dimensions 00403 i = dimension - 1; 00404 while (i >= 0) { 00405 loop_flag[i] += width[i]; 00406 if (loop_flag[i] > upperboundary[i] + HALF_Y_SIZE * width[i] - width[i] + EPSILON) { 00407 loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i]; 00408 i--; 00409 } 00410 else 00411 break; 00412 } 00413 } 00414 00415 // double integration 00416 std::vector<double> av(dimension, 0); 00417 std::vector<double> diff_av(dimension, 0); 00418 00419 std::vector<double> loop_flag_x(dimension, 0); 00420 std::vector<double> loop_flag_y(dimension, 0); 00421 for (i = 0; i < dimension; i++) { 00422 loop_flag_x[i] = lowerboundary[i]; 00423 loop_flag_y[i] = loop_flag_x[i] - HALF_Y_SIZE * width[i]; 00424 } 00425 00426 i = 0; 00427 while (i >= 0) { 00428 norm = 0; 00429 for (k = 0; k < dimension; k++) { 00430 av[k] = 0; 00431 diff_av[k] = 0; 00432 loop_flag_y[k] = loop_flag_x[k] - HALF_Y_SIZE * width[k]; 00433 } 00434 00435 j = 0; 00436 while (j >= 0) { 00437 norm += distribution_x_y.get_value(loop_flag_x, loop_flag_y); 00438 for (k = 0; k < dimension; k++) { 00439 if (sigma_square[k].get_value(loop_flag_y) > EPSILON || sigma_square[k].get_value(loop_flag_y) < -EPSILON) 00440 av[k] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * ( (loop_flag_x[k] + 0.5 * width[k]) - x_av[k].get_value(loop_flag_y)) / sigma_square[k].get_value(loop_flag_y); 00441 00442 diff_av[k] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * (loop_flag_x[k] - loop_flag_y[k]); 00443 } 00444 00445 // iterate over any dimensions 00446 j = dimension - 1; 00447 while (j >= 0) { 00448 loop_flag_y[j] += width[j]; 00449 if (loop_flag_y[j] > loop_flag_x[j] + HALF_Y_SIZE * width[j] - width[j] + EPSILON) { 00450 loop_flag_y[j] = loop_flag_x[j] - HALF_Y_SIZE * width[j]; 00451 j--; 00452 } 00453 else 00454 break; 00455 } 00456 } 00457 00458 std::vector<double> grad_temp(dimension, 0); 00459 for (k = 0; k < dimension; k++) { 00460 diff_av[k] /= (norm > 0 ? norm : 1); 00461 av[k] = proxy->boltzmann() * temperature * av[k] / (norm > 0 ? norm : 1); 00462 grad_temp[k] = av[k] - krestr[k] * diff_av[k]; 00463 } 00464 grad.set_value(loop_flag_x, grad_temp); 00465 count.set_value(loop_flag_x, norm); 00466 00467 // iterate over any dimensions 00468 i = dimension - 1; 00469 while (i >= 0) { 00470 loop_flag_x[i] += width[i]; 00471 if (loop_flag_x[i] > upperboundary[i] - width[i] + EPSILON) { 00472 loop_flag_x[i] = lowerboundary[i]; 00473 i--; 00474 } 00475 else 00476 break; 00477 } 00478 } 00479 } 00480 00481 00482 // calculate 1D pmf 00483 void calc_1D_pmf() 00484 { 00485 std::vector<double> last_position(1, 0); 00486 std::vector<double> position(1, 0); 00487 00488 double min = 0; 00489 double dG = 0; 00490 double i; 00491 00492 oneD_pmf.set_value(lowerboundary, 0); 00493 last_position = lowerboundary; 00494 for (i = lowerboundary[0] + width[0]; i < upperboundary[0] + EPSILON; i += width[0]) { 00495 position[0] = i + EPSILON; 00496 if (restart == false || input_count.get_value(last_position) == 0) { 00497 dG = oneD_pmf.get_value(last_position) + grad.get_value(last_position)[0] * width[0]; 00498 } 00499 else { 00500 dG = oneD_pmf.get_value(last_position) + ((grad.get_value(last_position)[0] * count.get_value(last_position) + input_grad.get_value(last_position)[0] * input_count.get_value(last_position)) / (count.get_value(last_position) + input_count.get_value(last_position))) * width[0]; 00501 } 00502 if (dG < min) 00503 min = dG; 00504 oneD_pmf.set_value(position, dG); 00505 last_position[0] = i + EPSILON; 00506 } 00507 00508 for (i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) { 00509 position[0] = i + EPSILON; 00510 oneD_pmf.set_value(position, oneD_pmf.get_value(position) - min); 00511 } 00512 } 00513 00514 // write 1D pmf 00515 void write_1D_pmf() { 00516 std::string pmf_filename = output_filename + ".UI.pmf"; 00517 00518 // only for colvars module! 00519 if (written_1D) cvm::backup_file(pmf_filename.c_str()); 00520 00521 std::ostream &ofile_pmf = cvm::proxy->output_stream(pmf_filename, 00522 "PMF file"); 00523 00524 std::vector<double> position(1, 0); 00525 for (double i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) { 00526 ofile_pmf << i << " "; 00527 position[0] = i + EPSILON; 00528 ofile_pmf << oneD_pmf.get_value(position) << std::endl; 00529 } 00530 cvm::proxy->close_output_stream(pmf_filename); 00531 00532 written_1D = true; 00533 } 00534 00535 // write heads of the output files 00536 void writehead(std::ostream& os) const { 00537 os << "# " << dimension << std::endl; 00538 for (int i = 0; i < dimension; i++) { 00539 os << "# " << lowerboundary[i] << " " << width[i] << " " << int((upperboundary[i] - lowerboundary[i]) / width[i] + EPSILON) << " " << 0 << std::endl; 00540 } 00541 os << std::endl; 00542 } 00543 00544 // write interal data, used for testing 00545 void write_interal_data() { 00546 std::string internal_filename = output_filename + ".UI.internal"; 00547 00548 std::ostream &ofile_internal = cvm::proxy->output_stream(internal_filename, 00549 "UI internal file"); 00550 00551 std::vector<double> loop_flag(dimension, 0); 00552 for (int i = 0; i < dimension; i++) { 00553 loop_flag[i] = lowerboundary[i]; 00554 } 00555 00556 int n = 0; 00557 while (n >= 0) { 00558 for (int j = 0; j < dimension; j++) { 00559 ofile_internal << loop_flag[j] + 0.5 * width[j] << " "; 00560 } 00561 00562 for (int k = 0; k < dimension; k++) { 00563 ofile_internal << grad.get_value(loop_flag)[k] << " "; 00564 } 00565 00566 std::vector<double> ii(dimension,0); 00567 for (double i = loop_flag[0] - 10; i < loop_flag[0] + 10 + EPSILON; i+= width[0]) { 00568 for (double j = loop_flag[1] - 10; j< loop_flag[1] + 10 + EPSILON; j+=width[1]) { 00569 ii[0] = i; 00570 ii[1] = j; 00571 ofile_internal << i <<" "<<j<<" "<< distribution_x_y.get_value(loop_flag,ii)<< " "; 00572 } 00573 } 00574 ofile_internal << std::endl; 00575 00576 // iterate over any dimensions 00577 n = dimension - 1; 00578 while (n >= 0) { 00579 loop_flag[n] += width[n]; 00580 if (loop_flag[n] > upperboundary[n] - width[n] + EPSILON) { 00581 loop_flag[n] = lowerboundary[n]; 00582 n--; 00583 } 00584 else 00585 break; 00586 } 00587 } 00588 cvm::proxy->close_output_stream(internal_filename.c_str()); 00589 } 00590 00591 // write output files 00592 void write_files() { 00593 std::string grad_filename = output_filename + ".UI.grad"; 00594 std::string hist_filename = output_filename + ".UI.hist.grad"; 00595 std::string count_filename = output_filename + ".UI.count"; 00596 00597 int i, j; 00598 // 00599 // only for colvars module! 00600 if (written) cvm::backup_file(grad_filename.c_str()); 00601 //if (written) cvm::backup_file(hist_filename.c_str()); 00602 if (written) cvm::backup_file(count_filename.c_str()); 00603 00604 std::ostream &ofile = cvm::proxy->output_stream(grad_filename, 00605 "gradient file"); 00606 std::ostream &ofile_hist = cvm::proxy->output_stream(hist_filename, 00607 "gradient history file"); 00608 std::ostream &ofile_count = cvm::proxy->output_stream(count_filename, 00609 "count file"); 00610 00611 writehead(ofile); 00612 writehead(ofile_hist); 00613 writehead(ofile_count); 00614 00615 if (dimension == 1) { 00616 calc_1D_pmf(); 00617 write_1D_pmf(); 00618 } 00619 00620 std::vector<double> loop_flag(dimension, 0); 00621 for (i = 0; i < dimension; i++) { 00622 loop_flag[i] = lowerboundary[i]; 00623 } 00624 00625 i = 0; 00626 while (i >= 0) { 00627 for (j = 0; j < dimension; j++) { 00628 ofile << loop_flag[j] + 0.5 * width[j] << " "; 00629 ofile_hist << loop_flag[j] + 0.5 * width[j] << " "; 00630 ofile_count << loop_flag[j] + 0.5 * width[j] << " "; 00631 } 00632 00633 if (restart == false) { 00634 for (j = 0; j < dimension; j++) { 00635 ofile << grad.get_value(loop_flag)[j] << " "; 00636 ofile_hist << grad.get_value(loop_flag)[j] << " "; 00637 } 00638 ofile << std::endl; 00639 ofile_hist << std::endl; 00640 ofile_count << count.get_value(loop_flag) << " " <<std::endl; 00641 } 00642 else { 00643 double final_grad = 0; 00644 for (j = 0; j < dimension; j++) { 00645 int total_count_temp = (count.get_value(loop_flag) + input_count.get_value(loop_flag)); 00646 if (input_count.get_value(loop_flag) == 0) 00647 final_grad = grad.get_value(loop_flag)[j]; 00648 else 00649 final_grad = ((grad.get_value(loop_flag)[j] * count.get_value(loop_flag) + input_grad.get_value(loop_flag)[j] * input_count.get_value(loop_flag)) / total_count_temp); 00650 ofile << final_grad << " "; 00651 ofile_hist << final_grad << " "; 00652 } 00653 ofile << std::endl; 00654 ofile_hist << std::endl; 00655 ofile_count << (count.get_value(loop_flag) + input_count.get_value(loop_flag)) << " " <<std::endl; 00656 } 00657 00658 // iterate over any dimensions 00659 i = dimension - 1; 00660 while (i >= 0) { 00661 loop_flag[i] += width[i]; 00662 if (loop_flag[i] > upperboundary[i] - width[i] + EPSILON) { 00663 loop_flag[i] = lowerboundary[i]; 00664 i--; 00665 ofile << std::endl; 00666 ofile_hist << std::endl; 00667 ofile_count << std::endl; 00668 } 00669 else 00670 break; 00671 } 00672 } 00673 cvm::proxy->close_output_stream(grad_filename.c_str()); 00674 // cvm::proxy->close_output_stream(hist_filename.c_str()); 00675 cvm::proxy->close_output_stream(count_filename.c_str()); 00676 00677 written = true; 00678 } 00679 00680 // read input files 00681 void read_inputfiles(const std::vector<std::string> filename) 00682 { 00683 char sharp; 00684 double nothing; 00685 int dimension_temp; 00686 int i, j, k, l, m; 00687 00688 colvarproxy *proxy = cvm::main()->proxy; 00689 std::vector<double> loop_bin_size(dimension, 0); 00690 std::vector<double> position_temp(dimension, 0); 00691 std::vector<double> grad_temp(dimension, 0); 00692 int count_temp = 0; 00693 for (i = 0; i < int(filename.size()); i++) { 00694 int size = 1 , size_temp = 0; 00695 00696 std::string count_filename = filename[i] + ".UI.count"; 00697 std::string grad_filename = filename[i] + ".UI.grad"; 00698 00699 std::istream &count_file = 00700 proxy->input_stream(count_filename, "count filename"); 00701 std::istream &grad_file = 00702 proxy->input_stream(grad_filename, "gradient filename"); 00703 00704 if (!count_file || !grad_file) { 00705 return; 00706 } 00707 00708 count_file >> sharp >> dimension_temp; 00709 grad_file >> sharp >> dimension_temp; 00710 00711 for (j = 0; j < dimension; j++) { 00712 count_file >> sharp >> nothing >> nothing >> size_temp >> nothing; 00713 grad_file >> sharp >> nothing >> nothing >> nothing >> nothing; 00714 size *= size_temp; 00715 } 00716 00717 for (j = 0; j < size; j++) { 00718 do { 00719 for (k = 0; k < dimension; k++) { 00720 count_file >> position_temp[k]; 00721 grad_file >> nothing; 00722 } 00723 00724 for (l = 0; l < dimension; l++) { 00725 grad_file >> grad_temp[l]; 00726 } 00727 count_file >> count_temp; 00728 } 00729 while (position_temp[i] < lowerboundary[i] - EPSILON || position_temp[i] > upperboundary[i] + EPSILON); 00730 00731 if (count_temp == 0) { 00732 continue; 00733 } 00734 00735 for (m = 0; m < dimension; m++) { 00736 grad_temp[m] = (grad_temp[m] * count_temp + input_grad.get_value(position_temp)[m] * input_count.get_value(position_temp)) / (count_temp + input_count.get_value(position_temp)); 00737 } 00738 input_grad.set_value(position_temp, grad_temp); 00739 input_count.increase_value(position_temp, count_temp); 00740 } 00741 00742 proxy->close_input_stream(count_filename); 00743 proxy->close_input_stream(grad_filename); 00744 } 00745 } 00746 }; 00747 } 00748 00749 #endif