colvarcomp.h

// -*- c++ -*-
 
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/Colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
 
#ifndef COLVARCOMP_H
#define COLVARCOMP_H
 
// Declaration of colvar::cvc base class and derived ones.
//
// Future cvc's could be declared on additional header files.
// After the declaration of a new derived class, its metric
// functions must be reimplemented as well.
// If the new cvc has no symmetry or periodicity,
// this can be done straightforwardly by using the macro:
// simple_scalar_dist_functions (derived_class)
 
#include <memory>
 
#include "colvarmodule.h"
#include "colvaratoms.h"
#include "colvar.h"
#include "colvar_geometricpath.h"
 
 
 
class colvar::cvc
  : public colvarparse, public colvardeps
{
public:
 
  std::string name;
 
  std::string function_type;
 
  std::string config_key;
 
  cvm::real sup_coeff;
  int       sup_np;
 
  cvm::real period;
 
  cvm::real wrap_center;
 
  cvc(std::string const &conf);
 
  int init_code = COLVARS_OK;
 
  int set_function_type(std::string const &type);
 
  virtual int init(std::string const &conf);
 
  virtual int init_dependencies();
 
  cvm::atom_group *parse_group(std::string const &conf,
                               char const *group_key,
                               bool optional = false);
 
  virtual int init_total_force_params(std::string const &conf);
 
  int setup();
 
  cvc();
 
  virtual ~cvc();
 
  static std::vector<feature *> cvc_features;
 
  virtual const std::vector<feature *> &features() const
  {
    return cvc_features;
  }
  virtual std::vector<feature *> &modify_features()
  {
    return cvc_features;
  }
  static void delete_features() {
    for (size_t i=0; i < cvc_features.size(); i++) {
      delete cvc_features[i];
    }
    cvc_features.clear();
  }
 
  virtual std::vector<std::vector<int> > get_atom_lists();
 
  virtual void read_data();
 
  virtual void calc_value() = 0;
 
  virtual void calc_gradients() {}
 
  void calc_fit_gradients();
 
  virtual void debug_gradients();
 
  virtual void collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_gradients);
 
  virtual void calc_force_invgrads();
 
  virtual void calc_Jacobian_derivative();
 
 
  colvarvalue const & value() const;
 
  colvarvalue const & total_force() const;
 
  colvarvalue const & 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 void wrap(colvarvalue &x_unwrapped) const;
 
  std::vector<cvm::atom_group *> atom_groups;
 
  void register_atom_group(cvm::atom_group *ag);
 
  virtual colvarvalue const *get_param_grad(std::string const &param_name);
 
  virtual int set_param(std::string const &param_name, void const *new_value);
 
  bool b_try_scalable;
 
  inline void set_value(colvarvalue const &new_value) {
    x = new_value;
  }
 
protected:
 
  int update_description();
 
  std::vector<std::string> function_types;
 
  colvarvalue x;
 
  colvarvalue x_old;
 
  colvarvalue ft;
 
  colvarvalue jd;
 
  void init_as_distance();
 
  void init_as_angle();
 
  void init_as_periodic_angle();
 
  void init_scalar_boundaries(cvm::real lb, cvm::real ub);
 
  colvarvalue lower_boundary;
 
  colvarvalue upper_boundary;
 
  cvm::real width;
};
 
 
inline colvarvalue const & colvar::cvc::value() const
{
  return x;
}
 
 
inline colvarvalue const & colvar::cvc::total_force() const
{
  return ft;
}
 
 
inline colvarvalue const & colvar::cvc::Jacobian_derivative() const
{
  return jd;
}
 
 
 
 
class colvar::distance
  : public colvar::cvc
{
protected:
  cvm::atom_group  *group1;
  cvm::atom_group  *group2;
  cvm::rvector     dist_v;
public:
  distance(std::string const &conf);
  distance();
  virtual ~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;
};
 
 
 
// \brief Colvar component: distance vector between centers of mass
// of two groups (\link colvarvalue::type_3vector \endlink type,
// range (-*:*)x(-*:*)x(-*:*))
class colvar::distance_vec
  : public colvar::distance
{
public:
  distance_vec(std::string const &conf);
  distance_vec();
  virtual ~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;
};
 
 
 
class colvar::distance_dir
  : public colvar::distance
{
public:
  distance_dir(std::string const &conf);
  distance_dir();
  virtual ~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;
};
 
 
 
class colvar::distance_z
  : public colvar::cvc
{
protected:
  cvm::atom_group  *main;
  cvm::atom_group  *ref1;
  cvm::atom_group  *ref2;
  cvm::rvector axis;
  cvm::real axis_norm;
  cvm::rvector     dist_v;
  bool fixed_axis;
public:
  distance_z(std::string const &conf);
  distance_z();
  virtual ~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 void wrap(colvarvalue &x_unwrapped) 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 ~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;
};
 
 
class colvar::polar_phi
  : public colvar::cvc
{
public:
  polar_phi(std::string const &conf);
  polar_phi();
  virtual ~polar_phi() {}
protected:
  cvm::atom_group  *atoms;
  cvm::real r, theta, phi;
public:
  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 void wrap(colvarvalue &x_unwrapped) const;
};
 
 
class colvar::polar_theta
  : public colvar::cvc
{
public:
  polar_theta(std::string const &conf);
  polar_theta();
  virtual ~polar_theta() {}
protected:
  cvm::atom_group  *atoms;
  cvm::real r, theta, phi;
public:
  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;
};
 
class colvar::distance_inv
  : public colvar::cvc
{
protected:
  cvm::atom_group  *group1;
  cvm::atom_group  *group2;
  int exponent;
public:
  distance_inv(std::string const &conf);
  virtual ~distance_inv() {}
  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;
};
 
 
 
class colvar::distance_pairs
  : public colvar::cvc
{
protected:
  cvm::atom_group  *group1;
  cvm::atom_group  *group2;
public:
  distance_pairs(std::string const &conf);
  distance_pairs();
  virtual ~distance_pairs() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
};
 
 
 
class colvar::dipole_magnitude
  : public colvar::cvc
{
protected:
  cvm::atom_group  *atoms;
  cvm::atom_pos dipoleV;
public:
  dipole_magnitude (std::string const &conf);
  dipole_magnitude (cvm::atom const &a1);
  dipole_magnitude();
  virtual ~dipole_magnitude() {}
  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;
};
 
 
 
class colvar::gyration
  : public colvar::cvc
{
protected:
  cvm::atom_group  *atoms;
public:
  gyration(std::string const &conf);
  virtual ~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;
};
 
 
 
class colvar::inertia
  : public colvar::gyration
{
public:
  inertia(std::string const &conf);
  inertia();
  virtual ~inertia() {}
  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;
};
 
 
 
class colvar::inertia_z
  : public colvar::inertia
{
protected:
  cvm::rvector axis;
public:
  inertia_z(std::string const &conf);
  inertia_z();
  virtual ~inertia_z() {}
  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;
};
 
 
 
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 ~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;
};
 
 
 
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);
  virtual ~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;
};
 
 
 
class colvar::dipole_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:
 
  dipole_angle (std::string const &conf);
  dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
  dipole_angle();
  virtual ~dipole_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;
};
 
 
 
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 ~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 void wrap(colvarvalue &x_unwrapped) const;
};
 
 
 
class colvar::coordnum
  : public colvar::cvc
{
protected:
  cvm::atom_group  *group1;
  cvm::atom_group  *group2;
  cvm::real     r0;
  cvm::rvector  r0_vec;
  bool b_anisotropic;
  int en;
  int ed;
 
  bool b_group2_center_only;
 
  cvm::real tolerance;
 
  int pairlist_freq;
 
  bool *pairlist;
 
public:
 
  coordnum(std::string const &conf);
  ~coordnum();
 
  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;
 
  enum {
    ef_null = 0,
    ef_gradients = 1,
    ef_anisotropic = (1<<8),
    ef_use_pairlist = (1<<9),
    ef_rebuild_pairlist = (1<<10)
  };
 
  template<int flags>
  static cvm::real switching_function(cvm::real const &r0,
                                      cvm::rvector const &r0_vec,
                                      int en,
                                      int ed,
                                      cvm::atom &A1,
                                      cvm::atom &A2,
                                      bool **pairlist_elem,
                                      cvm::real tolerance);
 
  template<int flags> int compute_coordnum();
 
  template<int flags> void main_loop(bool **pairlist_elem);
 
};
 
 
 
class colvar::selfcoordnum
  : public colvar::cvc
{
protected:
 
  cvm::atom_group  *group1;
  cvm::real     r0;
  int en;
  int ed;
  cvm::real tolerance;
  int pairlist_freq;
  bool *pairlist;
 
public:
 
  selfcoordnum(std::string const &conf);
  ~selfcoordnum();
  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;
 
  template<int flags> int compute_selfcoordnum();
};
 
 
 
class colvar::groupcoordnum
  : public colvar::distance
{
protected:
  cvm::real     r0;
  cvm::rvector  r0_vec;
  bool b_anisotropic;
  int en;
  int ed;
public:
  groupcoordnum(std::string const &conf);
  virtual ~groupcoordnum() {}
  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;
};
 
 
 
class colvar::h_bond
  : public colvar::cvc
{
protected:
  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;
};
 
 
 
// 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 ~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;
// };
 
 
 
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 collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_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;
};
 
 
 
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 collect_gradients(std::vector<int> const &atom_ids, std::vector<cvm::rvector> &atomic_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;
};
 
 
 
class colvar::orientation
  : public colvar::cvc
{
protected:
 
  cvm::atom_group  *          atoms;
  cvm::atom_pos              atoms_cog;
 
  std::vector<cvm::atom_pos> ref_pos;
 
  std::vector<cvm::atom_pos> shifted_pos;
 
  cvm::rotation              rot;
 
  cvm::quaternion            ref_quat;
 
  struct rotation_derivative_impl_;
  std::unique_ptr<rotation_derivative_impl_> rot_deriv_impl;
 
public:
 
  orientation(std::string const &conf);
  orientation();
  virtual int init(std::string const &conf);
  virtual ~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;
};
 
 
 
class colvar::orientation_angle
  : public colvar::orientation
{
public:
 
  orientation_angle(std::string const &conf);
  virtual int init(std::string const &conf);
  virtual ~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;
};
 
 
 
class colvar::orientation_proj
  : public colvar::orientation
{
public:
 
  orientation_proj(std::string const &conf);
  orientation_proj();
  virtual int init(std::string const &conf);
  virtual ~orientation_proj() {}
  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;
};
 
 
 
class colvar::tilt
  : public colvar::orientation
{
protected:
 
  cvm::rvector axis;
 
public:
 
  tilt(std::string const &conf);
  virtual int init(std::string const &conf);
  virtual ~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;
};
 
 
 
class colvar::spin_angle
  : public colvar::orientation
{
protected:
 
  cvm::rvector axis;
 
public:
 
  spin_angle(std::string const &conf);
  spin_angle();
  virtual int init(std::string const &conf);
  virtual ~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 void wrap(colvarvalue &x_unwrapped) const;
};
 
 
class colvar::euler_phi
  : public colvar::orientation
{
public:
  euler_phi(std::string const &conf);
  euler_phi();
  virtual int init(std::string const &conf);
  virtual ~euler_phi() {}
  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 void wrap(colvarvalue &x_unwrapped) const;
};
 
 
class colvar::euler_psi
  : public colvar::orientation
{
public:
  euler_psi(std::string const &conf);
  euler_psi();
  virtual int init(std::string const &conf);
  virtual ~euler_psi() {}
  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 void wrap(colvarvalue &x_unwrapped) const;
};
 
 
class colvar::euler_theta
  : public colvar::orientation
{
public:
  euler_theta(std::string const &conf);
  euler_theta();
  virtual int init(std::string const &conf);
  virtual ~euler_theta() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
  // theta angle is a scalar variable and not periodic
  // we need to override the virtual functions from orientation
  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;
};
 
 
class colvar::rmsd
  : public colvar::cvc
{
protected:
 
  cvm::atom_group  *atoms;
 
  std::vector<cvm::atom_pos>  ref_pos;
 
  size_t n_permutations;
 
  size_t best_perm_index;
public:
 
  rmsd(std::string const &conf);
  virtual ~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;
};
 
 
 
// \brief Colvar component: flat vector of Cartesian coordinates
// Mostly useful to compute scripted colvar values
class colvar::cartesian
  : public colvar::cvc
{
protected:
  cvm::atom_group  *atoms;
  std::vector<size_t> axes;
public:
  cartesian(std::string const &conf);
  cartesian();
  virtual ~cartesian() {}
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
};
 
 
// \brief Colvar component: alch_lambda
// To communicate value with back-end in lambda-dynamics
class colvar::alch_lambda
  : public colvar::cvc
{
protected:
  // No atom groups needed
public:
  alch_lambda(std::string const &conf);
  alch_lambda();
  virtual ~alch_lambda() {}
  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;
};
 
 
// \brief Colvar component: alch_Flambda
// To communicate force on lambda with back-end in lambda-dynamics
class colvar::alch_Flambda
  : public colvar::cvc
{
protected:
  // No atom groups needed
public:
  alch_Flambda(std::string const &conf);
  alch_Flambda();
  virtual ~alch_Flambda() {}
  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;
};
 
 
class colvar::componentDisabled
  : public colvar::cvc
{
public:
    componentDisabled(std::string const & /* conf */) {
        cvm::error("Error: this component is not enabled in the current build; please see https://colvars.github.io/README-c++11.html");
    }
    virtual ~componentDisabled() {}
    virtual void calc_value() {}
    virtual void calc_gradients() {}
    virtual void apply_force(colvarvalue const & /* force */) {}
};
 
 
 
class colvar::CartesianBasedPath
  : public colvar::cvc
{
protected:
    virtual void computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& result);
    virtual void computeDistanceToReferenceFrames(std::vector<cvm::real>& result);
    cvm::atom_group *atoms;
    bool has_user_defined_fitting;
    std::vector<std::vector<cvm::atom_pos>> reference_frames;
    std::vector<std::vector<cvm::atom_pos>> reference_fitting_frames;
    std::vector<cvm::atom_group*> comp_atoms;
    size_t total_reference_frames;
public:
    CartesianBasedPath(std::string const &conf);
    virtual ~CartesianBasedPath();
    virtual void calc_value() = 0;
    virtual void apply_force(colvarvalue const &force) = 0;
};
 
class colvar::gspath
  : public colvar::CartesianBasedPath, public GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::S>
{
private:
    // Optimal rotation for compute v3
    cvm::rotation rot_v3;
protected:
    virtual void prepareVectors();
    virtual void updateDistanceToReferenceFrames();
public:
    gspath(std::string const &conf);
    virtual ~gspath() {}
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
 
 
class colvar::gzpath
  : public colvar::CartesianBasedPath, public GeometricPathCV::GeometricPathBase<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::Z>
{
private:
    // Optimal rotation for compute v3, v4
    cvm::rotation rot_v3;
    cvm::rotation rot_v4;
protected:
    virtual void prepareVectors();
    virtual void updateDistanceToReferenceFrames();
public:
    gzpath(std::string const &conf);
    virtual ~gzpath() {}
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
class colvar::linearCombination
  : public colvar::cvc
{
protected:
    std::vector<colvar::cvc*> cv;
    bool use_explicit_gradients;
protected:
    cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const;
public:
    linearCombination(std::string const &conf);
    virtual ~linearCombination();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
 
class colvar::customColvar
  : public colvar::linearCombination
{
protected:
    bool use_custom_function;
#ifdef LEPTON
    std::vector<Lepton::CompiledExpression *> value_evaluators;
    std::vector<Lepton::CompiledExpression *> gradient_evaluators;
    std::vector<double *> value_eval_var_refs;
    std::vector<double *> grad_eval_var_refs;
    double dev_null;
#endif
public:
    customColvar(std::string const &conf);
    virtual ~customColvar();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
 
class colvar::CVBasedPath
  : public colvar::cvc
{
protected:
    std::vector<colvar::cvc*> cv;
    std::vector<std::vector<colvarvalue>> ref_cv;
    bool use_explicit_gradients;
    size_t total_reference_frames;
protected:
    virtual void computeDistanceToReferenceFrames(std::vector<cvm::real>& result);
    virtual void computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& result) const;
    cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const;
public:
    CVBasedPath(std::string const &conf);
    virtual ~CVBasedPath();
    virtual void calc_value() = 0;
    virtual void apply_force(colvarvalue const &force) = 0;
};
 
 
class colvar::gspathCV
  : public colvar::CVBasedPath, public GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::S>
{
protected:
    virtual void updateDistanceToReferenceFrames();
    virtual void prepareVectors();
public:
    gspathCV(std::string const &conf);
    virtual ~gspathCV();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
 
 
class colvar::gzpathCV
  : public colvar::CVBasedPath, public GeometricPathCV::GeometricPathBase<colvarvalue, cvm::real, GeometricPathCV::path_sz::Z>
{
protected:
    virtual void updateDistanceToReferenceFrames();
    virtual void prepareVectors();
public:
    gzpathCV(std::string const &conf);
    virtual ~gzpathCV();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
struct ArithmeticPathImpl;
 
class colvar::aspath
  : public colvar::CartesianBasedPath
{
private:
    std::unique_ptr<ArithmeticPathImpl> impl_;
    friend struct ArithmeticPathImpl;
public:
    aspath(std::string const &conf);
    virtual ~aspath();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
class colvar::azpath
  : public colvar::CartesianBasedPath
{
private:
    std::unique_ptr<ArithmeticPathImpl> impl_;
    friend struct ArithmeticPathImpl;
public:
    azpath(std::string const &conf);
    virtual ~azpath();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
class colvar::aspathCV
  : public colvar::CVBasedPath
{
private:
    std::unique_ptr<ArithmeticPathImpl> impl_;
    friend struct ArithmeticPathImpl;
public:
    aspathCV(std::string const &conf);
    virtual ~aspathCV();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
 
class colvar::azpathCV
  : public colvar::CVBasedPath
{
private:
    std::unique_ptr<ArithmeticPathImpl> impl_;
    friend struct ArithmeticPathImpl;
public:
    azpathCV(std::string const &conf);
    virtual ~azpathCV();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
// forward declaration
namespace neuralnetworkCV {
    class neuralNetworkCompute;
}
 
 
class colvar::neuralNetwork
  : public linearCombination
{
protected:
    std::unique_ptr<neuralnetworkCV::neuralNetworkCompute> nn;
    size_t m_output_index;
public:
    neuralNetwork(std::string const &conf);
    virtual ~neuralNetwork();
    virtual void calc_value();
    virtual void calc_gradients();
    virtual void apply_force(colvarvalue const &force);
};
 
 
// \brief Colvar component: total value of a scalar map
// (usually implemented as a grid by the simulation engine)
class colvar::map_total
  : public colvar::cvc
{
public:
 
  map_total();
  map_total(std::string const &conf);
  virtual ~map_total() {}
  virtual int init(std::string const &conf);
  virtual void calc_value();
  virtual void calc_gradients();
  virtual void apply_force(colvarvalue const &force);
 
protected:
 
  std::string volmap_name;
 
  int volmap_id;
 
  int volmap_index;
 
  cvm::atom_group *atoms;
 
  std::vector<cvm::real> atom_weights;
};
 
 
 
// metrics functions for cvc implementations
 
// 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)                              \
                                                                        \
                                                                        \
  cvm::real colvar::TYPE::dist2(colvarvalue const &x1,                  \
                                colvarvalue const &x2) const            \
  {                                                                     \
    return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
  }                                                                     \
                                                                        \
                                                                        \
  colvarvalue colvar::TYPE::dist2_lgrad(colvarvalue const &x1,          \
                                        colvarvalue const &x2) const    \
  {                                                                     \
    return 2.0 * (x1.real_value - x2.real_value);                       \
  }                                                                     \
                                                                        \
                                                                        \
  colvarvalue colvar::TYPE::dist2_rgrad(colvarvalue const &x1,          \
                                        colvarvalue const &x2) const    \
  {                                                                     \
    return this->dist2_lgrad(x2, x1);                                   \
  }                                                                     \
 
 
#endif