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

colvar_UIestimator.h

Go to the documentation of this file.
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

Generated on Fri Oct 11 02:43:15 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002