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;

  cvm::real wrap_center;

  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;

  virtual void wrap (colvarvalue &x) const;

  std::vector<cvm::atom_group *> atom_groups;

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;
  }
}

inline void colvar::cvc::wrap (colvarvalue &x) const
{
  return;
}


class colvar::distance
  : public colvar::cvc
{
protected:
  cvm::atom_group  group1;
  cvm::atom_group  group2;
  cvm::rvector     dist_v;
  bool b_no_PBC;
  bool b_1site_force;
public:
  distance (std::string const &conf, bool twogroups = true);
  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);
  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_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_no_PBC;
  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;
  virtual void wrap (colvarvalue &x) 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;
  cvm::rvector dxdr1, dxdr3;

  bool b_1site_force;
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 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::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;
  virtual void wrap (colvarvalue &x) 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::selfcoordnum
  : public colvar::distance
{
protected:
  cvm::real     r0;
  int en;
  int ed;
public:
  selfcoordnum (std::string const &conf);
  selfcoordnum();
  virtual inline ~selfcoordnum() {}
  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);
  
  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 ~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:

//   /// Alpha-helical reference phi value
//   cvm::real phi_ref;

//   /// Alpha-helical reference psi value
//   cvm::real psi_ref;

//   /// List of phi dihedral angles
//   std::vector<dihedral *> phi;

//   /// List of psi dihedral angles
//   std::vector<dihedral *> psi;

//   /// List of hydrogen bonds
//   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 ~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::dihedPC
  : public colvar::cvc
{
protected:

  std::vector<dihedral *> theta;
  std::vector<cvm::real> coeffs;

public:

  dihedPC (std::string const &conf);
  dihedPC();
  virtual  ~dihedPC();
  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::tilt
  : public colvar::orientation
{
protected:

  cvm::rvector axis;

public:

  tilt (std::string const &conf);
  tilt();
  virtual inline ~tilt() {}
  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::spin_angle
  : public colvar::orientation
{
protected:

  cvm::rvector axis;

public:

  spin_angle (std::string const &conf);
  spin_angle();
  virtual inline ~spin_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;
  virtual void wrap (colvarvalue &x) 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
{
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return diff * diff;
}

inline colvarvalue colvar::dihedral::dist2_lgrad (colvarvalue const &x1,
                                                  colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return 2.0 * diff;
}

inline colvarvalue colvar::dihedral::dist2_rgrad (colvarvalue const &x1,
                                                  colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return (-2.0) * diff;
}

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

inline void colvar::dihedral::wrap (colvarvalue &x) const
{
  if ((x.real_value - wrap_center) >= 180.0) {
    x.real_value -= 360.0;
    return;
  } 

  if ((x.real_value - wrap_center) < -180.0) {
    x.real_value += 360.0;
    return;
  } 

  return;
}

inline cvm::real colvar::spin_angle::dist2 (colvarvalue const &x1,
                                          colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return diff * diff;
}

inline colvarvalue colvar::spin_angle::dist2_lgrad (colvarvalue const &x1,
                                                  colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return 2.0 * diff;
}

inline colvarvalue colvar::spin_angle::dist2_rgrad (colvarvalue const &x1,
                                                  colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff));
  return (-2.0) * diff;
}

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

inline void colvar::spin_angle::wrap (colvarvalue &x) const
{
  if ((x.real_value - wrap_center) >= 180.0) {
    x.real_value -= 360.0;
    return;
  } 

  if ((x.real_value - wrap_center) < -180.0) {
    x.real_value += 360.0;
    return;
  } 

  return;
}

// simple definitions of the distance functions; these are useful only
// for optimization (the type check performed in the default
// colvarcomp 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 std::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)
  // NOTE: distance_z has explicit functions, see below 
  simple_scalar_dist_functions (distance_xy)
  simple_scalar_dist_functions (min_distance)
  simple_scalar_dist_functions (angle)
  simple_scalar_dist_functions (coordnum)
  simple_scalar_dist_functions (selfcoordnum)
  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 (tilt)
  simple_scalar_dist_functions (eigenvector)
  //  simple_scalar_dist_functions (alpha_dihedrals)
  simple_scalar_dist_functions (alpha_angles)
  simple_scalar_dist_functions (dihedPC)


// Projected distance
// Differences should always be wrapped around 0 (ignoring wrap_center)
inline cvm::real colvar::distance_z::dist2 (colvarvalue const &x1,
                                            colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  if (period != 0.0) {
    cvm::real shift = std::floor (diff/period + 0.5);
    diff -= shift * period;
  }
  return diff * diff;
}

inline colvarvalue colvar::distance_z::dist2_lgrad (colvarvalue const &x1,
                                                    colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  if (period != 0.0) {
    cvm::real shift = std::floor (diff/period + 0.5);
    diff -= shift * period;
  }
  return 2.0 * diff;
}

inline colvarvalue colvar::distance_z::dist2_rgrad (colvarvalue const &x1,
                                                    colvarvalue const &x2) const
{
  cvm::real diff = x1.real_value - x2.real_value;
  if (period != 0.0) {
    cvm::real shift = std::floor (diff/period + 0.5);
    diff -= shift * period;
  }
  return (-2.0) * diff;
}

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

inline void colvar::distance_z::wrap (colvarvalue &x) const
{
  if (! this->b_periodic) {
    // don't wrap if the period has not been set
    return;
  }

  cvm::real shift = std::floor ((x.real_value - wrap_center) / period + 0.5);
  x.real_value -= shift * period;
  return;
}


// distance between three dimensional vectors
//
// TODO apply PBC to distance_vec
// Note: differences should be centered around (0, 0, 0)!

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

inline colvarvalue colvar::distance_vec::dist2_lgrad (colvarvalue const &x1,
                                                      colvarvalue const &x2) const
{
  return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
}

inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1,
                                                      colvarvalue const &x2) const
{
  return 2.0 * cvm::position_distance(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;
}

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

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

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

inline cvm::real colvar::distance_dir::compare (colvarvalue const &x1,
                                                colvarvalue const &x2) const
{
  cvm::fatal_error ("Error: cannot compare() two distance directions.\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:

---