colvar.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 COLVAR_H
#define COLVAR_H
 
#include <functional>
#include <list>
#include <iosfwd>
#include <map>
#include <memory>
 
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvardeps.h"
 
#ifdef LEPTON
#include "Lepton.h" // for runtime custom expressions
#endif
 
 
class colvar : public colvarparse, public colvardeps {
 
public:
 
  std::string name;
 
  colvarvalue const & value() const;
 
  colvarvalue const & actual_value() const;
 
  colvarvalue const & run_ave() const;
 
  cvm::real const & force_constant() const;
 
  colvarvalue const & velocity() const;
 
  colvarvalue const & total_force() const;
 
  cvm::real width = 1.0;
 
  static std::vector<feature *> cv_features;
 
  virtual const std::vector<feature *> &features() const
  {
    return cv_features;
  }
  virtual std::vector<feature *> &modify_features()
  {
    return cv_features;
  }
  static void delete_features() {
    for (size_t i=0; i < cv_features.size(); i++) {
      delete cv_features[i];
    }
    cv_features.clear();
  }
 
  void do_feature_side_effects(int id);
 
  std::vector<colvarbias *> biases;
 
protected:
 
 
  /*
    extended:
    H = H_{0} + \sum_{i} 1/2*\lambda*(S_i(x(t))-s_i(t))^2 \\
    + \sum_{i} 1/2*m_i*(ds_i(t)/dt)^2 \\
    + \sum_{t'<t} W * exp(-1/2*\sum_{i} (s_i(t')-s_i(t))^2/(\delta{}s_i)^2) \\
    + \sum_{w} (\sum_{i}c_{w,i}s_i(t) - D_w)^M/(\Sigma_w)^M
 
    normal:
    H = H_{0} + \sum_{t'<t} W * exp(-1/2*\sum_{i} (S_i(x(t'))-S_i(x(t)))^2/(\delta{}S_i)^2) \\
    + \sum_{w} (\sum_{i}c_{w,i}S_i(t) - D_w)^M/(\Sigma_w)^M
 
    output:
    H = H_{0}   (only output S(x), no forces)
 
    Here:
    S(x(t)) = x
    s(t)    = x_ext
    DS = Ds = delta
  */
 
 
  colvarvalue x;
 
  // TODO: implement functionality to treat these
  // /// Vector of individual values from CVCs
  // colvarvalue x_cvc;
 
  // /// Jacobian matrix of individual values from CVCs
  // colvarvalue dx_cvc;
 
  colvarvalue x_reported;
 
  colvarvalue v_fdiff;
 
  inline colvarvalue fdiff_velocity(colvarvalue const &xold,
                                     colvarvalue const &xnew)
  {
    // using the gradient of the square distance to calculate the
    // velocity (non-scalar variables automatically taken into
    // account)
    cvm::real dt = cvm::dt();
    return ( ( (dt > 0.0) ? (1.0/dt) : 1.0 ) *
             0.5 * dist2_lgrad(xnew, xold) );
  }
 
  colvarvalue v_reported;
 
  // Options for extended_lagrangian
  colvarvalue x_ext;
  colvarvalue prev_x_ext;
  colvarvalue v_ext;
  colvarvalue prev_v_ext;
  cvm::real ext_mass = 0.0;
  cvm::real ext_force_k = 0.0;
  cvm::real ext_gamma = 0.0;
  cvm::real ext_sigma = 0.0;
 
  colvarvalue fr;
 
  colvarvalue fj;
 
  colvarvalue ft_reported;
 
public:
 
 
  colvarvalue fb;
 
  colvarvalue fb_actual;
 
  colvarvalue f;
 
  colvarvalue f_old;
 
  colvarvalue ft;
 
  cvm::real period = 0.0;
 
  cvm::real wrap_center = 0.0;
 
  bool expand_boundaries = false;
 
  colvarvalue lower_boundary;
 
  colvarvalue upper_boundary;
 
  bool periodic_boundaries() const;
 
  bool periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) const;
 
 
  colvar();
 
  int init(std::string const &conf);
 
  void define_component_types();
 
  int init_components(std::string const &conf);
 
  int init_custom_function(std::string const &conf);
 
  int init_grid_parameters(std::string const &conf);
 
  int check_grid_parameters();
 
  int parse_legacy_wall_params(std::string const &conf);
 
  int init_extended_Lagrangian(std::string const &conf);
 
  int init_output_flags(std::string const &conf);
 
  virtual int init_dependencies();
 
private:
 
  template <typename def_class_name>
  void add_component_type(char const *description, char const *config_key);
 
  int init_components_type(const std::string &conf, const char *config_key);
 
public:
 
  void setup();
 
  ~colvar();
 
 
  int calc();
 
  int end_of_step();
 
  int calc_cvcs(int first = 0, size_t num_cvcs = 0);
 
  int check_cvc_range(int first_cvc, size_t num_cvcs);
 
  int calc_cvc_values(int first, size_t num_cvcs);
  int calc_cvc_gradients(int first, size_t num_cvcs);
  int calc_cvc_total_force(int first, size_t num_cvcs);
  int calc_cvc_Jacobians(int first, size_t num_cvcs);
 
  int collect_cvc_data();
 
  int collect_cvc_values();
  int collect_cvc_gradients();
  int collect_cvc_total_forces();
  int collect_cvc_Jacobians();
  int calc_colvar_properties();
 
  inline colvarvalue const applied_force() const
  {
    if (is_enabled(f_cv_extended_Lagrangian)) {
      return fr;
    }
    return f;
  }
 
  void reset_bias_force();
 
  void add_bias_force(colvarvalue const &force);
 
  void add_bias_force_actual_value(colvarvalue const &force);
 
  cvm::real update_forces_energy();
 
  void update_extended_Lagrangian();
 
  void communicate_forces();
 
  int set_cvc_flags(std::vector<bool> const &flags);
 
  int update_cvc_flags();
 
  int update_cvc_config(std::vector<std::string> const &confs);
 
  int cvc_param_exists(std::string const &param_name);
 
  cvm::real get_cvc_param(std::string const &param_name);
 
  void const *get_cvc_param_ptr(std::string const &param_name);
 
  colvarvalue const *get_cvc_param_grad(std::string const &param_name);
 
  int set_cvc_param(std::string const &param_name, void const *new_value);
 
protected:
 
  size_t n_active_cvcs = 0;
 
  cvm::real active_cvc_square_norm = 0.0;
 
  void update_active_cvc_square_norm();
 
  cvm::step_number prev_timestep;
 
public:
 
  inline size_t num_dimensions() const
  {
    return value().size();
  }
 
  inline size_t num_cvcs() const
  {
    return cvcs.size();
  }
 
  inline size_t num_active_cvcs() const
  {
    return n_active_cvcs;
  }
 
  cvm::real dist2(colvarvalue const &x1,
                  colvarvalue const &x2) const;
 
  colvarvalue dist2_lgrad(colvarvalue const &x1,
                          colvarvalue const &x2) const;
 
  colvarvalue dist2_rgrad(colvarvalue const &x1,
                          colvarvalue const &x2) const;
 
  void wrap(colvarvalue &x_unwrapped) const;
 
  int parse_analysis(std::string const &conf);
 
  int analyze();
 
  std::istream & read_traj(std::istream &is);
 
  std::ostream & write_traj(std::ostream &os);
  std::ostream & write_traj_label(std::ostream &os);
 
  std::istream & read_state(std::istream &is);
 
  cvm::memory_stream & read_state(cvm::memory_stream &is);
 
  int check_matching_state(std::string const &state_conf);
 
  int set_state_params(std::string const &state_conf);
 
  std::string const get_state_params() const;
 
  std::ostream & write_state(std::ostream &os) const;
 
  cvm::memory_stream & write_state(cvm::memory_stream &os) const;
 
  int write_output_files();
 
protected:
 
  bool matching_state;
 
  colvarvalue            x_old;
 
  colvarvalue            x_restart;
 
  bool                   after_restart;
 
  std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
  std::list< std::list<colvarvalue> >::iterator acf_x_history_p, acf_v_history_p;
 
  std::list< std::list<colvarvalue> > x_history;
  std::list< std::list<colvarvalue> >::iterator x_history_p;
 
  std::string            acf_colvar_name;
  size_t                 acf_length;
  size_t                 acf_offset;
  size_t                 acf_stride;
  size_t                 acf_nframes;
  bool                   acf_normalize;
  std::vector<cvm::real> acf;
  std::string            acf_outfile;
 
  enum acf_type_e {
    acf_notset,
    acf_vel,
    acf_coor,
    acf_p2coor
  };
 
  acf_type_e             acf_type;
 
  void calc_vel_acf(std::list<colvarvalue> &v_history,
                    colvarvalue const      &v);
 
  void calc_coor_acf(std::list<colvarvalue> &x_history,
                     colvarvalue const      &x);
 
  void calc_p2coor_acf(std::list<colvarvalue> &x_history,
                       colvarvalue const      &x);
 
  int calc_acf();
  int write_acf(std::ostream &os);
 
  size_t         runave_length;
  size_t         runave_stride;
  std::string    runave_outfile;
  colvarvalue    runave;
  cvm::real      runave_variance = 0.0;
 
  int calc_runave();
 
  cvm::real kinetic_energy = 0.0;
  cvm::real potential_energy = 0.0;
 
public:
 
  // collective variable component base class
  class cvc;
 
  // list of available collective variable components
 
  // scalar colvar components
  class distance;
  class distance_z;
  class distance_xy;
  class polar_theta;
  class polar_phi;
  class distance_inv;
  class distance_pairs;
  class dipole_magnitude;
  class angle;
  class dipole_angle;
  class dihedral;
  class coordnum;
  class selfcoordnum;
  class groupcoordnum;
  class h_bond;
  class rmsd;
  class orientation_angle;
  class orientation_proj;
  class tilt;
  class spin_angle;
  class gyration;
  class inertia;
  class inertia_z;
  class eigenvector;
  class alpha_dihedrals;
  class alpha_angles;
  class dihedPC;
  class alch_lambda;
  class alch_Flambda;
  class CartesianBasedPath;
  class aspath;
  class azpath;
  class gspath;
  class gzpath;
  class linearCombination;
  class CVBasedPath;
  class gspathCV;
  class gzpathCV;
  class aspathCV;
  class azpathCV;
  class euler_phi;
  class euler_psi;
  class euler_theta;
  class neuralNetwork;
  class torchANN;
  class customColvar;
 
  // non-scalar components
  class distance_vec;
  class distance_dir;
  class cartesian;
  class orientation;
 
  // components that do not handle any atoms directly
  class map_total;
 
  static const std::map<std::string, std::function<colvar::cvc *()>> &get_global_cvc_map()
  {
    return global_cvc_map;
  }
 
  static bool compare_cvc(const colvar::cvc* const i, const colvar::cvc* const j);
 
protected:
 
  std::vector<std::shared_ptr<colvar::cvc>> cvcs;
 
  std::vector<bool> cvc_flags;
 
  void build_atom_list(void);
 
  std::string scripted_function;
 
  std::vector<const colvarvalue *> sorted_cvc_values;
 
#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
 
  static std::map<std::string, std::function<colvar::cvc *()>> global_cvc_map;
 
  static std::map<std::string, std::string> global_cvc_desc_map;
 
  std::vector<int> volmap_ids_;
 
public:
 
  std::vector<int> atom_ids;
 
  std::vector<cvm::rvector> atomic_gradients;
 
  virtual std::vector<std::vector<int> > get_atom_lists();
 
  std::vector<int> const &get_volmap_ids();
 
};
 
 
inline cvm::real const & colvar::force_constant() const
{
  return ext_force_k;
}
 
 
inline colvarvalue const & colvar::value() const
{
  return x_reported;
}
 
 
inline colvarvalue const & colvar::actual_value() const
{
  return x;
}
 
 
inline colvarvalue const & colvar::run_ave() const
{
  return runave;
}
 
 
inline colvarvalue const & colvar::velocity() const
{
  return v_reported;
}
 
 
inline colvarvalue const & colvar::total_force() const
{
  return ft_reported;
}
 
 
inline void colvar::add_bias_force(colvarvalue const &force)
{
  check_enabled(f_cv_gradient,
                std::string("applying a force to the variable \""+name+"\""));
  if (cvm::debug()) {
    cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
  }
  fb += force;
}
 
 
inline void colvar::add_bias_force_actual_value(colvarvalue const &force)
{
  if (cvm::debug()) {
    cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
  }
  fb_actual += force;
}
 
 
inline void colvar::reset_bias_force() {
  fb.type(value());
  fb.reset();
  fb_actual.type(value());
  fb_actual.reset();
}
 
 
namespace {
  // Tolerance parameter to decide when two boundaries coincide
  constexpr cvm::real colvar_boundaries_tol = 1.0e-10;
}
 
#endif