lammps/lib/colvars/colvarmodule.cpp

/// -*- c++ -*-

#include <sstream>
#include <string.h>

#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias.h"
#include "colvarbias_abf.h"
#include "colvarbias_alb.h"
#include "colvarbias_histogram.h"
#include "colvarbias_meta.h"
#include "colvarbias_restraint.h"
#include "colvarscript.h"


colvarmodule::colvarmodule(colvarproxy *proxy_in)
{
  // pointer to the proxy object
  if (proxy == NULL) {
    proxy = proxy_in;
    parse = new colvarparse();
  } else {
    // TODO relax this error to handle multiple molecules in VMD
    // once the module is not static anymore
    cvm::error("Error: trying to allocate the collective "
               "variable module twice.\n");
    return;
  }
  cvm::log(cvm::line_marker);
  cvm::log("Initializing the collective variables module, version "+
           cvm::to_str(COLVARS_VERSION)+".\n");
  cvm::log("Please cite Fiorin et al, Mol Phys 2013 in any publication "
           "based on this calculation.\n");

  if (proxy->smp_enabled() == COLVARS_OK) {
    cvm::log("SMP parallelism is available.\n");
  }

  // set initial default values

  // "it_restart" will be set by the input state file, if any;
  // "it" should be updated by the proxy
  colvarmodule::it = colvarmodule::it_restart = 0;
  colvarmodule::it_restart_from_state_file = true;

  colvarmodule::use_scripted_forces = false;

  colvarmodule::b_analysis = false;

  colvarmodule::debug_gradients_step_size = 1.0e-07;

  colvarmodule::rotation::monitor_crossings = false;
  colvarmodule::rotation::crossing_threshold = 1.0e-02;

  colvarmodule::cv_traj_freq = 100;
  colvarmodule::restart_out_freq = proxy->restart_frequency();

  // by default overwrite the existing trajectory file
  colvarmodule::cv_traj_append = false;
}


int colvarmodule::read_config_file(char const  *config_filename)
{
  cvm::log(cvm::line_marker);
  cvm::log("Reading new configuration from file \""+
           std::string(config_filename)+"\":\n");

  // open the configfile
  config_s.open(config_filename);
  if (!config_s.is_open()) {
    cvm::error("Error: in opening configuration file \""+
               std::string(config_filename)+"\".\n",
               FILE_ERROR);
    return COLVARS_ERROR;
  }

  // read the config file into a string
  std::string conf = "";
  std::string line;
  while (colvarparse::getline_nocomments(config_s, line)) {
    conf.append(line+"\n");
  }
  config_s.close();

  return parse_config(conf);
}


int colvarmodule::read_config_string(std::string const &config_str)
{
  cvm::log(cvm::line_marker);
  cvm::log("Reading new configuration:\n");
  std::istringstream config_s(config_str);

  // strip the comments away
  std::string conf = "";
  std::string line;
  while (colvarparse::getline_nocomments(config_s, line)) {
    conf.append(line+"\n");
  }
  return parse_config(conf);
}


int colvarmodule::parse_config(std::string &conf)
{
  // parse global options
  if (catch_input_errors(parse_global_params(conf))) {
    return get_error();
  }

  // parse the options for collective variables
  if (catch_input_errors(parse_colvars(conf))) {
    return get_error();
  }

  // parse the options for biases
  if (catch_input_errors(parse_biases(conf))) {
    return get_error();
  }

  // done parsing known keywords, check that all keywords found were valid ones
  if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) {
    return get_error();
  }

  cvm::log(cvm::line_marker);
  cvm::log("Collective variables module (re)initialized.\n");
  cvm::log(cvm::line_marker);

  if (cv_traj_os.is_open()) {
    // configuration might have changed, better redo the labels
    write_traj_label(cv_traj_os);
  }

  return get_error();
}


int colvarmodule::parse_global_params(std::string const &conf)
{
  std::string index_file_name;
  if (parse->get_keyval(conf, "indexFile", index_file_name)) {
    read_index_file(index_file_name.c_str());
  }

  parse->get_keyval(conf, "analysis", b_analysis, b_analysis);

  parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
                    debug_gradients_step_size,
                    colvarparse::parse_silent);

  parse->get_keyval(conf, "monitorEigenvalueCrossing",
                    colvarmodule::rotation::monitor_crossings,
                    colvarmodule::rotation::monitor_crossings,
                    colvarparse::parse_silent);
  parse->get_keyval(conf, "eigenvalueCrossingThreshold",
                    colvarmodule::rotation::crossing_threshold,
                    colvarmodule::rotation::crossing_threshold,
                    colvarparse::parse_silent);

  parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq);
  parse->get_keyval(conf, "colvarsRestartFrequency",
                    restart_out_freq, restart_out_freq);

  // if this is true when initializing, it means
  // we are continuing after a reset(): default to true
  parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append);

  parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false);

  parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true);

  if (use_scripted_forces && !proxy->force_script_defined) {
    cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::parse_colvars(std::string const &conf)
{
  if (cvm::debug())
    cvm::log("Initializing the collective variables.\n");

  std::string colvar_conf = "";
  size_t pos = 0;
  while (parse->key_lookup(conf, "colvar", colvar_conf, pos)) {

    if (colvar_conf.size()) {
      cvm::log(cvm::line_marker);
      cvm::increase_depth();
      colvars.push_back(new colvar(colvar_conf));
      if (cvm::get_error() ||
          ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) {
        cvm::log("Error while constructing colvar number " +
                 cvm::to_str(colvars.size()) + " : deleting.");
        delete colvars.back();  // the colvar destructor updates the colvars array
        return COLVARS_ERROR;
      }
      cvm::decrease_depth();
    } else {
      cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR);
      return COLVARS_ERROR;
    }
    cvm::decrease_depth();
    colvar_conf = "";
  }

  if (!colvars.size()) {
    cvm::log("Warning: no collective variables defined.\n");
  }

  if (colvars.size())
    cvm::log(cvm::line_marker);
  cvm::log("Collective variables initialized, "+
           cvm::to_str(colvars.size())+
           " in total.\n");

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


bool colvarmodule::check_new_bias(std::string &conf, char const *key)
{
  if (cvm::get_error() ||
      (biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
    cvm::log("Error while constructing bias number " +
             cvm::to_str(biases.size()) + " : deleting.\n");
    delete biases.back(); // the bias destructor updates the biases array
    return true;
  }
  return false;
}


template <class bias_type>
int colvarmodule::parse_biases_type(std::string const &conf,
                                    char const *keyword,
                                    size_t &bias_count)
{
  std::string bias_conf = "";
  size_t conf_saved_pos = 0;
  while (parse->key_lookup(conf, keyword, bias_conf, conf_saved_pos)) {
    if (bias_conf.size()) {
      cvm::log(cvm::line_marker);
      cvm::increase_depth();
      biases.push_back(new bias_type(bias_conf, keyword));
      if (cvm::check_new_bias(bias_conf, keyword)) {
        return COLVARS_ERROR;
      }
      cvm::decrease_depth();
      bias_count++;
    } else {
      cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n",
                 INPUT_ERROR);
      return COLVARS_ERROR;
    }
    bias_conf = "";
  }
  return COLVARS_OK;
}


int colvarmodule::parse_biases(std::string const &conf)
{
  if (cvm::debug())
    cvm::log("Initializing the collective variables biases.\n");

  /// initialize ABF instances
  parse_biases_type<colvarbias_abf>(conf, "abf", n_abf_biases);

  /// initialize adaptive linear biases
  parse_biases_type<colvarbias_alb>(conf, "ALB", n_rest_biases);

  /// initialize harmonic restraints
  parse_biases_type<colvarbias_restraint_harmonic>(conf, "harmonic", n_rest_biases);

  /// initialize histograms
  parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);

  /// initialize linear restraints
  parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases);

  /// initialize metadynamics instances
  parse_biases_type<colvarbias_meta>(conf, "metadynamics", n_meta_biases);

  if (use_scripted_forces) {
    cvm::log(cvm::line_marker);
    cvm::increase_depth();
    cvm::log("User forces script will be run at each bias update.");
    cvm::decrease_depth();
  }

  for (size_t i = 0; i < biases.size(); i++) {
    biases[i]->enable(cvm::deps::f_cvb_active);
    if (cvm::debug())
      biases[i]->print_state();
  }

  if (biases.size() || use_scripted_forces) {
    cvm::log(cvm::line_marker);
    cvm::log("Collective variables biases initialized, "+
             cvm::to_str(biases.size())+" in total.\n");
  } else {
    if (!use_scripted_forces) {
      cvm::log("No collective variables biases were defined.\n");
    }
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::catch_input_errors(int result)
{
  if (result != COLVARS_OK || get_error()) {
    set_error_bit(result);
    set_error_bit(INPUT_ERROR);
    parse->reset();
    return get_error();
  }
  return COLVARS_OK;
}


colvarbias * colvarmodule::bias_by_name(std::string const &name) {
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    if ((*bi)->name == name) {
      return (*bi);
    }
  }
  return NULL;
}


colvar *colvarmodule::colvar_by_name(std::string const &name) {
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    if ((*cvi)->name == name) {
      return (*cvi);
    }
  }
  return NULL;
}


int colvarmodule::change_configuration(std::string const &bias_name,
                                       std::string const &conf)
{
  // This is deprecated; supported strategy is to delete the bias
  // and parse the new config
  cvm::increase_depth();
  colvarbias *b;
  b = bias_by_name(bias_name);
  if (b == NULL) { cvm::error("Error: bias not found: " + bias_name); }
  b->change_configuration(conf);
  cvm::decrease_depth();
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


std::string colvarmodule::read_colvar(std::string const &name)
{
  cvm::increase_depth();
  colvar *c;
  std::stringstream ss;
  c = colvar_by_name(name);
  if (c == NULL) { cvm::fatal_error("Error: colvar not found: " + name); }
  ss << c->value();
  cvm::decrease_depth();
  return ss.str();
}

cvm::real colvarmodule::energy_difference(std::string const &bias_name,
                                          std::string const &conf)
{
  cvm::increase_depth();
  colvarbias *b;
  cvm::real energy_diff = 0.;
  b = bias_by_name(bias_name);
  if (b == NULL) { cvm::fatal_error("Error: bias not found: " + bias_name); }
  energy_diff = b->energy_difference(conf);
  cvm::decrease_depth();
  return energy_diff;
}

int colvarmodule::bias_current_bin(std::string const &bias_name)
{
  cvm::increase_depth();
  int ret;
  colvarbias *b = bias_by_name(bias_name);

  if (b != NULL) {
    ret = b->current_bin();
  } else {
    cvm::error("Error: bias not found.\n");
    ret = COLVARS_ERROR;
  }

  cvm::decrease_depth();
  return ret;
}

int colvarmodule::bias_bin_num(std::string const &bias_name)
{
  cvm::increase_depth();
  int ret;
  colvarbias *b = bias_by_name(bias_name);

  if (b != NULL) {
    ret = b->bin_num();
  } else {
    cvm::error("Error: bias not found.\n");
    ret = COLVARS_ERROR;
  }

  cvm::decrease_depth();
  return ret;
}

int colvarmodule::bias_bin_count(std::string const &bias_name, size_t bin_index)
{
  cvm::increase_depth();
  int ret;
  colvarbias *b = bias_by_name(bias_name);

  if (b != NULL) {
    ret = b->bin_count(bin_index);
  } else {
    cvm::error("Error: bias not found.\n");
    ret = COLVARS_ERROR;
  }

  cvm::decrease_depth();
  return ret;
}

int colvarmodule::bias_share(std::string const &bias_name)
{
  cvm::increase_depth();
  int ret;
  colvarbias *b = bias_by_name(bias_name);

  if (b != NULL) {
    b->replica_share();
    ret = COLVARS_OK;
  } else {
    cvm::error("Error: bias not found.\n");
    ret = COLVARS_ERROR;
  }

  cvm::decrease_depth();
  return ret;
}


int colvarmodule::calc()
{
  int error_code = COLVARS_OK;

  if (cvm::debug()) {
    cvm::log(cvm::line_marker);
    cvm::log("Collective variables module, step no. "+
             cvm::to_str(cvm::step_absolute())+"\n");
  }

  combine_errors(error_code, calc_colvars());
  // set biasing forces to zero before biases are calculated and summed over
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end(); cvi++) {
    (*cvi)->reset_bias_force();
  }
  combine_errors(error_code, calc_biases());
  combine_errors(error_code, update_colvar_forces());

  if (cvm::b_analysis) {
    combine_errors(error_code, analyze());
  }

  // write trajectory files, if needed
  if (cv_traj_freq && cv_traj_name.size()) {
    combine_errors(error_code, write_traj_files());
  }

  // write restart files, if needed
  if (restart_out_freq && restart_out_name.size()) {
    combine_errors(error_code, write_restart_files());
  }

  return error_code;
}


int colvarmodule::calc_colvars()
{
  if (cvm::debug())
    cvm::log("Calculating collective variables.\n");
  // calculate collective variables and their gradients

  int error_code = COLVARS_OK;
  std::vector<colvar *>::iterator cvi;

  // Determine which colvars are active at this time step
  for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
    (*cvi)->feature_states[cvm::deps::f_cv_active]->enabled = (step_absolute() % (*cvi)->get_time_step_factor() == 0);
  }

  // if SMP support is available, split up the work
  if (proxy->smp_enabled() == COLVARS_OK) {

    // first, calculate how much work (currently, how many active CVCs) each colvar has

    colvars_smp.resize(0);
    colvars_smp_items.resize(0);

    colvars_smp.reserve(colvars.size());
    colvars_smp_items.reserve(colvars.size());

    // set up a vector containing all components
    size_t num_colvar_items = 0;
    cvm::increase_depth();
    for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {

      if (!(*cvi)->is_enabled()) continue;
      combine_errors(error_code, (*cvi)->update_cvc_flags());

      size_t num_items = (*cvi)->num_active_cvcs();
      colvars_smp.reserve(colvars_smp.size() + num_items);
      colvars_smp_items.reserve(colvars_smp_items.size() + num_items);
      for (size_t icvc = 0; icvc < num_items; icvc++) {
        colvars_smp.push_back(*cvi);
        colvars_smp_items.push_back(icvc);
      }

      num_colvar_items += num_items;
    }
    cvm::decrease_depth();

    // calculate colvar components in parallel
    combine_errors(error_code, proxy->smp_colvars_loop());

    cvm::increase_depth();
    for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
      if (!(*cvi)->is_enabled()) continue;
      combine_errors(error_code, (*cvi)->collect_cvc_data());
    }
    cvm::decrease_depth();

  } else {

    // calculate colvars one at a time
    cvm::increase_depth();
    for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
      if (!(*cvi)->is_enabled()) continue;
      combine_errors(error_code, (*cvi)->calc());
      if (cvm::get_error()) {
        return COLVARS_ERROR;
      }
    }
    cvm::decrease_depth();
  }

  combine_errors(error_code, cvm::get_error());
  return error_code;
}


int colvarmodule::calc_biases()
{
  // update the biases and communicate their forces to the collective
  // variables
  if (cvm::debug() && biases.size())
    cvm::log("Updating collective variable biases.\n");

  std::vector<colvarbias *>::iterator bi;
  int error_code = COLVARS_OK;

  // if SMP support is available, split up the work
  if (proxy->smp_enabled() == COLVARS_OK) {

    if (use_scripted_forces && !scripting_after_biases) {
      // calculate biases and scripted forces in parallel
      combine_errors(error_code, proxy->smp_biases_script_loop());
    } else {
      // calculate biases in parallel
      combine_errors(error_code, proxy->smp_biases_loop());
    }

  } else {

    if (use_scripted_forces && !scripting_after_biases) {
      combine_errors(error_code, calc_scripted_forces());
    }

    cvm::increase_depth();
    for (bi = biases.begin(); bi != biases.end(); bi++) {
      combine_errors(error_code, (*bi)->update());
      if (cvm::get_error()) {
        return COLVARS_ERROR;
      }
    }
    cvm::decrease_depth();
  }

  cvm::real total_bias_energy = 0.0;
  for (bi = biases.begin(); bi != biases.end(); bi++) {
    total_bias_energy += (*bi)->get_energy();
  }

  proxy->add_energy(total_bias_energy);
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::update_colvar_forces()
{
  int error_code = COLVARS_OK;

  std::vector<colvar *>::iterator cvi;
  std::vector<colvarbias *>::iterator bi;

  // sum the forces from all biases for each collective variable
  if (cvm::debug() && biases.size())
    cvm::log("Collecting forces from all biases.\n");
  cvm::increase_depth();
  for (bi = biases.begin(); bi != biases.end(); bi++) {
    (*bi)->communicate_forces();
    if (cvm::get_error()) {
      return COLVARS_ERROR;
    }
  }
  cvm::decrease_depth();

  if (use_scripted_forces && scripting_after_biases) {
    combine_errors(error_code, calc_scripted_forces());
  }

  cvm::real total_colvar_energy = 0.0;
  // sum up the forces for each colvar, including wall forces
  // and integrate any internal
  // equation of motion (extended system)
  if (cvm::debug())
    cvm::log("Updating the internal degrees of freedom "
             "of colvars (if they have any).\n");
  cvm::increase_depth();
  for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
    // Here we call even inactive colvars, so they accumulate biasing forces
    // as well as update their extended-system dynamics
    total_colvar_energy += (*cvi)->update_forces_energy();
    if (cvm::get_error()) {
      return COLVARS_ERROR;
    }
  }
  cvm::decrease_depth();
  proxy->add_energy(total_colvar_energy);

  // make collective variables communicate their forces to their
  // coupled degrees of freedom (i.e. atoms)
  if (cvm::debug())
    cvm::log("Communicating forces from the colvars to the atoms.\n");
  cvm::increase_depth();
  for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
    if ((*cvi)->is_enabled(cvm::deps::f_cv_gradient)) {
      if (!(*cvi)->is_enabled()) continue;
      (*cvi)->communicate_forces();
      if (cvm::get_error()) {
        return COLVARS_ERROR;
      }
    }
  }
  cvm::decrease_depth();

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::calc_scripted_forces()
{
  // Run user force script, if provided,
  // potentially adding scripted forces to the colvars
  int res;
  res = proxy->run_force_callback();
  if (res == COLVARS_NOT_IMPLEMENTED) {
    cvm::error("Colvar forces scripts are not implemented.");
    return COLVARS_NOT_IMPLEMENTED;
  }
  if (res != COLVARS_OK) {
    cvm::error("Error running user colvar forces script");
    return COLVARS_ERROR;
  }
  return COLVARS_OK;
}


int colvarmodule::write_restart_files()
{
  if ( (cvm::step_relative() > 0) &&
       ((cvm::step_absolute() % restart_out_freq) == 0) ) {
    cvm::log("Writing the state file \""+
             restart_out_name+"\".\n");
    proxy->backup_file(restart_out_name.c_str());
    restart_out_os.open(restart_out_name.c_str());
    if (!restart_out_os.is_open() || !write_restart(restart_out_os))
      cvm::error("Error: in writing restart file.\n");
    restart_out_os.close();
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::write_traj_files()
{
  if (!cv_traj_os.is_open()) {
    open_traj_file(cv_traj_name);
  }

  // write labels in the traj file every 1000 lines and at first timestep
  if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
    write_traj_label(cv_traj_os);
  }

  if ((cvm::step_absolute() % cv_traj_freq) == 0) {
    write_traj(cv_traj_os);
  }

  if (restart_out_freq && cv_traj_os.is_open()) {
    // flush the trajectory file if we are at the restart frequency
    if ( (cvm::step_relative() > 0) &&
         ((cvm::step_absolute() % restart_out_freq) == 0) ) {
      cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+
               cv_traj_name+"\".\n");
      cv_traj_os.flush();
    }
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::analyze()
{
  if (cvm::debug()) {
    cvm::log("colvarmodule::analyze(), step = "+cvm::to_str(it)+".\n");
  }

  if (cvm::step_relative() == 0)
    cvm::log("Performing analysis.\n");

  // perform colvar-specific analysis
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    cvm::increase_depth();
    (*cvi)->analyze();
    cvm::decrease_depth();
  }

  // perform bias-specific analysis
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    cvm::increase_depth();
    (*bi)->analyze();
    cvm::decrease_depth();
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::setup()
{
  // loop over all components of all colvars to reset masses of all groups
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();  cvi++) {
    (*cvi)->setup();
  }
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


colvarmodule::~colvarmodule()
{
  if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) ||
      (proxy->smp_thread_id() == 0)) {
    reset();
    delete parse;
    proxy = NULL;
  }
}


int colvarmodule::reset()
{
  parse->reset();

  cvm::log("Resetting the Collective Variables Module.\n");

  // Iterate backwards because we are deleting the elements as we go
  for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
       bi != biases.rend();
       bi++) {
    delete *bi; // the bias destructor updates the biases array
  }
  biases.clear();

  // Iterate backwards because we are deleting the elements as we go
  for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
       cvi != colvars.rend();
       cvi++) {
    delete *cvi; // the colvar destructor updates the colvars array
  }
  colvars.clear();

  index_groups.clear();
  index_group_names.clear();

  if (cv_traj_os.is_open()) {
    // Do not close file here, as we might not be done with it yet.
    cv_traj_os.flush();
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::setup_input()
{
  // name of input state file
  restart_in_name = proxy->input_prefix().size() ?
    std::string(proxy->input_prefix()+".colvars.state") :
    std::string("") ;

  // read the restart configuration, if available
  if (restart_in_name.size()) {
    // read the restart file
    std::ifstream input_is(restart_in_name.c_str());
    if (!input_is.good()) {
      cvm::error("Error: in opening restart file \""+
                 std::string(restart_in_name)+"\".\n",
                 FILE_ERROR);
      return COLVARS_ERROR;
    } else {
      cvm::log("Restarting from file \""+restart_in_name+"\".\n");
      read_restart(input_is);
      if (cvm::get_error() != COLVARS_OK) {
        return COLVARS_ERROR;
      }
      cvm::log(cvm::line_marker);
    }
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


int colvarmodule::setup_output()
{
  int error_code = 0;

  // output state file (restart)
  restart_out_name = proxy->restart_output_prefix().size() ?
    std::string(proxy->restart_output_prefix()+".colvars.state") :
    std::string("");

  if (restart_out_name.size()) {
    cvm::log("The restart output state file will be \""+restart_out_name+"\".\n");
  }

  output_prefix = proxy->output_prefix();
  if (output_prefix.size()) {
    cvm::log("The final output state file will be \""+
             (output_prefix.size() ?
              std::string(output_prefix+".colvars.state") :
              std::string("colvars.state"))+"\".\n");
    // cvm::log (cvm::line_marker);
  }

  cv_traj_name =
    (output_prefix.size() ?
     std::string(output_prefix+".colvars.traj") :
     std::string(""));

  if (cv_traj_freq && cv_traj_name.size()) {
    combine_errors(error_code, open_traj_file(cv_traj_name));
  }

  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    combine_errors(error_code, (*bi)->setup_output());
  }

  if (error_code != COLVARS_OK || cvm::get_error()) {
    set_error_bit(FILE_ERROR);
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


std::istream & colvarmodule::read_restart(std::istream &is)
{
  {
    // read global restart information
    std::string restart_conf;
    if (is >> colvarparse::read_block("configuration", restart_conf)) {
      if (it_restart_from_state_file) {
        parse->get_keyval(restart_conf, "step",
                          it_restart, (size_t) 0,
                          colvarparse::parse_silent);
        it = it_restart;
      }
      std::string restart_version;
      parse->get_keyval(restart_conf, "version",
                        restart_version, std::string(""),
                        colvarparse::parse_silent);
      if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) {
        cvm::log("This state file was generated with version "+restart_version+"\n");
      }
    }
    is.clear();
    parse->clear_keyword_registry();
  }

  // colvars restart
  cvm::increase_depth();
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    if ( !((*cvi)->read_restart(is)) ) {
      cvm::error("Error: in reading restart configuration for collective variable \""+
                 (*cvi)->name+"\".\n",
                 INPUT_ERROR);
    }
  }

  // biases restart
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    if (!((*bi)->read_restart(is))) {
      cvm::error("Error: in reading restart configuration for bias \""+
                 (*bi)->name+"\".\n",
                 INPUT_ERROR);
    }
  }
  cvm::decrease_depth();

  return is;
}


int colvarmodule::backup_file(char const *filename)
{
  return proxy->backup_file(filename);
}


int colvarmodule::write_output_files()
{
  // if this is a simulation run (i.e. not a postprocessing), output data
  // must be written to be able to restart the simulation
  std::string const out_name =
    (output_prefix.size() ?
     std::string(output_prefix+".colvars.state") :
     std::string("colvars.state"));
  cvm::log("Saving collective variables state to \""+out_name+"\".\n");

  std::ostream * os = proxy->output_stream(out_name);
  os->setf(std::ios::scientific, std::ios::floatfield);
  this->write_restart(*os);
  proxy->close_output_stream(out_name);

  cvm::increase_depth();
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    (*cvi)->write_output_files();
  }
  cvm::decrease_depth();

  cvm::increase_depth();
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    (*bi)->write_output_files();
  }
  cvm::decrease_depth();

  if (cv_traj_os.is_open()) {
    // do not close to avoid problems with multiple NAMD runs
    cv_traj_os.flush();
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}



int colvarmodule::read_traj(char const *traj_filename,
                            long        traj_read_begin,
                            long        traj_read_end)
{
  cvm::log("Opening trajectory file \""+
           std::string(traj_filename)+"\".\n");
  std::ifstream traj_is(traj_filename);

  while (true) {
    while (true) {

      std::string line("");

      do {
        if (!colvarparse::getline_nocomments(traj_is, line)) {
          cvm::log("End of file \""+std::string(traj_filename)+
                   "\" reached, or corrupted file.\n");
          traj_is.close();
          return false;
        }
      } while (line.find_first_not_of(colvarparse::white_space) == std::string::npos);

      std::istringstream is(line);

      if (!(is >> it)) return false;

      if ( (it < traj_read_begin) ) {

        if ((it % 1000) == 0)
          std::cerr << "Skipping trajectory step " << it
                    << "                    \r";

        continue;

      } else {

        if ((it % 1000) == 0)
          std::cerr << "Reading from trajectory, step = " << it
                    << "                    \r";

        if ( (traj_read_end > traj_read_begin) &&
             (it > traj_read_end) ) {
          std::cerr << "\n";
          cvm::error("Reached the end of the trajectory, "
                     "read_end = "+cvm::to_str(traj_read_end)+"\n",
                     FILE_ERROR);
          return COLVARS_ERROR;
        }

        for (std::vector<colvar *>::iterator cvi = colvars.begin();
             cvi != colvars.end();
             cvi++) {
          if (!(*cvi)->read_traj(is)) {
            cvm::error("Error: in reading colvar \""+(*cvi)->name+
                       "\" from trajectory file \""+
                       std::string(traj_filename)+"\".\n",
                       FILE_ERROR);
            return COLVARS_ERROR;
          }
        }

        break;
      }
    }
  }
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}


std::ostream & colvarmodule::write_restart(std::ostream &os)
{
  os.setf(std::ios::scientific, std::ios::floatfield);
  os << "configuration {\n"
     << "  step " << std::setw(it_width)
     << it << "\n"
     << "  dt " << dt() << "\n"
     << "  version " << std::string(COLVARS_VERSION) << "\n"
     << "}\n\n";

  cvm::increase_depth();
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    (*cvi)->write_restart(os);
  }

  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    (*bi)->write_restart(os);
  }
  cvm::decrease_depth();

  return os;
}

int colvarmodule::open_traj_file(std::string const &file_name)
{
  if (cv_traj_os.is_open()) {
    return COLVARS_OK;
  }

  // (re)open trajectory file
  if (cv_traj_append) {
    cvm::log("Appending to colvar trajectory file \""+file_name+
             "\".\n");
    cv_traj_os.open(file_name.c_str(), std::ios::app);
  } else {
    cvm::log("Writing to colvar trajectory file \""+file_name+
             "\".\n");
    proxy->backup_file(file_name.c_str());
    cv_traj_os.open(file_name.c_str());
  }

  if (!cv_traj_os.is_open()) {
    cvm::error("Error: cannot write to file \""+file_name+"\".\n",
               FILE_ERROR);
  }

  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}

int colvarmodule::close_traj_file()
{
  if (cv_traj_os.is_open()) {
    cv_traj_os.close();
  }
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}

std::ostream & colvarmodule::write_traj_label(std::ostream &os)
{
  if (!os.good()) {
    cvm::error("Cannot write to trajectory file.");
    return os;
  }
  os.setf(std::ios::scientific, std::ios::floatfield);

  os << "# " << cvm::wrap_string("step", cvm::it_width-2)
     << " ";

  cvm::increase_depth();
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    (*cvi)->write_traj_label(os);
  }
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    (*bi)->write_traj_label(os);
  }
  os << "\n";
  if (cvm::debug()) {
    os.flush();
  }
  cvm::decrease_depth();
  return os;
}

std::ostream & colvarmodule::write_traj(std::ostream &os)
{
  os.setf(std::ios::scientific, std::ios::floatfield);

  os << std::setw(cvm::it_width) << it
     << " ";

  cvm::increase_depth();
  for (std::vector<colvar *>::iterator cvi = colvars.begin();
       cvi != colvars.end();
       cvi++) {
    (*cvi)->write_traj(os);
  }
  for (std::vector<colvarbias *>::iterator bi = biases.begin();
       bi != biases.end();
       bi++) {
    (*bi)->write_traj(os);
  }
  os << "\n";
  if (cvm::debug()) {
    os.flush();
  }
  cvm::decrease_depth();
  return os;
}


void cvm::log(std::string const &message)
{
  size_t const d = depth();
  if (d > 0)
    proxy->log((std::string(2*d, ' '))+message);
  else
    proxy->log(message);
}


void cvm::increase_depth()
{
  (depth())++;
}


void cvm::decrease_depth()
{
  if (depth() > 0) {
    (depth())--;
  }
}


size_t & cvm::depth()
{
  // NOTE: do not call log() or error() here, to avoid recursion
  size_t const nt = proxy->smp_num_threads();
  if (proxy->smp_enabled() == COLVARS_OK) {
    if (depth_v.size() != nt) {
      // update array of depths
      proxy->smp_lock();
      if (depth_v.size() > 0) { depth_s = depth_v[0]; }
      depth_v.clear();
      depth_v.assign(nt, depth_s);
      proxy->smp_unlock();
    }
    return depth_v[proxy->smp_thread_id()];
  }
  return depth_s;
}


void colvarmodule::set_error_bit(int code)
{
  if (code > 0) {
    cvm::fatal_error("Error: set_error_bit() received positive error code.\n");
    return;
  }
  proxy->smp_lock();
  errorCode = -1 * ((-1 * errorCode) | (-1 * code) | (-1 * COLVARS_ERROR));
  proxy->smp_unlock();
}

bool colvarmodule::get_error_bit(int code)
{
  return bool((-1 * errorCode) & (-1 * code));
}

void colvarmodule::combine_errors(int &target, int const code)
{
  if (code > 0 || target > 0) {
    cvm::fatal_error("Error: combine_errors() received positive error code.\n");
    return;
  }
  target = -1 * ((-1 * target) | (-1 * code));
}

void colvarmodule::clear_error()
{
  proxy->smp_lock();
  errorCode = COLVARS_OK;
  proxy->smp_unlock();
}


void cvm::error(std::string const &message, int code)
{
  set_error_bit(code);
  proxy->error(message);
}


void cvm::fatal_error(std::string const &message)
{
  // TODO once all non-fatal errors have been set to be handled by error(),
  // set DELETE_COLVARS here for VMD to handle it
  set_error_bit(FATAL_ERROR);
  proxy->fatal_error(message);
}


void cvm::exit(std::string const &message)
{
  proxy->exit(message);
}


int cvm::read_index_file(char const *filename)
{
  std::ifstream is(filename, std::ios::binary);
  if (!is.good()) {
    cvm::error("Error: in opening index file \""+
               std::string(filename)+"\".\n",
               FILE_ERROR);
  }

  while (is.good()) {
    char open, close;
    std::string group_name;
    if ( (is >> open) && (open == '[') &&
         (is >> group_name) &&
         (is >> close) && (close == ']') ) {
      for (std::list<std::string>::iterator names_i = index_group_names.begin();
           names_i != index_group_names.end();
           names_i++) {
        if (*names_i == group_name) {
          cvm::error("Error: the group name \""+group_name+
                     "\" appears in multiple index files.\n",
                     FILE_ERROR);
        }
      }
      cvm::index_group_names.push_back(group_name);
      cvm::index_groups.push_back(std::vector<int> ());
    } else {
      cvm::error("Error: in parsing index file \""+
                 std::string(filename)+"\".\n",
                 INPUT_ERROR);
    }

    int atom_number = 1;
    size_t pos = is.tellg();
    while ( (is >> atom_number) && (atom_number > 0) ) {
      (cvm::index_groups.back()).push_back(atom_number);
      pos = is.tellg();
    }
    is.clear();
    is.seekg(pos, std::ios::beg);
    std::string delim;
    if ( (is >> delim) && (delim == "[") ) {
      // new group
      is.clear();
      is.seekg(pos, std::ios::beg);
    } else {
      break;
    }
  }

  cvm::log("The following index groups were read from the index file \""+
           std::string(filename)+"\":\n");
  std::list<std::string>::iterator names_i = index_group_names.begin();
  std::list<std::vector<int> >::iterator lists_i = index_groups.begin();
  for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) {
    cvm::log("  "+(*names_i)+" ("+cvm::to_str(lists_i->size())+" atoms).\n");
  }
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}

int cvm::load_atoms(char const *file_name,
                    cvm::atom_group &atoms,
                    std::string const &pdb_field,
                    double const pdb_field_value)
{
  return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value);
}

int cvm::load_coords(char const *file_name,
                     std::vector<cvm::atom_pos> &pos,
                     const std::vector<int> &indices,
                     std::string const &pdb_field,
                     double const pdb_field_value)
{
  // Differentiate between PDB and XYZ files
  // for XYZ files, use CVM internal parser
  // otherwise call proxy function for PDB

  std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name);
  if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) {
    if ( pdb_field.size() > 0 ) {
      cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR);
      return COLVARS_ERROR;
    }
    return cvm::load_coords_xyz(file_name, pos, indices);
  } else {
    return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value);
  }
}


int cvm::load_coords_xyz(char const *filename,
                         std::vector<atom_pos> &pos,
                         const std::vector<int> &indices)
{
  std::ifstream xyz_is(filename);
  unsigned int natoms;
  char symbol[256];
  std::string line;

  if ( ! (xyz_is >> natoms) ) {
    cvm::error("Error: cannot parse XYZ file "
               + std::string(filename) + ".\n", INPUT_ERROR);
  }
  // skip comment line
  std::getline(xyz_is, line);
  std::getline(xyz_is, line);
  xyz_is.width(255);
  std::vector<atom_pos>::iterator pos_i = pos.begin();

  if (pos.size() != natoms) { // Use specified indices
    int next = 0; // indices are zero-based
    std::vector<int>::const_iterator index = indices.begin();
    for ( ; pos_i != pos.end() ; pos_i++, index++) {

      while ( next < *index ) {
        std::getline(xyz_is, line);
        next++;
      }
      xyz_is >> symbol;
      xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
    }
  } else {          // Use all positions
    for ( ; pos_i != pos.end() ; pos_i++) {
      xyz_is >> symbol;
      xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
    }
  }
  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}

// static pointers
std::vector<colvar *>     colvarmodule::colvars;
std::vector<colvarbias *> colvarmodule::biases;
size_t                    colvarmodule::n_abf_biases = 0;
size_t                    colvarmodule::n_rest_biases = 0;
size_t                    colvarmodule::n_histo_biases = 0;
size_t                    colvarmodule::n_meta_biases = 0;
colvarproxy              *colvarmodule::proxy = NULL;


// static runtime data
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
int       colvarmodule::errorCode = 0;
long      colvarmodule::it = 0;
long      colvarmodule::it_restart = 0;
size_t    colvarmodule::restart_out_freq = 0;
size_t    colvarmodule::cv_traj_freq = 0;
size_t    colvarmodule::depth_s = 0;
std::vector<size_t> colvarmodule::depth_v(0);
bool      colvarmodule::b_analysis = false;
std::list<std::string> colvarmodule::index_group_names;
std::list<std::vector<int> > colvarmodule::index_groups;
bool     colvarmodule::use_scripted_forces = false;
bool     colvarmodule::scripting_after_biases = true;

// file name prefixes
std::string colvarmodule::output_prefix = "";
std::string colvarmodule::restart_in_name = "";


// i/o constants
size_t const colvarmodule::it_width = 12;
size_t const colvarmodule::cv_prec  = 14;
size_t const colvarmodule::cv_width = 21;
size_t const colvarmodule::en_prec  = 14;
size_t const colvarmodule::en_width = 21;
std::string const colvarmodule::line_marker =
  "----------------------------------------------------------------------\n";