Homec4science
Diffusion lammps (lammm-master) Browse

lammps/lib/colvars/colvargrid.h
fee52f2fac2blammm-master

Last Change
Raw File

colvargrid.h
View Options

// -*- c++ -*-
#ifndef COLVARGRID_H
#define COLVARGRID_H
#include <iostream>
#include <iomanip>
#include <cmath>
#include "colvar.h"
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
/// \brief Grid of values of a function of several collective
/// variables \param T The data type
///
/// Only scalar colvars supported so far
template <class T> class colvar_grid : public colvarparse {
protected:
/// Number of dimensions
size_t nd;
/// Number of points along each dimension
std::vector<int> nx;
/// Cumulative number of points along each dimension
std::vector<int> nxc;
/// \brief Multiplicity of each datum (allow the binning of
/// non-scalar types)
size_t mult;
/// Total number of grid points
size_t nt;
/// Low-level array of values
std::vector<T> data;
/// Newly read data (used for count grids, when adding several grids read from disk)
std::vector<size_t> new_data;
/// Colvars collected in this grid
std::vector<colvar *> cv;
/// Do we request actual value (for extended-system colvars)?
std::vector<bool> actual_value;
/// Get the low-level index corresponding to an index
inline size_t address (std::vector<int> const &ix) const
{
size_t addr = 0;
for (size_t i = 0; i < nd; i++) {
addr += ix[i]*nxc[i];
if (cvm::debug()) {
if (ix[i] >= nx[i])
cvm::fatal_error ("Error: exceeding bounds in colvar_grid.\n");
}
}
return addr;
}
public:
/// Lower boundaries of the colvars in this grid
std::vector<colvarvalue> lower_boundaries;
/// Upper boundaries of the colvars in this grid
std::vector<colvarvalue> upper_boundaries;
/// Whether some colvars are periodic
std::vector<bool> periodic;
/// Whether some colvars have hard lower boundaries
std::vector<bool> hard_lower_boundaries;
/// Whether some colvars have hard upper boundaries
std::vector<bool> hard_upper_boundaries;
/// Widths of the colvars in this grid
std::vector<cvm::real> widths;
/// True if this is a count grid related to another grid of data
bool has_parent_data;
/// Whether this grid has been filled with data or is still empty
bool has_data;
/// Return the number of colvars
inline size_t number_of_colvars() const
{
return nd;
}
/// Return the number of points in the i-th direction, if provided, or
/// the total number
inline size_t number_of_points (int const icv = -1) const
{
if (icv < 0) {
return nt;
} else {
return nx[icv];
}
}
/// Get the sizes in each direction
inline std::vector<int> const & sizes() const
{
return nx;
}
/// Set the sizes in each direction
inline void set_sizes (std::vector<int> const &new_sizes)
{
nx = new_sizes;
}
/// Return the multiplicity of the type used
inline size_t multiplicity() const
{
return mult;
}
/// \brief Allocate data (allow initialization also after construction)
void create (std::vector<int> const &nx_i,
T const &t = T(),
size_t const &mult_i = 1)
{
mult = mult_i;
nd = nx_i.size();
nxc.resize (nd);
nx = nx_i;
nt = mult;
for (int i = nd-1; i >= 0; i--) {
if (nx_i[i] <= 0)
cvm::fatal_error ("Error: providing an invalid number of points, "+
cvm::to_str (nx_i[i])+".\n");
nxc[i] = nt;
nt *= nx[i];
}
data.reserve (nt);
data.assign (nt, t);
}
/// \brief Allocate data (allow initialization also after construction)
void create()
{
create (this->nx, T(), this->mult);
}
/// \brief Reset data (in case the grid is being reused)
void reset (T const &t = T())
{
data.assign (nt, t);
}
/// Default constructor
colvar_grid() : has_data (false)
{
save_delimiters = false;
nd = nt = 0;
}
/// Destructor
virtual ~colvar_grid()
{}
/// \brief "Almost copy-constructor": only copies configuration
/// parameters from another grid, but doesn't reallocate stuff;
/// create() must be called after that;
colvar_grid (colvar_grid<T> const &g) : has_data (false),
nd (g.nd),
mult (g.mult),
cv (g.cv),
lower_boundaries (g.lower_boundaries),
upper_boundaries (g.upper_boundaries),
hard_lower_boundaries (g.hard_lower_boundaries),
hard_upper_boundaries (g.hard_upper_boundaries),
periodic (g.periodic),
widths (g.widths),
actual_value (g.actual_value),
data()
{
save_delimiters = false;
}
/// \brief Constructor from explicit grid sizes \param nx_i Number
/// of grid points along each dimension \param t Initial value for
/// the function at each point (optional) \param mult_i Multiplicity
/// of each value
colvar_grid (std::vector<int> const &nx_i,
T const &t = T(),
size_t const &mult_i = 1) : has_data (false)
{
save_delimiters = false;
this->create (nx_i, t, mult_i);
}
/// \brief Constructor from a vector of colvars
colvar_grid (std::vector<colvar *> const &colvars,
T const &t = T(),
size_t const &mult_i = 1,
bool margin = false)
: cv (colvars), has_data (false)
{
save_delimiters = false;
std::vector<int> nx_i;
if (cvm::debug())
cvm::log ("Allocating a grid for "+cvm::to_str (colvars.size())+
" collective variables.\n");
for (size_t i = 0; i < cv.size(); i++) {
if (cv[i]->type() != colvarvalue::type_scalar) {
cvm::fatal_error ("Colvar grids can only be automatically "
"constructed for scalar variables. "
"ABF and histogram can not be used; metadynamics "
"can be used with useGrids disabled.\n");
}
if (cv[i]->width <= 0.0) {
cvm::fatal_error ("Tried to initialize a grid on a "
"variable with negative or zero width.\n");
}
if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) {
cvm::fatal_error ("Tried to initialize a grid on a "
"variable with undefined boundaries.\n");
}
widths.push_back (cv[i]->width);
hard_lower_boundaries.push_back (cv[i]->hard_lower_boundary);
hard_upper_boundaries.push_back (cv[i]->hard_upper_boundary);
periodic.push_back (cv[i]->periodic_boundaries());
// By default, get reported colvar value (for extended Lagrangian colvars)
actual_value.push_back (false);
// except if a colvar is specified twice in a row
// then the first instance is the actual value
if (i > 0 && cv[i-1] == cv[i]) {
actual_value[i-1] = true;
}
if (margin) {
if (periodic[i]) {
// Shift the grid by half the bin width (values at edges instead of center of bins)
lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
upper_boundaries.push_back (cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
} else {
// Make this grid larger by one bin width
lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
upper_boundaries.push_back (cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
}
} else {
lower_boundaries.push_back (cv[i]->lower_boundary);
upper_boundaries.push_back (cv[i]->upper_boundary);
}
{
cvm::real nbins = ( upper_boundaries[i].real_value -
lower_boundaries[i].real_value ) / widths[i];
int nbins_round = (int)(nbins+0.5);
if (std::fabs (nbins - cvm::real (nbins_round)) > 1.0E-10) {
cvm::log ("Warning: grid interval ("+
cvm::to_str (lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
cvm::to_str (upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
") is not commensurate to its bin width ("+
cvm::to_str (widths[i], cvm::cv_width, cvm::cv_prec)+").\n");
upper_boundaries[i].real_value = lower_boundaries[i].real_value +
(nbins_round * widths[i]);
}
if (cvm::debug())
cvm::log ("Number of points is "+cvm::to_str ((int) nbins_round)+
" for the colvar no. "+cvm::to_str (i+1)+".\n");
nx_i.push_back (nbins_round);
}
}
create (nx_i, t, mult_i);
}
/// Wrap an index vector around periodic boundary conditions
/// also checks validity of non-periodic indices
inline void wrap (std::vector<int> & ix) const
{
for (size_t i = 0; i < nd; i++) {
if (periodic[i]) {
ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
} else {
if (ix[i] < 0 || ix[i] >= nx[i])
cvm::fatal_error ("Trying to wrap illegal index vector (non-PBC): "
+ cvm::to_str (ix));
}
}
}
/// \brief Report the bin corresponding to the current value of variable i
inline int current_bin_scalar(int const i) const
{
return value_to_bin_scalar (actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
}
/// \brief Use the lower boundary and the width to report which bin
/// the provided value is in
inline int value_to_bin_scalar (colvarvalue const &value, const int i) const
{
return (int) std::floor ( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
}
/// \brief Same as the standard version, but uses another grid definition
inline int value_to_bin_scalar (colvarvalue const &value,
colvarvalue const &new_offset,
cvm::real const &new_width) const
{
return (int) std::floor ( (value.real_value - new_offset.real_value) / new_width );
}
/// \brief Use the two boundaries and the width to report the
/// central value corresponding to a bin index
inline colvarvalue bin_to_value_scalar (int const &i_bin, int const i) const
{
return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
}
/// \brief Same as the standard version, but uses different parameters
inline colvarvalue bin_to_value_scalar (int const &i_bin,
colvarvalue const &new_offset,
cvm::real const &new_width) const
{
return new_offset.real_value + new_width * (0.5 + i_bin);
}
/// Set the value at the point with index ix
inline void set_value (std::vector<int> const &ix,
T const &t,
size_t const &imult = 0)
{
data[this->address (ix)+imult] = t;
has_data = true;
}
/// \brief Get the binned value indexed by ix, or the first of them
/// if the multiplicity is larger than 1
inline T const & value (std::vector<int> const &ix,
size_t const &imult = 0) const
{
return data[this->address (ix) + imult];
}
/// \brief Add a constant to all elements (fast loop)
inline void add_constant (T const &t)
{
for (size_t i = 0; i < nt; i++)
data[i] += t;
has_data = true;
}
/// \brief Multiply all elements by a scalar constant (fast loop)
inline void multiply_constant (cvm::real const &a)
{
for (size_t i = 0; i < nt; i++)
data[i] *= a;
}
/// \brief Get the bin indices corresponding to the provided values of
/// the colvars
inline std::vector<int> const get_colvars_index (std::vector<colvarvalue> const &values) const
{
std::vector<int> index = new_index();
for (size_t i = 0; i < nd; i++) {
index[i] = value_to_bin_scalar (values[i], i);
}
return index;
}
/// \brief Get the bin indices corresponding to the current values
/// of the colvars
inline std::vector<int> const get_colvars_index() const
{
std::vector<int> index = new_index();
for (size_t i = 0; i < nd; i++) {
index[i] = current_bin_scalar (i);
}
return index;
}
/// \brief Get the minimal distance (in number of bins) from the
/// boundaries; a negative number is returned if the given point is
/// off-grid
inline cvm::real bin_distance_from_boundaries (std::vector<colvarvalue> const &values,
bool skip_hard_boundaries = false)
{
cvm::real minimum = 1.0E+16;
for (size_t i = 0; i < nd; i++) {
if (periodic[i]) continue;
cvm::real dl = std::sqrt (cv[i]->dist2 (values[i], lower_boundaries[i])) / widths[i];
cvm::real du = std::sqrt (cv[i]->dist2 (values[i], upper_boundaries[i])) / widths[i];
if (values[i].real_value < lower_boundaries[i])
dl *= -1.0;
if (values[i].real_value > upper_boundaries[i])
du *= -1.0;
if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum))
minimum = dl;
if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum))
minimum = du;
}
return minimum;
}
/// \brief Add data from another grid of the same type
///
/// Note: this function maps other_grid inside this one regardless
/// of whether it fits or not.
void map_grid (colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity())
cvm::fatal_error ("Error: trying to merge two grids with values of "
"different multiplicity.\n");
std::vector<colvarvalue> const &gb = this->lower_boundaries;
std::vector<cvm::real> const &gw = this->widths;
std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
std::vector<cvm::real> const &ogw = other_grid.widths;
std::vector<int> ix = this->new_index();
std::vector<int> oix = other_grid.new_index();
if (cvm::debug())
cvm::log ("Remapping grid...\n");
for ( ; this->index_ok (ix); this->incr (ix)) {
for (size_t i = 0; i < nd; i++) {
oix[i] =
value_to_bin_scalar (bin_to_value_scalar (ix[i], gb[i], gw[i]),
ogb[i],
ogw[i]);
}
if (! other_grid.index_ok (oix)) {
continue;
}
for (size_t im = 0; im < mult; im++) {
this->set_value (ix, other_grid.value (oix, im), im);
}
}
has_data = true;
if (cvm::debug())
cvm::log ("Remapping done.\n");
}
/// \brief Add data from another grid of the same type, AND
/// identical definition (boundaries, widths)
void add_grid (colvar_grid<T> const &other_grid,
cvm::real scale_factor = 1.0)
{
if (other_grid.multiplicity() != this->multiplicity())
cvm::fatal_error ("Error: trying to sum togetehr two grids with values of "
"different multiplicity.\n");
if (scale_factor != 1.0)
for (size_t i = 0; i < data.size(); i++) {
data[i] += scale_factor * other_grid.data[i];
}
else
// skip multiplication if possible
for (size_t i = 0; i < data.size(); i++) {
data[i] += other_grid.data[i];
}
has_data = true;
}
/// \brief Return the value suitable for output purposes (so that it
/// may be rescaled or manipulated without changing it permanently)
virtual inline T value_output (std::vector<int> const &ix,
size_t const &imult = 0)
{
return value (ix, imult);
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (the two may be different,
/// e.g. when using colvar_grid_count)
virtual inline void value_input (std::vector<int> const &ix,
T const &t,
size_t const &imult = 0,
bool add = false)
{
if ( add )
data[address (ix) + imult] += t;
else
data[address (ix) + imult] = t;
has_data = true;
}
// /// Get the pointer to the binned value indexed by ix
// inline T const *value_p (std::vector<int> const &ix)
// {
// return &(data[address (ix)]);
// }
/// \brief Get the index corresponding to the "first" bin, to be
/// used as the initial value for an index in looping
inline std::vector<int> const new_index() const
{
return std::vector<int> (nd, 0);
}
/// \brief Check that the index is within range in each of the
/// dimensions
inline bool index_ok (std::vector<int> const &ix) const
{
for (size_t i = 0; i < nd; i++) {
if ( (ix[i] < 0) || (ix[i] >= int (nx[i])) )
return false;
}
return true;
}
/// \brief Increment the index, in a way that will make it loop over
/// the whole nd-dimensional array
inline void incr (std::vector<int> &ix) const
{
for (int i = ix.size()-1; i >= 0; i--) {
ix[i]++;
if (ix[i] >= nx[i]) {
if (i > 0) {
ix[i] = 0;
continue;
} else {
// this is the last iteration, a non-valid index is being
// set for the outer index, which will be caught by
// index_ok()
ix[0] = nx[0];
return;
}
} else {
return;
}
}
}
/// \brief Write the grid parameters (number of colvars, boundaries, width and number of points)
std::ostream & write_params (std::ostream &os)
{
os << "grid_parameters {\n n_colvars " << nd << "\n";
os << " lower_boundaries ";
for (size_t i = 0; i < nd; i++)
os << " " << lower_boundaries[i];
os << "\n";
os << " upper_boundaries ";
for (size_t i = 0; i < nd; i++)
os << " " << upper_boundaries[i];
os << "\n";
os << " widths ";
for (size_t i = 0; i < nd; i++)
os << " " << widths[i];
os << "\n";
os << " sizes ";
for (size_t i = 0; i < nd; i++)
os << " " << nx[i];
os << "\n";
os << "}\n";
return os;
}
bool parse_params (std::string const &conf)
{
std::vector<int> old_nx = nx;
std::vector<colvarvalue> old_lb = lower_boundaries;
{
size_t nd_in = 0;
colvarparse::get_keyval (conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
if (nd_in != nd)
cvm::fatal_error ("Error: trying to read data for a grid "
"that contains a different number of colvars ("+
cvm::to_str (nd_in)+") than the grid defined "
"in the configuration file ("+cvm::to_str (nd)+
").\n");
}
colvarparse::get_keyval (conf, "lower_boundaries",
lower_boundaries, lower_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval (conf, "upper_boundaries",
upper_boundaries, upper_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval (conf, "widths", widths, widths, colvarparse::parse_silent);
colvarparse::get_keyval (conf, "sizes", nx, nx, colvarparse::parse_silent);
bool new_params = false;
for (size_t i = 0; i < nd; i++) {
if ( (old_nx[i] != nx[i]) ||
(std::sqrt (cv[i]->dist2 (old_lb[i],
lower_boundaries[i])) > 1.0E-10) ) {
new_params = true;
}
}
// reallocate the array in case the grid params have just changed
if (new_params) {
data.resize (0);
this->create (nx, T(), mult);
}
return true;
}
/// \brief Check that the grid information inside (boundaries,
/// widths, ...) is consistent with the current setting of the
/// colvars
void check_consistency()
{
for (size_t i = 0; i < nd; i++) {
if ( (std::sqrt (cv[i]->dist2 (cv[i]->lower_boundary,
lower_boundaries[i])) > 1.0E-10) ||
(std::sqrt (cv[i]->dist2 (cv[i]->upper_boundary,
upper_boundaries[i])) > 1.0E-10) ||
(std::sqrt (cv[i]->dist2 (cv[i]->width,
widths[i])) > 1.0E-10) ) {
cvm::fatal_error ("Error: restart information for a grid is "
"inconsistent with that of its colvars.\n");
}
}
}
/// \brief Check that the grid information inside (boundaries,
/// widths, ...) is consistent with the one of another grid
void check_consistency (colvar_grid<T> const &other_grid)
{
for (size_t i = 0; i < nd; i++) {
// we skip dist2(), because periodicities and the like should
// matter: boundaries should be EXACTLY the same (otherwise,
// map_grid() should be used)
if ( (std::fabs (other_grid.lower_boundaries[i] -
lower_boundaries[i]) > 1.0E-10) ||
(std::fabs (other_grid.upper_boundaries[i] -
upper_boundaries[i]) > 1.0E-10) ||
(std::fabs (other_grid.widths[i] -
widths[i]) > 1.0E-10) ||
(data.size() != other_grid.data.size()) ) {
cvm::fatal_error ("Error: inconsistency between "
"two grids that are supposed to be equal, "
"aside from the data stored.\n");
}
}
}
/// \brief Write the grid data without labels, as they are
/// represented in memory
/// \param buf_size Number of values per line
std::ostream & write_raw (std::ostream &os,
size_t const buf_size = 3)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
std::vector<int> ix = new_index();
size_t count = 0;
for ( ; index_ok (ix); incr (ix)) {
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw (w) << std::setprecision (p)
<< value_output (ix, imult);
if (((++count) % buf_size) == 0)
os << "\n";
}
}
// write a final newline only if buffer is not empty
if ((count % buf_size) != 0)
os << "\n";
return os;
}
/// \brief Read data written by colvar_grid::write_raw()
std::istream & read_raw (std::istream &is)
{
size_t const start_pos = is.tellg();
for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix)) {
for (size_t imult = 0; imult < mult; imult++) {
T new_value;
if (is >> new_value) {
value_input (ix, new_value, imult);
} else {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
}
}
has_data = true;
return is;
}
/// \brief To be called after colvar_grid::read_raw() returns an error
void read_raw_error()
{
cvm::fatal_error ("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n");
}
/// \brief Write the grid in a format which is both human readable
/// and suitable for visualization e.g. with gnuplot
void write_multicol (std::ostream &os)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
os << std::setw (2) << "# " << nd << "\n";
for (size_t i = 0; i < nd; i++) {
os << "# "
<< std::setw (10) << lower_boundaries[i]
<< std::setw (10) << widths[i]
<< std::setw (10) << nx[i] << " "
<< periodic[i] << "\n";
}
for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
if (ix.back() == 0) {
// if the last index is 0, add a new line to mark the new record
os << "\n";
}
for (size_t i = 0; i < nd; i++) {
os << " "
<< std::setw (w) << std::setprecision (p)
<< bin_to_value_scalar (ix[i], i);
}
os << " ";
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw (w) << std::setprecision (p)
<< value_output (ix, imult);
}
os << "\n";
}
}
/// \brief Read a grid written by colvar_grid::write_multicol()
/// Adding data if add is true, replacing if false
std::istream & read_multicol (std::istream &is, bool add = false)
{
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
std::string hash;
cvm::real lower, width, x;
size_t n, periodic;
bool remap;
std::vector<T> new_value;
std::vector<int> nx_read;
std::vector<int> bin;
if ( cv.size() != nd ) {
cvm::fatal_error ("Cannot read grid file: missing reference to colvars.");
}
if ( !(is >> hash) || (hash != "#") ) {
cvm::fatal_error ("Error reading grid at position "+
cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
}
is >> n;
if ( n != nd ) {
cvm::fatal_error ("Error reading grid: wrong number of collective variables.\n");
}
nx_read.resize (n);
bin.resize (n);
new_value.resize (mult);
if (this->has_parent_data && add) {
new_data.resize (data.size());
}
remap = false;
for (size_t i = 0; i < nd; i++ ) {
if ( !(is >> hash) || (hash != "#") ) {
cvm::fatal_error ("Error reading grid at position "+
cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
}
is >> lower >> width >> nx_read[i] >> periodic;
if ( (std::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) ||
(std::fabs (width - widths[i] ) > 1.0e-10) ||
(nx_read[i] != nx[i]) ) {
cvm::log ("Warning: reading from different grid definition (colvar "
+ cvm::to_str (i+1) + "); remapping data on new grid.\n");
remap = true;
}
}
if ( remap ) {
// re-grid data
while (is.good()) {
bool end_of_file = false;
for (size_t i = 0; i < nd; i++ ) {
if ( !(is >> x) ) end_of_file = true;
bin[i] = value_to_bin_scalar (x, i);
}
if (end_of_file) break;
for (size_t imult = 0; imult < mult; imult++) {
is >> new_value[imult];
}
if ( index_ok(bin) ) {
for (size_t imult = 0; imult < mult; imult++) {
value_input (bin, new_value[imult], imult, add);
}
}
}
} else {
// do not re-grid the data but assume the same grid is used
for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
for (size_t i = 0; i < nd; i++ ) {
is >> x;
}
for (size_t imult = 0; imult < mult; imult++) {
is >> new_value[imult];
value_input (ix, new_value[imult], imult, add);
}
}
}
has_data = true;
return is;
}
};
/// \brief Colvar_grid derived class to hold counters in discrete
/// n-dim colvar space
class colvar_grid_count : public colvar_grid<size_t>
{
public:
/// Default constructor
colvar_grid_count();
/// Destructor
virtual inline ~colvar_grid_count()
{}
/// Constructor
colvar_grid_count (std::vector<int> const &nx_i,
size_t const &def_count = 0);
/// Constructor from a vector of colvars
colvar_grid_count (std::vector<colvar *> &colvars,
size_t const &def_count = 0);
/// Increment the counter at given position
inline void incr_count (std::vector<int> const &ix)
{
++(data[this->address (ix)]);
}
/// \brief Get the binned count indexed by ix from the newly read data
inline size_t const & new_count (std::vector<int> const &ix,
size_t const &imult = 0)
{
return new_data[address (ix) + imult];
}
/// \brief Read the grid from a restart
std::istream & read_restart (std::istream &is);
/// \brief Write the grid to a restart
std::ostream & write_restart (std::ostream &os);
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual inline void value_input (std::vector<int> const &ix,
size_t const &t,
size_t const &imult = 0,
bool add = false)
{
if (add) {
data[address (ix)] += t;
if (this->has_parent_data) {
// save newly read data for inputting parent grid
new_data[address (ix)] = t;
}
} else {
data[address (ix)] = t;
}
has_data = true;
}
};
/// Class for accumulating a scalar function on a grid
class colvar_grid_scalar : public colvar_grid<cvm::real>
{
public:
/// \brief Provide the associated sample count by which each binned value
/// should be divided
colvar_grid_count *samples;
/// Default constructor
colvar_grid_scalar();
/// Copy constructor (needed because of the grad pointer)
colvar_grid_scalar (colvar_grid_scalar const &g);
/// Destructor
~colvar_grid_scalar();
/// Constructor from specific sizes arrays
colvar_grid_scalar (std::vector<int> const &nx_i);
/// Constructor from a vector of colvars
colvar_grid_scalar (std::vector<colvar *> &colvars,
bool margin = 0);
/// Accumulate the value
inline void acc_value (std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0)
{
// only legal value of imult here is 0
data[address (ix)] += new_value;
if (samples)
samples->incr_count (ix);
has_data = true;
}
/// Return the gradient of the scalar field from finite differences
inline const cvm::real * gradient_finite_diff ( const std::vector<int> &ix0 )
{
cvm::real A0, A1;
std::vector<int> ix;
if (nd != 2) cvm::fatal_error ("Finite differences available in dimension 2 only.");
for (int n = 0; n < nd; n++) {
ix = ix0;
A0 = data[address (ix)];
ix[n]++; wrap (ix);
A1 = data[address (ix)];
ix[1-n]++; wrap (ix);
A1 += data[address (ix)];
ix[n]--; wrap (ix);
A0 += data[address (ix)];
grad[n] = 0.5 * (A1 - A0) / widths[n];
}
return grad;
}
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual cvm::real value_output (std::vector<int> const &ix,
size_t const &imult = 0)
{
if (imult > 0)
cvm::fatal_error ("Error: trying to access a component "
"larger than 1 in a scalar data grid.\n");
if (samples)
return (samples->value (ix) > 0) ?
(data[address (ix)] / cvm::real (samples->value (ix))) :
0.0;
else
return data[address (ix)];
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual void value_input (std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0,
bool add = false)
{
if (imult > 0)
cvm::fatal_error ("Error: trying to access a component "
"larger than 1 in a scalar data grid.\n");
if (add) {
if (samples)
data[address (ix)] += new_value * samples->new_count (ix);
else
data[address (ix)] += new_value;
} else {
if (samples)
data[address (ix)] = new_value * samples->value (ix);
else
data[address (ix)] = new_value;
}
has_data = true;
}
/// \brief Read the grid from a restart
std::istream & read_restart (std::istream &is);
/// \brief Write the grid to a restart
std::ostream & write_restart (std::ostream &os);
/// \brief Return the highest value
inline cvm::real maximum_value()
{
cvm::real max = data[0];
for (size_t i = 0; i < nt; i++) {
if (data[i] > max) max = data[i];
}
return max;
}
/// \brief Return the lowest value
inline cvm::real minimum_value()
{
cvm::real min = data[0];
for (size_t i = 0; i < nt; i++) {
if (data[i] < min) min = data[i];
}
return min;
}
private:
// gradient
cvm::real * grad;
};
/// Class for accumulating the gradient of a scalar function on a grid
class colvar_grid_gradient : public colvar_grid<cvm::real>
{
public:
/// \brief Provide the sample count by which each binned value
/// should be divided
colvar_grid_count *samples;
/// Default constructor
colvar_grid_gradient();
/// Destructor
virtual inline ~colvar_grid_gradient()
{}
/// Constructor from specific sizes arrays
colvar_grid_gradient (std::vector<int> const &nx_i);
/// Constructor from a vector of colvars
colvar_grid_gradient (std::vector<colvar *> &colvars);
/// \brief Accumulate the gradient
inline void acc_grad (std::vector<int> const &ix, cvm::real const *grads) {
for (size_t imult = 0; imult < mult; imult++) {
data[address (ix) + imult] += grads[imult];
}
if (samples)
samples->incr_count (ix);
}
/// \brief Accumulate the gradient based on the force (i.e. sums the
/// opposite of the force)
inline void acc_force (std::vector<int> const &ix, cvm::real const *forces) {
for (size_t imult = 0; imult < mult; imult++) {
data[address (ix) + imult] -= forces[imult];
}
if (samples)
samples->incr_count (ix);
}
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual inline cvm::real value_output (std::vector<int> const &ix,
size_t const &imult = 0)
{
if (samples)
return (samples->value (ix) > 0) ?
(data[address (ix) + imult] / cvm::real (samples->value (ix))) :
0.0;
else
return data[address (ix) + imult];
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual inline void value_input (std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0,
bool add = false)
{
if (add) {
if (samples)
data[address (ix) + imult] += new_value * samples->new_count (ix);
else
data[address (ix) + imult] += new_value;
} else {
if (samples)
data[address (ix) + imult] = new_value * samples->value (ix);
else
data[address (ix) + imult] = new_value;
}
has_data = true;
}
/// \brief Read the grid from a restart
std::istream & read_restart (std::istream &is);
/// \brief Write the grid to a restart
std::ostream & write_restart (std::ostream &os);
/// Compute and return average value for a 1D gradient grid
inline cvm::real average()
{
size_t n = 0;
if (nd != 1 || nx[0] == 0) {
return 0.0;
}
cvm::real sum = 0.0;
std::vector<int> ix = new_index();
if (samples) {
for ( ; index_ok (ix); incr (ix)) {
if ( (n = samples->value (ix)) )
sum += value (ix) / n;
}
} else {
for ( ; index_ok (ix); incr (ix)) {
sum += value (ix);
}
}
return (sum / cvm::real (nx[0]));
}
/// \brief If the grid is 1-dimensional, integrate it and write the
/// integral to a file
void write_1D_integral (std::ostream &os);
};
#endif
c4science · Help