lammps/lib/colvars/colvargrid.cpp

/// -*- c++ -*-

#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
#include "colvargrid.h"


colvar_grid_count::colvar_grid_count()
  : colvar_grid<size_t>()
{
  mult = 1;
}

colvar_grid_count::colvar_grid_count(std::vector<int> const &nx_i,
                                     size_t const           &def_count)
  : colvar_grid<size_t>(nx_i, def_count, 1)
{}

colvar_grid_count::colvar_grid_count(std::vector<colvar *>  &colvars,
                                     size_t const           &def_count)
  : colvar_grid<size_t>(colvars, def_count, 1)
{}

std::istream & colvar_grid_count::read_restart(std::istream &is)
{
  size_t const start_pos = is.tellg();
  std::string key, conf;
  if ((is >> key) && (key == std::string("grid_parameters"))) {
    is.seekg(start_pos, std::ios::beg);
    is >> colvarparse::read_block("grid_parameters", conf);
    parse_params(conf);
  } else {
    cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
    is.seekg(start_pos, std::ios::beg);
  }
  read_raw(is);
  return is;
}

std::ostream & colvar_grid_count::write_restart(std::ostream &os)
{
  write_params(os);
  write_raw(os);
  return os;
}



colvar_grid_scalar::colvar_grid_scalar()
  : colvar_grid<cvm::real>(), samples(NULL), grad(NULL)
{}

colvar_grid_scalar::colvar_grid_scalar(colvar_grid_scalar const &g)
  : colvar_grid<cvm::real>(g), samples(NULL), grad(NULL)
{
  grad = new cvm::real[nd];
}

colvar_grid_scalar::colvar_grid_scalar(std::vector<int> const &nx_i)
  : colvar_grid<cvm::real>(nx_i, 0.0, 1), samples(NULL)
{
  grad = new cvm::real[nd];
}

colvar_grid_scalar::colvar_grid_scalar(std::vector<colvar *> &colvars, bool margin)
  : colvar_grid<cvm::real>(colvars, 0.0, 1, margin), samples(NULL)
{
  grad = new cvm::real[nd];
}

colvar_grid_scalar::~colvar_grid_scalar()
{
  if (grad) {
    delete [] grad;
    grad = NULL;
  }
}

std::istream & colvar_grid_scalar::read_restart(std::istream &is)
{
  size_t const start_pos = is.tellg();
  std::string key, conf;
  if ((is >> key) && (key == std::string("grid_parameters"))) {
    is.seekg(start_pos, std::ios::beg);
    is >> colvarparse::read_block("grid_parameters", conf);
    parse_params(conf);
  } else {
    cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
    is.seekg(start_pos, std::ios::beg);
  }
  read_raw(is);
  return is;
}

std::ostream & colvar_grid_scalar::write_restart(std::ostream &os)
{
  write_params(os);
  write_raw(os);
  return os;
}



colvar_grid_gradient::colvar_grid_gradient()
  : colvar_grid<cvm::real>(), samples(NULL)
{}

colvar_grid_gradient::colvar_grid_gradient(std::vector<int> const &nx_i)
  : colvar_grid<cvm::real>(nx_i, 0.0, nx_i.size()), samples(NULL)
{}

colvar_grid_gradient::colvar_grid_gradient(std::vector<colvar *> &colvars)
  : colvar_grid<cvm::real>(colvars, 0.0, colvars.size()), samples(NULL)
{}

std::istream & colvar_grid_gradient::read_restart(std::istream &is)
{
  size_t const start_pos = is.tellg();
  std::string key, conf;
  if ((is >> key) && (key == std::string("grid_parameters"))) {
    is.seekg(start_pos, std::ios::beg);
    is >> colvarparse::read_block("grid_parameters", conf);
    parse_params(conf);
  } else {
    cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
    is.seekg(start_pos, std::ios::beg);
  }
  read_raw(is);
  return is;
}

std::ostream & colvar_grid_gradient::write_restart(std::ostream &os)
{
  write_params(os);
  write_raw(os);
  return os;
}

void colvar_grid_gradient::write_1D_integral(std::ostream &os)
{
  cvm::real bin, min, integral;
  std::vector<cvm::real> int_vals;

  os << "#       xi            A(xi)\n";

  if ( cv.size() != 1 ) {
    cvm::fatal_error("Cannot write integral for multi-dimensional gradient grids.");
  }

  integral = 0.0;
  int_vals.push_back( 0.0 );
  bin = 0.0;
  min = 0.0;

  // correction for periodic colvars, so that the PMF is periodic
  cvm::real corr;
  if ( periodic[0] ) {
    corr = average();
  } else {
    corr = 0.0;
  }

  for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix), bin += 1.0 ) {

    if (samples) {
      size_t const samples_here = samples->value(ix);
      if (samples_here)
        integral += (value(ix) / cvm::real(samples_here) - corr) * cv[0]->width;
    } else {
      integral += (value(ix) - corr) * cv[0]->width;
    }

    if ( integral < min ) min = integral;
    int_vals.push_back( integral );
  }

  bin = 0.0;
  for ( int i = 0; i < nx[0]; i++, bin += 1.0 ) {
    os << std::setw(10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
       << std::setw(cvm::cv_width)
       << std::setprecision(cvm::cv_prec)
       << int_vals[i] - min << "\n";
  }

  os << std::setw(10) << cv[0]->lower_boundary.real_value + cv[0]->width * bin << " "
     << std::setw(cvm::cv_width)
     << std::setprecision(cvm::cv_prec)
     << int_vals[nx[0]] - min << "\n";

  return;
}