Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2017 of Haochuan Chen (excluding colvar_UIestimator.h)
3 : Copyright (c) 2017 of Haohao Fu (colvar_UIestimator.h)
4 :
5 : This program is free software: you can redistribute it and/or modify
6 : it under the terms of the GNU Lesser General Public License as published
7 : by the Free Software Foundation, either version 3 of the License, or
8 : (at your option) any later version.
9 :
10 : This program is distributed in the hope that it will be useful,
11 : but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : GNU Lesser General Public License for more details.
14 :
15 : You should have received a copy of the GNU Lesser General Public License
16 : along with this program. If not, see <http://www.gnu.org/licenses/>.
17 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
18 : #ifndef __PLUMED_drr_colvar_UIestimator_h
19 : #define __PLUMED_drr_colvar_UIestimator_h
20 : // The original code(https://github.com/fhh2626/colvars/blob/master/src/colvar_UIestimator.h) has been modified by Haochuan Chen.
21 : // Modifications:
22 : // 1. Disable colvars related code.
23 : // 2. Change boltzmann constant.
24 : // 3. Change output precision.
25 : // I(Haochuan Chen) don't know how to maintain this code and how it runs. If you are interested in it, please contact Haohao Fu.
26 :
27 : #include <cmath>
28 : #include <vector>
29 : #include <fstream>
30 : #include <string>
31 : #include <iomanip>
32 : #include <limits>
33 :
34 : #include <typeinfo>
35 :
36 : // only for colvar module!
37 : // when integrated into other code, just remove this line and "...cvm::backup_file(...)"
38 : // #include "colvarmodule.h"
39 :
40 : namespace PLMD {
41 : namespace drr {
42 :
43 : namespace UIestimator {
44 : const int Y_SIZE = 21;
45 : const int HALF_Y_SIZE = 10;
46 : const double BOLTZMANN = 0.0083144621;
47 : const int EXTENDED_X_SIZE = HALF_Y_SIZE;
48 :
49 : class n_matrix { // spare matrix, stores the distribution matrix of n(x,y)
50 : public:
51 21 : n_matrix() {}
52 9 : n_matrix(const std::vector<double> & lowerboundary_p, // lowerboundary of x
53 : const std::vector<double> & upperboundary_p, // upperboundary of
54 : const std::vector<double> & width_p, // width of x
55 : const int y_size) // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
56 9 : : lowerboundary(lowerboundary_p), upperboundary(upperboundary_p), width(width_p) {
57 9 : this->dimension = lowerboundary.size();
58 9 : this->y_size = y_size; // keep in mind the internal (spare) matrix is stored in diagonal form
59 9 : this->y_total_size = int(pow(y_size, dimension) + 0.000001);
60 :
61 : // the range of the matrix is [lowerboundary, upperboundary]
62 9 : x_total_size = 1;
63 24 : for (int i = 0; i < dimension; i++) {
64 15 : x_size.push_back(int((upperboundary[i] - lowerboundary[i]) / width[i] + 0.000001));
65 15 : x_total_size *= x_size[i];
66 : }
67 :
68 : // initialize the internal matrix
69 9 : matrix.reserve(x_total_size);
70 100 : for (int i = 0; i < x_total_size; i++) {
71 182 : matrix.push_back(std::vector<int>(y_total_size, 0));
72 : }
73 :
74 9 : temp.resize(dimension);
75 9 : }
76 :
77 171000 : int inline get_value(const std::vector<double> & x, const std::vector<double> & y) {
78 : //if (matrix[convert_x(x)][convert_y(x, y)]!=0)
79 : //{
80 : //std::cout<<convert_x(x)<<" "<<convert_y(x, y)<<" "<<x[0]<<" "<<x[1]<<" "<<y[0]<<" "<<y[1]<<" ";
81 : //std::cout<<matrix[convert_x(x)][convert_y(x, y)]<<"sadasfdasaaaaaaaa"<<std::endl;
82 : //}
83 171000 : return matrix[convert_x(x)][convert_y(x, y)];
84 : }
85 :
86 : void inline set_value(const std::vector<double> & x, const std::vector<double> & y, const int value) {
87 : matrix[convert_x(x)][convert_y(x,y)] = value;
88 : }
89 :
90 87 : void inline increase_value(const std::vector<double> & x, const std::vector<double> & y, const int value) {
91 87 : matrix[convert_x(x)][convert_y(x,y)] += value;
92 87 : }
93 :
94 : private:
95 : std::vector<double> lowerboundary;
96 : std::vector<double> upperboundary;
97 : std::vector<double> width;
98 : int dimension;
99 : std::vector<int> x_size; // the size of x in each dimension
100 : int x_total_size; // the size of x of the internal matrix
101 : int y_size; // the size of y in each dimension
102 : int y_total_size; // the size of y of the internal matrix
103 :
104 : std::vector<std::vector<int> > matrix; // the internal matrix
105 :
106 : std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource
107 :
108 171087 : int convert_x(const std::vector<double> & x) { // convert real x value to its interal index
109 510684 : for (int i = 0; i < dimension; i++) {
110 339597 : temp[i] = int((x[i] - lowerboundary[i]) / width[i] + 0.000001);
111 : }
112 :
113 : int index = 0;
114 510684 : for (int i = 0; i < dimension; i++) {
115 339597 : if (i + 1 < dimension) {
116 : int x_temp = 1;
117 337020 : for (int j = i + 1; j < dimension; j++) {
118 168510 : x_temp *= x_size[j];
119 : }
120 168510 : index += temp[i] * x_temp;
121 : } else {
122 171087 : index += temp[i];
123 : }
124 : }
125 171087 : return index;
126 : }
127 :
128 171087 : int convert_y(const std::vector<double> & x, const std::vector<double> & y) { // convert real y value to its interal index
129 510684 : for (int i = 0; i < dimension; i++) {
130 339597 : temp[i] = round((round(y[i] / width[i] + 0.000001) - round(x[i] / width[i] + 0.000001)) + (y_size - 1) / 2 + 0.000001);
131 : }
132 :
133 : int index = 0;
134 510684 : for (int i = 0; i < dimension; i++) {
135 339597 : if (i + 1 < dimension) {
136 168510 : index += temp[i] * int(pow(y_size, dimension - i - 1) + 0.000001);
137 : } else {
138 171087 : index += temp[i];
139 : }
140 : }
141 171087 : return index;
142 : }
143 :
144 1018791 : double round(double r) {
145 1018791 : return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
146 : }
147 : };
148 :
149 : // vector, store the sum_x, sum_x_square, count_y
150 : template <typename T>
151 : class n_vector {
152 : public:
153 126 : n_vector() {}
154 92 : n_vector(const std::vector<double> & lowerboundary, // lowerboundary of x
155 : const std::vector<double> & upperboundary, // upperboundary of
156 : const std::vector<double> & width_p, // width of x
157 : const int y_size, // size of y, for example, ysize=7, then when x=1, the distribution of y in [-2,4] is considered
158 : const T & default_value) // the default value of T
159 92 : :width(width_p) {
160 92 : this->dimension = lowerboundary.size();
161 :
162 92 : x_total_size = 1;
163 252 : for (int i = 0; i < dimension; i++) {
164 160 : this->lowerboundary.push_back(lowerboundary[i] - (y_size - 1) / 2 * width[i] - 0.000001);
165 160 : this->upperboundary.push_back(upperboundary[i] + (y_size - 1) / 2 * width[i] + 0.000001);
166 :
167 160 : x_size.push_back(int((this->upperboundary[i] - this->lowerboundary[i]) / this->width[i] + 0.000001));
168 160 : x_total_size *= x_size[i];
169 : }
170 :
171 : // initialize the internal vector
172 92 : vector.resize(x_total_size, default_value);
173 :
174 92 : temp.resize(dimension);
175 92 : }
176 :
177 288432 : T inline get_value(const std::vector<double> & x) {
178 288432 : return vector[convert_x(x)];
179 : }
180 :
181 27227 : void inline set_value(const std::vector<double> & x, const T & value) {
182 27227 : vector[convert_x(x)] = value;
183 27227 : }
184 :
185 409 : void inline increase_value(const std::vector<double> & x, const T & value) {
186 409 : vector[convert_x(x)] += value;
187 409 : }
188 : private:
189 : std::vector<double> lowerboundary;
190 : std::vector<double> upperboundary;
191 : std::vector<double> width;
192 : int dimension;
193 : std::vector<int> x_size; // the size of x in each dimension
194 : int x_total_size; // the size of x of the internal matrix
195 :
196 : std::vector<T> vector; // the internal vector
197 :
198 : std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource
199 :
200 316068 : int convert_x(const std::vector<double> & x) { // convert real x value to its interal index
201 943161 : for (int i = 0; i < dimension; i++) {
202 627093 : temp[i] = int((x[i] - lowerboundary[i]) / width[i] + 0.000001);
203 : }
204 :
205 : int index = 0;
206 943161 : for (int i = 0; i < dimension; i++) {
207 627093 : if (i + 1 < dimension) {
208 : int x_temp = 1;
209 622050 : for (int j = i + 1; j < dimension; j++) {
210 311025 : x_temp *= x_size[j];
211 : }
212 311025 : index += temp[i] * x_temp;
213 : } else {
214 316068 : index += temp[i];
215 : }
216 : }
217 316068 : return index;
218 : }
219 : };
220 :
221 : class UIestimator { // the implemension of UI estimator
222 : public:
223 12 : UIestimator() {}
224 :
225 : //called when (re)start an eabf simulation
226 9 : UIestimator(const std::vector<double>& lowerboundary_p,
227 : const std::vector<double>& upperboundary_p,
228 : const std::vector<double>& width_p,
229 : const std::vector<double>& krestr_p, // force constant in eABF
230 : const std::string& output_filename_p, // the prefix of output files
231 : const int output_freq_p,
232 : const bool restart_p, // whether restart from a .count and a .grad file
233 : const std::vector<std::string>& input_filename_p, // the prefixes of input files
234 : const double temperature_p)
235 9 : : lowerboundary(lowerboundary_p), upperboundary(upperboundary_p),
236 9 : width(width_p), krestr(krestr_p),
237 9 : output_filename(output_filename_p), output_freq(output_freq_p),
238 9 : restart(restart_p), input_filename(input_filename_p),
239 18 : temperature(temperature_p) {
240 :
241 9 : dimension = lowerboundary.size();
242 :
243 24 : for (int i = 0; i < dimension; i++) {
244 30 : sum_x.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
245 30 : sum_x_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
246 :
247 30 : x_av.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
248 30 : sigma_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
249 : }
250 :
251 9 : count_y = n_vector<int>(lowerboundary, upperboundary, width, Y_SIZE, 0);
252 9 : distribution_x_y = n_matrix(lowerboundary, upperboundary, width, Y_SIZE);
253 :
254 9 : grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0));
255 9 : count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0);
256 :
257 9 : written = false;
258 9 : written_1D = false;
259 :
260 9 : if (dimension == 1) {
261 3 : std::vector<double> upperboundary_temp = upperboundary;
262 3 : upperboundary_temp[0] = upperboundary[0] + width[0];
263 3 : oneD_pmf = n_vector<double>(lowerboundary, upperboundary_temp, width, 1, 0.0);
264 : }
265 :
266 9 : if (restart == true) {
267 1 : input_grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0));
268 1 : input_count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0);
269 :
270 : // initialize input_Grad and input_count
271 1 : std::vector<double> loop_flag(dimension, 0);
272 3 : for (int i = 0; i < dimension; i++) {
273 2 : loop_flag[i] = lowerboundary[i];
274 : }
275 : while (true) {
276 27 : for (int i = 0; i < dimension; i++) {
277 36 : input_grad.set_value(loop_flag, std::vector<double>(dimension,0));
278 : }
279 9 : input_count.set_value(loop_flag, 0);
280 :
281 : // iterate over any dimensions
282 9 : int i = dimension - 1;
283 : while (true) {
284 12 : loop_flag[i] += width[i];
285 12 : if (loop_flag[i] > upperboundary[i] - width[i] + 0.00001) {
286 4 : loop_flag[i] = lowerboundary[i];
287 4 : i--;
288 4 : if (i < 0) {
289 1 : goto INITIAL_LOOPEND;
290 : }
291 : } else {
292 : break;
293 : }
294 : }
295 : }
296 : INITIAL_LOOPEND:
297 1 : read_inputfiles(input_filename);
298 : }
299 9 : }
300 :
301 42 : ~UIestimator() {}
302 :
303 : // called from MD engine every step
304 87 : bool update(const int step, const std::vector<double> & x, std::vector<double> y) {
305 :
306 : //std::cout<<"weeeee: :"<<std::endl;
307 : //for (int i = 0; i < dimension; i++)
308 : //{
309 : // std::cout<<x[i]<<" "<<y[i]<<" ";
310 : //}
311 : //std::cout<<std::endl;
312 :
313 87 : if (step % output_freq == 0) {
314 21 : calc_pmf();
315 21 : write_files();
316 : //write_interal_data();
317 : }
318 :
319 246 : for (int i = 0; i < dimension; i++) {
320 : // for dihedral RC, it is possible that x = 179 and y = -179, should correct it
321 : // may have problem, need to fix
322 159 : if (x[i] > 150 && y[i] < -150) {
323 : //std::vector<double> x_temp(x);
324 : //x_temp[i] -= 360;
325 : //update(7, x_temp, y);
326 0 : y[i] += 360;
327 : }
328 159 : if (x[i] < -150 && y[i] > 150) {
329 : //std::vector<double> x_temp(x);
330 : //x_temp[i] += 360;
331 : //update(7, x_temp, y);
332 0 : y[i] -= 360;
333 : }
334 :
335 159 : if (x[i] < lowerboundary[i] - EXTENDED_X_SIZE * width[i] + 0.00001 || x[i] > upperboundary[i] + EXTENDED_X_SIZE * width[i] - 0.00001 \
336 159 : || y[i] - x[i] < -HALF_Y_SIZE * width[i] + 0.00001 || y[i] - x[i] > HALF_Y_SIZE * width[i] - 0.00001 \
337 318 : || y[i] - lowerboundary[i] < -HALF_Y_SIZE * width[i] + 0.00001 || y[i] - upperboundary[i] > HALF_Y_SIZE * width[i] - 0.00001) {
338 : return false;
339 : }
340 : }
341 :
342 : //for (int i = 0; i < dimension; i++)
343 : //{
344 : // std::cout<<x[i]<<" "<<y[i]<<" ";
345 : //}
346 : //std::cout<<std::endl;
347 :
348 246 : for (int i = 0; i < dimension; i++) {
349 159 : sum_x[i].increase_value(y, x[i]);
350 159 : sum_x_square[i].increase_value(y, x[i] * x[i]);
351 : }
352 87 : count_y.increase_value(y, 1);
353 :
354 246 : for (int i = 0; i < dimension; i++) {
355 : //if (x[i] < lowerboundary[i] + 0.000001 || x[i] > upperboundary[i] - 0.000001)
356 : // adapt colvars precision
357 159 : if (x[i] < lowerboundary[i] + 0.00001 || x[i] > upperboundary[i] - 0.00001) {
358 : return false;
359 : }
360 : }
361 87 : distribution_x_y.increase_value(x, y, 1);
362 :
363 87 : return true;
364 : }
365 :
366 : // update the output_filename
367 : void update_output_filename(const std::string& filename) {
368 87 : output_filename = filename;
369 : }
370 :
371 : private:
372 : std::vector<n_vector<double> > sum_x; // the sum of x in each y bin
373 : std::vector<n_vector<double> > sum_x_square; // the sum of x in each y bin
374 : n_vector<int> count_y; // the distribution of y
375 : n_matrix distribution_x_y; // the distribution of <x, y> pair
376 :
377 : int dimension;
378 :
379 : std::vector<double> lowerboundary;
380 : std::vector<double> upperboundary;
381 : std::vector<double> width;
382 : std::vector<double> krestr;
383 : std::string output_filename;
384 : int output_freq;
385 : bool restart;
386 : std::vector<std::string> input_filename;
387 : double temperature;
388 :
389 : n_vector<std::vector<double> > grad;
390 : n_vector<int> count;
391 :
392 : n_vector<double> oneD_pmf;
393 :
394 : n_vector<std::vector<double> > input_grad;
395 : n_vector<int> input_count;
396 :
397 : // used in double integration
398 : std::vector<n_vector<double> > x_av;
399 : std::vector<n_vector<double> > sigma_square;
400 :
401 : bool written;
402 : bool written_1D;
403 :
404 : // calculate gradients from the internal variables
405 21 : void calc_pmf() {
406 : int norm;
407 :
408 21 : std::vector<double> loop_flag(dimension, 0);
409 54 : for (int i = 0; i < dimension; i++) {
410 33 : loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i];
411 : }
412 :
413 : while (true) {
414 6783 : norm = count_y.get_value(loop_flag) > 0 ? count_y.get_value(loop_flag) : 1;
415 20106 : for (int i = 0; i < dimension; i++) {
416 13323 : x_av[i].set_value(loop_flag, sum_x[i].get_value(loop_flag) / norm);
417 13323 : sigma_square[i].set_value(loop_flag, sum_x_square[i].get_value(loop_flag) / norm - x_av[i].get_value(loop_flag) * x_av[i].get_value(loop_flag));
418 : }
419 :
420 : // iterate over any dimensions
421 6783 : int i = dimension - 1;
422 : while (true) {
423 7063 : loop_flag[i] += width[i];
424 7063 : if (loop_flag[i] > upperboundary[i] + HALF_Y_SIZE * width[i] - width[i] + 0.00001) {
425 301 : loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i];
426 301 : i--;
427 301 : if (i < 0) {
428 21 : goto LOOPEND;
429 : }
430 : } else {
431 : break;
432 : }
433 : }
434 : }
435 : LOOPEND:
436 :
437 : // double integration
438 21 : std::vector<double> av(dimension, 0);
439 21 : std::vector<double> diff_av(dimension, 0);
440 :
441 21 : std::vector<double> loop_flag_x(dimension, 0);
442 21 : std::vector<double> loop_flag_y(dimension, 0);
443 54 : for (int i = 0; i < dimension; i++) {
444 33 : loop_flag_x[i] = lowerboundary[i];
445 33 : loop_flag_y[i] = loop_flag_x[i] - HALF_Y_SIZE * width[i];
446 : }
447 :
448 : while (true) {
449 203 : norm = 0;
450 546 : for (int i = 0; i < dimension; i++) {
451 343 : av[i] = 0;
452 343 : diff_av[i] = 0;
453 343 : loop_flag_y[i] = loop_flag_x[i] - HALF_Y_SIZE * width[i];
454 : }
455 :
456 : while (true) {
457 : //std::cout<<"pppppppppppppppppppppp "<<loop_flag_x[0]<<" "<<loop_flag_x[1]<<" "<<loop_flag_y[0]<<" "<<loop_flag_y[1]<<std::endl;
458 57260 : norm += distribution_x_y.get_value(loop_flag_x, loop_flag_y);
459 170520 : for (int i = 0; i < dimension; i++) {
460 113260 : if (sigma_square[i].get_value(loop_flag_y) > 0.00001 || sigma_square[i].get_value(loop_flag_y) < -0.00001) {
461 480 : av[i] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * ( (loop_flag_x[i] + 0.5 * width[i]) - x_av[i].get_value(loop_flag_y)) / sigma_square[i].get_value(loop_flag_y);
462 : }
463 :
464 113260 : diff_av[i] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * (loop_flag_x[i] - loop_flag_y[i]);
465 : }
466 :
467 : // iterate over any dimensions
468 57260 : int i = dimension - 1;
469 : while (true) {
470 60060 : loop_flag_y[i] += width[i];
471 60060 : if (loop_flag_y[i] > loop_flag_x[i] + HALF_Y_SIZE * width[i] - width[i] + 0.00001) {
472 3003 : loop_flag_y[i] = loop_flag_x[i] - HALF_Y_SIZE * width[i];
473 3003 : i--;
474 3003 : if (i < 0) {
475 203 : goto LOOPEND2;
476 : }
477 : } else {
478 : break;
479 : }
480 : }
481 : }
482 : LOOPEND2:
483 :
484 203 : std::vector<double> grad_temp(dimension, 0);
485 546 : for (int i = 0; i < dimension; i++) {
486 343 : diff_av[i] /= (norm > 0 ? norm : 1);
487 343 : av[i] = BOLTZMANN * temperature * av[i] / (norm > 0 ? norm : 1);
488 343 : grad_temp[i] = av[i] - krestr[i] * diff_av[i];
489 : }
490 203 : grad.set_value(loop_flag_x, grad_temp);
491 203 : count.set_value(loop_flag_x, norm);
492 :
493 : // iterate over any dimensions
494 203 : int i = dimension - 1;
495 : while (true) {
496 243 : loop_flag_x[i] += width[i];
497 243 : if (loop_flag_x[i] > upperboundary[i] - width[i] + 0.00001) {
498 61 : loop_flag_x[i] = lowerboundary[i];
499 61 : i--;
500 61 : if (i < 0) {
501 21 : goto LOOPEND3;
502 : }
503 : } else {
504 : break;
505 : }
506 : }
507 : }
508 : LOOPEND3:
509 : ;
510 21 : }
511 :
512 :
513 : // calculate 1D pmf
514 9 : void calc_1D_pmf() {
515 9 : std::vector<double> last_position(1, 0);
516 9 : std::vector<double> position(1, 0);
517 :
518 : double min = 0;
519 9 : double dG = 0;
520 :
521 9 : oneD_pmf.set_value(lowerboundary, 0);
522 9 : last_position = lowerboundary;
523 72 : for (double i = lowerboundary[0] + width[0]; i < upperboundary[0] + 0.000001; i += width[0]) {
524 63 : position[0] = i + 0.000001;
525 63 : if (restart == false || input_count.get_value(last_position) == 0) {
526 63 : dG = oneD_pmf.get_value(last_position) + grad.get_value(last_position)[0] * width[0];
527 : } else {
528 0 : 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];
529 : }
530 63 : if (dG < min) {
531 : min = dG;
532 : }
533 63 : oneD_pmf.set_value(position, dG);
534 63 : last_position[0] = i + 0.000001;
535 : }
536 :
537 81 : for (double i = lowerboundary[0]; i < upperboundary[0] + 0.000001; i += width[0]) {
538 72 : position[0] = i + 0.000001;
539 72 : oneD_pmf.set_value(position, oneD_pmf.get_value(position) - min);
540 : }
541 9 : }
542 :
543 : // write 1D pmf
544 9 : void write_1D_pmf() {
545 9 : std::string pmf_filename = output_filename + ".UI.pmf";
546 :
547 : // only for colvars module!
548 : // if (written_1D) cvm::backup_file(pmf_filename.c_str());
549 :
550 9 : std::ofstream ofile_pmf(pmf_filename.c_str());
551 :
552 9 : std::vector<double> position(1, 0);
553 81 : for (double i = lowerboundary[0]; i < upperboundary[0] + 0.000001; i += width[0]) {
554 72 : ofile_pmf << i << " ";
555 72 : position[0] = i + 0.000001;
556 72 : ofile_pmf << oneD_pmf.get_value(position) << std::endl;
557 : }
558 9 : ofile_pmf.close();
559 :
560 9 : written_1D = true;
561 18 : }
562 :
563 : // write heads of the output files
564 63 : void writehead(std::ofstream& os) const {
565 63 : os << "# " << dimension << std::endl;
566 162 : for (int i = 0; i < dimension; i++) {
567 99 : os.precision(std::numeric_limits<double>::max_digits10);
568 198 : os << "# " << std::fixed << lowerboundary[i] << " " << width[i] << " " << int((upperboundary[i] - lowerboundary[i]) / width[i] + 0.000001) << " " << 0 << std::endl;
569 : }
570 : os << std::endl;
571 63 : }
572 :
573 : // write interal data, used for testing
574 : // void write_interal_data()
575 : // {
576 : // std::string internal_filaname = output_filename + ".UI.internal";
577 : //
578 : // std::ofstream ofile_internal(internal_filaname.c_str());
579 : //
580 : // std::vector<double> loop_flag(dimension, 0);
581 : // for (int i = 0; i < dimension; i++)
582 : // {
583 : // loop_flag[i] = lowerboundary[i];
584 : // }
585 : // while (true)
586 : // {
587 : // for (int i = 0; i < dimension; i++)
588 : // {
589 : // ofile_internal << loop_flag[i] + 0.5 * width[i] << " ";
590 : // }
591 : //
592 : // for (int i = 0; i < dimension; i++)
593 : // {
594 : // ofile_internal << grad.get_value(loop_flag)[i] << " ";
595 : // }
596 : //
597 : // std::vector<double> ii(dimension,0);
598 : // for (double i = loop_flag[0] - 10; i < loop_flag[0] + 10 + 0.00001; i+= width[0])
599 : // {
600 : // for (double j = loop_flag[1] - 10; j< loop_flag[1] + 10 +0.00001; j+=width[1])
601 : // {
602 : // ii[0] = i;
603 : // ii[1] = j;
604 : // ofile_internal << i <<" "<<j<<" "<< distribution_x_y.get_value(loop_flag,ii)<< " ";
605 : // }
606 : // }
607 : // ofile_internal << std::endl;
608 : //
609 : // // iterate over any dimensions
610 : // int i = dimension - 1;
611 : // while (true)
612 : // {
613 : // loop_flag[i] += width[i];
614 : // if (loop_flag[i] > upperboundary[i] - width[i] + 0.00001)
615 : // {
616 : // loop_flag[i] = lowerboundary[i];
617 : // i--;
618 : // if (i < 0)
619 : // goto LOOPEND5;
620 : // }
621 : // else
622 : // break;
623 : // }
624 : // }
625 : // LOOPEND5:
626 : // ofile_internal.close();
627 : // }
628 :
629 : // write output files
630 21 : void write_files() {
631 21 : std::string grad_filename = output_filename + ".UI.grad";
632 21 : std::string hist_filename = output_filename + ".UI.hist.grad";
633 21 : std::string count_filename = output_filename + ".UI.count";
634 :
635 : // only for colvars module!
636 : // if (written) cvm::backup_file(grad_filename.c_str());
637 : //if (written) cvm::backup_file(hist_filename.c_str());
638 : // if (written) cvm::backup_file(count_filename.c_str());
639 :
640 21 : std::ofstream ofile(grad_filename.c_str());
641 21 : std::ofstream ofile_hist(hist_filename.c_str(), std::ios::app);
642 21 : std::ofstream ofile_count(count_filename.c_str());
643 :
644 21 : writehead(ofile);
645 21 : writehead(ofile_hist);
646 21 : writehead(ofile_count);
647 :
648 21 : if (dimension == 1) {
649 9 : calc_1D_pmf();
650 9 : write_1D_pmf();
651 : }
652 :
653 21 : std::vector<double> loop_flag(dimension, 0);
654 54 : for (int i = 0; i < dimension; i++) {
655 33 : loop_flag[i] = lowerboundary[i];
656 : }
657 : while (true) {
658 546 : for (int i = 0; i < dimension; i++) {
659 343 : ofile << std::fixed << std::setprecision(6) << loop_flag[i] + 0.5 * width[i] << " ";
660 343 : ofile_hist << std::fixed << std::setprecision(6) << loop_flag[i] + 0.5 * width[i] << " ";
661 343 : ofile_count << std::fixed << std::setprecision(6) << loop_flag[i] + 0.5 * width[i] << " ";
662 : }
663 :
664 203 : if (restart == false) {
665 492 : for (int i = 0; i < dimension; i++) {
666 614 : ofile << std::fixed << std::setprecision(6) << grad.get_value(loop_flag)[i] << " ";
667 614 : ofile_hist << std::fixed << std::setprecision(6) << grad.get_value(loop_flag)[i] << " ";
668 : }
669 : ofile << std::endl;
670 : ofile_hist << std::endl;
671 185 : ofile_count << count.get_value(loop_flag) << " " <<std::endl;
672 : } else {
673 : double final_grad = 0;
674 54 : for (int i = 0; i < dimension; i++) {
675 36 : int total_count_temp = (count.get_value(loop_flag) + input_count.get_value(loop_flag));
676 36 : if (input_count.get_value(loop_flag) == 0) {
677 20 : final_grad = grad.get_value(loop_flag)[i];
678 : } else {
679 16 : final_grad = ((grad.get_value(loop_flag)[i] * count.get_value(loop_flag) + input_grad.get_value(loop_flag)[i] * input_count.get_value(loop_flag)) / total_count_temp);
680 : }
681 36 : ofile << std::fixed << std::setprecision(6) << final_grad << " ";
682 36 : ofile_hist << std::fixed << std::setprecision(6) << final_grad << " ";
683 : }
684 : ofile << std::endl;
685 : ofile_hist << std::endl;
686 18 : ofile_count << (count.get_value(loop_flag) + input_count.get_value(loop_flag)) << " " <<std::endl;
687 : }
688 :
689 : // iterate over any dimensions
690 203 : int i = dimension - 1;
691 : while (true) {
692 243 : loop_flag[i] += width[i];
693 243 : if (loop_flag[i] > upperboundary[i] - width[i] + 0.00001) {
694 61 : loop_flag[i] = lowerboundary[i];
695 61 : i--;
696 : ofile << std::endl;
697 : ofile_hist << std::endl;
698 : ofile_count << std::endl;
699 61 : if (i < 0) {
700 21 : goto LOOPEND4;
701 : }
702 : } else {
703 : break;
704 : }
705 : }
706 : }
707 : LOOPEND4:
708 21 : ofile.close();
709 21 : ofile_count.close();
710 21 : ofile_hist.close();
711 :
712 21 : written = true;
713 42 : }
714 :
715 : // read input files
716 1 : void read_inputfiles(const std::vector<std::string>& input_filename) {
717 : char sharp;
718 : double nothing;
719 : int dimension_temp;
720 :
721 1 : std::vector<double> loop_bin_size(dimension, 0);
722 1 : std::vector<double> position_temp(dimension, 0);
723 1 : std::vector<double> grad_temp(dimension, 0);
724 1 : int count_temp = 0;
725 2 : for (unsigned int i = 0; i < input_filename.size(); i++) {
726 1 : int size = 1, size_temp = 0;
727 :
728 1 : std::string count_filename = input_filename[i] + ".UI.count";
729 1 : std::string grad_filename = input_filename[i] + ".UI.grad";
730 :
731 1 : std::ifstream count_file(count_filename.c_str(), std::ios::in);
732 1 : std::ifstream grad_file(grad_filename.c_str(), std::ios::in);
733 :
734 1 : count_file >> sharp >> dimension_temp;
735 1 : grad_file >> sharp >> dimension_temp;
736 :
737 3 : for (int j = 0; j < dimension; j++) {
738 4 : count_file >> sharp >> nothing >> nothing >> size_temp >> nothing;
739 2 : grad_file >> sharp >> nothing >> nothing >> nothing >> nothing;
740 2 : size *= size_temp;
741 : }
742 :
743 10 : for (int j = 0; j < size; j++) {
744 : do {
745 27 : for (int k = 0; k < dimension; k++) {
746 18 : count_file >> position_temp[k];
747 : grad_file >> nothing;
748 : }
749 :
750 27 : for (int l = 0; l < dimension; l++) {
751 18 : grad_file >> grad_temp[l];
752 : }
753 9 : count_file >> count_temp;
754 9 : } while (position_temp[i] < lowerboundary[i] - 0.000001 || position_temp[i] > upperboundary[i] + 0.000001);
755 :
756 9 : if (count_temp == 0) {
757 5 : continue;
758 : }
759 :
760 12 : for (int m = 0; m < dimension; m++) {
761 8 : 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));
762 : }
763 4 : input_grad.set_value(position_temp, grad_temp);
764 4 : input_count.increase_value(position_temp, count_temp);
765 : }
766 :
767 1 : count_file.close();
768 1 : grad_file.close();
769 1 : }
770 1 : }
771 : };
772 : }
773 :
774 : }
775 : }
776 :
777 : #endif
|
