colvarmodule.h

Go to the documentation of this file.

// -*- 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 COLVARMODULE_H
#define COLVARMODULE_H
 
#include <cstdint>
 
#ifndef COLVARS_DEBUG
#define COLVARS_DEBUG false
#endif
 
#if defined(__FAST_MATH__)
// NOTE: This is used for fixing https://github.com/Colvars/colvars/issues/767
#define COLVARS_BOUNDED_INV_TRIGONOMETRIC_FUNC
#endif
 
#define COLVARS_USE_SOA
 
 
#include <cmath>
#include <iosfwd>
#include <string>
#include <vector>
 
#if defined(COLVARS_CUDA)
#include <cuda_runtime.h>
#endif
 
#if defined(COLVARS_HIP)
#include <hip/hip_runtime.h>
#define cudaHostAllocMapped hipHostMallocMapped
#define cudaHostAlloc hipHostMalloc
#define cudaFreeHost hipHostFree
#define cudaSuccess hipSuccess
#endif
 
class colvarparse;
class colvar;
class colvarbias;
class colvarproxy;
class colvarvalue;
 
 
class colvarmodule {
 
public:
 
  std::string version() const;
 
  int version_number() const;
 
  int patch_version_number() const;
 
#if ( defined(COLVARS_CUDA) || defined(COLVARS_HIP) )
  template <typename T>
  class CudaHostAllocator {
  public:
    using value_type = T;
 
    CudaHostAllocator() = default;
 
    template<typename U>
    constexpr CudaHostAllocator(const CudaHostAllocator<U>&) noexcept {}
 
    friend bool operator==(const CudaHostAllocator&, const CudaHostAllocator&) { return true; }
    friend bool operator!=(const CudaHostAllocator&, const CudaHostAllocator&) { return false; }
 
    T* allocate(size_t n) {
      T* ptr;
      if (cudaHostAlloc(&ptr, n * sizeof(T), cudaHostAllocMapped) != cudaSuccess) {
        throw std::bad_alloc();
      }
      return ptr;
    }
    void deallocate(T* ptr, size_t n) noexcept {
      cudaFreeHost(ptr);
    }
    template<typename U, typename... Args>
    void construct(U* p, Args&&... args) {
        new(p) U(std::forward<Args>(args)...);
    }
 
    template<typename U>
    void destroy(U* p) noexcept {
        p->~U();
    }
  };
#endif
 
private:
 
  int version_int = 0;
 
  int patch_version_int = 0;
 
public:
 
  typedef long long step_number;
 
  typedef  double    real;
 
 
  // Math functions
 
  static inline real integer_power(real const &x, int const n)
  {
    // Original code: math_special.h in LAMMPS
    double yy, ww;
    if (x == 0.0) return 0.0;
    int nn = (n > 0) ? n : -n;
    ww = x;
    for (yy = 1.0; nn != 0; nn >>= 1, ww *=ww) {
      if (nn & 1) yy *= ww;
    }
    return (n > 0) ? yy : 1.0/yy;
  }
 
  static inline real pow(real const &x, real const &y)
  {
    return ::pow(static_cast<double>(x), static_cast<double>(y));
  }
 
  static inline real floor(real const &x)
  {
    return ::floor(static_cast<double>(x));
  }
 
  static inline real fabs(real const &x)
  {
    return ::fabs(static_cast<double>(x));
  }
 
  static inline real sqrt(real const &x)
  {
    return ::sqrt(static_cast<double>(x));
  }
 
  static inline real sin(real const &x)
  {
    return ::sin(static_cast<double>(x));
  }
 
  static inline real cos(real const &x)
  {
    return ::cos(static_cast<double>(x));
  }
 
#ifndef PI
#define PI   3.14159265358979323846
#endif
#ifndef PI_2
#define PI_2 1.57079632679489661923
#endif
 
static inline real asin(real const &x)
{
#ifdef COLVARS_BOUNDED_INV_TRIGONOMETRIC_FUNC
    if (x <= -1.0) {
        return -PI_2;
    } else if (x >= 1.0) {
        return PI_2;
    } else {
        return ::asin(static_cast<double>(x));
    }
#else
    return ::asin(static_cast<double>(x));
#endif
}
 
static inline real acos(real const &x)
{
#ifdef COLVARS_BOUNDED_INV_TRIGONOMETRIC_FUNC
    if (x <= -1.0) {
        return PI;
    } else if (x >= 1.0) {
        return 0.0;
    } else {
        return ::acos(static_cast<double>(x));
    }
#else
    return ::acos(static_cast<double>(x));
#endif
}
 
  static inline real atan2(real const &x, real const &y)
  {
    return ::atan2(static_cast<double>(x), static_cast<double>(y));
  }
 
  static inline real exp(real const &x)
  {
    return ::exp(static_cast<double>(x));
  }
 
  static inline real logn(real const &x)
  {
    return ::log(static_cast<double>(x));
  }
 
  // Forward declarations
  class rvector;
  template <class T> class vector1d;
  template <class T> class matrix2d;
  class quaternion;
  class rotation;
 
  class usage;
  class memory_stream;
 
  typedef int residue_id;
 
  typedef rvector atom_pos;
 
  class rmatrix;
 
  // allow these classes to access protected data
  class atom;
  class atom_group;
  typedef std::vector<atom>::iterator       atom_iter;
  typedef std::vector<atom>::const_iterator atom_const_iter;
 
private:
 
  static int errorCode;
 
public:
 
  static void set_error_bits(int code);
 
  static bool get_error_bit(int code);
 
  static inline int get_error()
  {
    return errorCode;
  }
 
  static void clear_error();
 
  static step_number it;
  static step_number it_restart;
 
  static inline step_number step_relative()
  {
    return it - it_restart;
  }
 
  static inline step_number step_absolute()
  {
    return it;
  }
 
  bool binary_restart;
 
  static real debug_gradients_step_size;
 
private:
 
  std::string cvm_output_prefix;
 
public:
  static inline std::string &output_prefix()
  {
    colvarmodule *cv = colvarmodule::main();
    return cv->cvm_output_prefix;
  }
 
private:
 
  std::vector<colvar *> colvars;
 
  std::vector<colvar *> colvars_active;
 
  std::vector<colvar *> colvars_smp;
  std::vector<int> colvars_smp_items;
 
  std::vector<atom_group *> named_atom_groups_soa;
 
public:
 
  void register_named_atom_group_soa(atom_group *ag);
  void unregister_named_atom_group_soa(atom_group *ag);
 
  std::vector<colvar *> *variables();
 
  /* TODO: implement named CVCs
  static std::vector<cvc *>     cvcs;
  inline void register_cvc(cvc *p) {
    cvcs.push_back(p);
  }
  */
 
  std::vector<colvar *> *variables_active();
 
  std::vector<colvar *> *variables_active_smp();
 
  std::vector<int> *variables_active_smp_items();
 
  int calc_component_smp(int i);
 
  std::vector<colvarbias *> biases;
 
  real total_bias_energy;
 
private:
 
  void *num_biases_types_used_;
 
  std::vector<colvarbias *> biases_active_;
 
public:
 
  std::vector<colvarbias *> *biases_active();
 
  static inline bool debug()
  {
    return COLVARS_DEBUG;
  }
 
  size_t size() const;
 
  colvarmodule(colvarproxy *proxy);
 
private:
 
  colvarmodule();
 
public:
 
  ~colvarmodule();
 
  void set_initial_step(step_number it);
 
  int reset();
 
  int read_config_file(char const *config_file_name);
 
  int read_config_string(std::string const &conf);
 
  int parse_config(std::string &conf);
 
  std::string const & get_config() const;
 
  // Parse functions (setup internal data based on a string)
 
  static std::istream & getline(std::istream &is, std::string &line);
 
  int parse_global_params(std::string const &conf);
 
  int parse_colvars(std::string const &conf);
 
  int run_tcl_script(std::string const &filename);
 
  int parse_biases(std::string const &conf);
 
  int append_new_config(std::string const &conf);
 
  void config_changed();
 
private:
 
  std::string config_string;
 
  std::string extra_conf;
 
  template <class bias_type>
  int parse_biases_type(std::string const &conf, char const *keyword);
 
  bool check_new_bias(std::string &conf, char const *key);
 
  std::string source_Tcl_script;
 
public:
 
  size_t num_variables() const;
 
  size_t num_variables_feature(int feature_id) const;
 
  size_t num_biases() const;
 
  size_t num_biases_feature(int feature_id) const;
 
  size_t num_biases_type(std::string const &type) const;
 
  std::vector<std::string> const time_dependent_biases() const;
 
private:
  int catch_input_errors(int result);
 
public:
 
  // "Setup" functions (change internal data based on related data
  // from the proxy that may change during program execution)
  // No additional parsing is done within these functions
 
  int update_engine_parameters();
 
  int setup_input();
 
  int setup_output();
 
private:
 
  template <typename IST> IST & read_state_template_(IST &is);
 
  std::string default_input_state_file_;
 
  std::vector<unsigned char> input_state_buffer_;
 
public:
 
  std::istream & read_state(std::istream &is);
 
  memory_stream & read_state(memory_stream &is);
 
  int set_input_state_buffer(size_t n, unsigned char *buf);
 
  int set_input_state_buffer(std::vector<unsigned char> &buf);
 
  std::istream & read_objects_state(std::istream &is);
 
  memory_stream & read_objects_state(memory_stream &is);
 
  int print_total_forces_errning(bool warn_total_forces);
 
private:
  template <typename OST> OST &write_state_template_(OST &os);
 
public:
 
  std::ostream & write_state(std::ostream &os);
 
  memory_stream & write_state(memory_stream &os);
 
  int write_state_buffer(std::vector<unsigned char> &buffer);
 
  static std::string state_file_prefix(char const *filename);
 
  int open_traj_file(std::string const &file_name);
  int close_traj_file();
  std::ostream & write_traj(std::ostream &os);
  std::ostream & write_traj_label(std::ostream &os);
 
  int write_traj_files();
  int write_restart_file(std::string const &out_name);
  int write_output_files();
  static int backup_file(char const *filename);
 
  int write_restart_string(std::string &output);
 
  static colvarbias * bias_by_name(std::string const &name);
 
  static colvar * colvar_by_name(std::string const &name);
 
  static atom_group * atom_group_soa_by_name(std::string const& name);
 
  int change_configuration(std::string const &bias_name, std::string const &conf);
 
  std::string read_colvar(std::string const &name);
 
  real energy_difference(std::string const &bias_name, std::string const &conf);
 
  int calc();
 
  int calc_colvars();
 
  int calc_biases();
 
  int update_colvar_forces();
 
  int analyze();
 
  int end_of_step();
 
  int read_traj(char const *traj_filename,
                long        traj_read_begin,
                long        traj_read_end);
 
  static std::string to_str(char const *s);
 
  static std::string to_str(std::string const &s);
 
  static std::string to_str(bool x);
 
  static std::string to_str(int const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(size_t const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(long int const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(step_number const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(real const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(rvector const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(quaternion const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(colvarvalue const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(vector1d<real> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(matrix2d<real> const &x,
                            size_t width = 0, size_t prec = 0);
 
 
  static std::string to_str(std::vector<int> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<size_t> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<long int> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<real> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<rvector> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<quaternion> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<colvarvalue> const &x,
                            size_t width = 0, size_t prec = 0);
 
  static std::string to_str(std::vector<std::string> const &x,
                            size_t width = 0, size_t prec = 0);
 
#if ( defined(COLVARS_CUDA) || defined(COLVARS_HIP) )
  static std::string to_str(std::vector<rvector, CudaHostAllocator<rvector>> const &x,
                            size_t width = 0, size_t prec = 0);
  static std::string to_str(std::vector<real, CudaHostAllocator<real>> const &x,
                            size_t width = 0, size_t prec = 0);
#endif
 
 
  static std::string wrap_string(std::string const &s,
                                 size_t nchars);
 
  static size_t const it_width;
  static size_t const cv_prec;
  static size_t const cv_width;
  static size_t const en_prec;
  static size_t const en_width;
  static const char * const line_marker;
 
 
  // proxy functions
 
  static real dt();
 
  static void request_total_force();
 
  int cite_feature(std::string const &feature);
 
  std::string feature_report(int flag = 0);
 
  static void log(std::string const &message, int min_log_level = 10);
 
  static int error(std::string const &message, int code = -1);
 
private:
 
  static int log_level_;
 
public:
 
  static inline int log_level()
  {
    return log_level_;
  }
 
  static inline int log_init_messages()
  {
    return 1;
  }
 
  static inline int log_user_params()
  {
    return 2;
  }
 
  static inline int log_default_params()
  {
    return 3;
  }
 
  static inline int log_output_files()
  {
    return 4;
  }
 
  static inline int log_input_files()
  {
    return 5;
  }
 
  static rvector position_distance(atom_pos const &pos1,
                                   atom_pos const &pos2);
 
  std::vector<std::string> index_file_names;
 
  std::vector<std::string> index_group_names;
 
  std::vector<std::vector<int> *> index_groups;
 
  int read_index_file(char const *filename);
 
  int reset_index_groups();
 
  static int load_coords(char const *filename,
                         std::vector<rvector> *pos,
                         atom_group *atoms,
                         std::string const &pdb_field,
                         double pdb_field_value = 0.0);
 
  int load_coords_xyz(char const *filename,
                      std::vector<rvector> *pos,
                      atom_group *atoms,
                      bool keep_open = false);
 
  static size_t cv_traj_freq;
 
  static size_t restart_out_freq;
  std::string   restart_out_name;
 
  static real rand_gaussian();
 
protected:
 
  colvarparse *parse;
 
  std::string cv_traj_name;
 
  bool cv_traj_write_labels;
 
  std::string restart_version_str;
 
  int restart_version_int;
 
  size_t depth_s;
 
  std::vector<size_t> depth_v;
 
  int xyz_reader_use_count;
 
  usage *usage_;
 
  real max_gradient_error = 0.;
 
public:
 
  inline std::string restart_version() const
  {
    return restart_version_str;
  }
 
  inline int restart_version_number() const
  {
    return restart_version_int;
  }
 
  static size_t & depth();
 
  static void increase_depth();
 
  static void decrease_depth();
 
  static inline bool scripted_forces()
  {
    return use_scripted_forces;
  }
 
  static bool use_scripted_forces;
 
  static bool scripting_after_biases;
 
  int calc_scripted_forces();
 
  void record_gradient_error(real error) {
    if (error > max_gradient_error) max_gradient_error = error;
  }
 
  real get_max_gradient_error() {
    return max_gradient_error;
  }
 
  static colvarproxy *proxy;
 
  static colvarmodule *main();
 
};
 
 
typedef colvarmodule cvm;
 
 
std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
std::istream & operator >> (std::istream &is, cvm::rvector &v);
 
 
namespace {
  constexpr int32_t COLVARS_OK = 0;
  constexpr int32_t COLVARS_ERROR = 1;
  constexpr int32_t COLVARS_NOT_IMPLEMENTED = (1<<1);
  constexpr int32_t COLVARS_INPUT_ERROR     = (1<<2); // out of bounds or inconsistent input
  constexpr int32_t COLVARS_BUG_ERROR       = (1<<3); // Inconsistent state indicating bug
  constexpr int32_t COLVARS_FILE_ERROR      = (1<<4);
  constexpr int32_t COLVARS_MEMORY_ERROR    = (1<<5);
  constexpr int32_t COLVARS_NO_SUCH_FRAME   = (1<<6); // Cannot load the requested frame
}
 
 
#endif