colvarproxy_system.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 COLVARPROXY_SYSTEM_H
#define COLVARPROXY_SYSTEM_H
 
 
class colvarproxy_system {
 
public:
 
  colvarproxy_system();
 
  virtual ~colvarproxy_system();
 
  std::string units;
 
  virtual int set_unit_system(std::string const &units, bool check_only);
 
  inline cvm::real angstrom_to_internal(cvm::real l) const
  {
    return l * angstrom_value_;
  }
 
  inline cvm::real internal_to_angstrom(cvm::real l) const
  {
    return l / angstrom_value_;
  }
 
  inline cvm::real boltzmann() const
  {
    return boltzmann_;
  }
 
  inline cvm::real target_temperature() const
  {
    return target_temperature_;
  }
 
  virtual int set_target_temperature(cvm::real T);
 
  inline double dt() const
  {
    return timestep_;
  }
 
  virtual int set_integration_timestep(cvm::real dt);
 
  inline int time_step_factor() const
  {
    return time_step_factor_;
  }
 
  virtual int set_time_step_factor(int fact);
 
  virtual cvm::real rand_gaussian(void);
 
  virtual void add_energy(cvm::real energy);
 
  inline bool & use_internal_pbc() { return use_internal_pbc_; }
 
  virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
                                         cvm::atom_pos const &pos2) const;
 
  cvm::rvector position_distance_internal(cvm::atom_pos const &pos1,
                                          cvm::atom_pos const &pos2) const;
 
  static cvm::rvector position_distance_kernel(cvm::atom_pos const &pos1,
                                               cvm::atom_pos const &pos2,
                                               cvm::rvector const &a,
                                               cvm::rvector const &b,
                                               cvm::rvector const &c,
                                               cvm::rvector const &a_r,
                                               cvm::rvector const &b_r,
                                               cvm::rvector const &c_r,
                                               bool a_p = false,
                                               bool b_p = false,
                                               bool c_p = false);
 
  void update_pbc_lattice();
 
  void reset_pbc_lattice();
 
  virtual void request_total_force(bool yesno);
 
  virtual bool total_forces_enabled() const;
 
  virtual bool total_forces_same_step() const;
 
  virtual int get_molid(int &molid);
 
  virtual int get_alch_lambda(cvm::real* lambda);
 
  void set_alch_lambda(cvm::real lambda);
 
  virtual int send_alch_lambda();
 
  virtual int request_alch_energy_freq(int const /* freq */) {
    return COLVARS_OK;
  }
 
  virtual int get_dE_dlambda(cvm::real* dE_dlambda);
 
  virtual int apply_force_dE_dlambda(cvm::real* force);
 
  virtual int get_d2E_dlambda2(cvm::real* d2E_dlambda2);
 
  cvm::real indirect_lambda_biasing_force;
 
  virtual cvm::real get_accelMD_factor() const {
    cvm::error_static("Error: accessing the reweighting factor of accelerated MD  "
               "is not yet implemented in the MD engine.\n",
               COLVARS_NOT_IMPLEMENTED);
    return 1.0;
  }
  virtual bool accelMD_enabled() const {
    return false;
  }
 
protected:
 
  cvm::real cached_alch_lambda;
 
  bool cached_alch_lambda_changed;
 
  cvm::real boltzmann_;
 
  cvm::real target_temperature_;
 
  double timestep_;
 
  int time_step_factor_ = 1;
 
  cvm::real angstrom_value_;
 
  cvm::real kcal_mol_value_;
 
  bool total_force_requested;
 
  bool use_internal_pbc_ = false;
 
  enum Boundaries_type {
    boundaries_non_periodic,
    boundaries_pbc_ortho,
    boundaries_pbc_triclinic,
    boundaries_unsupported
  };
 
  Boundaries_type boundaries_type;
 
  cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
 
  cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
};
 
 
inline cvm::rvector colvarproxy_system::position_distance_internal(cvm::atom_pos const &pos1,
                                                                   cvm::atom_pos const &pos2) const
{
  if (boundaries_type == boundaries_unsupported) {
    cvm::error_static("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR);
    return cvm::rvector({0.0, 0.0, 0.0});
  }
 
  if (boundaries_type == boundaries_non_periodic) {
    return pos2 - pos1;
  }
 
  // Periodicity flags are hard-coded, because this is the only case supported so far other than
  // the two above
  return position_distance_kernel(pos1, pos2, unit_cell_x, unit_cell_y, unit_cell_z,
                                  reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z,
                                  true, true, true);
}
 
 
inline cvm::rvector colvarproxy_system::position_distance_kernel(cvm::atom_pos const &pos1,
                                                                 cvm::atom_pos const &pos2,
                                                                 cvm::rvector const &a,
                                                                 cvm::rvector const &b,
                                                                 cvm::rvector const &c,
                                                                 cvm::rvector const &a_r,
                                                                 cvm::rvector const &b_r,
                                                                 cvm::rvector const &c_r,
                                                                 bool a_p,
                                                                 bool b_p,
                                                                 bool c_p)
{
  cvm::rvector diff = (pos2 - pos1);
 
  cvm::real const x_shift = std::floor(a_r * diff + 0.5);
  cvm::real const y_shift = std::floor(b_r * diff + 0.5);
  cvm::real const z_shift = std::floor(c_r * diff + 0.5);
 
  if (a_p) {
    diff.x -= x_shift * a.x + y_shift * b.x + z_shift * c.x;
  }
 
  if (b_p) {
    diff.y -= x_shift * a.y + y_shift * b.y + z_shift * c.y;
  }
 
  if (c_p) {
    diff.z -= x_shift * a.z + y_shift * b.z + z_shift * c.z;
  }
 
  return diff;
}
 
 
#endif