Difference for src/colvarbias_restraint.C from version 1.3 to 1.4

version 1.3version 1.4
Line 1
Line 1
 /// -*- c++ -*- // -*- c++ -*-
  
 #include "colvarmodule.h" #include "colvarmodule.h"
 #include "colvarvalue.h" #include "colvarvalue.h"
 #include "colvarbias_restraint.h" #include "colvarbias_restraint.h"
  
  
 colvarbias_restraint::colvarbias_restraint(std::string const &conf, colvarbias_restraint::colvarbias_restraint(char const *key)
                                            char const *key)   : colvarbias(key)
   : colvarbias(conf, key), 
     target_nstages(0), 
     target_nsteps(0) 
 { {
  }
  
  
  int colvarbias_restraint::init(std::string const &conf)
  {
    colvarbias::init(conf);
  
    if (cvm::debug())
      cvm::log("Initializing a new restraint bias.\n");
  
    // TODO move these initializations to constructor and let get_keyval
    // only override existing values
    target_nstages = 0;
    target_nsteps = 0;
    force_k = 0.0;
  
   get_keyval(conf, "forceConstant", force_k, 1.0);   get_keyval(conf, "forceConstant", force_k, 1.0);
  
   {   {
Line 18
Line 31
     colvar_centers.resize(colvars.size());     colvar_centers.resize(colvars.size());
     colvar_centers_raw.resize(colvars.size());     colvar_centers_raw.resize(colvars.size());
     size_t i;     size_t i;
  
      enable(f_cvb_apply_force);
  
     for (i = 0; i < colvars.size(); i++) {     for (i = 0; i < colvars.size(); i++) {
       colvar_centers[i].type(colvars[i]->value());       colvar_centers[i].type(colvars[i]->value());
       colvar_centers_raw[i].type(colvars[i]->value());       colvar_centers_raw[i].type(colvars[i]->value());
Line 53
Line 69
  
     size_t i;     size_t i;
     if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {     if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
        if (colvar_centers.size() != colvars.size()) {
          cvm::error("Error: number of target centers does not match "
                     "that of collective variables.\n");
        }
       b_chg_centers = true;       b_chg_centers = true;
       for (i = 0; i < target_centers.size(); i++) {       for (i = 0; i < target_centers.size(); i++) {
         target_centers[i].apply_constraints();         target_centers[i].apply_constraints();
Line 115
Line 135
  
   if (cvm::debug())   if (cvm::debug())
     cvm::log("Done initializing a new restraint bias.\n");     cvm::log("Done initializing a new restraint bias.\n");
  
    return COLVARS_OK;
 } }
  
  
 colvarbias_restraint::~colvarbias_restraint() colvarbias_restraint::~colvarbias_restraint()
 { {
   if (cvm::n_rest_biases > 0)   if (cvm::n_rest_biases > 0)
Line 167
Line 190
 } }
  
  
 cvm::real colvarbias_restraint::update() int colvarbias_restraint::update()
 { {
   bias_energy = 0.0;   bias_energy = 0.0;
  
Line 333
Line 356
     cvm::log("Current forces for the restraint bias \""+     cvm::log("Current forces for the restraint bias \""+
               this->name+"\": "+cvm::to_str(colvar_forces)+".\n");               this->name+"\": "+cvm::to_str(colvar_forces)+".\n");
  
   return bias_energy;   return COLVARS_OK;
 } }
  
  
Line 498
Line 521
   return os;   return os;
 } }
  
 colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(std::string const &conf, char const *key) : 
   colvarbias_restraint(conf, key) { colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
    : colvarbias_restraint(key)
  {
  }
  
  
  int colvarbias_restraint_harmonic::init(std::string const &conf)
  {
    colvarbias_restraint::init(conf);
  
   for (size_t i = 0; i < colvars.size(); i++) {   for (size_t i = 0; i < colvars.size(); i++) {
     if (colvars[i]->width != 1.0)     if (colvars[i]->width != 1.0)
       cvm::log("The force constant for colvar \""+colvars[i]->name+       cvm::log("The force constant for colvar \""+colvars[i]->name+
Line 507
Line 539
                 cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+                 cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
                 " according to the specified width.\n");                 " according to the specified width.\n");
   }   }
    return COLVARS_OK;
 } }
  
 cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const 
  cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k,
                                                               colvar const *x,
                                                               colvarvalue const &xcenter) const
 { {
   return 0.5 * k * x->dist2(x->value(), xcenter);   return 0.5 * k * x->dist2(x->value(), xcenter);
 } }
  
 colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const 
  colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k,
                                                             colvar const *x,
                                                             colvarvalue const &xcenter) const
 { {
   return 0.5 * k * x->dist2_lgrad(x->value(), xcenter);   return 0.5 * k * x->dist2_lgrad(x->value(), xcenter);
 } }
  
 cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k, cvm::real dist_measure) const 
  cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k,
                                                               cvm::real dist_measure) const
 { {
   return k / (dist_measure * dist_measure);   return k / (dist_measure * dist_measure);
 } }
  
  
 colvarbias_restraint_linear::colvarbias_restraint_linear(std::string const &conf, char const *key) : 
   colvarbias_restraint(conf, key) { colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
    : colvarbias_restraint(key)
  {
  }
  
  
  int colvarbias_restraint_linear::init(std::string const &conf)
  {
    colvarbias_restraint::init(conf);
  
   for (size_t i = 0; i < colvars.size(); i++) {   for (size_t i = 0; i < colvars.size(); i++) {
     if (colvars[i]->width != 1.0)     if (colvars[i]->width != 1.0)
       cvm::log("The force constant for colvar \""+colvars[i]->name+       cvm::log("The force constant for colvar \""+colvars[i]->name+
Line 534
Line 584
                 cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+                 cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+
                 " according to the specified width.\n");                 " according to the specified width.\n");
   }   }
    return COLVARS_OK;
 } }
  
 cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const 
  cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k,
                                                             colvar const *x,
                                                             colvarvalue const &xcenter) const
 { {
   return k * (x->value() - xcenter);   return k * (x->value() - xcenter);
 } }
  
 colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k,  colvar* x,  const colvarvalue &xcenter) const 
  colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k,
                                                           colvar const *x,
                                                           colvarvalue const &xcenter) const
 { {
   return k;   return k;
 } }
  
 cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k, cvm::real dist_measure) const 
  cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k,
                                                             cvm::real dist_measure) const
 { {
   return k / dist_measure;   return k / dist_measure;
 } }
  
  
  
  colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key)
    : colvarbias(key)
  {
    lower_boundary = 0.0;
    upper_boundary = 0.0;
    width = 0.0;
    gaussian_width = 0.0;
  }
  
  
  int colvarbias_restraint_histogram::init(std::string const &conf)
  {
    colvarbias::init(conf);
  
    get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary);
    get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary);
    get_keyval(conf, "width", width, width);
  
    if (width <= 0.0) {
      cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
    }
  
    get_keyval(conf, "gaussianWidth", gaussian_width, 2.0 * width, colvarparse::parse_silent);
    get_keyval(conf, "gaussianSigma", gaussian_width, 2.0 * width);
  
    if (lower_boundary >= upper_boundary) {
      cvm::error("Error: the upper boundary, "+
                 cvm::to_str(upper_boundary)+
                 ", is not higher than the lower boundary, "+
                 cvm::to_str(lower_boundary)+".\n",
                 INPUT_ERROR);
    }
  
    cvm::real const nbins = (upper_boundary - lower_boundary) / width;
    int const nbins_round = (int)(nbins);
  
    if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
      cvm::log("Warning: grid interval ("+
               cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+
               cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+
               ") is not commensurate to its bin width ("+
               cvm::to_str(width, cvm::cv_width, cvm::cv_prec)+").\n");
    }
  
    p.resize(nbins_round);
    ref_p.resize(nbins_round);
    p_diff.resize(nbins_round);
  
    bool const inline_ref_p =
      get_keyval(conf, "refHistogram", ref_p.data_array(), ref_p.data_array());
    std::string ref_p_file;
    get_keyval(conf, "refHistogramFile", ref_p_file, std::string(""));
    if (ref_p_file.size()) {
      if (inline_ref_p) {
        cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n",
                   INPUT_ERROR);
      } else {
        std::ifstream is(ref_p_file.c_str());
        std::string data_s = "";
        std::string line;
        while (getline_nocomments(is, line)) {
          data_s.append(line+"\n");
        }
        if (data_s.size() == 0) {
          cvm::error("Error: file \""+ref_p_file+"\" empty or unreadable.\n", FILE_ERROR);
        }
        is.close();
        cvm::vector1d<cvm::real> data;
        if (data.from_simple_string(data_s) != 0) {
          cvm::error("Error: could not read histogram from file \""+ref_p_file+"\".\n");
        }
        if (data.size() == 2*ref_p.size()) {
          // file contains both x and p(x)
          size_t i;
          for (i = 0; i < ref_p.size(); i++) {
            ref_p[i] = data[2*i+1];
          }
        } else if (data.size() == ref_p.size()) {
          ref_p = data;
        } else {
          cvm::error("Error: file \""+ref_p_file+"\" contains a histogram of different length.\n",
                     INPUT_ERROR);
        }
      }
    }
    cvm::real const ref_integral = ref_p.sum() * width;
    if (std::fabs(ref_integral - 1.0) > 1.0e-03) {
      cvm::log("Reference distribution not normalized, normalizing to unity.\n");
      ref_p /= ref_integral;
    }
  
    get_keyval(conf, "writeHistogram", b_write_histogram, false);
    get_keyval(conf, "forceConstant", force_k, 1.0);
  
    return COLVARS_OK;
  }
  
  
  colvarbias_restraint_histogram::~colvarbias_restraint_histogram()
  {
    p.resize(0);
    ref_p.resize(0);
    p_diff.resize(0);
  }
  
  
  int colvarbias_restraint_histogram::update()
  {
    if (cvm::debug())
      cvm::log("Updating the histogram restraint bias \""+this->name+"\".\n");
  
    size_t vector_size = 0;
    size_t icv;
    for (icv = 0; icv < colvars.size(); icv++) {
      vector_size += colvars[icv]->value().size();
    }
  
    cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
  
    // calculate the histogram
    p.reset();
    for (icv = 0; icv < colvars.size(); icv++) {
      colvarvalue const &cv = colvars[icv]->value();
      if (cv.type() == colvarvalue::type_scalar) {
        cvm::real const cv_value = cv.real_value;
        size_t igrid;
        for (igrid = 0; igrid < p.size(); igrid++) {
          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
          p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                                      (2.0 * gaussian_width * gaussian_width));
        }
      } else if (cv.type() == colvarvalue::type_vector) {
        size_t idim;
        for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
          cvm::real const cv_value = cv.vector1d_value[idim];
          size_t igrid;
          for (igrid = 0; igrid < p.size(); igrid++) {
            cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
            p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                                        (2.0 * gaussian_width * gaussian_width));
          }
        }
      } else {
        // TODO
      }
    }
  
    cvm::real const force_k_cv = force_k * vector_size;
  
    // calculate the difference between current and reference
    p_diff = p - ref_p;
    bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
  
    // calculate the forces
    for (icv = 0; icv < colvars.size(); icv++) {
      colvarvalue const &cv = colvars[icv]->value();
      colvarvalue &cv_force = colvar_forces[icv];
      cv_force.type(cv);
      cv_force.reset();
  
      if (cv.type() == colvarvalue::type_scalar) {
        cvm::real const cv_value = cv.real_value;
        cvm::real &force = cv_force.real_value;
        size_t igrid;
        for (igrid = 0; igrid < p.size(); igrid++) {
          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
          force += force_k_cv * p_diff[igrid] *
            norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                            (2.0 * gaussian_width * gaussian_width)) *
            (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
        }
      } else if (cv.type() == colvarvalue::type_vector) {
        size_t idim;
        for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
          cvm::real const cv_value = cv.vector1d_value[idim];
          cvm::real &force = cv_force.vector1d_value[idim];
          size_t igrid;
          for (igrid = 0; igrid < p.size(); igrid++) {
            cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
            force += force_k_cv * p_diff[igrid] *
              norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
                              (2.0 * gaussian_width * gaussian_width)) *
              (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
          }
        }
      } else {
        // TODO
      }
    }
  
    return COLVARS_OK;
  }
  
  
  std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
  {
    if (b_write_histogram) {
      std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat");
      std::ofstream os(file_name.c_str());
      os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width)
         << "  " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3)
         << ")\n";
      size_t igrid;
      for (igrid = 0; igrid < p.size(); igrid++) {
        cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
        os << "  "
           << std::setprecision(cvm::cv_prec)
           << std::setw(cvm::cv_width)
           << x_grid
           << "  "
           << std::setprecision(cvm::cv_prec)
           << std::setw(cvm::cv_width)
           << p[igrid] << "\n";
      }
      os.close();
    }
    return os;
  }
  
  
  std::istream & colvarbias_restraint_histogram::read_restart(std::istream &is)
  {
    return is;
  }
  
  
  std::ostream & colvarbias_restraint_histogram::write_traj_label(std::ostream &os)
  {
    os << " ";
    if (b_output_energy) {
      os << " E_"
         << cvm::wrap_string(this->name, cvm::en_width-2);
    }
    return os;
  }
  
  
  std::ostream & colvarbias_restraint_histogram::write_traj(std::ostream &os)
  {
    os << " ";
    if (b_output_energy) {
      os << " "
         << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
         << bias_energy;
    }
    return os;
  }

Legend:
Removed in v.1.3 
changed lines
 Added in v.1.4


Made by using version 1.53 of cvs2html