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