diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 40f3c236ca..1c3d9c8215 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -1291,7 +1291,7 @@ std::ostream & colvar::write_traj_label (std::ostream & os)
   if (tasks[task_output_system_force]) {
 
     os << " fs_"
-       << cvm::wrap_string (this->name, this_cv_width-2);
+       << cvm::wrap_string (this->name, this_cv_width-3);
 
     if (tasks[task_extended_lagrangian]) {
       // restraint center
@@ -1302,7 +1302,7 @@ std::ostream & colvar::write_traj_label (std::ostream & os)
 
   if (tasks[task_output_applied_force]) {
     os << " fa_"
-       << cvm::wrap_string (this->name, this_cv_width-2);
+       << cvm::wrap_string (this->name, this_cv_width-3);
   }
 
   return os;
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index f3bfb9477f..10898f0fa9 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -105,9 +105,19 @@ void cvm::atom_group::parse (std::string const &conf,
   colvarparse::Parse_Mode mode = parse_normal;
 
   {
-    // get the atoms by numbers
+    //    std::vector<int> atom_indexes;
+    std::string numbers_conf = "";
+    size_t pos = 0;
     std::vector<int> atom_indexes;
-    while (get_keyval (group_conf, "atomNumbers", atom_indexes, atom_indexes, mode)) {
+    while (key_lookup (group_conf, "atomNumbers", numbers_conf, pos)) {
+      if (numbers_conf.size()) {
+        std::istringstream is (numbers_conf);
+        int ia;
+        while (is >> ia) {
+          atom_indexes.push_back (ia);
+        }
+      }
+
       if (atom_indexes.size()) {
         this->reserve (this->size()+atom_indexes.size());
         for (size_t i = 0; i < atom_indexes.size(); i++) {
@@ -116,6 +126,7 @@ void cvm::atom_group::parse (std::string const &conf,
       } else
         cvm::fatal_error ("Error: no numbers provided for \""
                           "atomNumbers\".\n");
+
       atom_indexes.clear();
     }
   }
@@ -381,7 +392,7 @@ void cvm::atom_group::parse (std::string const &conf,
   this->check_keywords (group_conf, key);
 
   cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
-            cvm::to_str (this->size())+" initialized: total mass = "+
+            cvm::to_str (this->size())+" atoms initialized: total mass = "+
             cvm::to_str (this->total_mass)+".\n");
 
   cvm::decrease_depth();
diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp
index 2ee494c58c..561bb4a688 100644
--- a/lib/colvars/colvarbias.cpp
+++ b/lib/colvars/colvarbias.cpp
@@ -159,7 +159,7 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
 
     if (target_nstages) {
       // This means that either numStages of lambdaSchedule has been provided
-      stage = -1;
+      stage = 0;
       restraint_FE = 0.0;
     }
 
@@ -236,8 +236,7 @@ cvm::real colvarbias_harmonic::update()
     cvm::log ("Updating the harmonic bias \""+this->name+"\".\n");
 
   // Setup first stage of staged variable force constant calculation
-  if (b_chg_force_k && target_nstages && stage == -1) {
-    stage = 0;
+  if (b_chg_force_k && target_nstages && cvm::step_absolute() == 0) {
     cvm::real lambda;
     if (lambda_schedule.size()) {
       lambda = lambda_schedule[0];
@@ -250,25 +249,6 @@ cvm::real colvarbias_harmonic::update()
         cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda));
     cvm::log ("Setting force constant to " + cvm::to_str (force_k));
   }
-  
-  // Force and energy calculation
-  for (size_t i = 0; i < colvars.size(); i++) {
-    colvar_forces[i] =
-      (-0.5) * force_k /
-      (colvars[i]->width * colvars[i]->width) *
-      colvars[i]->dist2_lgrad (colvars[i]->value(),
-                               colvar_centers[i]);
-    bias_energy += 0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
-              colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]);
-    if (cvm::debug())
-      cvm::log ("dist_grad["+cvm::to_str (i)+
-                "] = "+cvm::to_str (colvars[i]->dist2_lgrad (colvars[i]->value(),
-                               colvar_centers[i]))+"\n");
-  }
-
-  if (cvm::debug())
-    cvm::log ("Current forces for the harmonic bias \""+
-              this->name+"\": "+cvm::to_str (colvar_forces)+".\n");
 
   if (b_chg_centers) {
 
@@ -280,26 +260,21 @@ cvm::real colvarbias_harmonic::update()
       centers_incr.resize (colvars.size());
       for (size_t i = 0; i < colvars.size(); i++) {
         centers_incr[i].type (colvars[i]->type());
-        centers_incr[i] = (target_centers[i] - colvar_centers[i]) /
+        centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) /
           cvm::real ( target_nstages ? (target_nstages - stage) :
                                       (target_nsteps - cvm::step_absolute()));
       }
       if (cvm::debug())
         cvm::log ("Center increment for the harmonic bias \""+
-                  this->name+"\": "+cvm::to_str (centers_incr)+".\n");
+                  this->name+"\": "+cvm::to_str (centers_incr)+" at stage "+cvm::to_str (stage)+ ".\n");
+
     }
 
-    if (cvm::debug())
-      cvm::log ("Current centers for the harmonic bias \""+
-                this->name+"\": "+cvm::to_str (colvar_centers)+".\n");
-
     if (target_nstages) {
-      if (cvm::step_absolute() > 0
+      if ((cvm::step_relative() > 0)
             && (cvm::step_absolute() % target_nsteps) == 0
             && stage < target_nstages) {
 
-          // any per-stage calculation, e.g. free energy stuff
-          // should be done here
           for (size_t i = 0; i < colvars.size(); i++) {
             colvar_centers_raw[i] += centers_incr[i];
             colvar_centers[i] = colvar_centers_raw[i];
@@ -308,9 +283,10 @@ cvm::real colvarbias_harmonic::update()
           }
           stage++;
           cvm::log ("Moving restraint stage " + cvm::to_str(stage) +
-              " : setting centers to " + cvm::to_str (colvar_centers));
+              " : setting centers to " + cvm::to_str (colvar_centers) +
+              " at step " +  cvm::to_str (cvm::step_absolute()));
       }
-    } else if (cvm::step_absolute() < target_nsteps) {
+    } else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) {
       // move the restraint centers in the direction of the targets
       // (slow growth)
       for (size_t i = 0; i < colvars.size(); i++) {
@@ -320,6 +296,10 @@ cvm::real colvarbias_harmonic::update()
         colvar_centers[i].apply_constraints();
       }
     }
+
+    if (cvm::debug())
+      cvm::log ("Current centers for the harmonic bias \""+
+                this->name+"\": "+cvm::to_str (colvar_centers)+".\n");
   }
 
   if (b_chg_force_k) {
@@ -383,6 +363,26 @@ cvm::real colvarbias_harmonic::update()
 
   if (cvm::debug())
     cvm::log ("Done updating the harmonic bias \""+this->name+"\".\n");
+  
+  // Force and energy calculation
+  for (size_t i = 0; i < colvars.size(); i++) {
+    colvar_forces[i] =
+      (-0.5) * force_k /
+      (colvars[i]->width * colvars[i]->width) *
+      colvars[i]->dist2_lgrad (colvars[i]->value(),
+                               colvar_centers[i]);
+    bias_energy += 0.5 * force_k / (colvars[i]->width * colvars[i]->width) *
+              colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]);
+    if (cvm::debug())
+      cvm::log ("dist_grad["+cvm::to_str (i)+
+                "] = "+cvm::to_str (colvars[i]->dist2_lgrad (colvars[i]->value(),
+                               colvar_centers[i]))+"\n");
+  }
+
+  if (cvm::debug())
+    cvm::log ("Current forces for the harmonic bias \""+
+              this->name+"\": "+cvm::to_str (colvar_forces)+".\n");
+
   return bias_energy;
 }
 
@@ -423,7 +423,7 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is)
   }
 
   if (b_chg_centers) {
-    cvm::log ("Reading the updated restraint centers from the restart.\n");
+//    cvm::log ("Reading the updated restraint centers from the restart.\n");
     if (!get_keyval (conf, "centers", colvar_centers))
       cvm::fatal_error ("Error: restraint centers are missing from the restart.\n");
     if (!get_keyval (conf, "centers_raw", colvar_centers_raw))
@@ -431,13 +431,13 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is)
   }
 
   if (b_chg_force_k) {
-    cvm::log ("Reading the updated force constant from the restart.\n");
+//    cvm::log ("Reading the updated force constant from the restart.\n");
     if (!get_keyval (conf, "forceConstant", force_k))
       cvm::fatal_error ("Error: force constant is missing from the restart.\n");
   }
 
   if (target_nstages) {
-    cvm::log ("Reading current stage from the restart.\n");
+//    cvm::log ("Reading current stage from the restart.\n");
     if (!get_keyval (conf, "stage", stage))
       cvm::fatal_error ("Error: current stage is missing from the restart.\n");
   }
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 32d02893ce..a371cff605 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -754,7 +754,8 @@ colvar::rmsd::rmsd (std::string const &conf)
     cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
   } else {
     // Default: fit everything
-    cvm::log ("No explicit fitting parameters: enabling both centerReference and rotateReference to compute standard minimum RMSD.");
+    cvm::log ("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: "
+              "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
     atoms.b_center = true;
     atoms.b_rotate = true;
     // default case: reference positions for calculating the rmsd are also those used
@@ -936,7 +937,8 @@ colvar::eigenvector::eigenvector (std::string const &conf)
     cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
   } else {
     // default: fit everything
-    cvm::log ("Enabling centerReference and rotateReference (minimize RMSD before calculating the projection).\n");
+    cvm::log ("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: "
+              "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n");
     atoms.b_center = true;
     atoms.b_rotate = true;
     atoms.ref_pos = ref_pos;
@@ -996,34 +998,32 @@ colvar::eigenvector::eigenvector (std::string const &conf)
   eig_center *= 1.0 / atoms.size();
   cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n");
 
-  bool b_difference_vector;
+  bool b_difference_vector = false;
   get_keyval (conf, "differenceVector", b_difference_vector, false);
+
   if (b_difference_vector) {
 
-    cvm::log ("Subtracting the reference positions from the provided vector.\n");
-
     if (atoms.b_center) {
-      cvm::log ("centerReference is on: recentering the provided vector to the reference positions.\n");
       // both sets should be centered on the origin for fitting
       for (size_t i = 0; i < atoms.size(); i++) {
         eigenvec[i] -= eig_center;
         ref_pos[i]  -= ref_pos_center;
       }
     }
-
     if (atoms.b_rotate) {
-      cvm::log ("rotateReference is on: aligning the provided vector to the reference positions.\n");
       atoms.rot.calc_optimal_rotation (eigenvec, ref_pos);
       for (size_t i = 0; i < atoms.size(); i++) {
         eigenvec[i] = atoms.rot.rotate (eigenvec[i]);
-        eigenvec[i] -= ref_pos[i];
       }
     }
-
+    cvm::log ("\"differenceVector\" is on: subtracting the reference positions from the provided vector: v = v - x0.\n");
+    for (size_t i = 0; i < atoms.size(); i++) {
+      eigenvec[i] -= ref_pos[i];
+    }
     if (atoms.b_center) {
       // bring back the ref positions to where they were
       for (size_t i = 0; i < atoms.size(); i++) {
-        ref_pos[i] -= ref_pos_center;
+        ref_pos[i] += ref_pos_center;
       }
     }
 
@@ -1034,7 +1034,7 @@ colvar::eigenvector::eigenvector (std::string const &conf)
     }
   }
 
-  cvm::log ("The first three components (v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n");
+  // cvm::log ("The first three components (v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n");
 
   // for inverse gradients
   eigenvec_invnorm2 = 0.0;
@@ -1042,6 +1042,15 @@ colvar::eigenvector::eigenvector (std::string const &conf)
     eigenvec_invnorm2 += eigenvec[i].norm2();
   }
   eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
+
+  if (b_difference_vector) {
+    cvm::log ("\"differenceVector\" is on: normalizing the vector.\n");
+    for (size_t i = 0; i < atoms.size(); i++) {
+      eigenvec[i] *= eigenvec_invnorm2;
+    }
+  } else {
+    cvm::log ("The norm of the vector is |v| = "+cvm::to_str (eigenvec_invnorm2)+".\n");
+  }
 }
 
   
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 9cdd70539c..7449ae23f7 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -2,7 +2,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2013-01-30"
+#define COLVARS_VERSION "2013-04-17"
 #endif
 
 #ifndef COLVARS_DEBUG
diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp
index b1f4290549..e08746c547 100644
--- a/lib/colvars/colvarparse.cpp
+++ b/lib/colvars/colvarparse.cpp
@@ -102,12 +102,19 @@ size_t      colvarparse::dummy_pos = 0;
                                                                         \
     if (data.size()) {                                                  \
       std::istringstream is (data);                                     \
+      size_t data_count = 0;                                            \
       TYPE x (def_value);                                               \
-      if (is >> x)                                                      \
+      while (is >> x) {                                                 \
         value = x;                                                      \
-      else                                                              \
+        data_count++;                                                   \
+      }                                                                 \
+      if (data_count == 0)                                              \
         cvm::fatal_error ("Error: in parsing \""+                       \
                           std::string (key)+"\".\n");                   \
+      if (data_count > 1)                                               \
+        cvm::fatal_error ("Error: multiple values "                     \
+                          "are not allowed for keyword \""+             \
+                          std::string (key)+"\".\n");                   \
       if (parse_mode != parse_silent) {                                 \
         cvm::log ("# "+std::string (key)+" = "+                         \
                   cvm::to_str (value)+"\n");                            \
@@ -533,7 +540,7 @@ bool colvarparse::key_lookup (std::string const &conf,
         if (line[brace] == '}') brace_count--;
       }
 
-      // set data_begin afterthe opening brace
+      // set data_begin after the opening brace
       data_begin = line.find_first_of ('{', data_begin) + 1;
       data_begin = line.find_first_not_of (white_space,
                                            data_begin);