Collective Variables Module - Developer Documentation
Loading...
Searching...
No Matches
colvar_UIestimator.h
1// -*- Mode:c++; c-basic-offset: 4; -*-
2
3// This file is part of the Collective Variables module (Colvars).
4// The original version of Colvars and its updates are located at:
5// https://github.com/Colvars/colvars
6// Please update all Colvars source files before making any changes.
7// If you wish to distribute your changes, please submit them to the
8// Colvars repository at GitHub.
9
10#ifndef COLVAR_UIESTIMATOR_H
11#define COLVAR_UIESTIMATOR_H
12
13#include <cmath>
14#include <vector>
15#include <iostream>
16#include <fstream>
17#include <string>
18
19#include <typeinfo>
20
21// only for colvar module!
22// when integrated into other code, just remove this line and "...cvm::backup_file(...)"
23#include "colvarmodule.h"
24#include "colvarproxy.h"
25
26namespace UIestimator {
27 const int Y_SIZE = 21; // defines the range of extended CV with respect to a given CV
28 // For example, CV=10, width=1, Y_SIZE=21, then eCV=[0-20], having a size of 21
29 const int HALF_Y_SIZE = 10;
30 const int EXTENDED_X_SIZE = HALF_Y_SIZE;
31 const double EPSILON = 0.000001; // for comparison of float numbers
32
33 class n_matrix { // Stores the distribution matrix of n(x,y)
34
35 public:
36 n_matrix() {}
37 n_matrix(const std::vector<double> & lowerboundary_input, // lowerboundary of x
38 const std::vector<double> & upperboundary_input, // upperboundary of
39 const std::vector<double> & width_input, // width of x
40 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
41
42 int i;
43
44 this->lowerboundary = lowerboundary_input;
45 this->upperboundary = upperboundary_input;
46 this->width = width_input;
47 this->dimension = lowerboundary_input.size();
48 this->y_size = y_size_input; // keep in mind the internal (spare) matrix is stored in diagonal form
49 this->y_total_size = int(cvm::pow(double(y_size_input), double(dimension)) + EPSILON);
50
51 // the range of the matrix is [lowerboundary, upperboundary]
52 x_total_size = 1;
53 for (i = 0; i < dimension; i++) {
54 x_size.push_back(int((upperboundary_input[i] - lowerboundary_input[i]) / width_input[i] + EPSILON));
55 x_total_size *= x_size[i];
56 }
57
58 // initialize the internal matrix
59 matrix.reserve(x_total_size);
60 for (i = 0; i < x_total_size; i++) {
61 matrix.push_back(std::vector<int>(y_total_size, 0));
62 }
63
64 temp.resize(dimension);
65 }
66
67 int get_value(const std::vector<double> & x, const std::vector<double> & y) {
68 return matrix[convert_x(x)][convert_y(x, y)];
69 }
70
71 void set_value(const std::vector<double> & x, const std::vector<double> & y, const int value) {
72 matrix[convert_x(x)][convert_y(x,y)] = value;
73 }
74
75 void increase_value(const std::vector<double> & x, const std::vector<double> & y, const int value) {
76 matrix[convert_x(x)][convert_y(x,y)] += value;
77 }
78
79 private:
80 std::vector<double> lowerboundary;
81 std::vector<double> upperboundary;
82 std::vector<double> width;
83 int dimension;
84 std::vector<int> x_size; // the size of x in each dimension
85 int x_total_size; // the size of x of the internal matrix
86 int y_size; // the size of y in each dimension
87 int y_total_size; // the size of y of the internal matrix
88
89 std::vector<std::vector<int> > matrix; // the internal matrix
90
91 std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource
92
93 int convert_x(const std::vector<double> & x) { // convert real x value to its interal index
94
95 int i, j;
96
97 for (i = 0; i < dimension; i++) {
98 temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON);
99 }
100
101 int index = 0;
102 for (i = 0; i < dimension; i++) {
103 if (i + 1 < dimension) {
104 int x_temp = 1;
105 for (j = i + 1; j < dimension; j++)
106 x_temp *= x_size[j];
107 index += temp[i] * x_temp;
108 }
109 else
110 index += temp[i];
111 }
112 return index;
113 }
114
115 int convert_y(const std::vector<double> & x, const std::vector<double> & y) { // convert real y value to its interal index
116
117 int i;
118
119 for (i = 0; i < dimension; i++) {
120 temp[i] = int(round((round(y[i] / width[i] + EPSILON) - round(x[i] / width[i] + EPSILON)) + (y_size - 1) / 2 + EPSILON));
121 }
122
123 int index = 0;
124 for (i = 0; i < dimension; i++) {
125 if (i + 1 < dimension)
126 index += temp[i] * int(cvm::pow(double(y_size), double(dimension - i - 1)) + EPSILON);
127 else
128 index += temp[i];
129 }
130 return index;
131 }
132
133 double round(double r) {
134 return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5);
135 }
136 };
137
138 // vector, store the sum_x, sum_x_square, count_y
139 template <typename T>
140 class n_vector {
141
142 public:
143 n_vector() {}
144 n_vector(const std::vector<double> & lowerboundary_input, // lowerboundary of x
145 const std::vector<double> & upperboundary_input, // upperboundary of
146 const std::vector<double> & width_input, // width of x
147 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
148 const T & default_value) { // the default value of T
149
150 this->width = width_input;
151 this->dimension = lowerboundary_input.size();
152
153 x_total_size = 1;
154 for (int i = 0; i < dimension; i++) {
155 this->lowerboundary.push_back(lowerboundary_input[i] - (y_size_input - 1) / 2 * width_input[i] - EPSILON);
156 this->upperboundary.push_back(upperboundary_input[i] + (y_size_input - 1) / 2 * width_input[i] + EPSILON);
157
158 x_size.push_back(int((this->upperboundary[i] - this->lowerboundary[i]) / this->width[i] + EPSILON));
159 x_total_size *= x_size[i];
160 }
161
162 // initialize the internal vector
163 vector.resize(x_total_size, default_value);
164
165 temp.resize(dimension);
166 }
167
168 T & get_value(const std::vector<double> & x) {
169 return vector[convert_x(x)];
170 }
171
172 void set_value(const std::vector<double> & x, const T value) {
173 vector[convert_x(x)] = value;
174 }
175
176 void increase_value(const std::vector<double> & x, const T value) {
177 vector[convert_x(x)] += value;
178 }
179
180 private:
181 std::vector<double> lowerboundary;
182 std::vector<double> upperboundary;
183 std::vector<double> width;
184 int dimension;
185 std::vector<int> x_size; // the size of x in each dimension
186 int x_total_size; // the size of x of the internal matrix
187
188 std::vector<T> vector; // the internal vector
189
190 std::vector<int> temp; // this vector is used in convert_x and convert_y to save computational resource
191
192 int convert_x(const std::vector<double> & x) { // convert real x value to its interal index
193
194 int i, j;
195
196 for (i = 0; i < dimension; i++) {
197 temp[i] = int((x[i] - lowerboundary[i]) / width[i] + EPSILON);
198 }
199
200 int index = 0;
201 for (i = 0; i < dimension; i++) {
202 if (i + 1 < dimension) {
203 int x_temp = 1;
204 for (j = i + 1; j < dimension; j++)
205 x_temp *= x_size[j];
206 index += temp[i] * x_temp;
207 }
208 else
209 index += temp[i];
210 }
211 return index;
212 }
213 };
214
215 class UIestimator { // the implemension of UI estimator
216
217 public:
218 UIestimator() {}
219
220 //called when (re)start an eabf simulation
221 UIestimator(const std::vector<double> & lowerboundary_input,
222 const std::vector<double> & upperboundary_input,
223 const std::vector<double> & width_input,
224 const std::vector<double> & krestr_input, // force constant in eABF
225 const std::string & output_filename_input, // the prefix of output files
226 const int output_freq_input,
227 const bool restart_input, // whether restart from a .count and a .grad file
228 const std::vector<std::string> & input_filename_input, // the prefixes of input files
229 const double temperature_input) {
230
231 // initialize variables
232 this->lowerboundary = lowerboundary_input;
233 this->upperboundary = upperboundary_input;
234 this->width = width_input;
235 this->krestr = krestr_input;
236 this->output_filename = output_filename_input;
237 this->output_freq = output_freq_input;
238 this->restart = restart_input;
239 this->input_filename = input_filename_input;
240 this->temperature = temperature_input;
241
242 int i, j;
243
244 dimension = lowerboundary.size();
245
246 for (i = 0; i < dimension; i++) {
247 sum_x.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
248 sum_x_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
249
250 x_av.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
251 sigma_square.push_back(n_vector<double>(lowerboundary, upperboundary, width, Y_SIZE, 0.0));
252 }
253
254 count_y = n_vector<int>(lowerboundary, upperboundary, width, Y_SIZE, 0);
255 distribution_x_y = n_matrix(lowerboundary, upperboundary, width, Y_SIZE);
256
257 grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0));
258 count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0);
259
260 written = false;
261 written_1D = false;
262
263 if (dimension == 1) {
264 std::vector<double> upperboundary_temp = upperboundary;
265 upperboundary_temp[0] = upperboundary[0] + width[0];
266 oneD_pmf = n_vector<double>(lowerboundary, upperboundary_temp, width, 1, 0.0);
267 }
268
269 if (restart == true) {
270 input_grad = n_vector<std::vector<double> >(lowerboundary, upperboundary, width, 1, std::vector<double>(dimension, 0.0));
271 input_count = n_vector<int>(lowerboundary, upperboundary, width, 1, 0);
272
273 // initialize input_Grad and input_count
274 // the loop_flag is a n-dimensional vector, increae from lowerboundary to upperboundary when looping
275 std::vector<double> loop_flag(dimension, 0);
276 for (i = 0; i < dimension; i++) {
277 loop_flag[i] = lowerboundary[i];
278 }
279
280 i = 0;
281 while (i >= 0) {
282 for (j = 0; j < dimension; j++) {
283 input_grad.set_value(loop_flag, std::vector<double>(dimension,0));
284 }
285 input_count.set_value(loop_flag, 0);
286
287 // iterate over any dimensions
288 i = dimension - 1;
289 while (i >= 0) {
290 loop_flag[i] += width[i];
291 if (loop_flag[i] > upperboundary[i] - width[i] + EPSILON) {
292 loop_flag[i] = lowerboundary[i];
293 i--;
294 }
295 else
296 break;
297 }
298 }
299 read_inputfiles(input_filename);
300 }
301 }
302
303 ~UIestimator() {}
304
305 // called from MD engine every step
306 bool update(cvm::step_number /* step */,
307 std::vector<double> x, std::vector<double> y) {
308
309 int i;
310
311 for (i = 0; i < dimension; i++) {
312 // for dihedral RC, it is possible that x = 179 and y = -179, should correct it
313 // may have problem, need to fix
314 if (x[i] > 150 && y[i] < -150) {
315 y[i] += 360;
316 }
317 if (x[i] < -150 && y[i] > 150) {
318 y[i] -= 360;
319 }
320
321 if (x[i] < lowerboundary[i] - EXTENDED_X_SIZE * width[i] + EPSILON || x[i] > upperboundary[i] + EXTENDED_X_SIZE * width[i] - EPSILON \
322 || y[i] - x[i] < -HALF_Y_SIZE * width[i] + EPSILON || y[i] - x[i] > HALF_Y_SIZE * width[i] - EPSILON \
323 || y[i] - lowerboundary[i] < -HALF_Y_SIZE * width[i] + EPSILON || y[i] - upperboundary[i] > HALF_Y_SIZE * width[i] - EPSILON)
324 return false;
325 }
326
327 for (i = 0; i < dimension; i++) {
328 sum_x[i].increase_value(y, x[i]);
329 sum_x_square[i].increase_value(y, x[i] * x[i]);
330 }
331 count_y.increase_value(y, 1);
332
333 for (i = 0; i < dimension; i++) {
334 // adapt colvars precision
335 if (x[i] < lowerboundary[i] + EPSILON || x[i] > upperboundary[i] - EPSILON)
336 return false;
337 }
338 distribution_x_y.increase_value(x, y, 1);
339
340 return true;
341 }
342
343 // update the output_filename
344 void update_output_filename(const std::string& filename) {
345 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
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 public:
382 // calculate gradients from the internal variables
383 void calc_pmf() {
384 colvarproxy *proxy = cvm::main()->proxy;
385
386 int norm;
387 int i, j, k;
388
389 std::vector<double> loop_flag(dimension, 0);
390 for (i = 0; i < dimension; i++) {
391 loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i];
392 }
393
394 i = 0;
395 while (i >= 0) {
396 norm = count_y.get_value(loop_flag) > 0 ? count_y.get_value(loop_flag) : 1;
397 for (j = 0; j < dimension; j++) {
398 x_av[j].set_value(loop_flag, sum_x[j].get_value(loop_flag) / norm);
399 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));
400 }
401
402 // iterate over any dimensions
403 i = dimension - 1;
404 while (i >= 0) {
405 loop_flag[i] += width[i];
406 if (loop_flag[i] > upperboundary[i] + HALF_Y_SIZE * width[i] - width[i] + EPSILON) {
407 loop_flag[i] = lowerboundary[i] - HALF_Y_SIZE * width[i];
408 i--;
409 }
410 else
411 break;
412 }
413 }
414
415 // double integration
416 std::vector<double> av(dimension, 0);
417 std::vector<double> diff_av(dimension, 0);
418
419 std::vector<double> loop_flag_x(dimension, 0);
420 std::vector<double> loop_flag_y(dimension, 0);
421 for (i = 0; i < dimension; i++) {
422 loop_flag_x[i] = lowerboundary[i];
423 loop_flag_y[i] = loop_flag_x[i] - HALF_Y_SIZE * width[i];
424 }
425
426 i = 0;
427 while (i >= 0) {
428 norm = 0;
429 for (k = 0; k < dimension; k++) {
430 av[k] = 0;
431 diff_av[k] = 0;
432 loop_flag_y[k] = loop_flag_x[k] - HALF_Y_SIZE * width[k];
433 }
434
435 j = 0;
436 while (j >= 0) {
437 norm += distribution_x_y.get_value(loop_flag_x, loop_flag_y);
438 for (k = 0; k < dimension; k++) {
439 if (sigma_square[k].get_value(loop_flag_y) > EPSILON || sigma_square[k].get_value(loop_flag_y) < -EPSILON)
440 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);
441
442 diff_av[k] += distribution_x_y.get_value(loop_flag_x, loop_flag_y) * (loop_flag_x[k] - loop_flag_y[k]);
443 }
444
445 // iterate over any dimensions
446 j = dimension - 1;
447 while (j >= 0) {
448 loop_flag_y[j] += width[j];
449 if (loop_flag_y[j] > loop_flag_x[j] + HALF_Y_SIZE * width[j] - width[j] + EPSILON) {
450 loop_flag_y[j] = loop_flag_x[j] - HALF_Y_SIZE * width[j];
451 j--;
452 }
453 else
454 break;
455 }
456 }
457
458 std::vector<double> grad_temp(dimension, 0);
459 for (k = 0; k < dimension; k++) {
460 diff_av[k] /= (norm > 0 ? norm : 1);
461 av[k] = proxy->boltzmann() * temperature * av[k] / (norm > 0 ? norm : 1);
462 grad_temp[k] = av[k] - krestr[k] * diff_av[k];
463 }
464 grad.set_value(loop_flag_x, grad_temp);
465 count.set_value(loop_flag_x, norm);
466
467 // iterate over any dimensions
468 i = dimension - 1;
469 while (i >= 0) {
470 loop_flag_x[i] += width[i];
471 if (loop_flag_x[i] > upperboundary[i] - width[i] + EPSILON) {
472 loop_flag_x[i] = lowerboundary[i];
473 i--;
474 }
475 else
476 break;
477 }
478 }
479 }
480
481
482 // calculate 1D pmf
483 void calc_1D_pmf()
484 {
485 std::vector<double> last_position(1, 0);
486 std::vector<double> position(1, 0);
487
488 double min = 0;
489 double dG = 0;
490 double i;
491
492 oneD_pmf.set_value(lowerboundary, 0);
493 last_position = lowerboundary;
494 for (i = lowerboundary[0] + width[0]; i < upperboundary[0] + EPSILON; i += width[0]) {
495 position[0] = i + EPSILON;
496 if (restart == false || input_count.get_value(last_position) == 0) {
497 dG = oneD_pmf.get_value(last_position) + grad.get_value(last_position)[0] * width[0];
498 }
499 else {
500 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];
501 }
502 if (dG < min)
503 min = dG;
504 oneD_pmf.set_value(position, dG);
505 last_position[0] = i + EPSILON;
506 }
507
508 for (i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) {
509 position[0] = i + EPSILON;
510 oneD_pmf.set_value(position, oneD_pmf.get_value(position) - min);
511 }
512 }
513
514 // write 1D pmf
515 void write_1D_pmf() {
516 std::string pmf_filename = output_filename + ".UI.pmf";
517
518 // only for colvars module!
519 if (written_1D) cvm::backup_file(pmf_filename.c_str());
520
521 std::ostream &ofile_pmf = cvm::proxy->output_stream(pmf_filename,
522 "PMF file");
523
524 std::vector<double> position(1, 0);
525 for (double i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) {
526 ofile_pmf << i << " ";
527 position[0] = i + EPSILON;
528 ofile_pmf << oneD_pmf.get_value(position) << std::endl;
529 }
530 cvm::proxy->close_output_stream(pmf_filename);
531
532 written_1D = true;
533 }
534
535 // write heads of the output files
536 void writehead(std::ostream& os) const {
537 os << "# " << dimension << std::endl;
538 for (int i = 0; i < dimension; i++) {
539 os << "# " << lowerboundary[i] << " " << width[i] << " " << int((upperboundary[i] - lowerboundary[i]) / width[i] + EPSILON) << " " << 0 << std::endl;
540 }
541 os << std::endl;
542 }
543
544 // write interal data, used for testing
545 void write_interal_data() {
546 std::string internal_filename = output_filename + ".UI.internal";
547
548 std::ostream &ofile_internal = cvm::proxy->output_stream(internal_filename,
549 "UI internal file");
550
551 std::vector<double> loop_flag(dimension, 0);
552 for (int i = 0; i < dimension; i++) {
553 loop_flag[i] = lowerboundary[i];
554 }
555
556 int n = 0;
557 while (n >= 0) {
558 for (int j = 0; j < dimension; j++) {
559 ofile_internal << loop_flag[j] + 0.5 * width[j] << " ";
560 }
561
562 for (int k = 0; k < dimension; k++) {
563 ofile_internal << grad.get_value(loop_flag)[k] << " ";
564 }
565
566 std::vector<double> ii(dimension,0);
567 for (double i = loop_flag[0] - 10; i < loop_flag[0] + 10 + EPSILON; i+= width[0]) {
568 for (double j = loop_flag[1] - 10; j< loop_flag[1] + 10 + EPSILON; j+=width[1]) {
569 ii[0] = i;
570 ii[1] = j;
571 ofile_internal << i <<" "<<j<<" "<< distribution_x_y.get_value(loop_flag,ii)<< " ";
572 }
573 }
574 ofile_internal << std::endl;
575
576 // iterate over any dimensions
577 n = dimension - 1;
578 while (n >= 0) {
579 loop_flag[n] += width[n];
580 if (loop_flag[n] > upperboundary[n] - width[n] + EPSILON) {
581 loop_flag[n] = lowerboundary[n];
582 n--;
583 }
584 else
585 break;
586 }
587 }
588 cvm::proxy->close_output_stream(internal_filename.c_str());
589 }
590
591 // write output files
592 void write_files() {
593 std::string grad_filename = output_filename + ".UI.grad";
594 std::string hist_filename = output_filename + ".UI.hist.grad";
595 std::string count_filename = output_filename + ".UI.count";
596
597 int i, j;
598//
599 // only for colvars module!
600 if (written) cvm::backup_file(grad_filename.c_str());
601 //if (written) cvm::backup_file(hist_filename.c_str());
602 if (written) cvm::backup_file(count_filename.c_str());
603
604 std::ostream &ofile = cvm::proxy->output_stream(grad_filename,
605 "gradient file");
606 std::ostream &ofile_hist = cvm::proxy->output_stream(hist_filename,
607 "gradient history file");
608 std::ostream &ofile_count = cvm::proxy->output_stream(count_filename,
609 "count file");
610
611 writehead(ofile);
612 writehead(ofile_hist);
613 writehead(ofile_count);
614
615 if (dimension == 1) {
616 calc_1D_pmf();
617 write_1D_pmf();
618 }
619
620 std::vector<double> loop_flag(dimension, 0);
621 for (i = 0; i < dimension; i++) {
622 loop_flag[i] = lowerboundary[i];
623 }
624
625 i = 0;
626 while (i >= 0) {
627 for (j = 0; j < dimension; j++) {
628 ofile << loop_flag[j] + 0.5 * width[j] << " ";
629 ofile_hist << loop_flag[j] + 0.5 * width[j] << " ";
630 ofile_count << loop_flag[j] + 0.5 * width[j] << " ";
631 }
632
633 if (restart == false) {
634 for (j = 0; j < dimension; j++) {
635 ofile << grad.get_value(loop_flag)[j] << " ";
636 ofile_hist << grad.get_value(loop_flag)[j] << " ";
637 }
638 ofile << std::endl;
639 ofile_hist << std::endl;
640 ofile_count << count.get_value(loop_flag) << " " <<std::endl;
641 }
642 else {
643 double final_grad = 0;
644 for (j = 0; j < dimension; j++) {
645 int total_count_temp = (count.get_value(loop_flag) + input_count.get_value(loop_flag));
646 if (input_count.get_value(loop_flag) == 0)
647 final_grad = grad.get_value(loop_flag)[j];
648 else
649 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);
650 ofile << final_grad << " ";
651 ofile_hist << final_grad << " ";
652 }
653 ofile << std::endl;
654 ofile_hist << std::endl;
655 ofile_count << (count.get_value(loop_flag) + input_count.get_value(loop_flag)) << " " <<std::endl;
656 }
657
658 // iterate over any dimensions
659 i = dimension - 1;
660 while (i >= 0) {
661 loop_flag[i] += width[i];
662 if (loop_flag[i] > upperboundary[i] - width[i] + EPSILON) {
663 loop_flag[i] = lowerboundary[i];
664 i--;
665 ofile << std::endl;
666 ofile_hist << std::endl;
667 ofile_count << std::endl;
668 }
669 else
670 break;
671 }
672 }
673 cvm::proxy->close_output_stream(grad_filename.c_str());
674 // cvm::proxy->close_output_stream(hist_filename.c_str());
675 cvm::proxy->close_output_stream(count_filename.c_str());
676
677 written = true;
678 }
679
680 // read input files
681 void read_inputfiles(const std::vector<std::string> filename)
682 {
683 char sharp;
684 double nothing;
685 int dimension_temp;
686 int i, j, k, l, m;
687
688 colvarproxy *proxy = cvm::main()->proxy;
689 std::vector<double> loop_bin_size(dimension, 0);
690 std::vector<double> position_temp(dimension, 0);
691 std::vector<double> grad_temp(dimension, 0);
692 int count_temp = 0;
693 for (i = 0; i < int(filename.size()); i++) {
694 int size = 1 , size_temp = 0;
695
696 std::string count_filename = filename[i] + ".UI.count";
697 std::string grad_filename = filename[i] + ".UI.grad";
698
699 std::istream &count_file =
700 proxy->input_stream(count_filename, "count filename");
701 std::istream &grad_file =
702 proxy->input_stream(grad_filename, "gradient filename");
703
704 if (!count_file || !grad_file) {
705 return;
706 }
707
708 count_file >> sharp >> dimension_temp;
709 grad_file >> sharp >> dimension_temp;
710
711 for (j = 0; j < dimension; j++) {
712 count_file >> sharp >> nothing >> nothing >> size_temp >> nothing;
713 grad_file >> sharp >> nothing >> nothing >> nothing >> nothing;
714 size *= size_temp;
715 }
716
717 for (j = 0; j < size; j++) {
718 do {
719 for (k = 0; k < dimension; k++) {
720 count_file >> position_temp[k];
721 grad_file >> nothing;
722 }
723
724 for (l = 0; l < dimension; l++) {
725 grad_file >> grad_temp[l];
726 }
727 count_file >> count_temp;
728 }
729 while (position_temp[i] < lowerboundary[i] - EPSILON || position_temp[i] > upperboundary[i] + EPSILON);
730
731 if (count_temp == 0) {
732 continue;
733 }
734
735 for (m = 0; m < dimension; m++) {
736 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));
737 }
738 input_grad.set_value(position_temp, grad_temp);
739 input_count.increase_value(position_temp, count_temp);
740 }
741
742 proxy->close_input_stream(count_filename);
743 proxy->close_input_stream(grad_filename);
744 }
745 }
746 };
747}
748
749#endif
Definition: colvar_UIestimator.h:215
Definition: colvar_UIestimator.h:33
Definition: colvar_UIestimator.h:140
static real pow(real const &x, real const &y)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:115
static int backup_file(char const *filename)
Backup a file before writing it.
Definition: colvarmodule.cpp:1738
static colvarmodule * main()
Access the one instance of the Colvars module.
Definition: colvarmodule.cpp:174
static colvarproxy * proxy
Pointer to the proxy object, used to retrieve atomic data from the hosting program; it is static in o...
Definition: colvarmodule.h:860
long long step_number
Use a 64-bit integer to store the step number.
Definition: colvarmodule.h:92
virtual int close_output_stream(std::string const &output_name)
Closes the given output file/channel.
Definition: colvarproxy_io.cpp:410
virtual std::ostream & output_stream(std::string const &output_name, std::string const description)
Definition: colvarproxy_io.cpp:343
Definition: colvarproxy.h:573
Collective variables main module.
Colvars proxy classes.