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>
#include <unordered_map>
#include <memory>
#include <cstdio>
#include <atomic>
 
#include "colvars_version.h"
 
#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>
#include "colvar_gpu_support.h"
 
class colvarparse;
class colvar;
class colvarbias;
class colvarproxy;
class colvarvalue;
class colvardeps;
 
 
class colvarmodule {
 
public:
 
  std::string version() const;
 
  int version_number() const;
 
  int patch_version_number() const;
 
private:
 
  int version_int = 0;
 
  int patch_version_int = 0;
 
public:
 
  typedef long long step_number;
 
  typedef  double    real;
 
 
  // Math functions
 
  static COLVARS_HOST_DEVICE 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;
  }
 
  template <int n>
  static COLVARS_HOST_DEVICE inline real
  positive_integer_power(real const& x) {
    static_assert(n > 0, "n must be greater than zero in positive_integer_power.");
    real y = 1.0;
    for (int i = 0; i < n; ++i) {
      y *= x;
    }
    return y;
  }
 
  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 ceil(real const &x)
  {
    return ::ceil(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 tanh(real const &x)
  {
    return ::tanh(static_cast<double>(x));
  }
 
  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 system_boundary_conditions;
 
  class usage;
  class memory_stream;
 
  typedef int residue_id;
 
  typedef rvector atom_pos;
 
  class rmatrix;
 
  // NOTE: Just here for ensuring the compilation of GROMACS with Colvars
  struct atom;
  // class atom_group_base;
  class atom_group;
  // class atom_group_gpu;
  // typedef std::vector<atom>::iterator       atom_iter;
  // typedef std::vector<atom>::const_iterator atom_const_iter;
 
private:
 
  int errorCode = 0;
 
public:
 
  void set_error_bits(int code);
 
  bool get_error_bit(int code);
 
  inline int get_error()
  {
    return errorCode;
  }
 
  void clear_error();
 
  step_number it = 0;
  step_number it_restart = 0;
 
  inline step_number step_relative()
  {
    return it - it_restart;
  }
 
  inline step_number step_absolute()
  {
    return it;
  }
 
  bool binary_restart;
 
  static real debug_gradients_step_size;
 
private:
 
  std::string cvm_output_prefix;
 
public:
  inline std::string &output_prefix()
  {
    return this->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;
 
public:
 
  void register_named_atom_group(atom_group *ag);
  void unregister_named_atom_group(atom_group *ag);
 
  std::vector<colvar *> *variables();
 
  /* TODO: implement named CVCs
  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 constexpr debug()
  {
#if (defined(__HIP_DEVICE_COMPILE__)) || (defined(__CUDA_ARCH__))
    return false;
#else
    return COLVARS_DEBUG;
#endif
  }
 
  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);
 
  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();
  int backup_file(char const *filename);
 
  int write_restart_string(std::string &output);
 
  colvarbias * bias_by_name(std::string const &name);
 
  colvar * colvar_by_name(std::string const &name);
 
  atom_group * atom_group_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 proxy_gpu_buffers_reallocated_done();
 
  int analyze();
 
  int end_of_step();
 
  int read_traj(char const *traj_filename,
                long        traj_read_begin,
                long        traj_read_end);
 
  // In a first phase, keep to_str() static
 
  static std::string to_str(char const *s);
 
  static std::string to_str(const void* ptr);
 
  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, colvars_gpu::CudaHostAllocator<rvector>> const &x,
                            size_t width = 0, size_t prec = 0);
  static std::string to_str(std::vector<real, colvars_gpu::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);
 
  // i/o constants
 
  static constexpr size_t it_width = 12;
  static constexpr size_t cv_prec = 14;
  static constexpr size_t cv_width = 21;
  static constexpr size_t en_prec = 14;
  static constexpr size_t en_width = 21;
  static constexpr const char line_marker[] =
    "----------------------------------------------------------------------\n";
 
  // proxy functions
 
  real dt();
 
  void request_total_force();
 
  int cite_feature(std::string const &feature);
 
  std::string feature_report(int flag = 0);
 
  void log(std::string const &message, int min_log_level = 10);
 
  static void log_static(std::string const &message) {
    if (colvarmodule::main()) {
      colvarmodule::main()->log(message);
    } else {
      std::printf("colvars: %s\n", message.c_str());
    }
  }
 
  int error(std::string const &message, int code = -1);
 
  static int error_static(std::string const &message, int code = -1) {
    if (colvarmodule::main()) {
      code = colvarmodule::main()->error(message, code);
    } else {
      std::cerr << "colvars: " << message << std::endl;
      exit(-1);
    }
    return code;
  }
 
private:
 
  int log_level_ = 10;
 
public:
 
  inline int log_level()
  {
    return log_level_;
  }
 
  inline int log_init_messages()
  {
    return 1;
  }
 
  inline int log_user_params()
  {
    return 2;
  }
 
  inline int log_default_params()
  {
    return 3;
  }
 
  inline int log_output_files()
  {
    return 4;
  }
 
  inline int log_input_files()
  {
    return 5;
  }
 
  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();
 
  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);
 
  size_t cv_traj_freq = 0;
 
  size_t restart_out_freq = 0;
  std::string   restart_out_name;
 
  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_;
 
  std::atomic<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;
  }
 
  size_t & depth();
 
  void increase_depth();
 
  void decrease_depth();
 
  inline bool scripted_forces()
  {
    return use_scripted_forces;
  }
 
  bool use_scripted_forces = false;
 
  bool scripting_after_biases = true;
 
  int calc_scripted_forces();
 
  void record_gradient_error(real error) {
    real current = max_gradient_error.load(std::memory_order_relaxed);
    while (error > current &&
          !max_gradient_error.compare_exchange_weak(current, error, std::memory_order_relaxed)) {
    }
  }
 
  real get_max_gradient_error() {
    return max_gradient_error;
  }
 
  colvarproxy *proxy = nullptr;
  static colvarproxy *proxy_static;
 
  static colvarmodule *main();
 
#if defined (COLVARS_CUDA) || defined (COLVARS_HIP)
  template <typename T>
  using allocator_type = colvars_gpu::CudaHostAllocator<T>;
#else
  template <typename T>
  using allocator_type = std::allocator<T>;
#endif
  using ag_vector_real_t = std::vector<real, allocator_type<real>>;
 
};
 
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