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