colvar.h

Go to the documentation of this file.

// -*- c++ -*-

#ifndef COLVAR_H
#define COLVAR_H

#include <iostream>
#include <iomanip>
#include <cmath>

#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"



class colvar : public colvarparse {

public:

  std::string name;

  colvarvalue::Type type() const;

  colvarvalue const & value() const;

  colvarvalue const & velocity() const;

  colvarvalue const & system_force() const;


  cvm::real width;

  bool b_linear;

  bool b_inverse_gradients;

  bool b_Jacobian_force;

  enum task {
    task_gradients,
    task_fdiff_velocity,
    task_system_force,
    task_extended_lagrangian,
    task_Jacobian_force,
    task_report_Jacobian_force,
    task_output_value,
    task_output_velocity,
    task_output_applied_force,
    task_output_system_force,
    task_lower_boundary,
    task_upper_boundary,
    task_grid,
    task_lower_wall,
    task_upper_wall,
    task_runave,
    task_corrfunc,
    task_ntot
  };

  bool tasks[task_ntot];

protected:


  /*
    extended:
    H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\
    + \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\
    + \sum_{t'<t} W * exp (-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\
    + \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M

    normal:
    H = H_{0} + \sum_{t'<t} W * exp (-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\
    + \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M

    output:
    H = H_{0}   (only output S(x), no forces)

    Here:
    S(x(t)) = x
    s(t)    = xr
    DS = Ds = delta
  */


  colvarvalue x;

  colvarvalue x_reported;

  colvarvalue v_fdiff;

  inline colvarvalue fdiff_velocity (colvarvalue const &xold,
                                     colvarvalue const &xnew)
  {
    // using the gradient of the square distance to calculate the
    // velocity (non-scalar variables automatically taken into
    // account)
    return ( ( (cvm::dt > 0.0) ? (1.0/cvm::dt) : 1.0 ) *
             0.5 * dist2_lgrad (xnew, xold) );
  }

  colvarvalue v_reported;

  // Options for task_extended_lagrangian
  colvarvalue xr;
  colvarvalue vr;
  cvm::real ext_mass;
  cvm::real ext_force_k;
  colvarvalue fr;

  colvarvalue fj;

  colvarvalue ft_reported;

public:


  colvarvalue fb;

  colvarvalue f;

  colvarvalue ft;


  cvm::real period;

  bool b_periodic;


  bool expand_boundaries;

  colvarvalue lower_boundary;
  colvarvalue lower_wall;
  cvm::real   lower_wall_k;

  colvarvalue upper_boundary;
  colvarvalue upper_wall;
  cvm::real   upper_wall_k;

  // TODO: move everything below into the grid class;
  // Possibly: use vectors instead of single bin number / colvar values?

  int current_bin_scalar() const;

  int value_to_bin_scalar (colvarvalue const &value) const;

  int value_to_bin_scalar (colvarvalue const &value,
                           colvarvalue const &new_offset,
                           cvm::real const   &new_width) const;

  colvarvalue bin_to_value_scalar (int const &i_bin) const;

  colvarvalue bin_to_value_scalar (int const &i_bin,
                                   colvarvalue const &new_offset,
                                   cvm::real const &new_width) const;
  // End of TODO

  bool periodic_boundaries() const;

  bool periodic_boundaries (colvarvalue const &lb, colvarvalue const &ub) const;


  colvar (std::string const &conf);

  void enable (colvar::task const &t);

  void disable (colvar::task const &t);

  ~colvar();


  void calc();

  void calc_value (colvarvalue &xn);

  void reset_bias_force();

  void add_bias_force (colvarvalue const &force);

  void update();

  void communicate_forces();


  cvm::real dist2 (colvarvalue const &x1,
                   colvarvalue const &x2) const;

  colvarvalue dist2_lgrad (colvarvalue const &x1,
                           colvarvalue const &x2) const;

  colvarvalue dist2_rgrad (colvarvalue const &x1,
                           colvarvalue const &x2) const;

  cvm::real compare (colvarvalue const &x1,
                     colvarvalue const &x2) const;




  void parse_analysis (std::string const &conf);
  void analyse();


  std::istream & read_traj (std::istream &is);
  std::ostream & write_traj (std::ostream &os);
  std::ostream & write_traj_label (std::ostream &os);


  std::istream & read_restart (std::istream &is);
  std::ostream & write_restart (std::ostream &os);


protected:

  colvarvalue            x_old;

  std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
  std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p;

  std::list< std::list<colvarvalue> > x_history;
  std::list< std::list<colvarvalue> >::iterator x_history_p;

  std::string            acf_colvar_name;
  size_t                 acf_length;
  size_t                 acf_offset;
  size_t                 acf_stride;
  size_t                 acf_nframes;
  bool                   acf_normalize;
  std::vector<cvm::real> acf;
  std::string            acf_outfile;

  enum acf_type_e {
    acf_notset,
    acf_vel,
    acf_coor,
    acf_p2coor
  };

  acf_type_e             acf_type;

  void calc_vel_acf (std::list<colvarvalue> &v_history,
                     colvarvalue const      &v);

  void calc_coor_acf (std::list<colvarvalue> &x_history,
                      colvarvalue const      &x);

  void calc_p2coor_acf (std::list<colvarvalue> &x_history,
                        colvarvalue const      &x);

  void calc_acf();
  void write_acf (std::ostream &os);

  size_t         runave_length;
  size_t         runave_stride;
  std::ofstream  runave_os;
  colvarvalue    runave;
  cvm::real      runave_variance;

  void calc_runave();


public:


  // collective variable component base class
  class cvc;

  // currently available collective variable components

  // scalar colvar components
  class distance;
  class distance_z;
  class distance_xy;
  class min_distance;
  class angle;
  class dihedral;
  class coordnum;
  class h_bond;
  class rmsd;
  class logmsd;
  class orientation_angle;
  class gyration;
  class eigenvector;
  class alpha_dihedrals;
  class alpha_angles;

  // non-scalar components
//   class distance_vec;
  class distance_dir;
  class orientation;

protected:

  std::vector<cvc *> cvcs;

public:

  inline size_t n_components () const {
    return cvcs.size();
  }
};


inline colvar * cvm::colvar_p (std::string const &name)
{
  for (std::vector<colvar *>::iterator cvi = cvm::colvars.begin();
       cvi != cvm::colvars.end();
       cvi++) {
    if ((*cvi)->name == name) {
      return (*cvi);
    }
  }
  return NULL;
}


inline colvarvalue::Type colvar::type() const
{
  return x.type();
}


inline colvarvalue const & colvar::value() const
{
  return x_reported;
}


inline colvarvalue const & colvar::velocity() const
{
  return v_reported;
}


inline colvarvalue const & colvar::system_force() const
{
  return ft_reported;
}



inline int colvar::current_bin_scalar() const
{
  return this->value_to_bin_scalar (this->value());
}

inline int colvar::value_to_bin_scalar (colvarvalue const &value) const
{
  return (int) ::floor ( (value.real_value - lower_boundary.real_value) / width );
}

inline int colvar::value_to_bin_scalar (colvarvalue const &value,
                                        colvarvalue const &new_offset,
                                        cvm::real   const &new_width) const
{
  return (int) ::floor ( (value.real_value - new_offset.real_value) / new_width );
}


inline colvarvalue colvar::bin_to_value_scalar (int const &i_bin) const
{
  return lower_boundary.real_value + width * (0.5 + i_bin);
}

inline colvarvalue colvar::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);
}


inline void colvar::add_bias_force (colvarvalue const &force)
{
  fb += force;
}

inline void colvar::reset_bias_force() {
  fb.reset();
}

#endif

---