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