colvarcomp.h

Go to the documentation of this file.

#ifndef COLVARDEF_H
#define COLVARDEF_H

#include <fstream>
#include <cmath>


#include "colvarmodule.h"
#include "colvar.h"
#include "colvaratoms.h"



class colvar::cvc
  : public colvarparse
{
public:

  std::string name;

  std::string function_type;

  colvarvalue::Type type() const;

  cvm::real sup_coeff;
  int       sup_np;

  cvm::real period;

  bool b_periodic;

  cvc (std::string const &conf);

  void parse_group (std::string const &conf, 
                    char const *group_key,
                    cvm::atom_group &group,
                    bool optional = false);

  cvc();

  virtual ~cvc();

  bool b_debug_gradients;

  colvarvalue fdiff_change (cvm::atom_group &group);

  bool b_inverse_gradients;

  bool b_Jacobian_derivative;

  virtual void calc_value() = 0;

  virtual void calc_gradients() = 0;

  virtual void calc_force_invgrads();

  virtual void calc_Jacobian_derivative();


  virtual colvarvalue value() const;

  virtual colvarvalue system_force() const;

  virtual colvarvalue Jacobian_derivative() const;

  virtual void apply_force (colvarvalue const &cvforce) = 0;

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

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

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

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

protected:

  colvarvalue x;

  colvarvalue x_old;

  colvarvalue ft;

  colvarvalue jd;
};




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

inline colvarvalue colvar::cvc::value() const
{
  return x;
}

inline colvarvalue colvar::cvc::system_force() const
{
  return ft;
}

inline colvarvalue colvar::cvc::Jacobian_derivative() const
{
  return jd;
}


inline cvm::real colvar::cvc::dist2 (colvarvalue const &x1,
                                     colvarvalue const &x2) const
{
  return x1.dist2 (x2);
}

inline colvarvalue colvar::cvc::dist2_lgrad (colvarvalue const &x1,
                                             colvarvalue const &x2) const
{
  return x1.dist2_grad (x2);
}

inline colvarvalue colvar::cvc::dist2_rgrad (colvarvalue const &x1,
                                             colvarvalue const &x2) const
{
  return x2.dist2_grad (x1);
}

inline cvm::real colvar::cvc::compare (colvarvalue const &x1,
                                       colvarvalue const &x2) const
{
  if (this->type() == colvarvalue::type_scalar) {
    return cvm::real (x1 - x2);
  } else {
    cvm::fatal_error ("Error: you requested an operation which requires "
                      "comparison between two non-scalar values.\n");
    return 0.0;
  }
}



class colvar::distance
  : public colvar::cvc
{
protected:
  cvm::atom_group  group1;
  cvm::atom_group  group2;
  cvm::rvector     dist_v;
  bool b_1site_force;
public:
  distance (std::string const &conf);
  distance();
  virtual inline ~distance() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


// // \brief Colvar component: distance vector between centers of mass
// // of two groups (\link colvarvalue::type_vector \endlink type,
// // range (-*:*)x(-*:*)x(-*:*))
// class colvar::distance_vec
//   : public colvar::distance
// {
// public:
//   distance_vec (std::string const &conf);
//   distance_vec();
//   virtual inline ~distance_vec() {}
//   virtual void calc_value();
//   virtual void calc_gradients();
//   virtual void apply_force (colvarvalue const &force);
//   /// Redefined to handle the box periodicity
//   virtual cvm::real dist2 (colvarvalue const &x1,
//                            colvarvalue const &x2) const;
//   /// Redefined to handle the box periodicity
//   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
//                                    colvarvalue const &x2) const;
//   /// Redefined to handle the box periodicity
//   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
//                                    colvarvalue const &x2) const;
//   /// Redefined to handle the box periodicity
//   virtual cvm::real compare (colvarvalue const &x1,
//                              colvarvalue const &x2) const;
// };


class colvar::distance_dir
  : public colvar::distance
{
public:
  distance_dir (std::string const &conf);
  distance_dir();
  virtual inline ~distance_dir() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
//   virtual cvm::real dist2 (colvarvalue const &x1,
//                            colvarvalue const &x2) const;
//   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
//                                    colvarvalue const &x2) const;
//   virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
//                                    colvarvalue const &x2) const;
//   virtual cvm::real compare (colvarvalue const &x1,
//                              colvarvalue const &x2) const;
};


class colvar::distance_z
  : public colvar::cvc
{
protected:
  cvm::atom_group  main;
  cvm::atom_group  ref1;
  cvm::atom_group  ref2;
  bool b_1site_force;
  cvm::rvector axis;
  cvm::real axis_norm;
  cvm::rvector     dist_v;
  bool fixed_axis;
public:
  distance_z (std::string const &conf);
  distance_z();
  virtual inline ~distance_z() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::distance_xy
  : public colvar::distance_z
{
protected:
  cvm::rvector dist_v_ortho;
  cvm::rvector v12, v13;
public:
  distance_xy (std::string const &conf);
  distance_xy();
  virtual inline ~distance_xy() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::min_distance
  : public colvar::distance
{
protected:
  cvm::real smoothing;
public:
  min_distance (std::string const &conf);
  min_distance();
  virtual inline ~min_distance() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};




class colvar::gyration
  : public colvar::cvc
{
protected:
  cvm::atom_group atoms;
public:
  gyration (std::string const &conf);
  gyration();
  virtual inline ~gyration() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::eigenvector
  : public colvar::cvc
{
protected:

  cvm::atom_group             atoms;

  std::vector<cvm::atom_pos>  ref_pos;

  std::vector<cvm::rvector>   eigenvec;

  cvm::real                   eigenvec_invnorm2;

public:

  eigenvector (std::string const &conf);
  virtual inline ~eigenvector() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};



class colvar::angle
  : public colvar::cvc
{
protected:

  cvm::atom_group group1;
  cvm::atom_group group2;
  cvm::atom_group group3;

  cvm::rvector r21, r23;
  cvm::real r21l, r23l;

public:

  angle (std::string const &conf);
  angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
  angle();
  virtual inline ~angle() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::dihedral
  : public colvar::cvc
{
protected:

  cvm::atom_group group1;
  cvm::atom_group group2;
  cvm::atom_group group3;
  cvm::atom_group group4;
  cvm::rvector r12, r23, r34;

  bool b_1site_force;
  
public:

  dihedral (std::string  const &conf);
  dihedral (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4);
  dihedral();
  virtual inline ~dihedral() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);

  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::coordnum
  : public colvar::distance
{
protected:
  cvm::real     r0;
  cvm::rvector  r0_vec;
  bool b_anisotropic;
  int en;
  int ed;
  bool b_group2_center_only;
public:
  coordnum (std::string const &conf);
  coordnum();
  virtual inline ~coordnum() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
  template<bool b_gradients>
  static cvm::real switching_function (cvm::real const &r0,
                                       int const &exp_num, int const &exp_den,
                                       cvm::atom &A1, cvm::atom &A2);
  
  template<bool b_gradients>
  static cvm::real switching_function (cvm::rvector const &r0_vec,
                                       int const &exp_num, int const &exp_den,
                                       cvm::atom &A1, cvm::atom &A2);

  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::h_bond 
  : public colvar::cvc
{
protected:
  cvm::atom     acceptor, donor;
  cvm::real     r0;
  int en;
  int ed;
public:
  h_bond (std::string const &conf);
  h_bond (cvm::atom const &acceptor,
          cvm::atom const &donor,
          cvm::real r0, int en, int ed);
  h_bond();
  virtual inline ~h_bond() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);

  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::alpha_dihedrals
  : public colvar::cvc
{
protected:

  cvm::real phi_ref;

  cvm::real psi_ref;

  std::vector<dihedral *> phi;

  std::vector<dihedral *> psi;

  std::vector<h_bond *>   hb;

public:

  alpha_dihedrals (std::string const &conf);
  alpha_dihedrals();
  virtual inline ~alpha_dihedrals() {}
  virtual void calc_value();
  virtual void calc_gradients(); 
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::alpha_angles
  : public colvar::cvc
{
protected:

  cvm::real theta_ref;

  cvm::real theta_tol;

  std::vector<angle *> theta;

  std::vector<h_bond *>   hb;

  cvm::real hb_coeff;

public:

  alpha_angles (std::string const &conf);
  alpha_angles();
  virtual inline ~alpha_angles() {}
  void calc_value();
  void calc_gradients();
  void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::orientation
  : public colvar::cvc
{
protected:

  cvm::atom_group            atoms;
  cvm::atom_pos              atoms_cog;

  std::vector<cvm::atom_pos> ref_pos;

  cvm::rotation              rot;

  cvm::quaternion            ref_quat;

public:

  orientation (std::string const &conf);
  orientation();
  virtual inline ~orientation() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::orientation_angle
  : public colvar::orientation
{
public:

  orientation_angle (std::string const &conf);
  orientation_angle();
  virtual inline ~orientation_angle() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};



class colvar::rmsd
  : public colvar::orientation
{
protected:

  cvm::real                  ref_pos_sum2;

public:

  rmsd (std::string const &conf);
  virtual inline ~rmsd() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


class colvar::logmsd
  : public colvar::orientation
{
protected:

  cvm::real                 ref_pos_sum2;
  cvm::real                 MSD;

public:

  logmsd (std::string const &conf);
  virtual inline ~logmsd() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void calc_force_invgrads();
  virtual void calc_Jacobian_derivative();
  virtual void apply_force (colvarvalue const &force);
  virtual cvm::real dist2 (colvarvalue const &x1,
                           colvarvalue const &x2) const;
  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
                                   colvarvalue const &x2) const;
  virtual cvm::real compare (colvarvalue const &x1,
                             colvarvalue const &x2) const;
};


// metrics functions for cvc implementations with a periodicity

inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1,
                                          colvarvalue const &x2) const
{
  // This version is asymptotically accurate for small distances,
  // other choices are possible
  cvm::real const diff = x1.real_value - x2.real_value;
  return (180.0*180.0/PI/PI) * 2.0 * (1.0 - ::cos (diff * (PI/180.00)));
}

inline colvarvalue colvar::dihedral::dist2_lgrad (colvarvalue const &x1,
                                                  colvarvalue const &x2) const
{
  cvm::real const diff = x1.real_value - x2.real_value;
  return colvarvalue ((180.0/PI) * 2.0 * ::sin ( diff * (PI/180.0)));
}

inline colvarvalue colvar::dihedral::dist2_rgrad (colvarvalue const &x1,
                                                  colvarvalue const &x2) const
{
  cvm::real const diff = x1.real_value - x2.real_value;
  return colvarvalue ((180.0/PI) * (-2.0) * ::sin ( diff * (PI/180.0)));
}

inline cvm::real colvar::dihedral::compare (colvarvalue const &x1,
                                            colvarvalue const &x2) const
{
  return dist2_lgrad (x1, x2);
}


// simple definitions of the distance functions; they are useful for
// optimizing the calculation (the type check performed by the base
// class functions is skipped)


// definitions assuming the scalar type

#define simple_scalar_dist_functions(TYPE)                              \
                                                                        \
  inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1,          \
                                        colvarvalue const &x2) const    \
  {                                                                     \
    return ::pow (x1.real_value - x2.real_value, int (2));              \
  }                                                                     \
                                                                        \
  inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1,  \
                                                colvarvalue const &x2) const \
  {                                                                     \
    return 2.0 * (x1.real_value - x2.real_value);                       \
  }                                                                     \
                                                                        \
  inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1,  \
                                                colvarvalue const &x2) const \
  {                                                                     \
    return this->dist2_lgrad (x2, x1);                                  \
  }                                                                     \
                                                                        \
  inline cvm::real colvar::TYPE::compare (colvarvalue const &x1,        \
                                          colvarvalue const &x2) const  \
  {                                                                     \
    return this->dist2_lgrad (x1, x2);                                  \
  }                                                                     \
                                                                        \

  simple_scalar_dist_functions (distance)
  simple_scalar_dist_functions (distance_xy)
  simple_scalar_dist_functions (distance_z)
  simple_scalar_dist_functions (min_distance)
  simple_scalar_dist_functions (angle)
  simple_scalar_dist_functions (coordnum)
  simple_scalar_dist_functions (h_bond)
  simple_scalar_dist_functions (gyration)
  simple_scalar_dist_functions (rmsd)
  simple_scalar_dist_functions (logmsd)
  simple_scalar_dist_functions (orientation_angle)
  simple_scalar_dist_functions (eigenvector)
  simple_scalar_dist_functions (alpha_dihedrals)
  simple_scalar_dist_functions (alpha_angles)



// generic definitions


// distance between three dimensional vectors

// inline cvm::real colvar::distance_vec::dist2 (colvarvalue const &x1,
//                                               colvarvalue const &x2) const
// {
//   return (x1.rvector_value - x2.rvector_value).norm2();
// }

// inline colvarvalue colvar::distance_vec::dist2_lgrad (colvarvalue const &x1,
//                                                       colvarvalue const &x2) const
// {
//   return (x1.rvector_value - x2.rvector_value);
// }

// inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1,
//                                                       colvarvalue const &x2) const
// {
//   return (x2.rvector_value - x1.rvector_value);
// }

// inline cvm::real colvar::distance_vec::compare (colvarvalue const &x1,
//                                                 colvarvalue const &x2) const
// {
//   cvm::fatal_error ("Error: cannot compare() two distance vectors.\n");
//   return 0.0;
// }


// distance between quaternions

inline cvm::real colvar::orientation::dist2 (colvarvalue const &x1,
                                             colvarvalue const &x2) const
{
  return x1.quaternion_value.dist2 (x2);
}

inline colvarvalue colvar::orientation::dist2_lgrad (colvarvalue const &x1,
                                                     colvarvalue const &x2) const
{
  return x1.quaternion_value.dist2_grad (x2);
}

inline colvarvalue colvar::orientation::dist2_rgrad (colvarvalue const &x1,
                                                     colvarvalue const &x2) const
{
  return x2.quaternion_value.dist2_grad (x1);
}

inline cvm::real colvar::orientation::compare (colvarvalue const &x1,
                                               colvarvalue const &x2) const
{
  cvm::fatal_error ("Error: cannot compare() two quaternions.\n");
  return 0.0;
}





#endif


// Emacs
// Local Variables:
// mode: C++
// End:

---