git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14703 f3b2605a-c512-4ea7-a41b-209d697bcdaa

lib/colvars/colvar.cpp

@@ -128,6 +128,7 @@ colvar::colvar(std::string const &conf)
                         "number",           "selfCoordNum",   selfcoordnum);
 
   initialize_components("angle",            "angle",          angle);
+  initialize_components("dipole angle",     "dipoleAngle",    dipole_angle);
   initialize_components("dihedral",         "dihedral",       dihedral);
 
   initialize_components("hydrogen bond",    "hBond",          h_bond);
@@ -520,7 +521,8 @@ void colvar::build_atom_list(void)
   for (size_t i = 0; i < cvcs.size(); i++) {
     for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
       for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
-        temp_id_list.push_back(cvcs[i]->atom_groups[j]->at(k).id);
+        cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
+        temp_id_list.push_back(ag[k].id);
       }
     }
   }
@@ -830,8 +832,9 @@ void colvar::setup() {
   for (size_t i = 0; i < cvcs.size(); i++) {
     for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
       cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
-      atoms.read_positions();
+      atoms.setup();
       atoms.reset_mass(name,i,ig);
+      atoms.read_positions();
     }
   }
 }
@@ -861,37 +864,40 @@ colvar::~colvar()
 
 void colvar::calc()
 {
-  size_t i, ig;
   if (cvm::debug())
     cvm::log("Calculating colvar \""+this->name+"\".\n");
+  update_cvc_flags();
+  calc_cvcs();
+  calc_colvar_properties();
+}
+
+
+int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
+{
+  size_t i, ig;
+  size_t cvc_count;
+  size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs();
 
   // prepare atom groups for calculation
   if (cvm::debug())
     cvm::log("Collecting data from atom groups.\n");
 
-  // Update the enabled/disabled status of cvcs if necessary
+  if (first_cvc >= cvcs.size()) {
-  if (cvc_flags.size()) {
+    cvm::error("Error: trying to address a component outside the "
-    bool any = false;
+               "range defined for colvar \""+name+"\".\n", BUG_ERROR);
-    for (i = 0; i < cvcs.size(); i++) {
+    return BUG_ERROR;
-      cvcs[i]->b_enabled = cvc_flags[i];
-      any = any || cvc_flags[i];
-    }
-    if (!any) {
-      cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
-      return;
-    }
-    cvc_flags.resize(0);
   }
 
-  for (i = 0; i < cvcs.size(); i++) {
+  for (i = first_cvc, cvc_count = 0;
+       (i < cvcs.size()) && (cvc_count < cvc_max_count);
+       i++) {
     if (!cvcs[i]->b_enabled) continue;
+    cvc_count++;
     for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
       cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
       atoms.reset_atoms_data();
       atoms.read_positions();
-      if (atoms.b_center || atoms.b_rotate) {
+      atoms.calc_required_properties();
-        atoms.calc_apply_roto_translation();
-      }
       // each atom group will take care of its own ref_pos_group, if defined
     }
   }
@@ -906,7 +912,11 @@ void colvar::calc()
 //   }
 
   if (tasks[task_system_force]) {
-    for (i = 0; i < cvcs.size(); i++) {
+    for (i = first_cvc, cvc_count = 0;
+         (i < cvcs.size()) && (cvc_count < cvc_max_count);
+         i++) {
+      if (!cvcs[i]->b_enabled) continue;
+      cvc_count++;
       for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
         cvcs[i]->atom_groups[ig]->read_system_forces();
       }
@@ -919,9 +929,12 @@ void colvar::calc()
     cvm::log("Calculating colvar components.\n");
   x.reset();
 
-  // First, update component values
+  // First, calculate component values
-  for (i = 0; i < cvcs.size(); i++) {
+  for (i = first_cvc, cvc_count = 0;
+       (i < cvcs.size()) && (cvc_count < cvc_max_count);
+       i++) {
     if (!cvcs[i]->b_enabled) continue;
+    cvc_count++;
     cvm::increase_depth();
     (cvcs[i])->calc_value();
     cvm::decrease_depth();
@@ -932,22 +945,25 @@ void colvar::calc()
                 cvm::cv_width, cvm::cv_prec)+".\n");
   }
 
-  // Then combine them appropriately
+  // Then combine them appropriately, using either a scripted function or a polynomial
   if (tasks[task_scripted]) {
     // cvcs combined by user script
     int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x);
     if (res == COLVARS_NOT_IMPLEMENTED) {
       cvm::error("Scripted colvars are not implemented.");
-      return;
+      return COLVARS_NOT_IMPLEMENTED;
     }
     if (res != COLVARS_OK) {
       cvm::error("Error running scripted colvar");
-      return;
+      return COLVARS_OK;
     }
   } else if (x.type() == colvarvalue::type_scalar) {
     // polynomial combination allowed
-    for (i = 0; i < cvcs.size(); i++) {
+    for (i = first_cvc, cvc_count = 0;
+        (i < cvcs.size()) && (cvc_count < cvc_max_count);
+        i++) {
       if (!cvcs[i]->b_enabled) continue;
+      cvc_count++;
       x += (cvcs[i])->sup_coeff *
       ( ((cvcs[i])->sup_np != 1) ?
         std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
@@ -955,8 +971,11 @@ void colvar::calc()
     }
   } else {
     // only linear combination allowed
-    for (i = 0; i < cvcs.size(); i++) {
+    for (i = first_cvc, cvc_count = 0;
+        (i < cvcs.size()) && (cvc_count < cvc_max_count);
+        i++) {
       if (!cvcs[i]->b_enabled) continue;
+      cvc_count++;
       x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
     }
   }
@@ -970,9 +989,12 @@ void colvar::calc()
     if (cvm::debug())
       cvm::log("Calculating gradients of colvar \""+this->name+"\".\n");
 
-    for (i = 0; i < cvcs.size(); i++) {
     // calculate the gradients of each component
+    for (i = first_cvc, cvc_count = 0;
+        (i < cvcs.size()) && (cvc_count < cvc_max_count);
+        i++) {
       if (!cvcs[i]->b_enabled) continue;
+      cvc_count++;
       cvm::increase_depth();
       (cvcs[i])->calc_gradients();
       // if requested, propagate (via chain rule) the gradients above
@@ -991,40 +1013,44 @@ void colvar::calc()
 
       if (tasks[task_scripted]) {
         cvm::error("Collecting atomic gradients is not implemented for "
-          "scripted colvars.");
+                   "scripted colvars.", COLVARS_NOT_IMPLEMENTED);
-        return;
+        return COLVARS_NOT_IMPLEMENTED;
       }
 
       // Collect the atomic gradients inside colvar object
       for (unsigned int a = 0; a < atomic_gradients.size(); a++) {
         atomic_gradients[a].reset();
       }
-      for (i = 0; i < cvcs.size(); i++) {
+      for (i = first_cvc, cvc_count = 0;
+          (i < cvcs.size()) && (cvc_count < cvc_max_count);
+          i++) {
         if (!cvcs[i]->b_enabled) continue;
+        cvc_count++;
         // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx
         cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) *
           std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1);
 
         for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) {
 
+          cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]);
+
           // If necessary, apply inverse rotation to get atomic
           // gradient in the laboratory frame
           if (cvcs[i]->atom_groups[j]->b_rotate) {
             cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse();
 
             for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
-              int a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+              size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
-                  cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
+                                          ag[k].id) - atom_ids.begin();
-              atomic_gradients[a] += coeff *
+              atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad);
-                rot_inv.rotate(cvcs[i]->atom_groups[j]->at(k).grad);
             }
 
           } else {
 
             for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
-              int a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
+              size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(),
-                  cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
+                                          ag[k].id) - atom_ids.begin();
-              atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad;
+              atomic_gradients[a] += coeff * ag[k].grad;
             }
           }
         }
@@ -1040,8 +1066,8 @@ void colvar::calc()
       // TODO see if this could reasonably be done in a generic way
       // from generic inverse gradients
       cvm::error("System force is not implemented for "
-        "scripted colvars.");
+                 "scripted colvars.", COLVARS_NOT_IMPLEMENTED);
-      return;
+      return COLVARS_NOT_IMPLEMENTED;
     }
     if (cvm::debug())
       cvm::log("Calculating system force of colvar \""+this->name+"\".\n");
@@ -1054,7 +1080,11 @@ void colvar::calc()
 
     if (cvm::step_relative() > 0) {
       // get from the cvcs the system forces from the PREVIOUS step
-      for (i = 0; i < cvcs.size(); i++) {
+      for (i = first_cvc, cvc_count = 0;
+          (i < cvcs.size()) && (cvc_count < cvc_max_count);
+          i++) {
+        if (!cvcs[i]->b_enabled) continue;
+        cvc_count++;
         (cvcs[i])->calc_force_invgrads();
         // linear combination is assumed
         cvm::increase_depth();
@@ -1072,7 +1102,12 @@ void colvar::calc()
     if (cvm::debug())
       cvm::log("Done calculating system force of colvar \""+this->name+"\".\n");
   }
+  return COLVARS_OK;
+}
 
+
+int colvar::calc_colvar_properties()
+{
   if (tasks[task_fdiff_velocity]) {
     // calculate the velocity by finite differences
     if (cvm::step_relative() == 0)
@@ -1109,6 +1144,8 @@ void colvar::calc()
 
   if (cvm::debug())
     cvm::log("Done calculating colvar \""+this->name+"\".\n");
+
+  return COLVARS_OK;
 }
 
 
@@ -1243,6 +1280,7 @@ void colvar::communicate_forces()
 
   if (tasks[task_scripted]) {
     std::vector<cvm::matrix2d<cvm::real> > func_grads;
+    func_grads.reserve(cvcs.size());
     for (i = 0; i < cvcs.size(); i++) {
       if (!cvcs[i]->b_enabled) continue;
       func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
@@ -1252,7 +1290,7 @@ void colvar::communicate_forces()
 
     if (res != COLVARS_OK) {
       if (res == COLVARS_NOT_IMPLEMENTED) {
-        cvm::error("Colvar gradient scripts are not implemented.");
+        cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED);
       } else {
         cvm::error("Error running colvar gradient script");
       }
@@ -1296,8 +1334,8 @@ void colvar::communicate_forces()
 }
 
 
-int colvar::set_cvc_flags(std::vector<bool> const &flags) {
+int colvar::set_cvc_flags(std::vector<bool> const &flags)
-
+{
   if (flags.size() != cvcs.size()) {
     cvm::error("ERROR: Wrong number of CVC flags provided.");
     return COLVARS_ERROR;
@@ -1309,6 +1347,36 @@ int colvar::set_cvc_flags(std::vector<bool> const &flags) {
 }
 
 
+int colvar::update_cvc_flags()
+{
+  size_t i;
+  // Update the enabled/disabled status of cvcs if necessary
+  if (cvc_flags.size()) {
+    bool any = false;
+    for (i = 0; i < cvcs.size(); i++) {
+      cvcs[i]->b_enabled = cvc_flags[i];
+      any = any || cvc_flags[i];
+    }
+    if (!any) {
+      cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n");
+      return COLVARS_ERROR;
+    }
+    cvc_flags.resize(0);
+  }
+  return COLVARS_OK;
+}
+
+
+int colvar::num_active_cvcs() const
+{
+  int result = 0;
+  for (size_t i = 0; i < cvcs.size(); i++) {
+    if (cvcs[i]->b_enabled) result++;
+  }
+  return result;
+}
+
+
 // ******************** METRIC FUNCTIONS ********************
 // Use the metrics defined by \link cvc \endlink objects
 

lib/colvars/colvar.h

@@ -310,8 +310,11 @@ public:
   /// quantities
   void calc();
 
-  /// \brief Calculate just the value, and store it in the argument
+  /// \brief Calculate a given subset of colvar components (CVCs) (default: all CVCs)
-  void calc_value(colvarvalue &xn);
+  int calc_cvcs(int first = 0, size_t num_cvcs = 0);
+
+  /// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
+  int calc_colvar_properties();
 
   /// Set the total biasing force to zero
   void reset_bias_force();
@@ -332,6 +335,12 @@ public:
   /// \brief Enables and disables individual CVCs based on flags
   int set_cvc_flags(std::vector<bool> const &flags);
 
+  /// \brief Updates the flags in the CVC objects
+  int update_cvc_flags();
+
+  /// \brief Return the number of CVC objects with an active flag
+  int num_active_cvcs() const;
+
   /// \brief Use the internal metrics (as from \link cvc
   /// \endlink objects) to calculate square distances and gradients
   ///
@@ -486,6 +495,7 @@ public:
   class distance_inv;
   class distance_pairs;
   class angle;
+  class dipole_angle;
   class dihedral;
   class coordnum;
   class selfcoordnum;

lib/colvars/colvaratoms.h

@@ -6,43 +6,44 @@
 #include "colvarmodule.h"
 #include "colvarparse.h"
 
+
 /// \brief Stores numeric id, mass and all mutable data for an atom,
 /// mostly used by a \link cvc \endlink
 ///
-/// This class may be used (although not necessarily) to keep atomic
+/// This class may be used to keep atomic data such as id, mass,
-/// data (id, mass, position and collective variable derivatives)
+/// position and collective variable derivatives) altogether.
-/// altogether.  There may be multiple instances with identical
+/// There may be multiple instances with identical
 /// numeric id, all acting independently: forces communicated through
 /// these instances will be summed together.
-///
+
-/// Read/write operations depend on the underlying code: hence, some
-/// member functions are defined in colvarproxy_xxx.h.
 class colvarmodule::atom {
 
 protected:
 
-  /// \brief Index in the list of atoms involved by the colvars (\b
+  /// Index in the colvarproxy arrays (\b NOT in the global topology!)
-  /// NOT in the global topology!)
   int index;
 
 public:
 
-  /// Internal identifier (zero-based)
+  /// Identifier for the MD program (0-based)
   int             id;
 
   /// Mass
   cvm::real       mass;
 
+  /// Charge
+  cvm::real       charge;
+
   /// \brief Current position (copied from the program, can be
-  /// manipulated)
+  /// modified if necessary)
   cvm::atom_pos   pos;
 
   /// \brief Current velocity (copied from the program, can be
-  /// manipulated)
+  /// modified if necessary)
   cvm::rvector    vel;
 
   /// \brief System force at the previous step (copied from the
-  /// program, can be manipulated)
+  /// program, can be modified if necessary)
   cvm::rvector    system_force;
 
   /// \brief Gradient of a scalar collective variable with respect
@@ -57,13 +58,13 @@ public:
   /// implementation
   cvm::rvector   grad;
 
-  /// \brief Default constructor, setting index and id to invalid numbers
+  /// \brief Default constructor (sets index and id both to -1)
-  atom() : index(-1), id(-1) { reset_data(); }
+  atom();
 
   /// \brief Initialize an atom for collective variable calculation
   /// and get its internal identifier \param atom_number Atom index in
   /// the system topology (starting from 1)
-  atom(int const &atom_number);
+  atom(int atom_number);
 
   /// \brief Initialize an atom for collective variable calculation
   /// and get its internal identifier \param residue Residue number
@@ -72,7 +73,7 @@ public:
   /// type of topologies, may not be required
   atom(cvm::residue_id const &residue,
        std::string const     &atom_name,
-        std::string const     &segment_id = std::string(""));
+       std::string const     &segment_id);
 
   /// Copy constructor
   atom(atom const &a);
@@ -80,57 +81,176 @@ public:
   /// Destructor
   ~atom();
 
-  /// Set non-constant data (everything except id and mass) to zero
+  /// Set mutable data (everything except id and mass) to zero; update mass
-  inline void reset_data() {
+  inline void reset_data()
-    pos = atom_pos(0.0);
+  {
-    vel = grad = system_force = rvector(0.0);
+    pos = cvm::atom_pos(0.0);
+    vel = grad = system_force = cvm::rvector(0.0);
+  }
+
+  /// Get the latest value of the mass
+  inline void update_mass()
+  {
+    mass = (cvm::proxy)->get_atom_mass(index);
+  }
+
+  /// Get the latest value of the charge
+  inline void update_charge()
+  {
+    charge = (cvm::proxy)->get_atom_charge(index);
   }
 
   /// Get the current position
-  void read_position();
+  inline void read_position()
+  {
+    pos = (cvm::proxy)->get_atom_position(index);
+  }
 
   /// Get the current velocity
-  void read_velocity();
+  inline void read_velocity()
+  {
+    vel = (cvm::proxy)->get_atom_velocity(index);
+  }
 
   /// Get the system force
-  void read_system_force();
+  inline void read_system_force()
+  {
+    system_force = (cvm::proxy)->get_atom_system_force(index);
+  }
 
   /// \brief Apply a force to the atom
   ///
-  /// The force will be used later by the MD integrator, the
+  /// Note: the force is not applied instantly, but will be used later
-  /// collective variables module does not integrate equations of
+  /// by the MD integrator (the colvars module does not integrate
-  /// motion.  Multiple calls to this function by either the same
+  /// equations of motion.
+  ///
+  /// Multiple calls to this function by either the same
   /// \link atom \endlink object or different objects with identical
-  /// \link id \endlink, will all add to the existing MD force.
+  /// \link id \endlink will all be added together.
-  void apply_force(cvm::rvector const &new_force);
+  inline void apply_force(cvm::rvector const &new_force) const
+  {
+    (cvm::proxy)->apply_atom_force(index, new_force);
+  }
 };
 
 
 
-
 /// \brief Group of \link atom \endlink objects, mostly used by a
-/// \link cvc \endlink
+/// \link cvc \endlink object to gather all atomic data
-///
-/// This class inherits from \link colvarparse \endlink and from
-/// std::vector<colvarmodule::atom>, and hence all functions and
-/// operators (including the bracket operator, group[i]) can be used
-/// on an \link atom_group \endlink object.  It can be initialized as
-/// a vector, or by parsing a keyword in the configuration.
 class colvarmodule::atom_group
-  : public std::vector<cvm::atom>,
+  : public colvarparse
-    public colvarparse
 {
 public:
-  // Note: all members here are kept public, to allow any
-  // object accessing and manipulating them
 
+  /// \brief Initialize the group by looking up its configuration
+  /// string in conf and parsing it; this is actually done by parse(),
+  /// which is a member function so that a group can be initialized
+  /// also after construction
+  atom_group(std::string const &conf,
+             char const        *key);
+
+  /// \brief Keyword used to define the group
+  // TODO Make this field part of the data structures that link a group to a CVC
+  std::string key;
+
+  /// \brief Set default values for common flags
+  int init();
+
+  /// \brief Update data required to calculate cvc's
+  int setup();
+
+  /// \brief Initialize the group by looking up its configuration
+  /// string in conf and parsing it
+  int parse(std::string const &conf);
+
+  int add_atom_numbers(std::string const &numbers_conf);
+  int add_index_group(std::string const &index_group_name);
+  int add_atom_numbers_range(std::string const &range_conf);
+  int add_atom_name_residue_range(std::string const &psf_segid,
+                                  std::string const &range_conf);
+  int parse_fitting_options(std::string const &group_conf);
+
+  /// \brief Initialize the group after a (temporary) vector of atoms
+  atom_group(std::vector<cvm::atom> const &atoms_in);
+
+  /// \brief Add an atom object to this group
+  int add_atom(cvm::atom const &a);
+
+  /// \brief Add an atom ID to this group (the actual atomicdata will be not be handled by the group)
+  int add_atom_id(int aid);
+
+  /// \brief Remove an atom object from this group
+  int remove_atom(cvm::atom_iter ai);
+
+  /// \brief Re-initialize the total mass of a group.
+  /// This is needed in case the hosting MD code has an option to
+  /// change atom masses after their initialization.
+  void reset_mass(std::string &name, int i, int j);
+
+  /// \brief Whether or not the properties of this group will be computed in parallel
+  bool b_scalable;
+
+  /// \brief Default constructor
+  atom_group();
+
+  /// \brief Destructor
+  ~atom_group();
+
+protected:
+
+  /// \brief Array of atom objects
+  std::vector<cvm::atom> atoms;
+
+  /// \brief Array of atom identifiers for the MD program (0-based)
+  std::vector<int> atoms_ids;
+
+  /// \brief Dummy atom position
+  cvm::atom_pos dummy_atom_pos;
+
+  /// \brief Index in the colvarproxy arrays (if the group is scalable)
+  int index;
+
+public:
+
+  inline cvm::atom & operator [] (size_t const i)
+  {
+    return atoms[i];
+  }
+
+  inline cvm::atom const & operator [] (size_t const i) const
+  {
+    return atoms[i];
+  }
+
+  inline cvm::atom_iter begin()
+  {
+    return atoms.begin();
+  }
+
+  inline cvm::atom_const_iter begin() const
+  {
+    return atoms.begin();
+  }
+
+  inline cvm::atom_iter end()
+  {
+    return atoms.end();
+  }
+
+  inline cvm::atom_const_iter end() const
+  {
+    return atoms.end();
+  }
+
+  inline size_t size() const
+  {
+    return atoms.size();
+  }
 
   /// \brief If this option is on, this group merely acts as a wrapper
   /// for a fixed position; any calls to atoms within or to
   /// functions that return disaggregated data will fail
   bool b_dummy;
-  /// \brief dummy atom position
-  cvm::atom_pos dummy_atom_pos;
 
   /// Sorted list of zero-based (internal) atom ids
   /// (populated on-demand by create_sorted_ids)
@@ -177,56 +297,35 @@ public:
 
   /// Total mass of the atom group
   cvm::real total_mass;
+  void update_total_mass();
+
+  /// Total charge of the atom group
+  cvm::real total_charge;
+  void update_total_charge();
 
   /// \brief Don't apply any force on this group (use its coordinates
   /// only to calculate a colvar)
   bool        noforce;
 
-
+  /// \brief Get the current positions
-  /// \brief Initialize the group by looking up its configuration
-  /// string in conf and parsing it; this is actually done by parse(),
-  /// which is a member function so that a group can be initialized
-  /// also after construction
-  atom_group(std::string const &conf,
-              char const        *key);
-
-  /// \brief Initialize the group by looking up its configuration
-  /// string in conf and parsing it
-  int parse(std::string const &conf,
-              char const        *key);
-
-  /// \brief Initialize the group after a temporary vector of atoms
-  atom_group(std::vector<cvm::atom> const &atoms);
-
-  /// \brief Add an atom to this group
-  void add_atom(cvm::atom const &a);
-
-  /// \brief Re-initialize the total mass of a group.
-  /// This is needed in case the hosting MD code has an option to
-  /// change atom masses after their initialization.
-  void reset_mass(std::string &name, int i, int j);
-
-  /// \brief Default constructor
-  atom_group();
-
-  /// \brief Destructor
-  ~atom_group();
-
-  /// \brief Get the current positions; if b_center or b_rotate are
-  /// true, calc_apply_roto_translation() will be called too
   void read_positions();
 
   /// \brief (Re)calculate the optimal roto-translation
   void calc_apply_roto_translation();
 
-  /// \brief Save center of geometry fo ref positions,
+  /// \brief Save aside the center of geometry of the reference positions,
-  /// then subtract it
+  /// then subtract it from them
+  ///
+  /// In this way it will be possible to use ref_pos also for the
+  /// rotational fit.
+  /// This is called either by atom_group::parse or by CVCs that assign
+  /// reference positions (eg. RMSD, eigenvector).
   void center_ref_pos();
 
   /// \brief Move all positions
   void apply_translation(cvm::rvector const &t);
 
-  /// \brief Rotate all positions
+  /// \brief Rotate all positions around the center of geometry
   void apply_rotation(cvm::rotation const &q);
 
 
@@ -243,36 +342,68 @@ public:
   /// Call reset_data() for each atom
   inline void reset_atoms_data()
   {
-    for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++)
+    for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++)
       ai->reset_data();
     if (ref_pos_group)
       ref_pos_group->reset_atoms_data();
   }
 
+  /// \brief Recompute all mutable quantities that are required to compute CVCs
+  int calc_required_properties();
+
   /// \brief Return a copy of the current atom positions
   std::vector<cvm::atom_pos> positions() const;
 
+  /// \brief Calculate the center of geometry of the atomic positions, assuming
+  /// that they are already pbc-wrapped
+  int calc_center_of_geometry();
+private:
+  /// \brief Center of geometry
+  cvm::atom_pos cog;
+public:
+  /// \brief Return the center of geometry of the atomic positions
+  inline cvm::atom_pos center_of_geometry() const
+  {
+    return cog;
+  }
+
+  /// \brief Calculate the center of mass of the atomic positions, assuming that
+  /// they are already pbc-wrapped
+  int calc_center_of_mass();
+private:
+  /// \brief Center of mass
+  cvm::atom_pos com;
+  /// \brief The derivative of a scalar variable with respect to the COM
+  // TODO for scalable calculations of more complex variables (e.g. rotation),
+  // use a colvarvalue of vectors to hold the entire derivative
+  cvm::rvector scalar_com_gradient;
+public:
+  /// \brief Return the center of mass of the atomic positions
+  inline cvm::atom_pos center_of_mass() const
+  {
+    return com;
+  }
+
   /// \brief Return a copy of the current atom positions, shifted by a constant vector
   std::vector<cvm::atom_pos> positions_shifted(cvm::rvector const &shift) const;
 
-  /// \brief Return the center of geometry of the positions, assuming
-  /// that coordinates are already pbc-wrapped
-  cvm::atom_pos center_of_geometry() const;
-
-  /// \brief Return the center of mass of the positions, assuming that
-  /// coordinates are already pbc-wrapped
-  cvm::atom_pos center_of_mass() const;
-
-  /// \brief Atom positions at the previous step
-  std::vector<cvm::atom_pos> old_pos;
-
-
   /// \brief Return a copy of the current atom velocities
   std::vector<cvm::rvector> velocities() const;
 
+  ///\brief Calculate the dipole of the atom group around the specified center
+  int calc_dipole(cvm::atom_pos const &com);
+private:
+  cvm::rvector dip;
+public:
+  ///\brief Return the (previously calculated) dipole of the atom group
+  inline cvm::rvector dipole() const
+  {
+    return dip;
+  }
 
   /// \brief Return a copy of the system forces
   std::vector<cvm::rvector> system_forces() const;
+
   /// \brief Return a copy of the aggregated total force on the group
   cvm::rvector system_force() const;
 

lib/colvars/colvarcomp.cpp

@@ -51,7 +51,9 @@ void colvar::cvc::parse_group(std::string const &conf,
                               bool optional)
 {
   if (key_lookup(conf, group_key)) {
-    if (group.parse(conf, group_key) != COLVARS_OK) {
+    // TODO turn on scalable flag for group objects in cvc init function
+    group.key = group_key;
+    if (group.parse(conf) != COLVARS_OK) {
       cvm::error("Error parsing definition for atom group \""+
                          std::string(group_key)+"\".\n");
       return;
@@ -130,10 +132,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group)
       group.read_positions();
       // change one coordinate
       group[ia].pos[id] += cvm::debug_gradients_step_size;
-      // (re)do the fit (if defined)
+      group.calc_required_properties();
-      if (group.b_center || group.b_rotate) {
-        group.calc_apply_roto_translation();
-      }
       calc_value();
       cvm::real x_1 = x.real_value;
       if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];
@@ -167,7 +166,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group)
 //         ref.read_positions();
 //         // change one coordinate
 //         ref[ia].pos[id] += cvm::debug_gradients_step_size;
-//         group.calc_apply_roto_translation();
+//         group.update();
 //         calc_value();
 //         cvm::real const x_1 = x.real_value;
 //         cvm::log("refPosGroup atom "+cvm::to_str(ia)+", component "+cvm::to_str (id)+":\n");

lib/colvars/colvarcomp.h

@@ -670,6 +670,49 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+/// \brief Colvar component: angle between the dipole of a molecule and an axis
+/// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI])
+class colvar::dipole_angle
+  : public colvar::cvc
+{
+protected:
+
+  /// Dipole atom group
+  cvm::atom_group group1;
+  /// Atom group
+  cvm::atom_group group2;
+  /// Atom group
+  cvm::atom_group group3;
+
+  /// Inter site vectors
+  cvm::rvector r21, r23;
+  /// Inter site vector norms
+  cvm::real r21l, r23l;
+  /// Derivatives wrt group centers of mass
+  cvm::rvector dxdr1, dxdr3;
+
+  /// Compute system force on first site only to avoid unwanted
+  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
+  /// (or to allow dummy atoms)
+  bool b_1site_force;
+public:
+
+  /// Initialize by parsing the configuration
+  dipole_angle (std::string const &conf);
+  /// \brief Initialize the three groups after three atoms
+  dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3);
+  dipole_angle();
+  virtual inline ~dipole_angle() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force (colvarvalue const &force);
+  virtual cvm::real dist2 (colvarvalue const &x1,
+                           colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
+                                   colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
+                                   colvarvalue const &x2) const;
+};
 
 /// \brief Colvar component: dihedral between the centers of mass of
 /// four groups (colvarvalue::type_scalar type, range [-PI:PI])
@@ -1190,6 +1233,7 @@ simple_scalar_dist_functions(distance)
 simple_scalar_dist_functions(distance_xy)
 simple_scalar_dist_functions(distance_inv)
 simple_scalar_dist_functions(angle)
+simple_scalar_dist_functions(dipole_angle)
 simple_scalar_dist_functions(coordnum)
 simple_scalar_dist_functions(selfcoordnum)
 simple_scalar_dist_functions(h_bond)

lib/colvars/colvarcomp_angles.cpp

@@ -95,6 +95,12 @@ void colvar::angle::calc_gradients()
     group3[i].grad = (group3[i].mass/group3.total_mass) *
       (dxdr3);
   }
+
+  if (b_debug_gradients) {
+    debug_gradients(group1);
+    debug_gradients(group2);
+    debug_gradients(group3);
+  }
 }
 
 void colvar::angle::calc_force_invgrads()
@@ -143,6 +149,116 @@ void colvar::angle::apply_force(colvarvalue const &force)
 
 
 
+colvar::dipole_angle::dipole_angle(std::string const &conf)
+  : cvc(conf)
+{
+  function_type = "dipole_angle";
+  parse_group(conf, "group1", group1);
+  parse_group(conf, "group2", group2);
+  parse_group(conf, "group3", group3);
+
+  atom_groups.push_back(&group1);
+  atom_groups.push_back(&group2);
+  atom_groups.push_back(&group3);
+  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
+    cvm::log("Computing system force on group 1 only");
+  }
+  x.type(colvarvalue::type_scalar);
+}
+
+
+colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
+                      cvm::atom const &a2,
+                      cvm::atom const &a3)
+  : group1(std::vector<cvm::atom>(1, a1)),
+    group2(std::vector<cvm::atom>(1, a2)),
+    group3(std::vector<cvm::atom>(1, a3))
+{
+  function_type = "dipole_angle";
+  b_1site_force = false;
+  atom_groups.push_back(&group1);
+  atom_groups.push_back(&group2);
+  atom_groups.push_back(&group3);
+
+  x.type(colvarvalue::type_scalar);
+}
+
+
+colvar::dipole_angle::dipole_angle()
+{
+  function_type = "dipole_angle";
+  x.type(colvarvalue::type_scalar);
+}
+
+
+void colvar::dipole_angle::calc_value()
+{
+  cvm::atom_pos const g1_pos = group1.center_of_mass();
+  cvm::atom_pos const g2_pos = group2.center_of_mass();
+  cvm::atom_pos const g3_pos = group3.center_of_mass();
+
+  group1.calc_dipole(g1_pos);
+
+  r21 = group1.dipole();
+  r21l = r21.norm();
+  r23  = cvm::position_distance(g2_pos, g3_pos);
+  r23l = r23.norm();
+
+  cvm::real     const cos_theta = (r21*r23)/(r21l*r23l);
+
+  x.real_value = (180.0/PI) * std::acos(cos_theta);
+}
+
+//to be implemented
+//void colvar::dipole_angle::calc_force_invgrads(){}
+//void colvar::dipole_angle::calc_Jacobian_derivative(){}
+
+void colvar::dipole_angle::calc_gradients()
+{
+  cvm::real const cos_theta = (r21*r23)/(r21l*r23l);
+  cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta);
+
+  dxdr1 = (180.0/PI) * dxdcos *
+  (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l );
+
+  dxdr3 =  (180.0/PI) * dxdcos *
+    (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l );
+
+  //this auxiliar variables are to avoid numerical errors inside "for"
+  double aux1 = group1.total_charge/group1.total_mass;
+  // double aux2 = group2.total_charge/group2.total_mass;
+  // double aux3 = group3.total_charge/group3.total_mass;
+
+  size_t i;
+  for (i = 0; i < group1.size(); i++) {
+    group1[i].grad =(group1[i].charge + (-1)* group1[i].mass * aux1) * (dxdr1);
+  }
+
+  for (i = 0; i < group2.size(); i++) {
+    group2[i].grad = (group2[i].mass/group2.total_mass)* dxdr3 * (-1.0);
+  }
+
+  for (i = 0; i < group3.size(); i++) {
+    group3[i].grad =(group3[i].mass/group3.total_mass) * (dxdr3);
+  }
+}
+
+
+void colvar::dipole_angle::apply_force(colvarvalue const &force)
+{
+  if (!group1.noforce)
+    group1.apply_colvar_force(force.real_value);
+
+  if (!group2.noforce)
+    group2.apply_colvar_force(force.real_value);
+
+  if (!group3.noforce)
+    group3.apply_colvar_force(force.real_value);
+}
+
+
+
+
 colvar::dihedral::dihedral(std::string const &conf)
   : cvc(conf)
 {

lib/colvars/colvarcomp_rotations.cpp

@@ -88,9 +88,6 @@ colvar::orientation::orientation()
 
 void colvar::orientation::calc_value()
 {
-  // atoms.reset_atoms_data();
-  // atoms.read_positions();
-
   atoms_cog = atoms.center_of_geometry();
 
   rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
@@ -145,9 +142,6 @@ colvar::orientation_angle::orientation_angle()
 
 void colvar::orientation_angle::calc_value()
 {
-  // atoms.reset_atoms_data();
-  // atoms.read_positions();
-
   atoms_cog = atoms.center_of_geometry();
 
   rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
@@ -261,9 +255,6 @@ colvar::tilt::tilt()
 
 void colvar::tilt::calc_value()
 {
-  // atoms.reset_atoms_data();
-  // atoms.read_positions();
-
   atoms_cog = atoms.center_of_geometry();
 
   rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));
@@ -331,9 +322,6 @@ colvar::spin_angle::spin_angle()
 
 void colvar::spin_angle::calc_value()
 {
-  // atoms.reset_atoms_data();
-  // atoms.read_positions();
-
   atoms_cog = atoms.center_of_geometry();
 
   rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog));

lib/colvars/colvargrid.cpp

@@ -127,7 +127,6 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os)
 
   integral = 0.0;
   int_vals.push_back( 0.0 );
-  bin = 0.0;
   min = 0.0;
 
   // correction for periodic colvars, so that the PMF is periodic
@@ -138,7 +137,7 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os)
     corr = 0.0;
   }
 
-  for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix), bin += 1.0 ) {
+  for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
 
     if (samples) {
       size_t const samples_here = samples->value(ix);

lib/colvars/colvarmodule.cpp

@@ -100,38 +100,24 @@ int colvarmodule::read_config_string(std::string const &config_str)
 
 int colvarmodule::parse_config(std::string &conf)
 {
-  int error_code = 0;
-
   // parse global options
-  error_code |= parse_global_params(conf);
+  if (catch_input_errors(parse_global_params(conf))) {
-
+    return get_error();
-  if (error_code != COLVARS_OK || cvm::get_error()) {
-    set_error_bits(INPUT_ERROR);
-    return COLVARS_ERROR;
   }
 
   // parse the options for collective variables
-  error_code |= parse_colvars(conf);
+  if (catch_input_errors(parse_colvars(conf))) {
-
+    return get_error();
-  if (error_code != COLVARS_OK || cvm::get_error()) {
-    set_error_bits(INPUT_ERROR);
-    return COLVARS_ERROR;
   }
 
   // parse the options for biases
-  error_code |= parse_biases(conf);
+  if (catch_input_errors(parse_biases(conf))) {
-
+    return get_error();
-  if (error_code != COLVARS_OK || cvm::get_error()) {
-    set_error_bits(INPUT_ERROR);
-    return COLVARS_ERROR;
   }
 
   // done parsing known keywords, check that all keywords found were valid ones
-  error_code |= parse->check_keywords(conf, "colvarmodule");
+  if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) {
-
+    return get_error();
-  if (error_code != COLVARS_OK || cvm::get_error()) {
-    set_error_bits(INPUT_ERROR);
-    return COLVARS_ERROR;
   }
 
   cvm::log(cvm::line_marker);
@@ -143,7 +129,7 @@ int colvarmodule::parse_config(std::string &conf)
     write_traj_label(cv_traj_os);
   }
 
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+  return get_error();
 }
 
 
@@ -181,7 +167,7 @@ int colvarmodule::parse_global_params(std::string const &conf)
                     colvarparse::parse_silent);
 
   if (use_scripted_forces && !proxy->force_script_defined) {
-    cvm::fatal_error("User script for scripted colvar forces not found.");
+    cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
   }
 
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@@ -315,6 +301,17 @@ int colvarmodule::parse_biases(std::string const &conf)
 }
 
 
+int colvarmodule::catch_input_errors(int result)
+{
+  if (result != COLVARS_OK || get_error()) {
+    set_error_bits(result | 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();
@@ -652,7 +649,6 @@ int colvarmodule::setup()
     (*cvi)->setup();
   }
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
-
 }
 
 colvarmodule::~colvarmodule()
@@ -664,6 +660,8 @@ colvarmodule::~colvarmodule()
 
 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<colvar *>::reverse_iterator cvi = colvars.rbegin();
@@ -1072,6 +1070,8 @@ void cvm::error(std::string const &message, int code)
 
 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_bits(FATAL_ERROR);
   proxy->fatal_error(message);
 }
@@ -1143,7 +1143,7 @@ int cvm::read_index_file(char const *filename)
 }
 
 int cvm::load_atoms(char const *file_name,
-                    std::vector<cvm::atom> &atoms,
+                    cvm::atom_group &atoms,
                     std::string const &pdb_field,
                     double const pdb_field_value)
 {

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2015-10-22"
+#define COLVARS_VERSION "2016-02-28"
 #endif
 
 #ifndef COLVARS_DEBUG
@@ -31,6 +31,7 @@
 #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>
@@ -109,6 +110,7 @@ public:
   static inline void set_error_bits(int code)
   {
     errorCode |= code;
+    errorCode |= COLVARS_ERROR;
   }
   static inline int get_error()
   {
@@ -227,6 +229,12 @@ public:
   /// on error, delete new bias
   bool check_new_bias(std::string &conf, char const *key);
 
+private:
+  /// Useful wrapper to interrupt parsing if any error occurs
+  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
@@ -420,7 +428,7 @@ public:
   /// \param pdb_field (optiona) if "filename" is a PDB file, use this
   /// field to determine which are the atoms to be set
   static int load_atoms(char const *filename,
-                          std::vector<atom> &atoms,
+                        atom_group &atoms,
                         std::string const &pdb_field,
                         double const pdb_field_value = 0.0);
 

lib/colvars/colvarparse.cpp

@@ -110,8 +110,8 @@ size_t      colvarparse::dummy_pos = 0;
         data_count++;                                                   \
       }                                                                 \
       if (data_count == 0)                                              \
-        cvm::fatal_error("Error: in parsing \""+                        \
+        cvm::error("Error: in parsing \""+                              \
-                         std::string(key)+"\".\n");                     \
+                   std::string(key)+"\".\n", INPUT_ERROR);              \
       if (data_count > 1) {                                             \
         cvm::error("Error: multiple values "                            \
                    "are not allowed for keyword \""+                    \
@@ -283,8 +283,8 @@ bool colvarparse::get_keyval(std::string const &conf,
                 (data == std::string("false")) ) {
       value = false;
     } else
-      cvm::fatal_error("Error: boolean values only are allowed "
+      cvm::error("Error: boolean values only are allowed "
-                       "for \""+std::string(key)+"\".\n");
+                 "for \""+std::string(key)+"\".\n", INPUT_ERROR);
     if (parse_mode != parse_silent) {
       cvm::log("# "+std::string(key)+" = "+
                (value ? "on" : "off")+"\n");
@@ -398,10 +398,9 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
       }
     }
     if (!found_keyword) {
-      cvm::log("Error: keyword \""+uk+"\" is not supported, "
+      cvm::error("Error: keyword \""+uk+"\" is not supported, "
-               "or not recognized in this context.\n");
+                 "or not recognized in this context.\n", INPUT_ERROR);
-      cvm::set_error_bits(INPUT_ERROR);
+      return INPUT_ERROR;
-      return COLVARS_ERROR;
     }
   }
 
@@ -527,10 +526,10 @@ bool colvarparse::key_lookup(std::string const &conf,
         while (brace == std::string::npos) {
           // add a new line
           if (line_end >= conf.size()) {
-            cvm::fatal_error("Parse error: reached the end while "
+            cvm::error("Parse error: reached the end while "
                        "looking for closing brace; until now "
                        "the following was parsed: \"\n"+
-                             line+"\".\n");
+                       line+"\".\n", INPUT_ERROR);
             return false;
           }
           size_t const old_end = line.size();

lib/colvars/colvarparse.h

@@ -147,6 +147,7 @@ public:
   /// \brief Use this after parsing a config string (note that check_keywords() calls it already)
   void clear_keyword_registry();
 
+  inline void reset() { clear_keyword_registry(); }
 
   /// \brief Return a lowercased copy of the string
   static inline std::string to_lower_cppstr(std::string const &in)

lib/colvars/colvarproxy.h

@@ -9,16 +9,13 @@
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 
-// return values for the frame() routine
-#define COLVARS_NO_SUCH_FRAME -1
-
 // forward declarations
 class colvarscript;
 
 /// \brief Interface between the collective variables module and
-/// the simulation or analysis program.
+/// the simulation or analysis program (NAMD, VMD, LAMMPS...).
-/// This is the base class: each program is supported by a derived class.
+/// This is the base class: each interfaced program is supported by a derived class.
-/// Only pure virtual functions ("= 0") must be reimplemented in a new interface.
+/// Only pure virtual functions ("= 0") must be reimplemented to ensure baseline functionality.
 
 class colvarproxy {
 
@@ -31,13 +28,29 @@ public:
   inline colvarproxy() : script(NULL) {}
 
   /// Default destructor
-  virtual inline ~colvarproxy() {}
+  virtual ~colvarproxy() {}
 
-  /// (Re)initialize member data after construction
+  /// (Re)initialize required member data after construction
-  virtual void setup() {}
+  virtual int setup()
+  {
+    return COLVARS_OK;
+  }
 
+  /// \brief Update data required by the colvars module (e.g. cache atom positions)
+  ///
+  /// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these
+  virtual int update_input()
+  {
+    return COLVARS_OK;
+  }
 
-  // **************** SYSTEM-WIDE PHYSICAL QUANTITIES ****************
+  /// \brief Update data based from the results of a module update (e.g. send forces)
+  virtual int update_output()
+  {
+    return COLVARS_OK;
+  }
+
+  // **************** SIMULATION PARAMETERS ****************
 
   /// \brief Value of the unit for atomic coordinates with respect to
   /// angstroms (used by some variables for hard-coded default values)
@@ -62,6 +75,45 @@ public:
   // return 0 on success, -1 on failure
   virtual int frame(int) { return COLVARS_NOT_IMPLEMENTED; }
 
+  /// \brief Prefix to be used for input files (restarts, not
+  /// configuration)
+  std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
+
+  inline std::string input_prefix()
+  {
+    return input_prefix_str;
+  }
+
+  /// \brief Prefix to be used for output restart files
+  inline std::string restart_output_prefix()
+  {
+    return restart_output_prefix_str;
+  }
+
+  /// \brief Prefix to be used for output files (final system
+  /// configuration)
+  inline std::string output_prefix()
+  {
+    return output_prefix_str;
+  }
+
+  /// \brief Restarts will be written each time this number of steps has passed
+  virtual size_t restart_frequency()
+  {
+    return 0;
+  }
+
+protected:
+
+  /// \brief Currently opened output files: by default, these are ofstream objects.
+  /// Allows redefinition to implement different output mechanisms
+  std::list<std::ostream *> output_files;
+  /// \brief Identifiers for output_stream objects: by default, these are the names of the files
+  std::list<std::string>    output_stream_names;
+
+public:
+
+  // **************** MULTIPLE REPLICAS COMMUNICATION ****************
 
   // Replica exchange commands:
 
@@ -87,68 +139,6 @@ public:
     return COLVARS_NOT_IMPLEMENTED;
   }
 
-  // **************** SIMULATION PARAMETERS ****************
-
-  /// \brief Prefix to be used for input files (restarts, not
-  /// configuration)
-  std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
-
-  inline std::string input_prefix()
-  {
-    return input_prefix_str;
-  }
-
-  /// \brief Prefix to be used for output restart files
-  inline std::string restart_output_prefix()
-  {
-    return restart_output_prefix_str;
-  }
-
-  /// \brief Prefix to be used for output files (final system
-  /// configuration)
-  inline std::string output_prefix()
-  {
-    return output_prefix_str;
-  }
-
-  /// \brief Restarts will be written each time this number of steps has passed
-  virtual size_t restart_frequency() = 0;
-
-
-
-  // **************** ACCESS ATOMIC DATA ****************
-
-  /// Pass restraint energy value for current timestep to MD engine
-  virtual void add_energy(cvm::real energy) = 0;
-
-  /// Tell the proxy whether system forces are needed (may not always be available)
-  virtual void request_system_force(bool yesno) = 0;
-
-
-
-  // **************** PERIODIC BOUNDARY CONDITIONS ****************
-
-  /// \brief Get the PBC-aware distance vector between two positions
-  virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
-                                         cvm::atom_pos const &pos2) = 0;
-
-  /// \brief Get the PBC-aware square distance between two positions;
-  /// may be implemented independently from position_distance() for optimization purposes
-  virtual cvm::real position_dist2(cvm::atom_pos const &pos1,
-                                   cvm::atom_pos const &pos2);
-
-  /// \brief Get the closest periodic image to a reference position
-  /// \param pos The position to look for the closest periodic image
-  /// \param ref_pos The reference position
-  virtual void select_closest_image(cvm::atom_pos &pos,
-                                    cvm::atom_pos const &ref_pos) = 0;
-
-  /// \brief Perform select_closest_image() on a set of atomic positions
-  ///
-  /// After that, distance vectors can then be calculated directly,
-  /// without using position_distance()
-  void select_closest_images(std::vector<cvm::atom_pos> &pos,
-                             cvm::atom_pos const &ref_pos);
 
   // **************** SCRIPTING INTERFACE ****************
 
@@ -175,6 +165,7 @@ public:
                                            std::vector<cvm::matrix2d<cvm::real> > &gradient)
   { return COLVARS_NOT_IMPLEMENTED; }
 
+
   // **************** INPUT/OUTPUT ****************
 
   /// Print a message to the main log
@@ -187,36 +178,11 @@ public:
   virtual void fatal_error(std::string const &message) = 0;
 
   /// Print a message to the main log and exit normally
-  virtual void exit(std::string const &message) = 0;
+  virtual void exit(std::string const &message)
-
+  {
-  /// \brief Read atom identifiers from a file \param filename name of
+    cvm::error("Error: exiting without error is not implemented, returning error code.\n",
-  /// the file (usually a PDB) \param atoms array to which atoms read
+               COLVARS_NOT_IMPLEMENTED);
-  /// from "filename" will be appended \param pdb_field (optiona) if
+  }
-  /// "filename" is a PDB file, use this field to determine which are
-  /// the atoms to be set
-  virtual int load_atoms(char const *filename,
-                         std::vector<cvm::atom> &atoms,
-                         std::string const &pdb_field,
-                         double const pdb_field_value = 0.0) = 0;
-
-  /// \brief Load the coordinates for a group of atoms from a file
-  /// (usually a PDB); if "pos" is already allocated, the number of its
-  /// elements must match the number of atoms in "filename"
-  virtual int load_coords(char const *filename,
-                          std::vector<cvm::atom_pos> &pos,
-                          const std::vector<int> &indices,
-                          std::string const &pdb_field,
-                          double const pdb_field_value = 0.0) = 0;
-
-protected:
-
-  /// \brief Open output files: by default, these are ofstream objects.
-  /// Allows redefinition to implement different output mechanisms
-  std::list<std::ostream *> output_files;
-  /// \brief Identifiers for output_stream objects: by default, these are the names of the files
-  std::list<std::string>    output_stream_names;
-
-public:
 
   // TODO the following definitions may be moved to a .cpp file
 
@@ -254,6 +220,8 @@ public:
         return COLVARS_OK;
       }
     }
+    cvm::error("Error: trying to close an output file or stream that wasn't open.\n",
+               BUG_ERROR);
     return COLVARS_ERROR;
   }
 
@@ -262,10 +230,48 @@ public:
   {
     return COLVARS_NOT_IMPLEMENTED;
   }
-};
 
 
-inline void colvarproxy::select_closest_images(std::vector<cvm::atom_pos> &pos,
+
+  // **************** ACCESS SYSTEM DATA ****************
+
+  /// Pass restraint energy value for current timestep to MD engine
+  virtual void add_energy(cvm::real energy) = 0;
+
+  /// Tell the proxy whether system forces are needed (may not always be available)
+  virtual void request_system_force(bool yesno)
+  {
+    if (yesno == true)
+      cvm::error("Error: system forces are currently not implemented.\n",
+                 COLVARS_NOT_IMPLEMENTED);
+  }
+
+  /// \brief Get the PBC-aware distance vector between two positions
+  virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
+                                         cvm::atom_pos const &pos2) = 0;
+
+  /// \brief Get the PBC-aware square distance between two positions;
+  /// may need to be reimplemented independently from position_distance() for optimization purposes
+  virtual cvm::real position_dist2(cvm::atom_pos const &pos1,
+                                   cvm::atom_pos const &pos2)
+  {
+    return (position_distance(pos1, pos2)).norm2();
+  }
+
+  /// \brief Get the closest periodic image to a reference position
+  /// \param pos The position to look for the closest periodic image
+  /// \param ref_pos The reference position
+  virtual void select_closest_image(cvm::atom_pos &pos,
+                                    cvm::atom_pos const &ref_pos)
+  {
+    pos = position_distance(ref_pos, pos) + ref_pos;
+  }
+
+  /// \brief Perform select_closest_image() on a set of atomic positions
+  ///
+  /// After that, distance vectors can then be calculated directly,
+  /// without using position_distance()
+  void select_closest_images(std::vector<cvm::atom_pos> &pos,
                              cvm::atom_pos const &ref_pos)
   {
     for (std::vector<cvm::atom_pos>::iterator pi = pos.begin();
@@ -274,10 +280,281 @@ inline void colvarproxy::select_closest_images(std::vector<cvm::atom_pos> &pos,
     }
   }
 
-inline cvm::real colvarproxy::position_dist2(cvm::atom_pos const &pos1,
+
-                                             cvm::atom_pos const &pos2)
+  // **************** ACCESS ATOMIC DATA ****************
+protected:
+
+  /// \brief Array of 0-based integers used to uniquely associate atoms
+  /// within the host program
+  std::vector<int>          atoms_ids;
+  /// \brief Keep track of how many times each atom is used by a separate colvar object
+  std::vector<size_t>       atoms_ncopies;
+  /// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
+  std::vector<cvm::real>    atoms_masses;
+  /// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS)
+  std::vector<cvm::real>    atoms_charges;
+  /// \brief Current three-dimensional positions of the atoms
+  std::vector<cvm::rvector> atoms_positions;
+  /// \brief Most recent total forces on each atom
+  std::vector<cvm::rvector> atoms_total_forces;
+  /// \brief Most recent forces applied by external potentials onto each atom
+  std::vector<cvm::rvector> atoms_applied_forces;
+  /// \brief Forces applied from colvars, to be communicated to the MD integrator
+  std::vector<cvm::rvector> atoms_new_colvar_forces;
+
+  /// Used by all init_atom() functions: create a slot for an atom not requested yet
+  inline int add_atom_slot(int atom_id)
   {
-  return (position_distance(pos1, pos2)).norm2();
+    atoms_ids.push_back(atom_id);
+    atoms_ncopies.push_back(1);
+    atoms_masses.push_back(1.0);
+    atoms_charges.push_back(0.0);
+    atoms_positions.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    atoms_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    atoms_applied_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    atoms_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    return (atoms_ids.size() - 1);
   }
 
+public:
+
+  /// Prepare this atom for collective variables calculation, selecting it by numeric index (1-based)
+  virtual int init_atom(int atom_number) = 0;
+
+  /// Check that this atom number is valid, but do not initialize the corresponding atom yet
+  virtual int check_atom_id(int atom_number) = 0;
+
+  /// Select this atom for collective variables calculation, using name and residue number.
+  /// Not all programs support this: leave this function as is in those cases.
+  virtual int init_atom(cvm::residue_id const &residue,
+                        std::string const     &atom_name,
+                        std::string const     &segment_id)
+  {
+    cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n",
+               COLVARS_NOT_IMPLEMENTED);
+    return COLVARS_NOT_IMPLEMENTED;
+  }
+
+  /// Check that this atom is valid, but do not initialize it yet
+  virtual int check_atom_id(cvm::residue_id const &residue,
+                            std::string const     &atom_name,
+                            std::string const     &segment_id)
+  {
+    cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n",
+               COLVARS_NOT_IMPLEMENTED);
+    return COLVARS_NOT_IMPLEMENTED;
+  }
+
+  /// \brief Used by the atom class destructor: rather than deleting the array slot
+  /// (costly) set the corresponding atoms_ncopies to zero
+  virtual void clear_atom(int index)
+  {
+    if (((size_t) index) >= atoms_ids.size()) {
+      cvm::error("Error: trying to disable an atom that was not previously requested.\n",
+                 INPUT_ERROR);
+    }
+    if (atoms_ncopies[index] > 0) {
+      atoms_ncopies[index] -= 1;
+    }
+  }
+
+  /// Get the numeric ID of the given atom (for the program)
+  inline int get_atom_id(int index) const
+  {
+    return atoms_ids[index];
+  }
+
+  /// Get the mass of the given atom
+  inline cvm::real get_atom_mass(int index) const
+  {
+    return atoms_masses[index];
+  }
+
+  /// Get the charge of the given atom
+  inline cvm::real get_atom_charge(int index) const
+  {
+    return atoms_charges[index];
+  }
+
+  /// Read the current position of the given atom
+  inline cvm::rvector get_atom_position(int index) const
+  {
+    return atoms_positions[index];
+  }
+
+  /// Read the current total system force of the given atom
+  inline cvm::rvector get_atom_system_force(int index) const
+  {
+    return atoms_total_forces[index] - atoms_applied_forces[index];
+  }
+
+  /// Request that this force is applied to the given atom
+  inline void apply_atom_force(int index, cvm::rvector const &new_force)
+  {
+    atoms_new_colvar_forces[index] += new_force;
+  }
+
+  /// Read the current velocity of the given atom
+  virtual cvm::rvector get_atom_velocity(int index)
+  {
+    cvm::error("Error: reading the current velocity of an atom is not yet implemented.\n",
+               COLVARS_NOT_IMPLEMENTED);
+    return cvm::rvector(0.0);
+  }
+
+  // useful functions for data management outside this class
+  inline std::vector<int> *modify_atom_ids() { return &atoms_ids; }
+  inline std::vector<cvm::real> *modify_atom_masses() { return &atoms_masses; }
+  inline std::vector<cvm::real> *modify_atom_charges() { return &atoms_charges; }
+  inline std::vector<cvm::rvector> *modify_atom_positions() { return &atoms_positions; }
+  inline std::vector<cvm::rvector> *modify_atom_total_forces() { return &atoms_total_forces; }
+  inline std::vector<cvm::rvector> *modify_atom_applied_forces() { return &atoms_applied_forces; }
+  inline std::vector<cvm::rvector> *modify_atom_new_colvar_forces() { return &atoms_new_colvar_forces; }
+
+  /// \brief Read atom identifiers from a file \param filename name of
+  /// the file (usually a PDB) \param atoms array to which atoms read
+  /// from "filename" will be appended \param pdb_field (optiona) if
+  /// "filename" is a PDB file, use this field to determine which are
+  /// the atoms to be set
+  virtual int load_atoms(char const *filename,
+                         cvm::atom_group &atoms,
+                         std::string const &pdb_field,
+                         double const pdb_field_value = 0.0)
+  {
+    cvm::error("Error: loading atom identifiers from a file is currently not implemented.\n",
+               COLVARS_NOT_IMPLEMENTED);
+    return COLVARS_NOT_IMPLEMENTED;
+  }
+
+  /// \brief Load the coordinates for a group of atoms from a file
+  /// (usually a PDB); if "pos" is already allocated, the number of its
+  /// elements must match the number of atoms in "filename"
+  virtual int load_coords(char const *filename,
+                          std::vector<cvm::atom_pos> &pos,
+                          const std::vector<int> &indices,
+                          std::string const &pdb_field,
+                          double const pdb_field_value = 0.0)
+  {
+    cvm::error("Error: loading atomic coordinates from a file is currently not implemented.\n");
+    return COLVARS_NOT_IMPLEMENTED;
+  }
+
+  // **************** ACCESS GROUP DATA ****************
+
+protected:
+
+  /// \brief Array of 0-based integers used to uniquely associate atom groups
+  /// within the host program
+  std::vector<int>          atom_groups_ids;
+  /// \brief Keep track of how many times each group is used by a separate cvc
+  std::vector<size_t>       atom_groups_ncopies;
+  /// \brief Total masses of the atom groups
+  std::vector<cvm::real>    atom_groups_masses;
+  /// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS)
+  std::vector<cvm::real>    atom_groups_charges;
+  /// \brief Current centers of mass of the atom groups
+  std::vector<cvm::rvector> atom_groups_coms;
+  /// \brief Most recently updated total forces on the com of each group
+  std::vector<cvm::rvector> atom_groups_total_forces;
+  /// \brief Most recent forces applied by external potentials onto each group
+  std::vector<cvm::rvector> atom_groups_applied_forces;
+  /// \brief Forces applied from colvars, to be communicated to the MD integrator
+  std::vector<cvm::rvector> atom_groups_new_colvar_forces;
+
+  /// TODO Add here containers of handles to cvc objects that are computed in parallel
+
+public:
+
+  /// \brief Whether this proxy implementation has capability for scalable groups
+  virtual bool has_scalable_groups() const
+  {
+    return false;
+  }
+
+  /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet
+  // TODO Add a handle to cvc objects
+  inline int add_atom_group_slot(int atom_group_id)
+  {
+    atom_groups_ids.push_back(atom_group_id);
+    atom_groups_masses.push_back(1.0);
+    atom_groups_charges.push_back(0.0);
+    atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    atom_groups_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    atom_groups_applied_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    atom_groups_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0));
+    return (atom_groups_ids.size() - 1);
+  }
+
+  /// Prepare this group for collective variables calculation, selecting atoms by internal ids (0-based)
+  virtual int init_atom_group(std::vector<int> const &atoms_ids) // TODO Add a handle to cvc objects
+  {
+    cvm::error("Error: initializing a group outside of the colvars module is currently not supported.\n",
+               COLVARS_NOT_IMPLEMENTED);
+    return COLVARS_NOT_IMPLEMENTED;
+  }
+
+  /// \brief Used by the atom_group class destructor
+  virtual void clear_atom_group(int index)
+  {
+    if (((size_t) index) >= atom_groups_ids.size()) {
+      cvm::error("Error: trying to disable an atom group that was not previously requested.\n",
+                 INPUT_ERROR);
+    }
+
+    atom_groups_ids.erase(atom_groups_ids.begin()+index);
+    atom_groups_masses.erase(atom_groups_masses.begin()+index);
+    atom_groups_charges.erase(atom_groups_charges.begin()+index);
+    atom_groups_coms.erase(atom_groups_coms.begin()+index);
+    atom_groups_total_forces.erase(atom_groups_total_forces.begin()+index);
+    atom_groups_applied_forces.erase(atom_groups_applied_forces.begin()+index);
+    atom_groups_new_colvar_forces.erase(atom_groups_new_colvar_forces.begin()+index);
+  }
+
+  /// Get the numeric ID of the given atom group (for the MD program)
+  inline cvm::real get_atom_group_id(int index) const
+  {
+    return atom_groups_ids[index];
+  }
+
+  /// Get the mass of the given atom group
+  inline cvm::real get_atom_group_mass(int index) const
+  {
+    return atom_groups_masses[index];
+  }
+
+  /// Get the charge of the given atom group
+  inline cvm::real get_atom_group_charge(int index) const
+  {
+    return atom_groups_charges[index];
+  }
+
+  /// Read the current position of the center of mass given atom group
+  inline cvm::rvector get_atom_group_com(int index) const
+  {
+    return atom_groups_coms[index];
+  }
+
+  /// Read the current total system force of the given atom group
+  inline cvm::rvector get_atom_group_system_force(int index) const
+  {
+    return atom_groups_total_forces[index] - atom_groups_applied_forces[index];
+  }
+
+  /// Request that this force is applied to the given atom group
+  inline void apply_atom_group_force(int index, cvm::rvector const &new_force)
+  {
+    atom_groups_new_colvar_forces[index] += new_force;
+  }
+
+  /// Read the current velocity of the given atom group
+  virtual cvm::rvector get_atom_group_velocity(int index)
+  {
+    cvm::error("Error: reading the current velocity of an atom group is not yet implemented.\n",
+               COLVARS_NOT_IMPLEMENTED);
+    return cvm::rvector(0.0);
+  }
+
+};
+
+
 #endif

lib/colvars/colvarscript.cpp

@@ -26,11 +26,13 @@ int colvarscript::run(int argc, char const *argv[]) {
 
   if (argc < 2) {
     result = help_string();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   std::string cmd = argv[1];
 
+  int error_code = COLVARS_OK;
+
   if (cmd == "colvar") {
     return proc_colvar(argc-1, &(argv[1]));
   }
@@ -41,13 +43,13 @@ int colvarscript::run(int argc, char const *argv[]) {
 
   if (cmd == "version") {
     result = COLVARS_VERSION;
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (cmd == "reset") {
     /// Delete every child object
     colvars->reset();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (cmd == "delete") {
@@ -55,12 +57,17 @@ int colvarscript::run(int argc, char const *argv[]) {
     // Note: the delete bit may be ignored by some backends
     // it is mostly useful in VMD
     colvars->set_error_bits(DELETE_COLVARS);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (cmd == "update") {
-    colvars->calc();
+    error_code |= proxy->update_input();
-    return COLVARSCRIPT_OK;
+    error_code |= colvars->calc();
+    error_code |= proxy->update_output();
+    if (error_code) {
+      result += "Error updating the colvars module.\n";
+    }
+    return error_code;
   }
 
   if (cmd == "list") {
@@ -70,14 +77,14 @@ int colvarscript::run(int argc, char const *argv[]) {
            ++cvi) {
         result += (cvi == colvars->colvars.begin() ? "" : " ") + (*cvi)->name;
       }
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     } else if (argc == 3 && !strcmp(argv[2], "biases")) {
       for (std::vector<colvarbias *>::iterator bi = colvars->biases.begin();
            bi != colvars->biases.end();
            ++bi) {
         result += (bi == colvars->biases.begin() ? "" : " ") + (*bi)->name;
       }
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     } else {
       result = "Wrong arguments to command \"list\"\n" + help_string();
       return COLVARSCRIPT_ERROR;
@@ -91,7 +98,7 @@ int colvarscript::run(int argc, char const *argv[]) {
       return COLVARSCRIPT_ERROR;
     }
     if (colvars->read_config_file(argv[2]) == COLVARS_OK) {
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     } else {
       result = "Error parsing configuration file";
       return COLVARSCRIPT_ERROR;
@@ -106,7 +113,7 @@ int colvarscript::run(int argc, char const *argv[]) {
     }
     std::string conf = argv[2];
     if (colvars->read_config_string(conf) == COLVARS_OK) {
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     } else {
       result = "Error parsing configuration string";
       return COLVARSCRIPT_ERROR;
@@ -121,7 +128,7 @@ int colvarscript::run(int argc, char const *argv[]) {
     }
     proxy->input_prefix_str = argv[2];
     if (colvars->setup_input() == COLVARS_OK) {
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     } else {
       result = "Error loading state file";
       return COLVARSCRIPT_ERROR;
@@ -138,7 +145,7 @@ int colvarscript::run(int argc, char const *argv[]) {
     int error = 0;
     error |= colvars->setup_output();
     error |= colvars->write_output_files();
-    return error ? COLVARSCRIPT_ERROR : COLVARSCRIPT_OK;
+    return error ? COLVARSCRIPT_ERROR : COLVARS_OK;
   }
 
   /// Print the values that would go on colvars.traj
@@ -146,13 +153,13 @@ int colvarscript::run(int argc, char const *argv[]) {
     std::ostringstream os;
     colvars->write_traj_label(os);
     result = os.str();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
   if (cmd == "printframe") {
     std::ostringstream os;
     colvars->write_traj(os);
     result = os.str();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (cmd == "frame") {
@@ -160,7 +167,7 @@ int colvarscript::run(int argc, char const *argv[]) {
       int f = proxy->frame();
       if (f >= 0) {
         result = cvm::to_str(f);
-        return COLVARSCRIPT_OK;
+        return COLVARS_OK;
       } else {
         result = "Frame number is not available";
         return COLVARSCRIPT_ERROR;
@@ -171,7 +178,7 @@ int colvarscript::run(int argc, char const *argv[]) {
       long int f = proxy->frame(strtol(argv[2], NULL, 10));
       colvars->it = proxy->frame();
       result = cvm::to_str(f);
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     } else {
       result = "Wrong arguments to command \"frame\"\n" + help_string();
       return COLVARSCRIPT_ERROR;
@@ -199,24 +206,24 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
 
   if (subcmd == "value") {
     result = (cv->value()).to_simple_string();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "width") {
     result = cvm::to_str(cv->width, 0, cvm::cv_prec);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "type") {
     result = cv->value().type_desc(cv->value().value_type);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "update") {
     cv->calc();
     cv->update();
     result = (cv->value()).to_simple_string();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "delete") {
@@ -228,12 +235,12 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
     delete cv;
     // TODO this could be done by the destructors
     colvars->write_traj_label(colvars->cv_traj_os);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "getconfig") {
     result = cv->get_config();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "addforce") {
@@ -253,7 +260,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
     }
     cv->add_bias_force(force);
     result = force.to_simple_string();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "cvcflags") {
@@ -276,7 +283,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
       return COLVARSCRIPT_ERROR;
     }
     result = "0";
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   result = "Syntax error\n" + help_string();
@@ -301,25 +308,25 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
 
   if (subcmd == "energy") {
     result = cvm::to_str(b->get_energy());
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "update") {
     b->update();
     result = cvm::to_str(b->get_energy());
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "getconfig") {
     result = b->get_config();
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   // Subcommands for MW ABF
   if (subcmd == "bin") {
     int r = b->current_bin();
     result = cvm::to_str(r);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "binnum") {
@@ -329,7 +336,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
       return COLVARSCRIPT_ERROR;
     }
     result = cvm::to_str(r);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (subcmd == "share") {
@@ -339,7 +346,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
       return COLVARSCRIPT_ERROR;
     }
     result = cvm::to_str(r);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
   // End commands for MW ABF
 
@@ -348,7 +355,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
     delete b;
     // TODO this could be done by the destructors
     colvars->write_traj_label(colvars->cv_traj_os);
-    return COLVARSCRIPT_OK;
+    return COLVARS_OK;
   }
 
   if (argc >= 4) {
@@ -360,7 +367,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
         return COLVARSCRIPT_ERROR;
       }
       result = cvm::to_str(b->bin_count(index));
-      return COLVARSCRIPT_OK;
+      return COLVARS_OK;
     }
 
     result = "Syntax error\n" + help_string();

lib/colvars/colvarscript.h

@@ -9,6 +9,7 @@
 #include "colvarbias.h"
 #include "colvarproxy.h"
 
+// TODO merge these into colvarmodule.h
 #define COLVARSCRIPT_ERROR -1
 #define COLVARSCRIPT_OK 0