colvarmodule.h

Go to the documentation of this file.

// -*- c++ -*-

#ifndef COLVARMODULE_H
#define COLVARMODULE_H

#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2017-03-13"
#endif

#ifndef COLVARS_DEBUG
#define COLVARS_DEBUG false
#endif


#define COLVARS_OK 0
#define COLVARS_ERROR   1
#define COLVARS_NOT_IMPLEMENTED (1<<1)
#define INPUT_ERROR     (1<<2) // out of bounds or inconsistent input
#define BUG_ERROR       (1<<3) // Inconsistent state indicating bug
#define FILE_ERROR      (1<<4)
#define MEMORY_ERROR    (1<<5)
#define FATAL_ERROR     (1<<6) // Should be set, or not, together with other bits
#define DELETE_COLVARS  (1<<7) // Instruct the caller to delete cvm
#define COLVARS_NO_SUCH_FRAME (1<<8) // Cannot load the requested frame

#include <iostream>
#include <iomanip>
#include <string>
#include <cstring>
#include <sstream>
#include <fstream>
#include <cmath>
#include <vector>
#include <list>

#ifdef NAMD_VERSION
// use Lustre-friendly wrapper to POSIX write()
#include "fstream_namd.h"
#endif

class colvarparse;
class colvar;
class colvarbias;
class colvarproxy;
class colvarscript;


class colvarmodule {

private:

  colvarmodule();

public:

  friend class colvarproxy;
  // TODO colvarscript should be unaware of colvarmodule's internals
  friend class colvarscript;

  typedef  double    real;
  typedef  int       residue_id;

  class rvector;
  template <class T> class vector1d;
  template <class T> class matrix2d;
  class quaternion;
  class rotation;

  typedef rvector atom_pos;

  class rmatrix;

  // allow these classes to access protected data
  class atom;
  class atom_group;
  friend class atom;
  friend class atom_group;
  typedef std::vector<atom>::iterator       atom_iter;
  typedef std::vector<atom>::const_iterator atom_const_iter;

protected:

  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 long it;
  static long it_restart;

  static inline long step_relative()
  {
    return it - it_restart;
  }

  static inline long step_absolute()
  {
    return it;
  }

  bool it_restart_from_state_file;

  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;

public:

  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();

  std::vector<colvarbias *> biases;

  real total_bias_energy;

private:

  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);

  ~colvarmodule();

  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);

  // 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 parse_biases(std::string const &conf);

  int append_new_config(std::string const &conf);

private:

  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);

public:

  static int num_biases_feature(int feature_id);

  static int num_biases_type(std::string const &type);

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 setup();

  int setup_input();

  int setup_output();

#ifdef NAMD_VERSION
  typedef ofstream_namd ofstream;
#else
  typedef std::ofstream ofstream;
#endif

  std::istream & read_restart(std::istream &is);
  std::ostream & write_restart(std::ostream &os);

  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_files();
  int write_output_files();
  static int backup_file(char const *filename);

  static colvarbias * bias_by_name(std::string const &name);

  static colvar * colvar_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 bias_bin_num(std::string const &bias_name);
  int bias_current_bin(std::string const &bias_name);
  int bias_bin_count(std::string const &bias_name, size_t bin_index);
  int bias_share(std::string const &bias_name);

  int calc();

  int calc_colvars();

  int calc_biases();

  int update_colvar_forces();

  int analyze();

  int read_traj(char const *traj_filename,
                long        traj_read_begin,
                long        traj_read_end);

  template<typename T> static std::string to_str(T const &x,
                                                  size_t const &width = 0,
                                                  size_t const &prec = 0);
  template<typename T> static std::string to_str(std::vector<T> const &x,
                                                  size_t const &width = 0,
                                                  size_t const &prec = 0);

  static inline std::string wrap_string(std::string const &s,
                                         size_t const &nchars)
  {
    if (!s.size())
      return std::string(nchars, ' ');
    else
      return ( (s.size() <= size_t(nchars)) ?
               (s+std::string(nchars-s.size(), ' ')) :
               (std::string(s, 0, 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 unit_angstrom();

  static real boltzmann();

  static real temperature();

  static real dt();

  static void request_total_force();

  static void log(std::string const &message);

  static void fatal_error(std::string const &message);

  static void error(std::string const &message, int code = COLVARS_ERROR);

  static void exit(std::string const &message);

  // Replica exchange commands.
  static bool replica_enabled();
  static int replica_index();
  static int replica_num();
  static void replica_comm_barrier();
  static int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
  static int replica_comm_send(char* msg_data, int msg_len, int dest_rep);

  static rvector position_distance(atom_pos const &pos1,
                                    atom_pos const &pos2);


  static real position_dist2(atom_pos const &pos1,
                              atom_pos const &pos2);

  static void select_closest_image(atom_pos &pos,
                                    atom_pos const &ref_pos);

  static void select_closest_images(std::vector<atom_pos> &pos,
                                     atom_pos const &ref_pos);


  std::list<std::string> index_group_names;

  std::list<std::vector<int> > index_groups;

  int read_index_file(char const *filename);


  static int load_atoms(char const *filename,
                        atom_group &atoms,
                        std::string const &pdb_field,
                        double const pdb_field_value = 0.0);

  static int load_coords(char const *filename,
                         std::vector<atom_pos> &pos,
                         const std::vector<int> &indices,
                         std::string const &pdb_field,
                         double const pdb_field_value = 0.0);

  static int load_coords_xyz(char const *filename,
                              std::vector<atom_pos> &pos,
                              const std::vector<int> &indices);

  static size_t cv_traj_freq;

  static bool   b_analysis;

  static size_t restart_out_freq;
  std::string   restart_out_name;

  static real rand_gaussian(void);

protected:

  std::ifstream config_s;

  colvarparse *parse;

  std::string cv_traj_name;

  colvarmodule::ofstream cv_traj_os;

  bool cv_traj_append;

  colvarmodule::ofstream restart_out_os;

private:

  size_t depth_s;

  std::vector<size_t> depth_v;

public:

  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();

  static colvarproxy *proxy;

  static colvarmodule *main();

};


typedef colvarmodule cvm;


#include "colvartypes.h"


std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
std::istream & operator >> (std::istream &is, cvm::rvector &v);


template<typename T> std::string cvm::to_str(T const &x,
                                              size_t const &width,
                                              size_t const &prec) {
  std::ostringstream os;
  if (width) os.width(width);
  if (prec) {
    os.setf(std::ios::scientific, std::ios::floatfield);
    os.precision(prec);
  }
  os << x;
  return os.str();
}

template<typename T> std::string cvm::to_str(std::vector<T> const &x,
                                              size_t const &width,
                                              size_t const &prec) {
  if (!x.size()) return std::string("");
  std::ostringstream os;
  if (prec) {
    os.setf(std::ios::scientific, std::ios::floatfield);
  }
  os << "{ ";
  if (width) os.width(width);
  if (prec) os.precision(prec);
  os << x[0];
  for (size_t i = 1; i < x.size(); i++) {
    os << ", ";
    if (width) os.width(width);
    if (prec) os.precision(prec);
    os << x[i];
  }
  os << " }";
  return os.str();
}


#include "colvarproxy.h"


inline cvm::real cvm::unit_angstrom()
{
  return proxy->unit_angstrom();
}

inline cvm::real cvm::boltzmann()
{
  return proxy->boltzmann();
}

inline cvm::real cvm::temperature()
{
  return proxy->temperature();
}

inline cvm::real cvm::dt()
{
  return proxy->dt();
}

// Replica exchange commands
inline bool cvm::replica_enabled() {
  return proxy->replica_enabled();
}
inline int cvm::replica_index() {
  return proxy->replica_index();
}
inline int cvm::replica_num() {
  return proxy->replica_num();
}
inline void cvm::replica_comm_barrier() {
  return proxy->replica_comm_barrier();
}
inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
  return proxy->replica_comm_recv(msg_data,buf_len,src_rep);
}
inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
  return proxy->replica_comm_send(msg_data,msg_len,dest_rep);
}


inline void cvm::request_total_force()
{
  proxy->request_total_force(true);
}

inline void cvm::select_closest_image(atom_pos &pos,
                                       atom_pos const &ref_pos)
{
  proxy->select_closest_image(pos, ref_pos);
}

inline void cvm::select_closest_images(std::vector<atom_pos> &pos,
                                        atom_pos const &ref_pos)
{
  proxy->select_closest_images(pos, ref_pos);
}

inline cvm::rvector cvm::position_distance(atom_pos const &pos1,
                                            atom_pos const &pos2)
{
  return proxy->position_distance(pos1, pos2);
}

inline cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
                                      cvm::atom_pos const &pos2)
{
  return proxy->position_dist2(pos1, pos2);
}

inline cvm::real cvm::rand_gaussian(void)
{
  return proxy->rand_gaussian();
}

#endif