Collective Variables Module - Developer Documentation
Loading...
Searching...
No Matches
colvargrid_def.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
11
12#ifndef COLVARGRID_DEF_H
13#define COLVARGRID_DEF_H
14
15#include <iostream>
16#include <iomanip>
17
18#include "colvarmodule.h"
19#include "colvarproxy.h"
20#include "colvar.h"
21#include "colvargrid.h"
22#include "colvars_memstream.h"
23
24
25template <class T>
26colvar_grid<T>::colvar_grid(std::string const &filename, size_t mult_i)
27{
28std::istream &is = cvm::main()->proxy->input_stream(filename, "multicol grid file");
29if (!is) {
30 return;
31}
32
33// Data in the header: nColvars, then for each
34// xiMin, dXi, nPoints, periodic flag
35
36std::string hash;
37size_t i;
38
39if ( !(is >> hash) || (hash != "#") ) {
40 cvm::error("Error reading grid at position "+
41 cvm::to_str(static_cast<size_t>(is.tellg()))+
42 " in stream(read \"" + hash + "\")\n");
43 return;
44}
45
46is >> nd;
47mult = (mult_i == 0) ? nd : mult_i;
48
49std::vector<cvm::real> lower_in(nd), widths_in(nd);
50std::vector<int> nx_in(nd);
51std::vector<int> periodic_in(nd);
52
53for (i = 0; i < nd; i++ ) {
54 if ( !(is >> hash) || (hash != "#") ) {
55 cvm::error("Error reading grid at position "+
56 cvm::to_str(static_cast<size_t>(is.tellg()))+
57 " in stream(read \"" + hash + "\")\n");
58 return;
59 }
60
61 is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i];
62}
63
64this->setup(nx_in, 0., mult);
65
66widths = widths_in;
67
68for (i = 0; i < nd; i++ ) {
69 lower_boundaries.push_back(colvarvalue(lower_in[i]));
70 periodic.push_back(static_cast<bool>(periodic_in[i]));
71}
72
73// Reset the istream for read_multicol, which expects the whole file
74is.clear();
75is.seekg(0);
76read_multicol(is);
78}
79
80
81template <class T, class IST> IST &read_restart_template_(colvar_grid<T> &g, IST &is)
82{
83 auto const start_pos = is.tellg();
84 std::string conf;
85 if ((is >> colvarparse::read_block("grid_parameters", &conf)) &&
86 (g.parse_params(conf, colvarparse::parse_restart) == COLVARS_OK) && g.read_raw(is)) {
87 return is;
88 }
89 auto const error_pos = is.tellg();
90 is.clear();
91 is.seekg(start_pos);
92 is.setstate(std::ios::failbit);
93 cvm::error("Error: in reading grid state from stream at position " + cvm::to_str(error_pos) +
94 "\n",
95 COLVARS_INPUT_ERROR);
96 return is;
97}
98
99
100template <class T> std::istream &colvar_grid<T>::read_restart(std::istream &is)
101{
102 return read_restart_template_<T, std::istream>(*this, is);
103}
104
105
107{
108 return read_restart_template_<T, cvm::memory_stream>(*this, is);
109}
110
111
112template <class T> std::ostream &colvar_grid<T>::write_restart(std::ostream &os)
113{
114 os << "grid_parameters {\n" << get_state_params() << "}\n";
115 write_raw(os);
116 return os;
117}
118
119
121{
122 os << std::string("grid_parameters") << get_state_params();
123 write_raw(os);
124 return os;
125}
126
127
128template <class T, class IST> IST &read_raw_template_(colvar_grid<T> &g, IST &is)
129{
130 auto const start_pos = is.tellg();
131
132 for (std::vector<int> ix = g.new_index(); g.index_ok(ix); g.incr(ix)) {
133 for (size_t imult = 0; imult < g.mult; imult++) {
134 T new_value;
135 if (is >> new_value) {
136 g.value_input(ix, new_value, imult);
137 } else {
138 is.clear();
139 is.seekg(start_pos);
140 is.setstate(std::ios::failbit);
142 "Error: failed to read all of the grid points from file. Possible explanations: grid "
143 "parameters in the configuration (lowerBoundary, upperBoundary, width) are different "
144 "from those in the file, or the file is corrupt/incomplete.\n",
145 COLVARS_INPUT_ERROR);
146 return is;
147 }
148 }
149 }
150
151 g.has_data = true;
152 return is;
153}
154
155
156template <class T> std::istream &colvar_grid<T>::read_raw(std::istream &is)
157{
158 return read_raw_template_<T, std::istream>(*this, is);
159}
160
161
163{
164 return read_raw_template_<T, cvm::memory_stream>(*this, is);
165}
166
167
168template <class T>
169std::ostream &colvar_grid<T>::write_raw(std::ostream &os, size_t const buf_size) const
170{
171 auto const w = os.width();
172 auto const p = os.precision();
173
174 size_t count = 0;
175 for (auto ix = new_index(); index_ok(ix); incr(ix)) {
176 for (size_t imult = 0; imult < mult; imult++) {
177 os << " " << std::setw(w) << std::setprecision(p) << value_output(ix, imult);
178 if (((++count) % buf_size) == 0)
179 os << "\n";
180 }
181 }
182 // write a final newline only if buffer is not empty
183 if ((count % buf_size) != 0)
184 os << "\n";
185
186 return os;
187}
188
189
190template <class T>
192{
193 for (auto ix = new_index(); index_ok(ix); incr(ix)) {
194 for (size_t imult = 0; imult < mult; imult++) {
195 os << value_output(ix, imult);
196 }
197 }
198 return os;
199}
200
201
202template <class T> std::string colvar_grid<T>::get_state_params() const
203{
204 std::ostringstream os;
205 size_t i;
206 os << " n_colvars " << nd << "\n";
207
208 os << " lower_boundaries ";
209 for (i = 0; i < nd; i++)
210 os << " " << lower_boundaries[i];
211 os << "\n";
212
213 os << " upper_boundaries ";
214 for (i = 0; i < nd; i++)
215 os << " " << upper_boundaries[i];
216 os << "\n";
217
218 os << " widths ";
219 for (i = 0; i < nd; i++)
220 os << " " << widths[i];
221 os << "\n";
222
223 os << " sizes ";
224 for (i = 0; i < nd; i++)
225 os << " " << nx[i];
226 os << "\n";
227
228 return os.str();
229}
230
231
232template <class T> int colvar_grid<T>::parse_params(std::string const &conf,
233 colvarparse::Parse_Mode const parse_mode)
234{
235 if (cvm::debug())
236 cvm::log("Reading grid configuration from string.\n");
237
238 std::vector<int> old_nx = nx;
239 std::vector<colvarvalue> old_lb = lower_boundaries;
240 std::vector<colvarvalue> old_ub = upper_boundaries;
241 std::vector<cvm::real> old_w = widths;
242
243 {
244 size_t nd_in = 0;
245 // this is only used in state files
246 colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
247 if (nd_in != nd) {
248 cvm::error("Error: trying to read data for a grid "
249 "that contains a different number of colvars ("+
250 cvm::to_str(nd_in)+") than the grid defined "
251 "in the configuration file("+cvm::to_str(nd)+
252 ").\n");
253 return COLVARS_ERROR;
254 }
255 }
256
257 // underscore keywords are used in state file
258 colvarparse::get_keyval(conf, "lower_boundaries",
259 lower_boundaries, lower_boundaries, colvarparse::parse_silent);
260 colvarparse::get_keyval(conf, "upper_boundaries",
261 upper_boundaries, upper_boundaries, colvarparse::parse_silent);
262 // plural form is used in state file
263 colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent);
264
265 // camel case keywords are used in config file
266 colvarparse::get_keyval(conf, "lowerBoundary",
267 lower_boundaries, lower_boundaries, parse_mode);
268 colvarparse::get_keyval(conf, "upperBoundary",
269 upper_boundaries, upper_boundaries, parse_mode);
270
271 colvarparse::get_keyval(conf, "width", widths, widths, parse_mode);
272
273 // only used in state file
274 colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);
275
276 if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
277
278 if (! use_actual_value.size()) use_actual_value.assign(nd, false);
279 if (! periodic.size()) periodic.assign(nd, false);
280 if (! widths.size()) widths.assign(nd, 1.0);
281
282 cvm::real eps = 1.e-10;
283
284 bool new_params = false;
285 if (old_nx.size()) {
286 for (size_t i = 0; i < nd; i++) {
287 if (old_nx[i] != nx[i] ||
288 cvm::sqrt(cv[i]->dist2(old_lb[i], lower_boundaries[i])) > eps ||
289 cvm::sqrt(cv[i]->dist2(old_ub[i], upper_boundaries[i])) > eps ||
290 cvm::fabs(old_w[i] - widths[i]) > eps) {
291 new_params = true;
292 }
293 }
294 } else {
295 new_params = true;
296 }
297
298 // reallocate the array in case the grid params have just changed
299 if (new_params) {
300 init_from_boundaries();
301 // data.clear(); // no longer needed: setup calls clear()
302 return this->setup(nx, T(), mult);
303 }
304
305 return COLVARS_OK;
306}
307
308
309template <class T>
310std::istream & colvar_grid<T>::read_multicol(std::istream &is, bool add)
311{
312 // Data in the header: nColvars, then for each
313 // xiMin, dXi, nPoints, periodic flag
314
315 std::string hash;
316 cvm::real lower, width, x;
317 size_t n, periodic_flag;
318 bool remap;
319 std::vector<T> new_value;
320 std::vector<int> nx_read;
321 std::vector<int> bin;
322
323 if ( cv.size() > 0 && cv.size() != nd ) {
324 cvm::error("Cannot read grid file: number of variables in file differs from number referenced by grid.\n");
325 return is;
326 }
327
328 if ( !(is >> hash) || (hash != "#") ) {
329 cvm::error("Error reading grid at position "+
330 cvm::to_str(static_cast<size_t>(is.tellg()))+
331 " in stream(read \"" + hash + "\")\n", COLVARS_INPUT_ERROR);
332 return is;
333 }
334
335 is >> n;
336 if ( n != nd ) {
337 cvm::error("Error reading grid: wrong number of collective variables.\n");
338 return is;
339 }
340
341 nx_read.resize(n);
342 bin.resize(n);
343 new_value.resize(mult);
344
345 if (this->has_parent_data && add) {
346 new_data.resize(data.size());
347 }
348
349 remap = false;
350 for (size_t i = 0; i < nd; i++ ) {
351 if ( !(is >> hash) || (hash != "#") ) {
352 cvm::error("Error reading grid at position "+
353 cvm::to_str(static_cast<size_t>(is.tellg()))+
354 " in stream(read \"" + hash + "\")\n");
355 return is;
356 }
357
358 is >> lower >> width >> nx_read[i] >> periodic_flag;
359
360
361 if ( (cvm::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
362 (cvm::fabs(width - widths[i] ) > 1.0e-10) ||
363 (nx_read[i] != nx[i]) ) {
364 cvm::log("Warning: reading from different grid definition (colvar "
365 + cvm::to_str(i+1) + "); remapping data on new grid.\n");
366 remap = true;
367 }
368 }
369
370 if ( remap ) {
371 // re-grid data
372 while (is.good()) {
373 bool end_of_file = false;
374
375 for (size_t i = 0; i < nd; i++ ) {
376 if ( !(is >> x) ) end_of_file = true;
377 bin[i] = value_to_bin_scalar(x, i);
378 // if x is out of bounds and we are using PBC, wrap it
379 // Ignore out of bounds points in non-PBC
380 wrap_detect_edge(bin);
381 }
382 if (end_of_file) break;
383
384 for (size_t imult = 0; imult < mult; imult++) {
385 is >> new_value[imult];
386 }
387
388 if ( index_ok(bin) ) {
389 for (size_t imult = 0; imult < mult; imult++) {
390 value_input(bin, new_value[imult], imult, add);
391 }
392 }
393 }
394 } else {
395 // do not re-grid the data but assume the same grid is used
396 for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix) ) {
397 for (size_t i = 0; i < nd; i++ ) {
398 is >> x;
399 }
400 for (size_t imult = 0; imult < mult; imult++) {
401 is >> new_value[imult];
402 value_input(ix, new_value[imult], imult, add);
403 }
404 }
405 }
406 has_data = true;
407 return is;
408}
409
410
411template <class T>
412int colvar_grid<T>::read_multicol(std::string const &filename,
413 std::string description,
414 bool add)
415{
416 std::istream &is = cvm::main()->proxy->input_stream(filename, description);
417 if (!is) {
418 return COLVARS_FILE_ERROR;
419 }
420 if (colvar_grid<T>::read_multicol(is, add)) {
421 cvm::main()->proxy->close_input_stream(filename);
422 return COLVARS_OK;
423 }
424 return COLVARS_FILE_ERROR;
425}
426
427
428template <class T>
429std::ostream & colvar_grid<T>::write_multicol(std::ostream &os) const
430{
431 // Save the output formats
432 std::ios_base::fmtflags prev_flags(os.flags());
433
434 // Data in the header: nColvars, then for each
435 // xiMin, dXi, nPoints, periodic
436
437 os << std::setw(2) << "# " << nd << "\n";
438 // Write the floating numbers in full precision
439 os.setf(std::ios::scientific, std::ios::floatfield);
440 for (size_t i = 0; i < nd; i++) {
441 os << "# "
442 << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << lower_boundaries[i] << " "
443 << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << widths[i] << " "
444 << std::setw(10) << nx[i] << " "
445 << periodic[i] << "\n";
446 }
447
448 for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix) ) {
449
450 if (ix.back() == 0) {
451 // if the last index is 0, add a new line to mark the new record
452 os << "\n";
453 }
454
455 for (size_t i = 0; i < nd; i++) {
456 os << " "
457 << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec)
458 << bin_to_value_scalar(ix[i], i);
459 }
460 os << " ";
461 for (size_t imult = 0; imult < mult; imult++) {
462 os << " "
463 << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec)
464 << value_output(ix, imult);
465 }
466 os << "\n";
467 }
468
469 // Restore the output formats
470 os.flags(prev_flags);
471
472 return os;
473}
474
475
476template <class T>
477int colvar_grid<T>::write_multicol(std::string const &filename,
478 std::string description) const
479{
480 int error_code = COLVARS_OK;
481 std::ostream &os = cvm::main()->proxy->output_stream(filename, description);
482 if (!os) {
483 return COLVARS_FILE_ERROR;
484 }
485 error_code |= colvar_grid<T>::write_multicol(os) ? COLVARS_OK :
486 COLVARS_FILE_ERROR;
487 cvm::main()->proxy->close_output_stream(filename);
488 return error_code;
489}
490
491
492template <class T>
493std::ostream & colvar_grid<T>::write_opendx(std::ostream &os) const
494{
495 // write the header
496 os << "object 1 class gridpositions counts";
497 size_t icv;
498 for (icv = 0; icv < num_variables(); icv++) {
499 os << " " << number_of_points(icv);
500 }
501 os << "\n";
502
503 os << "origin";
504 for (icv = 0; icv < num_variables(); icv++) {
505 os << " " << (lower_boundaries[icv].real_value + 0.5 * widths[icv]);
506 }
507 os << "\n";
508
509 for (icv = 0; icv < num_variables(); icv++) {
510 os << "delta";
511 for (size_t icv2 = 0; icv2 < num_variables(); icv2++) {
512 if (icv == icv2) os << " " << widths[icv];
513 else os << " " << 0.0;
514 }
515 os << "\n";
516 }
517
518 os << "object 2 class gridconnections counts";
519 for (icv = 0; icv < num_variables(); icv++) {
520 os << " " << number_of_points(icv);
521 }
522 os << "\n";
523
524 os << "object 3 class array type double rank 0 items "
525 << number_of_points() << " data follows\n";
526
527 write_raw(os);
528
529 os << "object \"collective variables scalar field\" class field\n";
530 return os;
531}
532
533
534template <class T>
535int colvar_grid<T>::write_opendx(std::string const &filename,
536 std::string description) const
537{
538 int error_code = COLVARS_OK;
539 std::ostream &os = cvm::main()->proxy->output_stream(filename, description);
540 if (!os) {
541 return COLVARS_FILE_ERROR;
542 }
543 error_code |= colvar_grid<T>::write_opendx(os) ? COLVARS_OK :
544 COLVARS_FILE_ERROR;
545 cvm::main()->proxy->close_output_stream(filename);
546 return error_code;
547}
548
549#endif
Grid of values of a function of several collective variables.
Definition: colvargrid.h:52
std::istream & read_restart(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid_def.h:100
std::string get_state_params() const
Write the current grid parameters to a string.
Definition: colvargrid_def.h:202
std::ostream & write_raw(std::ostream &os, size_t const buf_size=3) const
Definition: colvargrid_def.h:169
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:854
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:872
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:833
size_t mult
Multiplicity of each datum (allow the binning of non-scalar types such as atomic gradients)
Definition: colvargrid.h:59
bool index_ok(std::vector< int > const &ix) const
Check that the index is within range in each of the dimensions.
Definition: colvargrid.h:861
bool has_data
Whether this grid has been filled with data or is still empty.
Definition: colvargrid.h:106
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:232
std::ostream & write_opendx(std::ostream &os) const
Write the grid data without labels, as they are represented in memory.
Definition: colvargrid_def.h:493
std::istream & read_raw(std::istream &is)
Read all grid parameters and data from a formatted stream.
Definition: colvargrid_def.h:156
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:429
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:310
colvar_grid()
Default constructor.
Definition: colvargrid.h:213
std::ostream & write_restart(std::ostream &os)
Write all grid parameters and data to a formatted stream.
Definition: colvargrid_def.h:112
double real
Defining an abstract real number allows to switch precision.
Definition: colvarmodule.h:100
static void log(std::string const &message, int min_log_level=10)
Definition: colvarmodule.cpp:1993
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:2070
static size_t const cv_prec
Number of digits to represent a collective variables value(s)
Definition: colvarmodule.h:699
static colvarmodule * main()
Access the one instance of the Colvars module.
Definition: colvarmodule.cpp:175
static real sqrt(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:138
static size_t const cv_width
Number of characters to represent a collective variables value(s)
Definition: colvarmodule.h:701
static bool debug()
Whether debug output should be enabled (compile-time option)
Definition: colvarmodule.h:365
static real fabs(real const &x)
Reimplemented to work around MS compiler issues.
Definition: colvarmodule.h:132
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:895
static std::string to_str(char const *s)
Convert to string for output purposes.
Definition: colvarmodule.cpp:2416
Definition: colvarparse.h:276
Parse_Mode
How a keyword is parsed in a string.
Definition: colvarparse.h:53
@ parse_restart
The call is being executed from a read_restart() function.
Definition: colvarparse.h:70
@ parse_silent
Do not print the keyword.
Definition: colvarparse.h:63
virtual int close_output_stream(std::string const &output_name)
Closes the given output file/channel.
Definition: colvarproxy_io.cpp:475
std::istream & input_stream(std::string const &input_name, std::string const description="file/channel", bool error_on_fail=true)
Definition: colvarproxy_io.cpp:265
virtual std::ostream & output_stream(std::string const &output_name, std::string const description)
Definition: colvarproxy_io.cpp:408
int close_input_stream(std::string const &input_name)
Closes the given input stream.
Definition: colvarproxy_io.cpp:339
Value of a collective variable: this is a metatype which can be set at runtime. By default it is set ...
Definition: colvarvalue.h:43
Definition: colvars_memstream.h:30
Collective variables main module.
Colvars proxy classes.