Collective Variables Module - Developer Documentation
Loading...
Searching...
No Matches
colvargrid.h
1// -*- c++ -*-
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 COLVARGRID_H
11#define COLVARGRID_H
12
13#include <iosfwd>
14#include <memory>
15
16#include "colvar.h"
17#include "colvarmodule.h"
18#include "colvarvalue.h"
19#include "colvarparse.h"
20
21
26template <class T> class colvar_grid : public colvarparse {
27
28 //protected:
29public: // TODO create accessors for these after all instantiations work
30
32 size_t nd;
33
35 std::vector<int> nx;
36
38 std::vector<int> nxc;
39
42 size_t mult;
43
45 size_t nt;
46
48 std::vector<T> data;
49
51 std::vector<size_t> new_data;
52
54 std::vector<colvar *> cv;
55
57 std::vector<bool> use_actual_value;
58
60 inline size_t address(std::vector<int> const &ix) const
61 {
62 size_t addr = 0;
63 for (size_t i = 0; i < nd; i++) {
64 addr += ix[i]*static_cast<size_t>(nxc[i]);
65 if (cvm::debug()) {
66 if (ix[i] >= nx[i]) {
67 cvm::error("Error: exceeding bounds in colvar_grid.\n", COLVARS_BUG_ERROR);
68 return 0;
69 }
70 }
71 }
72 return addr;
73 }
74
75public:
76
78 std::vector<colvarvalue> lower_boundaries;
79
81 std::vector<colvarvalue> upper_boundaries;
82
84 std::vector<bool> periodic;
85
87 std::vector<bool> hard_lower_boundaries;
88
90 std::vector<bool> hard_upper_boundaries;
91
93 std::vector<cvm::real> widths;
94
97
100
102 inline size_t num_variables() const
103 {
104 return nd;
105 }
106
108 inline std::vector<int> const &number_of_points_vec() const
109 {
110 return nx;
111 }
112
115 inline size_t number_of_points(int const icv = -1) const
116 {
117 if (icv < 0) {
118 return nt;
119 } else {
120 return nx[icv];
121 }
122 }
123
125 inline std::vector<int> const & sizes() const
126 {
127 return nx;
128 }
129
131 inline void set_sizes(std::vector<int> const &new_sizes)
132 {
133 nx = new_sizes;
134 }
135
137 inline size_t multiplicity() const
138 {
139 return mult;
140 }
141
143 inline void request_actual_value(bool b = true)
144 {
145 size_t i;
146 for (i = 0; i < use_actual_value.size(); i++) {
147 use_actual_value[i] = b;
148 }
149 }
150
152 int setup(std::vector<int> const &nx_i,
153 T const &t = T(),
154 size_t const &mult_i = 1)
155 {
156 if (cvm::debug()) {
157 cvm::log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+
158 ", dimensions = "+cvm::to_str(nx_i)+".\n");
159 }
160
161 mult = mult_i;
162
163 data.clear();
164
165 nx = nx_i;
166 nd = nx.size();
167
168 nxc.resize(nd);
169
170 // setup dimensions
171 nt = mult;
172 for (int i = nd-1; i >= 0; i--) {
173 if (nx[i] <= 0) {
174 cvm::error("Error: providing an invalid number of grid points, "+
175 cvm::to_str(nx[i])+".\n", COLVARS_BUG_ERROR);
176 return COLVARS_ERROR;
177 }
178 nxc[i] = nt;
179 nt *= nx[i];
180 }
181
182 if (cvm::debug()) {
183 cvm::log("Total number of grid elements = "+cvm::to_str(nt)+".\n");
184 }
185
186 data.reserve(nt);
187 data.assign(nt, t);
188
189 return COLVARS_OK;
190 }
191
193 int setup()
194 {
195 return setup(this->nx, T(), this->mult);
196 }
197
199 void reset(T const &t = T())
200 {
201 data.assign(nt, t);
202 }
203
204
207 {
208 nd = nt = 0;
209 mult = 1;
210 has_parent_data = false;
211 this->setup();
212 }
213
215 virtual ~colvar_grid()
216 {}
217
222 nd(g.nd),
223 nx(g.nx),
224 mult(g.mult),
225 data(),
226 cv(g.cv),
233 widths(g.widths),
234 has_parent_data(false),
235 has_data(false)
236 {}
237
242 colvar_grid(std::vector<int> const &nx_i,
243 T const &t = T(),
244 size_t mult_i = 1)
245 : has_parent_data(false), has_data(false)
246 {
247 this->setup(nx_i, t, mult_i);
248 }
249
253 colvar_grid(std::vector<colvar *> const &colvars,
254 T const &t = T(),
255 size_t mult_i = 1,
256 bool add_extra_bin = false)
257 : has_parent_data(false), has_data(false)
258 {
259 (void) t;
260 this->init_from_colvars(colvars, mult_i, add_extra_bin);
261 }
262
263 int init_from_colvars(std::vector<colvar *> const &colvars,
264 size_t mult_i = 1,
265 bool add_extra_bin = false)
266 {
267 if (cvm::debug()) {
268 cvm::log("Reading grid configuration from collective variables.\n");
269 }
270
271 cv = colvars;
272 nd = colvars.size();
273 mult = mult_i;
274
275 size_t i;
276
277 if (cvm::debug()) {
278 cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+
279 " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");
280 }
281
282 for (i = 0; i < cv.size(); i++) {
283
284 if (cv[i]->value().type() != colvarvalue::type_scalar) {
285 cvm::error("Colvar grids can only be automatically "
286 "constructed for scalar variables. "
287 "ABF and histogram can not be used; metadynamics "
288 "can be used with useGrids disabled.\n", COLVARS_INPUT_ERROR);
289 return COLVARS_ERROR;
290 }
291
292 if (cv[i]->width <= 0.0) {
293 cvm::error("Tried to initialize a grid on a "
294 "variable with negative or zero width.\n", COLVARS_INPUT_ERROR);
295 return COLVARS_ERROR;
296 }
297
298 widths.push_back(cv[i]->width);
301 periodic.push_back(cv[i]->periodic_boundaries());
302
303 // By default, get reported colvar value (for extended Lagrangian colvars)
304 use_actual_value.push_back(false);
305
306 // except if a colvar is specified twice in a row
307 // then the first instance is the actual value
308 // For histograms of extended-system coordinates
309 if (i > 0 && cv[i-1] == cv[i]) {
310 use_actual_value[i-1] = true;
311 }
312
313 if (add_extra_bin) {
314 if (periodic[i]) {
315 // Shift the grid by half the bin width (values at edges instead of center of bins)
316 lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
317 upper_boundaries.push_back(cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
318 } else {
319 // Make this grid larger by one bin width
320 lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
321 upper_boundaries.push_back(cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
322 }
323 } else {
324 lower_boundaries.push_back(cv[i]->lower_boundary);
325 upper_boundaries.push_back(cv[i]->upper_boundary);
326 }
327 }
328
329 this->init_from_boundaries();
330 return this->setup();
331 }
332
333 int init_from_boundaries()
334 {
335 if (cvm::debug()) {
336 cvm::log("Configuring grid dimensions from colvars boundaries.\n");
337 }
338
339 // these will have to be updated
340 nx.clear();
341 nxc.clear();
342 nt = 0;
343
344 for (size_t i = 0; i < lower_boundaries.size(); i++) {
345 // Re-compute periodicity using current grid boundaries
346 periodic[i] = cv[i]->periodic_boundaries(lower_boundaries[i].real_value,
347 upper_boundaries[i].real_value);
348
349 cvm::real nbins = ( upper_boundaries[i].real_value -
350 lower_boundaries[i].real_value ) / widths[i];
351 int nbins_round = (int)(nbins+0.5);
352
353 if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
354 cvm::log("Warning: grid interval("+
357 ") is not commensurate to its bin width("+
359 upper_boundaries[i].real_value = lower_boundaries[i].real_value +
360 (nbins_round * widths[i]);
361 }
362
363 if (cvm::debug())
364 cvm::log("Number of points is "+cvm::to_str((int) nbins_round)+
365 " for the colvar no. "+cvm::to_str(i+1)+".\n");
366
367 nx.push_back(nbins_round);
368 }
369
370 return COLVARS_OK;
371 }
372
375 inline void wrap(std::vector<int> & ix) const
376 {
377 for (size_t i = 0; i < nd; i++) {
378 if (periodic[i]) {
379 ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined)
380 } else {
381 if (ix[i] < 0 || ix[i] >= nx[i]) {
382 cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: "
383 + cvm::to_str(ix), COLVARS_BUG_ERROR);
384 return;
385 }
386 }
387 }
388 }
389
392 inline bool wrap_detect_edge(std::vector<int> & ix) const
393 {
394 bool edge = false;
395 for (size_t i = 0; i < nd; i++) {
396 if (periodic[i]) {
397 ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined)
398 } else if (ix[i] < 0 || ix[i] >= nx[i]) {
399 edge = true;
400 }
401 }
402 return edge;
403 }
404
407 inline bool wrap_to_edge(std::vector<int> & ix, std::vector<int> & edge_bin) const
408 {
409 bool edge = false;
410 edge_bin = ix;
411 for (size_t i = 0; i < nd; i++) {
412 if (periodic[i]) {
413 ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined)
414 edge_bin[i] = ix[i];
415 } else if (ix[i] < 0) {
416 edge = true;
417 edge_bin[i] = 0;
418 } else if (ix[i] >= nx[i]) {
419 edge = true;
420 edge_bin[i] = nx[i] - 1;
421 }
422 }
423 return edge;
424 }
425
426
428 inline int current_bin_scalar(int const i) const
429 {
430 return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
431 }
432
435 inline int current_bin_flat_bound() const
436 {
437 std::vector<int> index = new_index();
438 for (size_t i = 0; i < nd; i++) {
439 index[i] = current_bin_scalar_bound(i);
440 }
441 return address(index);
442 }
443
446 inline int current_bin_scalar_bound(int const i) const
447 {
448 return value_to_bin_scalar_bound(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
449 }
450
452 inline int current_bin_scalar(int const i, int const iv) const
453 {
455 cv[i]->actual_value().vector1d_value[iv] :
456 cv[i]->value().vector1d_value[iv], i);
457 }
458
461 inline int value_to_bin_scalar(colvarvalue const &value, const int i) const
462 {
463 return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
464 }
465
467 inline cvm::real current_bin_scalar_fraction(int const i) const
468 {
469 return value_to_bin_scalar_fraction(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
470 }
471
474 inline cvm::real value_to_bin_scalar_fraction(colvarvalue const &value, const int i) const
475 {
476 cvm::real x = (value.real_value - lower_boundaries[i].real_value) / widths[i];
477 return x - cvm::floor(x);
478 }
479
482 inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
483 {
484 int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
485
486 // Wrap bins for periodic dimensions before truncating
487 if (periodic[i]) bin_index %= nx[i];
488 if (bin_index < 0) bin_index=0;
489 if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
490 return (int) bin_index;
491 }
492
494 inline int value_to_bin_scalar(colvarvalue const &value,
495 colvarvalue const &new_offset,
496 cvm::real const &new_width) const
497 {
498 return (int) cvm::floor( (value.real_value - new_offset.real_value) / new_width );
499 }
500
503 inline colvarvalue bin_to_value_scalar(int const &i_bin, int const i) const
504 {
505 return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
506 }
507
509 inline colvarvalue bin_to_value_scalar(int const &i_bin,
510 colvarvalue const &new_offset,
511 cvm::real const &new_width) const
512 {
513 return new_offset.real_value + new_width * (0.5 + i_bin);
514 }
515
517 inline void set_value(std::vector<int> const &ix,
518 T const &t,
519 size_t const &imult = 0)
520 {
521 data[this->address(ix)+imult] = t;
522 has_data = true;
523 }
524
526 inline void set_value(size_t i, T const &t)
527 {
528 data[i] = t;
529 }
530
532 inline T get_value(size_t i) const
533 {
534 return data[i];
535 }
536
537
542 void delta_grid(colvar_grid<T> const &other_grid)
543 {
544
545 if (other_grid.multiplicity() != this->multiplicity()) {
546 cvm::error("Error: trying to subtract two grids with "
547 "different multiplicity.\n");
548 return;
549 }
550
551 if (other_grid.data.size() != this->data.size()) {
552 cvm::error("Error: trying to subtract two grids with "
553 "different size.\n");
554 return;
555 }
556
557 for (size_t i = 0; i < data.size(); i++) {
558 data[i] = other_grid.data[i] - data[i];
559 }
560 has_data = true;
561 }
562
563
567 void copy_grid(colvar_grid<T> const &other_grid)
568 {
569 if (other_grid.multiplicity() != this->multiplicity()) {
570 cvm::error("Error: trying to copy two grids with "
571 "different multiplicity.\n");
572 return;
573 }
574
575 if (other_grid.data.size() != this->data.size()) {
576 cvm::error("Error: trying to copy two grids with "
577 "different size.\n");
578 return;
579 }
580
581
582 for (size_t i = 0; i < data.size(); i++) {
583 data[i] = other_grid.data[i];
584 }
585 has_data = true;
586 }
587
590 void raw_data_out(T* out_data) const
591 {
592 for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i];
593 }
594 void raw_data_out(std::vector<T>& out_data) const
595 {
596 out_data = data;
597 }
599 void raw_data_in(const T* in_data)
600 {
601 for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i];
602 has_data = true;
603 }
604 void raw_data_in(const std::vector<T>& in_data)
605 {
606 data = in_data;
607 has_data = true;
608 }
610 size_t raw_data_num() const { return data.size(); }
611
612
615 inline T const & value(std::vector<int> const &ix,
616 size_t const &imult = 0) const
617 {
618 return data[this->address(ix) + imult];
619 }
620
622 inline T const & value(size_t i) const
623 {
624 return data[i];
625 }
626
628 inline void add_constant(T const &t)
629 {
630 for (size_t i = 0; i < nt; i++)
631 data[i] += t;
632 has_data = true;
633 }
634
636 inline void multiply_constant(cvm::real const &a)
637 {
638 for (size_t i = 0; i < nt; i++)
639 data[i] *= a;
640 }
641
643 inline void remove_small_values(cvm::real const &a)
644 {
645 for (size_t i = 0; i < nt; i++)
646 if(data[i]<a) data[i] = a;
647 }
648
649
652 inline std::vector<int> const get_colvars_index(std::vector<colvarvalue> const &values) const
653 {
654 std::vector<int> index = new_index();
655 for (size_t i = 0; i < nd; i++) {
656 index[i] = value_to_bin_scalar(values[i], i);
657 }
658 return index;
659 }
660
663 inline std::vector<int> const get_colvars_index() const
664 {
665 std::vector<int> index = new_index();
666 for (size_t i = 0; i < nd; i++) {
667 index[i] = current_bin_scalar(i);
668 }
669 return index;
670 }
671
674 inline std::vector<int> const get_colvars_index_bound() const
675 {
676 std::vector<int> index = new_index();
677 for (size_t i = 0; i < nd; i++) {
678 index[i] = current_bin_scalar_bound(i);
679 }
680 return index;
681 }
682
686 inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
687 bool skip_hard_boundaries = false)
688 {
689 cvm::real minimum = 1.0E+16;
690 for (size_t i = 0; i < nd; i++) {
691
692 if (periodic[i]) continue;
693
694 cvm::real dl = cvm::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i];
695 cvm::real du = cvm::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i];
696
697 if (values[i].real_value < lower_boundaries[i])
698 dl *= -1.0;
699 if (values[i].real_value > upper_boundaries[i])
700 du *= -1.0;
701
702 if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum))
703 minimum = dl;
704 if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum))
705 minimum = du;
706 }
707
708 return minimum;
709 }
710
711
716 void map_grid(colvar_grid<T> const &other_grid)
717 {
718 if (other_grid.multiplicity() != this->multiplicity()) {
719 cvm::error("Error: trying to merge two grids with values of "
720 "different multiplicity.\n");
721 return;
722 }
723
724 std::vector<colvarvalue> const &gb = this->lower_boundaries;
725 std::vector<cvm::real> const &gw = this->widths;
726 std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
727 std::vector<cvm::real> const &ogw = other_grid.widths;
728
729 std::vector<int> ix = this->new_index();
730 std::vector<int> oix = other_grid.new_index();
731
732 if (cvm::debug())
733 cvm::log("Remapping grid...\n");
734 for ( ; this->index_ok(ix); this->incr(ix)) {
735
736 for (size_t i = 0; i < nd; i++) {
737 oix[i] =
738 value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
739 ogb[i],
740 ogw[i]);
741 }
742
743 if (! other_grid.index_ok(oix)) {
744 continue;
745 }
746
747 for (size_t im = 0; im < mult; im++) {
748 this->set_value(ix, other_grid.value(oix, im), im);
749 }
750 }
751
752 has_data = true;
753 if (cvm::debug())
754 cvm::log("Remapping done.\n");
755 }
756
759 void add_grid(colvar_grid<T> const &other_grid,
760 cvm::real scale_factor = 1.0)
761 {
762 if (other_grid.multiplicity() != this->multiplicity()) {
763 cvm::error("Error: trying to sum togetehr two grids with values of "
764 "different multiplicity.\n");
765 return;
766 }
767 if (scale_factor != 1.0)
768 for (size_t i = 0; i < data.size(); i++) {
769 data[i] += static_cast<T>(scale_factor * other_grid.data[i]);
770 }
771 else
772 // skip multiplication if possible
773 for (size_t i = 0; i < data.size(); i++) {
774 data[i] += other_grid.data[i];
775 }
776 has_data = true;
777 }
778
781 virtual T value_output(std::vector<int> const &ix,
782 size_t const &imult = 0) const
783 {
784 return value(ix, imult);
785 }
786
790 virtual void value_input(std::vector<int> const &ix,
791 T const &t,
792 size_t const &imult = 0,
793 bool add = false)
794 {
795 if ( add )
796 data[address(ix) + imult] += t;
797 else
798 data[address(ix) + imult] = t;
799 has_data = true;
800 }
801
802
803 // /// Get the pointer to the binned value indexed by ix
804 // inline T const *value_p (std::vector<int> const &ix)
805 // {
806 // return &(data[address (ix)]);
807 // }
808
811 inline std::vector<int> const new_index() const
812 {
813 return std::vector<int> (nd, 0);
814 }
815
818 inline bool index_ok(std::vector<int> const &ix) const
819 {
820 for (size_t i = 0; i < nd; i++) {
821 if ( (ix[i] < 0) || (ix[i] >= int(nx[i])) )
822 return false;
823 }
824 return true;
825 }
826
829 inline void incr(std::vector<int> &ix) const
830 {
831 for (int i = ix.size()-1; i >= 0; i--) {
832
833 ix[i]++;
834
835 if (ix[i] >= nx[i]) {
836
837 if (i > 0) {
838 ix[i] = 0;
839 continue;
840 } else {
841 // this is the last iteration, a non-valid index is being
842 // set for the outer index, which will be caught by
843 // index_ok()
844 ix[0] = nx[0];
845 return;
846 }
847 } else {
848 return;
849 }
850 }
851 }
852
854 std::string get_state_params() const;
855
857 int parse_params(std::string const &conf,
859
864 {
865 for (size_t i = 0; i < nd; i++) {
866 if ( (cvm::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
867 lower_boundaries[i])) > 1.0E-10) ||
868 (cvm::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
869 upper_boundaries[i])) > 1.0E-10) ||
870 (cvm::sqrt(cv[i]->dist2(cv[i]->width,
871 widths[i])) > 1.0E-10) ) {
872 cvm::error("Error: restart information for a grid is "
873 "inconsistent with that of its colvars.\n");
874 return;
875 }
876 }
877 }
878
879
882 void check_consistency(colvar_grid<T> const &other_grid)
883 {
884 for (size_t i = 0; i < nd; i++) {
885 // we skip dist2(), because periodicities and the like should
886 // matter: boundaries should be EXACTLY the same (otherwise,
887 // map_grid() should be used)
888 if ( (cvm::fabs(other_grid.lower_boundaries[i] -
889 lower_boundaries[i]) > 1.0E-10) ||
890 (cvm::fabs(other_grid.upper_boundaries[i] -
891 upper_boundaries[i]) > 1.0E-10) ||
892 (cvm::fabs(other_grid.widths[i] -
893 widths[i]) > 1.0E-10) ||
894 (data.size() != other_grid.data.size()) ) {
895 cvm::error("Error: inconsistency between "
896 "two grids that are supposed to be equal, "
897 "aside from the data stored.\n");
898 return;
899 }
900 }
901 }
902
904 std::istream & read_restart(std::istream &is);
905
908
910 std::ostream & write_restart(std::ostream &os);
911
914
916 std::istream &read_raw(std::istream &is);
917
920
924 std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
925
929 cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
930
932 std::istream & read_multicol(std::istream &is, bool add = false);
933
935 int read_multicol(std::string const &filename,
936 std::string description = "grid file",
937 bool add = false);
938
940 std::ostream & write_multicol(std::ostream &os) const;
941
943 int write_multicol(std::string const &filename,
944 std::string description = "grid file") const;
945
947 std::ostream & write_opendx(std::ostream &os) const;
948
950 int write_opendx(std::string const &filename,
951 std::string description = "grid file") const;
952};
953
954
955
958class colvar_grid_count : public colvar_grid<size_t>
959{
960public:
961
964
967 {}
968
970 colvar_grid_count(std::vector<int> const &nx_i,
971 size_t const &def_count = 0);
972
974 colvar_grid_count(std::vector<colvar *> &colvars,
975 size_t const &def_count = 0,
976 bool add_extra_bin = false);
977
979 inline void incr_count(std::vector<int> const &ix)
980 {
981 ++(data[this->address(ix)]);
982 }
983
985 inline size_t const & new_value(std::vector<int> const &ix)
986 {
987 return new_data[address(ix)];
988 }
989
991 std::string get_state_params() const;
992
994 int parse_params(std::string const &conf,
996
998 std::istream & read_restart(std::istream &is);
999
1002
1004 std::ostream & write_restart(std::ostream &os);
1005
1008
1010 std::istream &read_raw(std::istream &is);
1011
1014
1018 std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
1019
1023 cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
1024
1026 std::istream & read_multicol(std::istream &is, bool add = false);
1027
1029 int read_multicol(std::string const &filename,
1030 std::string description = "grid file",
1031 bool add = false);
1032
1034 std::ostream & write_multicol(std::ostream &os) const;
1035
1037 int write_multicol(std::string const &filename,
1038 std::string description = "grid file") const;
1039
1041 std::ostream & write_opendx(std::ostream &os) const;
1042
1044 int write_opendx(std::string const &filename,
1045 std::string description = "grid file") const;
1046
1048 virtual void value_input(std::vector<int> const &ix,
1049 size_t const &t,
1050 size_t const &imult = 0,
1051 bool add = false)
1052 {
1053 (void) imult;
1054 if (add) {
1055 data[address(ix)] += t;
1056 if (this->has_parent_data) {
1057 // save newly read data for inputting parent grid
1058 new_data[address(ix)] = t;
1059 }
1060 } else {
1061 data[address(ix)] = t;
1062 }
1063 has_data = true;
1064 }
1065
1068 inline int local_sample_count(int radius)
1069 {
1070 std::vector<int> ix0 = new_index();
1071 std::vector<int> ix = new_index();
1072
1073 for (size_t i = 0; i < nd; i++) {
1074 ix0[i] = current_bin_scalar_bound(i);
1075 }
1076 if (radius < 1) {
1077 // Simple case: no averaging
1078 if (index_ok(ix0))
1079 return value(ix0);
1080 else
1081 return 0;
1082 }
1083 size_t count = 0;
1084 size_t nbins = 0;
1085 int i, j, k;
1086 bool edge;
1087 ix = ix0;
1088 // Treat each dimension separately to simplify code
1089 switch (nd)
1090 {
1091 case 1:
1092 for (i = -radius; i <= radius; i++) {
1093 ix[0] = ix0[0] + i;
1094 edge = wrap_detect_edge(ix);
1095 if (!edge) {
1096 nbins++;
1097 count += value(ix);
1098 }
1099 }
1100 break;
1101 case 2:
1102 for (i = -radius; i <= radius; i++) {
1103 ix[0] = ix0[0] + i;
1104 for (j = -radius; j <= radius; j++) {
1105 ix[1] = ix0[1] + j;
1106 edge = wrap_detect_edge(ix);
1107 if (!edge) {
1108 nbins++;
1109 count += value(ix);
1110 }
1111 }
1112 }
1113 break;
1114 case 3:
1115 for (i = -radius; i <= radius; i++) {
1116 ix[0] = ix0[0] + i;
1117 for (j = -radius; j <= radius; j++) {
1118 ix[1] = ix0[1] + j;
1119 for (k = -radius; k <= radius; k++) {
1120 ix[2] = ix0[2] + k;
1121 edge = wrap_detect_edge(ix);
1122 if (!edge) {
1123 nbins++;
1124 count += value(ix);
1125 }
1126 }
1127 }
1128 }
1129 break;
1130 default:
1131 cvm::error("Error: local_sample_count is not implemented for grids of dimension > 3", COLVARS_NOT_IMPLEMENTED);
1132 break;
1133 }
1134
1135 if (nbins)
1136 // Integer division - an error on the order of 1 doesn't matter
1137 return count / nbins;
1138 else
1139 return 0.0;
1140 }
1141
1142
1146 inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
1147 int n = 0, int offset = 0)
1148 {
1149 cvm::real A0, A1, A2;
1150 std::vector<int> ix = ix0;
1151
1152 // TODO this can be rewritten more concisely with wrap_edge()
1153 if (periodic[n]) {
1154 ix[n]--; wrap(ix);
1155 A0 = value(ix) + offset;
1156 ix = ix0;
1157 ix[n]++; wrap(ix);
1158 A1 = value(ix) + offset;
1159 if (A0 * A1 == 0) {
1160 return 0.; // can't handle empty bins
1161 } else {
1162 return (cvm::logn(A1) - cvm::logn(A0))
1163 / (widths[n] * 2.);
1164 }
1165 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
1166 ix[n]--;
1167 A0 = value(ix) + offset;
1168 ix = ix0;
1169 ix[n]++;
1170 A1 = value(ix) + offset;
1171 if (A0 * A1 == 0) {
1172 return 0.; // can't handle empty bins
1173 } else {
1174 return (cvm::logn(A1) - cvm::logn(A0))
1175 / (widths[n] * 2.);
1176 }
1177 } else {
1178 // edge: use 2nd order derivative
1179 int increment = (ix[n] == 0 ? 1 : -1);
1180 // move right from left edge, or the other way around
1181 A0 = value(ix) + offset;
1182 ix[n] += increment; A1 = value(ix) + offset;
1183 ix[n] += increment; A2 = value(ix) + offset;
1184 if (A0 * A1 * A2 == 0) {
1185 return 0.; // can't handle empty bins
1186 } else {
1187 return (-1.5 * cvm::logn(A0) + 2. * cvm::logn(A1)
1188 - 0.5 * cvm::logn(A2)) * increment / widths[n];
1189 }
1190 }
1191 }
1192
1193
1197 inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
1198 int n = 0)
1199 {
1200 cvm::real A0, A1, A2;
1201 std::vector<int> ix = ix0;
1202
1203 // FIXME this can be rewritten more concisely with wrap_edge()
1204 if (periodic[n]) {
1205 ix[n]--; wrap(ix);
1206 A0 = value(ix);
1207 ix = ix0;
1208 ix[n]++; wrap(ix);
1209 A1 = value(ix);
1210 if (A0 * A1 == 0) {
1211 return 0.; // can't handle empty bins
1212 } else {
1213 return (A1 - A0) / (widths[n] * 2.);
1214 }
1215 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
1216 ix[n]--;
1217 A0 = value(ix);
1218 ix = ix0;
1219 ix[n]++;
1220 A1 = value(ix);
1221 if (A0 * A1 == 0) {
1222 return 0.; // can't handle empty bins
1223 } else {
1224 return (A1 - A0) / (widths[n] * 2.);
1225 }
1226 } else {
1227 // edge: use 2nd order derivative
1228 int increment = (ix[n] == 0 ? 1 : -1);
1229 // move right from left edge, or the other way around
1230 A0 = value(ix);
1231 ix[n] += increment; A1 = value(ix);
1232 ix[n] += increment; A2 = value(ix);
1233 return (-1.5 * A0 + 2. * A1
1234 - 0.5 * A2) * increment / widths[n];
1235 }
1236 }
1237};
1238
1239
1241class colvar_grid_scalar : public colvar_grid<cvm::real>
1242{
1243public:
1244
1248
1251
1254
1256 virtual ~colvar_grid_scalar();
1257
1259 colvar_grid_scalar(std::vector<int> const &nx_i);
1260
1262 colvar_grid_scalar(std::vector<colvar *> &colvars,
1263 bool add_extra_bin = false);
1264
1266 inline void acc_value(std::vector<int> const &ix,
1267 cvm::real const &new_value,
1268 size_t const &imult = 0)
1269 {
1270 (void) imult;
1271 // only legal value of imult here is 0
1272 data[address(ix)] += new_value;
1273 if (samples)
1274 samples->incr_count(ix);
1275 has_data = true;
1276 }
1277
1279 std::string get_state_params() const;
1280
1282 int parse_params(std::string const &conf,
1284
1286 std::istream & read_restart(std::istream &is);
1287
1290
1292 std::ostream & write_restart(std::ostream &os);
1293
1296
1298 std::istream &read_raw(std::istream &is);
1299
1302
1306 std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
1307
1311 cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
1312
1314 std::istream & read_multicol(std::istream &is, bool add = false);
1315
1317 int read_multicol(std::string const &filename,
1318 std::string description = "grid file",
1319 bool add = false);
1320
1322 std::ostream & write_multicol(std::ostream &os) const;
1323
1325 int write_multicol(std::string const &filename,
1326 std::string description = "grid file") const;
1327
1329 std::ostream & write_opendx(std::ostream &os) const;
1330
1332 int write_opendx(std::string const &filename,
1333 std::string description = "grid file") const;
1334
1339 inline void vector_gradient_finite_diff( const std::vector<int> &ix0, std::vector<cvm::real> &grad)
1340 {
1341 cvm::real A0, A1;
1342 std::vector<int> ix;
1343 size_t i, j, k, n;
1344
1345 if (nd == 2) {
1346 for (n = 0; n < 2; n++) {
1347 ix = ix0;
1348 A0 = value(ix);
1349 ix[n]++; wrap(ix);
1350 A1 = value(ix);
1351 ix[1-n]++; wrap(ix);
1352 A1 += value(ix);
1353 ix[n]--; wrap(ix);
1354 A0 += value(ix);
1355 grad[n] = 0.5 * (A1 - A0) / widths[n];
1356 }
1357 } else if (nd == 3) {
1358
1359 cvm::real p[8]; // potential values within cube, indexed in binary (4 i + 2 j + k)
1360 ix = ix0;
1361 int index = 0;
1362 for (i = 0; i<2; i++) {
1363 ix[1] = ix0[1];
1364 for (j = 0; j<2; j++) {
1365 ix[2] = ix0[2];
1366 for (k = 0; k<2; k++) {
1367 wrap(ix);
1368 p[index++] = value(ix);
1369 ix[2]++;
1370 }
1371 ix[1]++;
1372 }
1373 ix[0]++;
1374 }
1375
1376 // The following would be easier to read using binary literals
1377 // 100 101 110 111 000 001 010 011
1378 grad[0] = 0.25 * ((p[4] + p[5] + p[6] + p[7]) - (p[0] + p[1] + p[2] + p[3])) / widths[0];
1379 // 010 011 110 111 000 001 100 101
1380 grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[1];
1381 // 001 011 101 111 000 010 100 110
1382 grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[2];
1383 } else {
1384 cvm::error("Finite differences available in dimension 2 and 3 only.");
1385 }
1386 }
1387
1388
1392 inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
1393 int n = 0, int offset = 0)
1394 {
1395 cvm::real A0, A1, A2;
1396 std::vector<int> ix = ix0;
1397
1398 // TODO this can be rewritten more concisely with wrap_edge()
1399 if (periodic[n]) {
1400 ix[n]--; wrap(ix);
1401 A0 = value(ix) + offset;
1402 ix = ix0;
1403 ix[n]++; wrap(ix);
1404 A1 = value(ix) + offset;
1405 if (A0 * A1 == 0) {
1406 return 0.; // can't handle empty bins
1407 } else {
1408 return (cvm::logn(A1) - cvm::logn(A0))
1409 / (widths[n] * 2.);
1410 }
1411 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
1412 ix[n]--;
1413 A0 = value(ix) + offset;
1414 ix = ix0;
1415 ix[n]++;
1416 A1 = value(ix) + offset;
1417 if (A0 * A1 == 0) {
1418 return 0.; // can't handle empty bins
1419 } else {
1420 return (cvm::logn(A1) - cvm::logn(A0))
1421 / (widths[n] * 2.);
1422 }
1423 } else {
1424 // edge: use 2nd order derivative
1425 int increment = (ix[n] == 0 ? 1 : -1);
1426 // move right from left edge, or the other way around
1427 A0 = value(ix) + offset;
1428 ix[n] += increment; A1 = value(ix) + offset;
1429 ix[n] += increment; A2 = value(ix) + offset;
1430 if (A0 * A1 * A2 == 0) {
1431 return 0.; // can't handle empty bins
1432 } else {
1433 return (-1.5 * cvm::logn(A0) + 2. * cvm::logn(A1)
1434 - 0.5 * cvm::logn(A2)) * increment / widths[n];
1435 }
1436 }
1437 }
1438
1439
1443 inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
1444 int n = 0)
1445 {
1446 cvm::real A0, A1, A2;
1447 std::vector<int> ix = ix0;
1448
1449 // FIXME this can be rewritten more concisely with wrap_edge()
1450 if (periodic[n]) {
1451 ix[n]--; wrap(ix);
1452 A0 = value(ix);
1453 ix = ix0;
1454 ix[n]++; wrap(ix);
1455 A1 = value(ix);
1456 if (A0 * A1 == 0) {
1457 return 0.; // can't handle empty bins
1458 } else {
1459 return (A1 - A0) / (widths[n] * 2.);
1460 }
1461 } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
1462 ix[n]--;
1463 A0 = value(ix);
1464 ix = ix0;
1465 ix[n]++;
1466 A1 = value(ix);
1467 if (A0 * A1 == 0) {
1468 return 0.; // can't handle empty bins
1469 } else {
1470 return cvm::real(A1 - A0) / (widths[n] * 2.);
1471 }
1472 } else {
1473 // edge: use 2nd order derivative
1474 int increment = (ix[n] == 0 ? 1 : -1);
1475 // move right from left edge, or the other way around
1476 A0 = value(ix);
1477 ix[n] += increment; A1 = value(ix);
1478 ix[n] += increment; A2 = value(ix);
1479 return (-1.5 * cvm::real(A0) + 2. * cvm::real(A1)
1480 - 0.5 * cvm::real(A2)) * increment / widths[n];
1481 }
1482 }
1483
1484
1487 virtual inline cvm::real value_output(std::vector<int> const &ix,
1488 size_t const &imult = 0) const override
1489 {
1490 int s;
1491 if (imult > 0) {
1492 cvm::error("Error: trying to access a component "
1493 "larger than 1 in a scalar data grid.\n");
1494 return 0.;
1495 }
1496 if (samples) {
1497 return ( (s = samples->value(ix)) > 0) ?
1498 (data[address(ix) + imult] / cvm::real(s)) :
1499 0.0;
1500 } else {
1501 return data[address(ix) + imult];
1502 }
1503 }
1504
1506 virtual void value_input(std::vector<int> const &ix,
1507 cvm::real const &new_value,
1508 size_t const &imult = 0,
1509 bool add = false) override
1510 {
1511 if (imult > 0) {
1512 cvm::error("Error: trying to access a component "
1513 "larger than 1 in a scalar data grid.\n");
1514 return;
1515 }
1516 if (add) {
1517 if (samples)
1518 data[address(ix)] += new_value * samples->new_value(ix);
1519 else
1520 data[address(ix)] += new_value;
1521 } else {
1522 if (samples)
1523 data[address(ix)] = new_value * samples->value(ix);
1524 else
1525 data[address(ix)] = new_value;
1526 }
1527 has_data = true;
1528 }
1529
1531 cvm::real maximum_value() const;
1532
1534 cvm::real minimum_value() const;
1535
1538
1540 cvm::real integral() const;
1541
1544 cvm::real entropy() const;
1545
1548 cvm::real grid_rmsd(colvar_grid_scalar const &other_grid) const;
1549};
1550
1551
1552
1554class colvar_grid_gradient : public colvar_grid<cvm::real>
1555{
1556public:
1557
1560 std::shared_ptr<colvar_grid_count> samples;
1561
1564
1567 {}
1568
1570 colvar_grid_gradient(std::vector<int> const &nx_i);
1571
1573 colvar_grid_gradient(std::vector<colvar *> &colvars);
1574
1576 colvar_grid_gradient(std::string &filename);
1577
1579 colvar_grid_gradient(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_count> samples_in);
1580
1583 int min_samples;
1584
1586 std::string get_state_params() const;
1587
1589 int parse_params(std::string const &conf,
1591
1593 std::istream & read_restart(std::istream &is);
1594
1597
1599 std::ostream & write_restart(std::ostream &os);
1600
1603
1605 std::istream &read_raw(std::istream &is);
1606
1609
1613 std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
1614
1618 cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
1619
1621 std::istream & read_multicol(std::istream &is, bool add = false);
1622
1624 int read_multicol(std::string const &filename,
1625 std::string description = "grid file",
1626 bool add = false);
1627
1629 std::ostream & write_multicol(std::ostream &os) const;
1630
1632 int write_multicol(std::string const &filename,
1633 std::string description = "grid file") const;
1634
1636 std::ostream & write_opendx(std::ostream &os) const;
1637
1639 int write_opendx(std::string const &filename,
1640 std::string description = "grid file") const;
1641
1643 inline void vector_value(std::vector<int> const &ix, std::vector<cvm::real> &v) const
1644 {
1645 cvm::real const * p = &value(ix);
1646 if (samples) {
1647 int count = samples->value(ix);
1648 if (count) {
1649 cvm::real invcount = 1.0 / count;
1650 for (size_t i = 0; i < mult; i++) {
1651 v[i] = invcount * p[i];
1652 }
1653 } else {
1654 for (size_t i = 0; i < mult; i++) {
1655 v[i] = 0.0;
1656 }
1657 }
1658 } else {
1659 for (size_t i = 0; i < mult; i++) {
1660 v[i] = p[i];
1661 }
1662 }
1663 }
1664
1665
1667 inline void acc_value(std::vector<int> const &ix, std::vector<colvarvalue> const &values) {
1668 for (size_t imult = 0; imult < mult; imult++) {
1669 data[address(ix) + imult] += values[imult].real_value;
1670 }
1671 if (samples)
1672 samples->incr_count(ix);
1673 }
1674
1677 inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
1678 for (size_t imult = 0; imult < mult; imult++) {
1679 data[address(ix) + imult] -= forces[imult];
1680 }
1681 if (samples)
1682 samples->incr_count(ix);
1683 }
1684
1687 virtual cvm::real value_output(std::vector<int> const &ix,
1688 size_t const &imult = 0) const override
1689 {
1690 int s;
1691 if (samples) {
1692 return ( (s = samples->value(ix)) > 0) ?
1693 (data[address(ix) + imult] / cvm::real(s)) :
1694 0.0;
1695 } else {
1696 return data[address(ix) + imult];
1697 }
1698 }
1699
1703 {
1704 cvm::real fact;
1705 if ( weight <= min_samples ) {
1706 fact = 0.0;
1707 } else if ( weight < full_samples ) {
1708 fact = (weight - min_samples) / (weight * cvm::real(full_samples - min_samples));
1709 } else {
1710 fact = 1.0 / weight;
1711 }
1712 return fact;
1713 }
1714
1715
1720 virtual inline cvm::real value_output_smoothed(std::vector<int> const &ix, bool smoothed = true)
1721 {
1722 cvm::real weight, fact;
1723
1724 if (samples) {
1725 weight = cvm::real(samples->value(ix));
1726 } else {
1727 weight = 1.;
1728 }
1729
1730 if (smoothed) {
1731 fact = smooth_inverse_weight(weight);
1732 } else {
1733 fact = weight > 0. ? 1. / weight : 0.;
1734 }
1735
1736 return fact * data[address(ix)];
1737 }
1738
1742 inline void vector_value_smoothed(std::vector<int> const &ix, cvm::real *grad, bool smoothed = true)
1743 {
1744 cvm::real weight, fact;
1745
1746 if (samples) {
1747 weight = cvm::real(samples->value(ix));
1748 } else {
1749 weight = 1.;
1750 }
1751
1752 if (smoothed) {
1753 fact = smooth_inverse_weight(weight);
1754 } else {
1755 fact = weight > 0. ? 1. / weight : 0.;
1756 }
1757
1758 cvm::real *p = &(data[address(ix)]);
1759
1760 // Appease Clang analyzer, which likes to assume that mult is zero
1761 #ifdef __clang_analyzer__
1762 assert(mult > 0);
1763 #endif
1764
1765 for (size_t imult = 0; imult < mult; imult++) {
1766 grad[imult] = fact * p[imult];
1767 }
1768 }
1769
1773 virtual void value_input(std::vector<int> const &ix,
1774 cvm::real const &new_value,
1775 size_t const &imult = 0,
1776 bool add = false) override
1777 {
1778 if (add) {
1779 if (samples)
1780 data[address(ix) + imult] += new_value * samples->new_value(ix);
1781 else
1782 data[address(ix) + imult] += new_value;
1783 } else {
1784 if (samples)
1785 data[address(ix) + imult] = new_value * samples->value(ix);
1786 else
1787 data[address(ix) + imult] = new_value;
1788 }
1789 has_data = true;
1790 }
1791
1792
1794 inline cvm::real average(bool smoothed = false)
1795 {
1796 if (nd != 1 || nx[0] == 0) {
1797 return 0.0;
1798 }
1799
1800 cvm::real sum = 0.0;
1801 for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
1802 sum += value_output_smoothed(ix, smoothed);
1803 }
1804
1805 return (sum / cvm::real(nx[0]));
1806 }
1807
1810 cvm::real grid_rmsd(colvar_grid_gradient const &other_grid) const;
1811
1814 void write_1D_integral(std::ostream &os);
1815
1816};
1817
1818
1819
1821
1823{
1824 public:
1825
1827
1828 virtual ~integrate_potential()
1829 {}
1830
1832 integrate_potential(std::vector<colvar *> &colvars, std::shared_ptr<colvar_grid_gradient> gradients);
1833
1836
1838 int integrate(const int itmax, const cvm::real & tol, cvm::real & err, bool verbose = true);
1839
1842 void update_div_neighbors(const std::vector<int> &ix);
1843
1846 void update_div_local(const std::vector<int> &ix);
1847
1850 void set_div();
1851
1854 inline void set_zero_minimum() {
1855 add_constant(-1.0 * minimum_value());
1856 }
1857
1860
1861
1862 protected:
1863
1864 // Reference to gradient grid
1865 std::shared_ptr<colvar_grid_gradient> gradients;
1866
1868 std::vector<cvm::real> divergence;
1869
1870// std::vector<cvm::real> inv_lap_diag; // Inverse of the diagonal of the Laplacian; for conditioning
1871
1875 void get_grad(cvm::real * g, std::vector<int> &ix);
1876
1878 void nr_linbcg_sym(const std::vector<cvm::real> &b, std::vector<cvm::real> &x,
1879 const cvm::real &tol, const int itmax, int &iter, cvm::real &err);
1880
1882 cvm::real l2norm(const std::vector<cvm::real> &x);
1883
1885 void atimes(const std::vector<cvm::real> &x, std::vector<cvm::real> &r);
1886
1887// /// Inversion of preconditioner matrix
1888// void asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x);
1889};
1890
1891#endif
1892
Colvar_grid derived class to hold counters in discrete n-dim colvar space.
Definition: colvargrid.h:959
colvar_grid_count()
Default constructor.
Definition: colvargrid.cpp:21
int local_sample_count(int radius)
Return the average number of samples in a given "radius" around current bin Really a hypercube of len...
Definition: colvargrid.h:1068
std::istream & read_raw(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid.cpp:69
std::string get_state_params() const
Write the current grid parameters to a string.
Definition: colvargrid.cpp:38
cvm::real gradient_finite_diff(const std::vector< int > &ix0, int n=0)
Return the gradient of discrete count from finite differences on the same grid for dimension n (colva...
Definition: colvargrid.h:1197
std::ostream & write_restart(std::ostream &os)
Write all grid parameters and data to a formatted stream.
Definition: colvargrid.cpp:59
std::ostream & write_opendx(std::ostream &os) const
Write the grid data without labels, as they are represented in memory.
Definition: colvargrid.cpp:113
std::ostream & write_raw(std::ostream &os, size_t const buf_size=3) const
Definition: colvargrid.cpp:79
int parse_params(std::string const &conf, colvarparse::Parse_Mode const parse_mode=colvarparse::parse_normal)
Read new grid parameters from a string.
Definition: colvargrid.cpp:43
virtual void value_input(std::vector< int > const &ix, size_t const &t, size_t const &imult=0, bool add=false)
Enter or add a value, but also handle parent grid.
Definition: colvargrid.h:1048
std::istream & read_restart(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid.cpp:49
std::istream & read_multicol(std::istream &is, bool add=false)
Read a grid written by write_multicol(), incrementin if data is true.
Definition: colvargrid.cpp:90
std::ostream & write_multicol(std::ostream &os) const
Write grid in a format which is both human-readable and gnuplot-friendly.
Definition: colvargrid.cpp:102
cvm::real log_gradient_finite_diff(const std::vector< int > &ix0, int n=0, int offset=0)
Return the log-gradient from finite differences on the same grid for dimension n (colvar_grid_count)
Definition: colvargrid.h:1146
size_t const & new_value(std::vector< int > const &ix)
Get the binned count indexed by ix from the newly read data.
Definition: colvargrid.h:985
virtual ~colvar_grid_count()
Destructor.
Definition: colvargrid.h:966
void incr_count(std::vector< int > const &ix)
Increment the counter at given position.
Definition: colvargrid.h:979
Class for accumulating the gradient of a scalar function on a grid.
Definition: colvargrid.h:1555
std::string get_state_params() const
Write the current grid parameters to a string.
Definition: colvargrid.cpp:418
cvm::real smooth_inverse_weight(cvm::real weight)
Definition: colvargrid.h:1702
void vector_value_smoothed(std::vector< int > const &ix, cvm::real *grad, bool smoothed=true)
Obtain the vector value of the function at ix divided by its number of samples (if the count grid is ...
Definition: colvargrid.h:1742
virtual cvm::real value_output_smoothed(std::vector< int > const &ix, bool smoothed=true)
Return the scalar value of the function at ix divided by its number of samples (if the count grid is ...
Definition: colvargrid.h:1720
void acc_value(std::vector< int > const &ix, std::vector< colvarvalue > const &values)
Accumulate the value.
Definition: colvargrid.h:1667
cvm::real average(bool smoothed=false)
Compute and return average value for a 1D gradient grid.
Definition: colvargrid.h:1794
std::shared_ptr< colvar_grid_count > samples
Provide the sample count by which each binned value should be divided.
Definition: colvargrid.h:1560
std::istream & read_multicol(std::istream &is, bool add=false)
Read a grid written by write_multicol(), incrementin if data is true.
Definition: colvargrid.cpp:470
colvar_grid_gradient()
Default constructor.
Definition: colvargrid.cpp:332
std::ostream & write_restart(std::ostream &os)
Write all grid parameters and data to a formatted stream.
Definition: colvargrid.cpp:439
void write_1D_integral(std::ostream &os)
If the grid is 1-dimensional, integrate it and write the integral to a file (DEPRECATED by the integr...
Definition: colvargrid.cpp:505
int parse_params(std::string const &conf, colvarparse::Parse_Mode const parse_mode=colvarparse::parse_normal)
Read new grid parameters from a string.
Definition: colvargrid.cpp:423
virtual cvm::real value_output(std::vector< int > const &ix, size_t const &imult=0) const override
Return the value of the function at ix divided by its number of samples (if the count grid is defined...
Definition: colvargrid.h:1687
std::istream & read_raw(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid.cpp:449
virtual ~colvar_grid_gradient()
Destructor.
Definition: colvargrid.h:1566
virtual void value_input(std::vector< int > const &ix, cvm::real const &new_value, size_t const &imult=0, bool add=false) override
Get the value from a formatted output and transform it into the internal representation (it may have ...
Definition: colvargrid.h:1773
void acc_force(std::vector< int > const &ix, cvm::real const *forces)
Accumulate the gradient based on the force (i.e. sums the opposite of the force)
Definition: colvargrid.h:1677
int full_samples
Parameters for smoothing data with low sampling.
Definition: colvargrid.h:1582
std::ostream & write_raw(std::ostream &os, size_t const buf_size=3) const
Definition: colvargrid.cpp:459
std::ostream & write_multicol(std::ostream &os) const
Write grid in a format which is both human-readable and gnuplot-friendly.
Definition: colvargrid.cpp:482
std::istream & read_restart(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid.cpp:429
cvm::real grid_rmsd(colvar_grid_gradient const &other_grid) const
Return the RMSD between this grid's data and another one Grids must have the same dimensions.
Definition: colvargrid.cpp:562
void vector_value(std::vector< int > const &ix, std::vector< cvm::real > &v) const
Get a vector with the binned value(s) indexed by ix, normalized if applicable.
Definition: colvargrid.h:1643
std::ostream & write_opendx(std::ostream &os) const
Write the grid data without labels, as they are represented in memory.
Definition: colvargrid.cpp:493
Class for accumulating a scalar function on a grid.
Definition: colvargrid.h:1242
void acc_value(std::vector< int > const &ix, cvm::real const &new_value, size_t const &imult=0)
Accumulate the value.
Definition: colvargrid.h:1266
cvm::real minimum_value() const
Return the lowest value.
Definition: colvargrid.cpp:246
cvm::real gradient_finite_diff(const std::vector< int > &ix0, int n=0)
Return the gradient of discrete count from finite differences on the same grid for dimension n (colva...
Definition: colvargrid.h:1443
std::istream & read_restart(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid.cpp:160
std::string get_state_params() const
Write the current grid parameters to a string.
Definition: colvargrid.cpp:149
virtual ~colvar_grid_scalar()
Destructor.
Definition: colvargrid.cpp:145
cvm::real minimum_pos_value() const
Return the lowest positive value.
Definition: colvargrid.cpp:255
virtual cvm::real value_output(std::vector< int > const &ix, size_t const &imult=0) const override
Return the value of the function at ix divided by its number of samples (if the count grid is defined...
Definition: colvargrid.h:1487
cvm::real grid_rmsd(colvar_grid_scalar const &other_grid) const
Return the RMSD between this grid's data and another one Grids must have the same dimensions.
Definition: colvargrid.cpp:302
std::ostream & write_opendx(std::ostream &os) const
Write the grid data without labels, as they are represented in memory.
Definition: colvargrid.cpp:224
std::ostream & write_restart(std::ostream &os)
Write all grid parameters and data to a formatted stream.
Definition: colvargrid.cpp:170
std::ostream & write_multicol(std::ostream &os) const
Write grid in a format which is both human-readable and gnuplot-friendly.
Definition: colvargrid.cpp:213
std::istream & read_multicol(std::istream &is, bool add=false)
Read a grid written by write_multicol(), incrementin if data is true.
Definition: colvargrid.cpp:201
cvm::real integral() const
Calculates the integral of the map (uses widths if they are defined)
Definition: colvargrid.cpp:271
int parse_params(std::string const &conf, colvarparse::Parse_Mode const parse_mode=colvarparse::parse_normal)
Read new grid parameters from a string.
Definition: colvargrid.cpp:154
std::ostream & write_raw(std::ostream &os, size_t const buf_size=3) const
Definition: colvargrid.cpp:190
virtual void value_input(std::vector< int > const &ix, cvm::real const &new_value, size_t const &imult=0, bool add=false) override
Enter or add value but also deal with count grid.
Definition: colvargrid.h:1506
cvm::real entropy() const
Assuming that the map is a normalized probability density, calculates the entropy (uses widths if the...
Definition: colvargrid.cpp:285
colvar_grid_scalar()
Default constructor.
Definition: colvargrid.cpp:126
cvm::real maximum_value() const
Return the highest value.
Definition: colvargrid.cpp:236
cvm::real log_gradient_finite_diff(const std::vector< int > &ix0, int n=0, int offset=0)
Return the log-gradient from finite differences on the same grid for dimension n (colvar_grid_scalar)
Definition: colvargrid.h:1392
colvar_grid_count * samples
Provide the associated sample count by which each binned value should be divided.
Definition: colvargrid.h:1247
void vector_gradient_finite_diff(const std::vector< int > &ix0, std::vector< cvm::real > &grad)
Return the gradient of the scalar field from finite differences Input coordinates are those of gradie...
Definition: colvargrid.h:1339
std::istream & read_raw(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid.cpp:180
Grid of values of a function of several collective variables.
Definition: colvargrid.h:26
std::vector< colvarvalue > upper_boundaries
Upper boundaries of the colvars in this grid.
Definition: colvargrid.h:81
int setup()
Allocate data (allow initialization also after construction)
Definition: colvargrid.h:193
std::istream & read_restart(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid_def.h:44
void set_value(size_t i, T const &t)
Set the value at the point with linear address i (for speed)
Definition: colvargrid.h:526
void check_consistency()
Check that the grid information inside (boundaries, widths, ...) is consistent with the current setti...
Definition: colvargrid.h:863
std::string get_state_params() const
Write the current grid parameters to a string.
Definition: colvargrid_def.h:146
std::ostream & write_raw(std::ostream &os, size_t const buf_size=3) const
Definition: colvargrid_def.h:113
int current_bin_scalar_bound(int const i) const
Report the bin corresponding to the current value of variable i and assign first or last bin if out o...
Definition: colvargrid.h:446
std::vector< int > const new_index() const
Get the index corresponding to the "first" bin, to be used as the initial value for an index in loopi...
Definition: colvargrid.h:811
void add_constant(T const &t)
Add a constant to all elements (fast loop)
Definition: colvargrid.h:628
std::vector< size_t > new_data
Newly read data (used for count grids, when adding several grids read from disk)
Definition: colvargrid.h:51
cvm::real value_to_bin_scalar_fraction(colvarvalue const &value, const int i) const
Use the lower boundary and the width to report the fraction of bin beyond value_to_bin_scalar() that ...
Definition: colvargrid.h:474
void incr(std::vector< int > &ix) const
Increment the index, in a way that will make it loop over the whole nd-dimensional array.
Definition: colvargrid.h:829
size_t number_of_points(int const icv=-1) const
Definition: colvargrid.h:115
std::vector< bool > hard_upper_boundaries
Whether some colvars have hard upper boundaries.
Definition: colvargrid.h:90
virtual void value_input(std::vector< int > const &ix, T const &t, size_t const &imult=0, bool add=false)
Get the value from a formatted output and transform it into the internal representation (the two may ...
Definition: colvargrid.h:790
T const & value(size_t i) const
Get the binned value indexed by linear address i.
Definition: colvargrid.h:622
std::vector< int > const get_colvars_index_bound() const
Get the bin indices corresponding to the provided values of the colvars and assign first or last bin ...
Definition: colvargrid.h:674
int current_bin_flat_bound() const
Report the flattened bin address corresponding to the current value of all variables and assign first...
Definition: colvargrid.h:435
size_t mult
Multiplicity of each datum (allow the binning of non-scalar types such as atomic gradients)
Definition: colvargrid.h:42
size_t multiplicity() const
Return the multiplicity of the type used.
Definition: colvargrid.h:137
std::vector< int > const & sizes() const
Get the sizes in each direction.
Definition: colvargrid.h:125
void set_value(std::vector< int > const &ix, T const &t, size_t const &imult=0)
Set the value at the point with index ix.
Definition: colvargrid.h:517
colvarvalue bin_to_value_scalar(int const &i_bin, int const i) const
Use the two boundaries and the width to report the central value corresponding to a bin index.
Definition: colvargrid.h:503
std::vector< T > data
Low-level array of values.
Definition: colvargrid.h:48
void multiply_constant(cvm::real const &a)
Multiply all elements by a scalar constant (fast loop)
Definition: colvargrid.h:636
std::vector< int > const & number_of_points_vec() const
Return the numbers of points in all dimensions.
Definition: colvargrid.h:108
bool index_ok(std::vector< int > const &ix) const
Check that the index is within range in each of the dimensions.
Definition: colvargrid.h:818
cvm::real current_bin_scalar_fraction(int const i) const
Report the fraction of bin beyond current_bin_scalar()
Definition: colvargrid.h:467
bool has_data
Whether this grid has been filled with data or is still empty.
Definition: colvargrid.h:99
void add_grid(colvar_grid< T > const &other_grid, cvm::real scale_factor=1.0)
Add data from another grid of the same type, AND identical definition (boundaries,...
Definition: colvargrid.h:759
std::vector< int > const get_colvars_index(std::vector< colvarvalue > const &values) const
Get the bin indices corresponding to the provided values of the colvars.
Definition: colvargrid.h:652
void check_consistency(colvar_grid< T > const &other_grid)
Check that the grid information inside (boundaries, widths, ...) is consistent with that of another g...
Definition: colvargrid.h:882
colvar_grid(std::vector< int > const &nx_i, T const &t=T(), size_t mult_i=1)
Constructor from explicit grid sizes.
Definition: colvargrid.h:242
int parse_params(std::string const &conf, colvarparse::Parse_Mode const parse_mode=colvarparse::parse_normal)
Read new grid parameters from a string.
Definition: colvargrid_def.h:176
size_t nt
Total number of grid points.
Definition: colvargrid.h:45
T get_value(size_t i) const
Get the value at the point with linear address i (for speed)
Definition: colvargrid.h:532
void reset(T const &t=T())
Reset data (in case the grid is being reused)
Definition: colvargrid.h:199
size_t num_variables() const
Return the number of colvar objects.
Definition: colvargrid.h:102
int setup(std::vector< int > const &nx_i, T const &t=T(), size_t const &mult_i=1)
Allocate data.
Definition: colvargrid.h:152
colvar_grid(std::vector< colvar * > const &colvars, T const &t=T(), size_t mult_i=1, bool add_extra_bin=false)
Constructor from a vector of colvars.
Definition: colvargrid.h:253
colvar_grid(colvar_grid< T > const &g)
"Almost copy-constructor": only copies configuration parameters from another grid,...
Definition: colvargrid.h:221
std::vector< cvm::real > widths
Widths of the colvars in this grid.
Definition: colvargrid.h:93
bool has_parent_data
True if this is a count grid related to another grid of data.
Definition: colvargrid.h:96
std::ostream & write_opendx(std::ostream &os) const
Write the grid data without labels, as they are represented in memory.
Definition: colvargrid_def.h:435
bool wrap_to_edge(std::vector< int > &ix, std::vector< int > &edge_bin) const
Definition: colvargrid.h:407
std::vector< int > const get_colvars_index() const
Get the bin indices corresponding to the current values of the colvars.
Definition: colvargrid.h:663
int current_bin_scalar(int const i, int const iv) const
Report the bin corresponding to the current value of item iv in variable i.
Definition: colvargrid.h:452
void delta_grid(colvar_grid< T > const &other_grid)
Get the change from this to other_grid and store the result in this. this_grid := other_grid - this_g...
Definition: colvargrid.h:542
void request_actual_value(bool b=true)
Request grid to use actual values of extended coords.
Definition: colvargrid.h:143
cvm::real bin_distance_from_boundaries(std::vector< colvarvalue > const &values, bool skip_hard_boundaries=false)
Get the minimal distance (in number of bins) from the boundaries; a negative number is returned if th...
Definition: colvargrid.h:686
std::vector< int > nxc
Cumulative number of points along each dimension.
Definition: colvargrid.h:38
int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
Use the lower boundary and the width to report which bin the provided value is in and assign first or...
Definition: colvargrid.h:482
int current_bin_scalar(int const i) const
Report the bin corresponding to the current value of variable i.
Definition: colvargrid.h:428
std::istream & read_raw(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid_def.h:100
std::vector< bool > hard_lower_boundaries
Whether some colvars have hard lower boundaries.
Definition: colvargrid.h:87
T const & value(std::vector< int > const &ix, size_t const &imult=0) const
Get the binned value indexed by ix, or the first of them if the multiplicity is larger than 1.
Definition: colvargrid.h:615
std::vector< bool > use_actual_value
Do we request actual value (for extended-system colvars)?
Definition: colvargrid.h:57
size_t raw_data_num() const
Size of the data as they are represented in memory.
Definition: colvargrid.h:610
size_t address(std::vector< int > const &ix) const
Get the low-level index corresponding to an index.
Definition: colvargrid.h:60
void remove_small_values(cvm::real const &a)
Assign values that are smaller than scalar constant the latter value (fast loop)
Definition: colvargrid.h:643
std::vector< int > nx
Number of points along each dimension.
Definition: colvargrid.h:35
std::ostream & write_multicol(std::ostream &os) const
Write grid in a format which is both human-readable and gnuplot-friendly.
Definition: colvargrid_def.h:371
virtual ~colvar_grid()
Destructor.
Definition: colvargrid.h:215
int value_to_bin_scalar(colvarvalue const &value, const int i) const
Use the lower boundary and the width to report which bin the provided value is in.
Definition: colvargrid.h:461
size_t nd
Number of dimensions.
Definition: colvargrid.h:32
virtual T value_output(std::vector< int > const &ix, size_t const &imult=0) const
Return the value suitable for output purposes (so that it may be rescaled or manipulated without chan...
Definition: colvargrid.h:781
void set_sizes(std::vector< int > const &new_sizes)
Set the sizes in each direction.
Definition: colvargrid.h:131
void copy_grid(colvar_grid< T > const &other_grid)
Copy data from another grid of the same type, AND identical definition (boundaries,...
Definition: colvargrid.h:567
bool wrap_detect_edge(std::vector< int > &ix) const
Definition: colvargrid.h:392
void raw_data_out(T *out_data) const
Extract the grid data as they are represented in memory. Put the results in "out_data".
Definition: colvargrid.h:590
void map_grid(colvar_grid< T > const &other_grid)
Add data from another grid of the same type.
Definition: colvargrid.h:716
std::istream & read_multicol(std::istream &is, bool add=false)
Read a grid written by write_multicol(), incrementing if add is true.
Definition: colvargrid_def.h:252
std::vector< bool > periodic
Whether some colvars are periodic.
Definition: colvargrid.h:84
int value_to_bin_scalar(colvarvalue const &value, colvarvalue const &new_offset, cvm::real const &new_width) const
Same as the standard version, but uses another grid definition.
Definition: colvargrid.h:494
std::vector< colvarvalue > lower_boundaries
Lower boundaries of the colvars in this grid.
Definition: colvargrid.h:78
void wrap(std::vector< int > &ix) const
Definition: colvargrid.h:375
colvarvalue bin_to_value_scalar(int const &i_bin, colvarvalue const &new_offset, cvm::real const &new_width) const
Same as the standard version, but uses different parameters.
Definition: colvargrid.h:509
colvar_grid()
Default constructor.
Definition: colvargrid.h:206
std::vector< colvar * > cv
Colvars collected in this grid.
Definition: colvargrid.h:54
std::ostream & write_restart(std::ostream &os)
Write all grid parameters and data to a formatted stream.
Definition: colvargrid_def.h:56
void raw_data_in(const T *in_data)
Input the data as they are represented in memory.
Definition: colvargrid.h:599
@ f_cv_hard_lower_boundary
The lower boundary is not defined from user's choice.
Definition: colvardeps.h:314
@ f_cv_hard_upper_boundary
The upper boundary is not defined from user's choice.
Definition: colvardeps.h:316
static real logn(real const &x)
Definition: colvarmodule.h:177
double real
Defining an abstract real number allows to switch precision.
Definition: colvarmodule.h:95
static real floor(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:121
static void log(std::string const &message, int min_log_level=10)
Definition: colvarmodule.cpp:1955
static int error(std::string const &message, int code=-1)
Print a message to the main log and set global error code.
Definition: colvarmodule.cpp:2032
static size_t const cv_prec
Number of digits to represent a collective variables value(s)
Definition: colvarmodule.h:664
static real sqrt(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:133
static size_t const cv_width
Number of characters to represent a collective variables value(s)
Definition: colvarmodule.h:666
static bool debug()
Whether debug output should be enabled (compile-time option)
Definition: colvarmodule.h:330
static real fabs(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:127
static std::string to_str(char const *s)
Convert to string for output purposes.
Definition: colvarmodule.cpp:2378
Base class containing parsing functions; all objects which need to parse input inherit from this.
Definition: colvarparse.h:27
Parse_Mode
How a keyword is parsed in a string.
Definition: colvarparse.h:53
@ parse_normal
Alias for old default behavior (should be phased out)
Definition: colvarparse.h:72
Value of a collective variable: this is a metatype which can be set at runtime. By default it is set ...
Definition: colvarvalue.h:43
@ type_scalar
Scalar number, implemented as colvarmodule::real (default)
Definition: colvarvalue.h:56
cvm::real real_value
Real data member.
Definition: colvarvalue.h:77
Definition: colvars_memstream.h:30
Integrate (1D, 2D or 3D) gradients.
Definition: colvargrid.h:1823
void set_zero_minimum()
Add constant to potential so that its minimum value is zero Useful e.g. for output.
Definition: colvargrid.h:1854
void update_div_local(const std::vector< int > &ix)
Update matrix containing divergence and boundary conditions called by update_div_neighbors and by col...
Definition: colvargrid.cpp:747
std::vector< cvm::real > divergence
Array holding divergence + boundary terms (modified Neumann) if not periodic.
Definition: colvargrid.h:1868
void set_div()
Set matrix containing divergence and boundary conditions based on complete gradient grid.
Definition: colvargrid.cpp:684
void update_div_neighbors(const std::vector< int > &ix)
Update matrix containing divergence and boundary conditions based on new gradient point value,...
Definition: colvargrid.cpp:693
cvm::real l2norm(const std::vector< cvm::real > &x)
l2 norm of a vector
Definition: colvargrid.cpp:1268
bool b_smoothed
Flag requesting the use of a smoothed version of the gradient (default: false)
Definition: colvargrid.h:1859
int integrate(const int itmax, const cvm::real &tol, cvm::real &err, bool verbose=true)
Calculate potential from divergence (in 2D); return number of steps.
Definition: colvargrid.cpp:644
void get_grad(cvm::real *g, std::vector< int > &ix)
Definition: colvargrid.cpp:731
void nr_linbcg_sym(const std::vector< cvm::real > &b, std::vector< cvm::real > &x, const cvm::real &tol, const int itmax, int &iter, cvm::real &err)
Solve linear system based on CG, valid for symmetric matrices only.
Definition: colvargrid.cpp:1215
void atimes(const std::vector< cvm::real > &x, std::vector< cvm::real > &r)
Multiplication by sparse matrix representing Lagrangian (or its transpose)
Definition: colvargrid.cpp:799
Collective variables main module.
Parsing functions for collective variables.