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

lib/colvars/colvar.cpp

@@ -22,7 +22,7 @@ colvar::colvar (std::string const &conf)
        cvi < cvm::colvars.end();
        cvi++) {
     if ((*cvi)->name == this->name)
-      cvm::fatal_error ("Error: this colvar has the same name, \""+this->name+
+      cvm::fatal_error ("Error: this colvar cannot have the same name, \""+this->name+
                         "\", as another colvar.\n");
   }
 
@@ -98,7 +98,8 @@ colvar::colvar (std::string const &conf)
                          "on an axis",       "distanceZ",      distance_z);
   initialize_components ("distance projection "
                          "on a plane",       "distanceXY",     distance_xy);
-  initialize_components ("minimum distance", "minDistance",    min_distance);
+  initialize_components ("average distance weighted by inverse sixth power",
+                         "distance6", distance6);
 
   initialize_components ("coordination "
                          "number",           "coordNum",       coordnum);
@@ -128,6 +129,10 @@ colvar::colvar (std::string const &conf)
 
   initialize_components ("radius of "
                          "gyration",         "gyration",       gyration);
+  initialize_components ("moment of "
+                         "inertia",          "inertia",        inertia);
+  initialize_components ("moment of inertia around an axis",
+                                             "inertia_z",      inertia_z);
   initialize_components ("eigenvector",      "eigenvector",    eigenvector);
 
   if (!cvcs.size())
@@ -249,6 +254,13 @@ colvar::colvar (std::string const &conf)
     }
   }
 
+  if (tasks[task_lower_boundary]) {
+    get_keyval (conf, "hardLowerBoundary", hard_lower_boundary, false);
+  }
+  if (tasks[task_upper_boundary]) {
+    get_keyval (conf, "hardUpperBoundary", hard_upper_boundary, false);
+  }
+
   // consistency checks for boundaries and walls
   if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) {
     if (lower_boundary >= upper_boundary) {

lib/colvars/colvar.h

@@ -51,6 +51,9 @@ public:
   /// \brief Current value (previously obtained from calc() or read_traj())
   colvarvalue const & value() const;
 
+  /// \brief Current actual value (not extended DOF)
+  colvarvalue const & actual_value() const;
+
   /// \brief Current velocity (previously obtained from calc() or read_traj())
   colvarvalue const & velocity() const;
 
@@ -257,6 +260,8 @@ public:
   colvarvalue lower_wall;
   /// \brief Force constant for the lower boundary potential (|x-xb|^2)
   cvm::real   lower_wall_k;
+  /// \brief Whether this colvar has a hard lower boundary
+  bool        hard_lower_boundary;
 
   /// \brief Location of the upper boundary
   colvarvalue upper_boundary;
@@ -264,6 +269,8 @@ public:
   colvarvalue upper_wall;
   /// \brief Force constant for the upper boundary potential (|x-xb|^2)
   cvm::real   upper_wall_k;
+  /// \brief Whether this colvar has a hard upper boundary
+  bool        hard_upper_boundary;
 
   /// \brief Is the interval defined by the two boundaries periodic?
   bool periodic_boundaries() const;
@@ -464,7 +471,7 @@ public:
   class distance;
   class distance_z;
   class distance_xy;
-  class min_distance;
+  class distance6;
   class angle;
   class dihedral;
   class coordnum;
@@ -476,6 +483,8 @@ public:
   class tilt;
   class spin_angle;
   class gyration;
+  class inertia;
+  class inertia_z;
   class eigenvector;
   class alpha_dihedrals;
   class alpha_angles;
@@ -535,6 +544,12 @@ inline colvarvalue const & colvar::value() const
 }
 
 
+inline colvarvalue const & colvar::actual_value() const
+{
+  return x;
+}
+
+
 inline colvarvalue const & colvar::velocity() const
 {
   return v_reported;

lib/colvars/colvarbias.cpp

@@ -26,6 +26,14 @@ colvarbias::colvarbias (std::string const &conf, char const *key)
 
   get_keyval (conf, "name", name, key_str+cvm::to_str (rank));
 
+  for (std::vector<colvarbias *>::iterator bi = cvm::biases.begin();
+       bi != cvm::biases.end();
+       bi++) {
+    if ((*bi)->name == this->name)
+      cvm::fatal_error ("Error: this bias cannot have the same name, \""+this->name+
+                        "\", of another bias.\n");
+  }
+
   // lookup the associated colvars
   std::vector<std::string> colvars_str;
   if (get_keyval (conf, "colvars", colvars_str)) {
@@ -78,7 +86,8 @@ void colvarbias::communicate_forces()
 colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
                                           char const *key)
   : colvarbias (conf, key), 
-    target_nsteps (0)
+    target_nsteps (0),
+    target_nstages (0)
 {
   get_keyval (conf, "forceConstant", force_k, 1.0);
   for (size_t i = 0; i < colvars.size(); i++) {
@@ -112,7 +121,6 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
 
   if (get_keyval (conf, "targetCenters", target_centers, colvar_centers)) {
     b_chg_centers = true;
-    target_nstages = 0;
     for (size_t i = 0; i < target_centers.size(); i++) {
       target_centers[i].apply_constraints();
     }
@@ -134,8 +142,6 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
     if (lambda_schedule.size()) {
       // There is one more lambda-point than stages
       target_nstages = lambda_schedule.size() - 1;
-    } else {
-      target_nstages = 0;
     }
   } else {
     b_chg_force_k = false;

lib/colvars/colvarbias.h

@@ -146,10 +146,10 @@ protected:
 
   /// \brief Number of stages over which to perform the change
   /// If zero, perform a continuous change
-  size_t target_nstages;
+  int target_nstages;
 
   /// \brief Number of current stage of the perturbation
-  size_t stage;
+  int stage;
   
   /// \brief Intermediate quantity to compute the restraint free energy
   /// (in TI, would be the accumulating FE derivative)

lib/colvars/colvarbias_meta.cpp

@@ -222,6 +222,17 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
                         traj_file_name + "\".\n");
   }
 
+  // for well-tempered metadynamics
+  get_keyval (conf, "wellTempered", well_tempered, false);
+  get_keyval (conf, "biasTemperature", bias_temperature, -1.0);
+  if ((bias_temperature == -1.0) && well_tempered) {
+    cvm::fatal_error ("Error: biasTemperature is not set.\n");
+  }
+  if (well_tempered) {
+    cvm::log("Well-tempered metadynamics is used.\n");
+    cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n");
+  }
+  
   if (cvm::debug())
     cvm::log ("Done initializing the metadynamics bias \""+this->name+"\""+
               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
@@ -271,7 +282,7 @@ colvarbias_meta::create_hill (colvarbias_meta::hill const &h)
     // also add it to the list of hills that are off-grid, which may
     // need to be computed analytically when the colvar returns
     // off-grid
-    cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers);
+    cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers, true);
     if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
       hills_off_grid.push_back (h);
     }
@@ -441,11 +452,25 @@ cvm::real colvarbias_meta::update()
     switch (comm) {
 
     case single_replica:
-      create_hill (hill (hill_weight, colvars, hill_width));
+      if (well_tempered) {
+        std::vector<int> curr_bin = hills_energy->get_colvars_index();
+        cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin);
+        cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
+        create_hill (hill ((hill_weight*exp_weight), colvars, hill_width));
+      } else { 
+        create_hill (hill (hill_weight, colvars, hill_width));
+      } 
       break;
 
     case multiple_replicas:
-      create_hill (hill (hill_weight, colvars, hill_width, replica_id));
+      if (well_tempered) {
+        std::vector<int> curr_bin = hills_energy->get_colvars_index();
+        cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin);
+        cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
+        create_hill (hill ((hill_weight*exp_weight), colvars, hill_width, replica_id));
+      } else { 
+        create_hill (hill (hill_weight, colvars, hill_width, replica_id));
+      } 
       if (replica_hills_os.good()) {
         replica_hills_os << hills.back();
       } else {
@@ -699,7 +724,8 @@ void colvarbias_meta::calc_hills_force (size_t const &i,
 void colvarbias_meta::project_hills (colvarbias_meta::hill_iter  h_first,
                                      colvarbias_meta::hill_iter  h_last,
                                      colvar_grid_scalar         *he,
-                                     colvar_grid_gradient       *hg)
+                                     colvar_grid_gradient       *hg,
+                                     bool print_progress)
 {
   if (cvm::debug())
     cvm::log ("Metadynamics bias \""+this->name+"\""+
@@ -716,12 +742,15 @@ void colvarbias_meta::project_hills (colvarbias_meta::hill_iter  h_first,
   cvm::real hills_energy_here = 0.0;
   std::vector<colvarvalue> hills_forces_here (colvars.size(), 0.0);
 
+  size_t count = 0;
+  size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1)));
+
   if (hg != NULL) {
 
     // loop over the points of the grid
     for ( ;
           (he->index_ok (he_ix)) && (hg->index_ok (hg_ix));
-          ) {
+          count++) {
 
       for (size_t i = 0; i < colvars.size(); i++) {
         colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
@@ -741,6 +770,18 @@ void colvarbias_meta::project_hills (colvarbias_meta::hill_iter  h_first,
 
       he->incr (he_ix);
       hg->incr (hg_ix);
+
+      if ((count % print_frequency) == 0) {
+        if (print_progress) {
+          cvm::real const progress = cvm::real (count) / cvm::real (hg->number_of_points());
+          std::ostringstream os;
+          os.setf (std::ios::fixed, std::ios::floatfield);
+          os << std::setw (6) << std::setprecision (2)
+             << 100.0 * progress
+             << "% done.";
+          cvm::log (os.str());
+        }
+      }
     }
 
   } else {
@@ -758,9 +799,26 @@ void colvarbias_meta::project_hills (colvarbias_meta::hill_iter  h_first,
       he->acc_value (he_ix, hills_energy_here);
 
       he->incr (he_ix);
+
+      count++;
+      if ((count % print_frequency) == 0) {
+        if (print_progress) {
+          cvm::real const progress = cvm::real (count) / cvm::real (he->number_of_points());
+          std::ostringstream os;
+          os.setf (std::ios::fixed, std::ios::floatfield);
+          os << std::setw (6) << std::setprecision (2)
+             << 100.0 * progress
+             << "% done.";
+          cvm::log (os.str());
+        }
+      }
     }
   }
 
+  if (print_progress) {
+    cvm::log ("100.00% done.");
+  }
+
   if (! keep_hills) {
     hills.erase (hills.begin(), hills.end());
   }
@@ -774,7 +832,7 @@ void colvarbias_meta::recount_hills_off_grid (colvarbias_meta::hill_iter  h_firs
   hills_off_grid.clear();
 
   for (hill_iter h = h_first; h != h_last; h++) {
-    cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers);
+    cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers, true);
     if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
       hills_off_grid.push_back (*h);
     }
@@ -1175,8 +1233,8 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
       is.clear();
       is.seekg (hills_energy_gradients_pos, std::ios::beg);
       grids_from_restart_file = false;
-      if (hills_energy_backup == NULL) {
+      if (!rebin_grids) { 
-        if (!rebin_grids) 
+        if (hills_energy_backup == NULL)
           cvm::fatal_error ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
                             this->name+"\""+
                             ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
@@ -1252,7 +1310,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
       cvm::log ("Rebinning the energy and forces grids from "+
                 cvm::to_str (hills.size())+" hills (this may take a while)...\n");
       project_hills (hills.begin(), hills.end(),
-                     new_hills_energy, new_hills_energy_gradients);
+                     new_hills_energy, new_hills_energy_gradients, true);
       cvm::log ("rebinning done.\n");
 
     } else {
@@ -1379,7 +1437,7 @@ std::istream & colvarbias_meta::read_hill (std::istream &is)
     // add this also to the list of hills that are off-grid, which will
     // be computed analytically
     cvm::real const min_dist =
-      hills_energy->bin_distance_from_boundaries ((hills.back()).centers);
+      hills_energy->bin_distance_from_boundaries ((hills.back()).centers, true);
     if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
       hills_off_grid.push_back (hills.back());
     }
@@ -1478,6 +1536,10 @@ void colvarbias_meta::write_pmf()
     cvm::real const max = pmf->maximum_value();
     pmf->add_constant (-1.0 * max);
     pmf->multiply_constant (-1.0);
+    if (well_tempered) {
+      cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
+      pmf->multiply_constant (well_temper_scale);
+    }
     std::string const fes_file_name (fes_file_name_prefix +
                                      ((comm != single_replica) ? ".partial" : "") +
                                      (dump_fes_save ?
@@ -1499,6 +1561,10 @@ void colvarbias_meta::write_pmf()
     cvm::real const max = pmf->maximum_value();
     pmf->add_constant (-1.0 * max);
     pmf->multiply_constant (-1.0);
+    if (well_tempered) {
+      cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
+      pmf->multiply_constant (well_temper_scale);
+    }
     std::string const fes_file_name (fes_file_name_prefix +
                                      (dump_fes_save ?
                                       "."+cvm::to_str (cvm::step_absolute()) : "") +

lib/colvars/colvarbias_meta.h

@@ -147,6 +147,12 @@ protected:
   /// time steps, appending the step number to each file
   bool       dump_fes_save;
 
+  /// \brief Whether to use well-tempered metadynamics 
+  bool       well_tempered; 
+
+  /// \brief Biasing temperature in well-tempered metadynamics 
+  cvm::real  bias_temperature;
+
   /// \brief Try to read the restart information by allocating new
   /// grids before replacing the current ones (used e.g. in
   /// multiple_replicas)
@@ -160,7 +166,8 @@ protected:
 
   /// \brief Project the selected hills onto grids
   void project_hills (hill_iter h_first, hill_iter h_last,
-                      colvar_grid_scalar *ge, colvar_grid_gradient *gf);
+                      colvar_grid_scalar *ge, colvar_grid_gradient *gf,
+                      bool print_progress = false);
 
 
   // Multiple Replicas variables and functions

lib/colvars/colvarcomp.h

@@ -473,18 +473,18 @@ public:
 };
 
 
-/// \brief Colvar component: projection of the distance vector along
+/// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
-/// a fixed axis (colvarvalue::type_scalar type, range (-*:*))
+/// as in NMR refinements (colvarvalue::type_scalar type, range (0:*))
-class colvar::min_distance
+class colvar::distance6
   : public colvar::distance
 {
 protected:
   /// Components of the distance vector orthogonal to the axis
   cvm::real smoothing;
 public:
-  min_distance (std::string const &conf);
+  distance6 (std::string const &conf);
-  min_distance();
+  distance6();
-  virtual inline ~min_distance() {}
+  virtual inline ~distance6() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
@@ -500,7 +500,6 @@ public:
 
 
 
-
 /// \brief Colvar component: Radius of gyration of an atom group
 /// (colvarvalue::type_scalar type, range [0:*))
 class colvar::gyration
@@ -530,6 +529,60 @@ public:
 };
 
 
+/// \brief Colvar component: moment of inertia of an atom group
+/// (colvarvalue::type_scalar type, range [0:*))
+class colvar::inertia
+  : public colvar::cvc
+{
+protected:
+  /// Atoms involved
+  cvm::atom_group atoms;
+public:
+  /// Constructor
+  inertia (std::string const &conf);
+  inertia();
+  virtual inline ~inertia() {}
+  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;
+  virtual cvm::real compare (colvarvalue const &x1,
+                             colvarvalue const &x2) const;
+};
+
+
+/// \brief Colvar component: moment of inertia of an atom group
+/// around a user-defined axis (colvarvalue::type_scalar type, range [0:*))
+class colvar::inertia_z
+  : public colvar::inertia
+{
+protected:
+  /// Vector on which the inertia tensor is projected
+  cvm::rvector axis;
+public:
+  /// Constructor
+  inertia_z (std::string const &conf);
+  inertia_z();
+  virtual inline ~inertia_z() {}
+  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;
+  virtual cvm::real compare (colvarvalue const &x1,
+                             colvarvalue const &x2) const;
+};
+
+
 /// \brief Colvar component: projection of 3N coordinates onto an
 /// eigenvector (colvarvalue::type_scalar type, range (-*:*))
 class colvar::eigenvector
@@ -1223,8 +1276,7 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const
   inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1,          \
                                         colvarvalue const &x2) const    \
   {                                                                     \
-    const cvm::real tmp = x1.real_value - x2.real_value;                \
+    return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
-    return tmp*tmp;                                                     \
   }                                                                     \
                                                                         \
   inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1,  \
@@ -1249,12 +1301,14 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const
   simple_scalar_dist_functions (distance)
   // NOTE: distance_z has explicit functions, see below 
   simple_scalar_dist_functions (distance_xy)
-  simple_scalar_dist_functions (min_distance)
+  simple_scalar_dist_functions (distance6)
   simple_scalar_dist_functions (angle)
   simple_scalar_dist_functions (coordnum)
   simple_scalar_dist_functions (selfcoordnum)
   simple_scalar_dist_functions (h_bond)
   simple_scalar_dist_functions (gyration)
+  simple_scalar_dist_functions (inertia)
+  simple_scalar_dist_functions (inertia_z)
   simple_scalar_dist_functions (rmsd)
   simple_scalar_dist_functions (logmsd)
   simple_scalar_dist_functions (orientation_angle)

lib/colvars/colvarcomp_distances.cpp

@@ -89,10 +89,10 @@ void colvar::distance::calc_Jacobian_derivative()
 void colvar::distance::apply_force (colvarvalue const &force)
 {
   if (!group1.noforce)
-    group1.apply_colvar_force (force);
+    group1.apply_colvar_force (force.real_value);
 
   if (!group2.noforce)
-    group2.apply_colvar_force (force);
+    group2.apply_colvar_force (force.real_value);
 }
 
 
@@ -394,83 +394,6 @@ void colvar::distance_xy::apply_force (colvarvalue const &force)
 
 
 
-colvar::min_distance::min_distance (std::string const &conf)
-  : distance (conf)
-{
-  function_type = "min_distance";
-  x.type (colvarvalue::type_scalar);
-
-  get_keyval (conf, "smoothing", smoothing, (1.0 * cvm::unit_angstrom()));
-}
-
-colvar::min_distance::min_distance()
-  : distance()
-{
-  function_type = "min_distance";
-  x.type (colvarvalue::type_scalar);
-}
-
-void colvar::min_distance::calc_value()
-{
-  group1.reset_atoms_data();
-  group2.reset_atoms_data();
-
-  group1.read_positions();
-  group2.read_positions();
-
-  x.real_value = 0.0;
-
-  bool zero_dist = false;
-
-  for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
-    for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
-      cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
-      cvm::real const d = dv.norm();
-      if (d > 0.0)
-        x.real_value += std::exp (smoothing / d);
-      else
-        zero_dist = true;
-    }
-  }
-
-  x.real_value = zero_dist ? 0.0 : smoothing/(std::log (x.real_value));
-}
-
-void colvar::min_distance::calc_gradients()
-{
-  if (x.real_value > 0.0) {
-    cvm::real const sum = std::exp (smoothing/x.real_value);
-    cvm::real const dxdsum = -1.0 *
-      (x.real_value/smoothing) * (x.real_value/smoothing) *
-      (1.0 / sum);
-
-    for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
-      for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
-        cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
-        cvm::real const d = dv.norm();
-        if (d > 0.0) {
-          cvm::rvector const dvu = dv / dv.norm();
-          ai1->grad += dxdsum * std::exp (smoothing / d) *
-            smoothing * (-1.0/(d*d)) * (-1.0) * dvu;
-          ai2->grad += dxdsum * std::exp (smoothing / d) *
-            smoothing * (-1.0/(d*d)) * dvu;
-        }
-      }
-    }
-  }
-}
-
-void colvar::min_distance::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);
-}
-
-
-
 colvar::distance_dir::distance_dir (std::string const &conf)
   : distance (conf)
 {
@@ -528,6 +451,76 @@ void colvar::distance_dir::apply_force (colvarvalue const &force)
 
 
 
+colvar::distance6::distance6 (std::string const &conf)
+  : distance (conf)
+{
+  function_type = "distance6";
+  b_inverse_gradients = false;
+  b_Jacobian_derivative = false;
+  x.type (colvarvalue::type_scalar);
+}
+
+colvar::distance6::distance6()
+{
+  function_type = "distance6";
+  b_inverse_gradients = false;
+  b_Jacobian_derivative = false;
+  b_1site_force = false;
+  x.type (colvarvalue::type_scalar);
+}
+
+void colvar::distance6::calc_value()
+{
+  group1.reset_atoms_data();
+  group2.reset_atoms_data();
+
+  group1.read_positions();
+  group2.read_positions();
+
+  x.real_value = 0.0;
+  if (b_no_PBC) {
+    for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
+      for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+        cvm::rvector const dv = ai2->pos - ai1->pos;
+        cvm::real const d2 = dv.norm2();
+        x.real_value += 1.0/(d2*d2*d2);
+        cvm::rvector const dsumddv = -6.0/(d2*d2*d2*d2) * dv;
+        ai1->grad += -1.0 * dsumddv;
+        ai2->grad +=        dsumddv;
+      }
+    }
+  } else {
+    for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
+      for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+        cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
+        cvm::real const d2 = dv.norm2();
+        x.real_value += 1.0/(d2*d2*d2);
+        cvm::rvector const dsumddv = -6.0/(d2*d2*d2*d2) * dv;
+        ai1->grad += -1.0 * dsumddv;
+        ai2->grad +=        dsumddv;
+      }
+    }
+  }
+
+  x.real_value = std::pow (x.real_value, -1.0/6.0);
+}
+
+void colvar::distance6::calc_gradients()
+{
+}
+
+void colvar::distance6::apply_force (colvarvalue const &force)
+{
+  cvm::real const dxdsum = (-1.0/6.0) * std::pow (x.real_value, -7.0/6.0);
+
+  if (!group1.noforce)
+    group1.apply_colvar_force (dxdsum * force.real_value);
+
+  if (!group2.noforce)
+    group2.apply_colvar_force (dxdsum * force.real_value);
+}
+
+
 
 colvar::gyration::gyration (std::string const &conf)
   : cvc (conf)
@@ -558,9 +551,9 @@ void colvar::gyration::calc_value()
 
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
-    x.real_value += (ai->mass/atoms.total_mass) * (ai->pos).norm2();
+    x.real_value += (ai->pos).norm2();
   }
-  x.real_value = std::sqrt (x.real_value);
+  x.real_value = std::sqrt (x.real_value / cvm::real (atoms.size()));
 }
 
 
@@ -600,6 +593,105 @@ void colvar::gyration::apply_force (colvarvalue const &force)
 
 
 
+colvar::inertia::inertia (std::string const &conf)
+  : cvc (conf)
+{
+  function_type = "inertia";
+  parse_group (conf, "atoms", atoms);
+  atom_groups.push_back (&atoms);
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::inertia::inertia()
+{
+  function_type = "inertia";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::inertia::calc_value()
+{
+  atoms.reset_atoms_data();
+  atoms.read_positions();
+  atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
+
+  x.real_value = 0.0;
+  for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
+    x.real_value += ai->mass * (ai->pos).norm2();
+  }
+  x.real_value = std::sqrt (x.real_value / atoms.total_mass);
+}
+
+
+void colvar::inertia::calc_gradients()
+{
+  cvm::real const drdx = 1.0/(atoms.total_mass * x.real_value);
+  for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
+    ai->grad = drdx * ai->mass * ai->pos;
+  }
+}
+
+
+void colvar::inertia::apply_force (colvarvalue const &force)
+{
+  if (!atoms.noforce)
+    atoms.apply_colvar_force (force.real_value);
+}
+
+
+colvar::inertia_z::inertia_z (std::string const &conf)
+  : inertia (conf)
+{
+  function_type = "inertia_z";
+  if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
+    if (axis.norm2() == 0.0)
+      cvm::fatal_error ("Axis vector is zero!");
+    axis = axis.unit();
+  }
+  x.type (colvarvalue::type_scalar);
+}
+
+
+colvar::inertia_z::inertia_z()
+{
+  function_type = "inertia_z";
+  x.type (colvarvalue::type_scalar);
+}
+
+
+void colvar::inertia_z::calc_value()
+{
+  atoms.reset_atoms_data();
+  atoms.read_positions();
+  atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
+
+  x.real_value = 0.0;
+  for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
+    cvm::real const iprod = ai->pos * axis;
+    x.real_value += ai->mass * iprod * iprod;
+  }
+  x.real_value = std::sqrt (x.real_value / atoms.total_mass);
+}
+
+
+void colvar::inertia_z::calc_gradients()
+{
+  cvm::real const drdx = 1.0/(atoms.total_mass * x.real_value);
+  for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
+    ai->grad = drdx * ai->mass * (ai->pos * axis) * axis;
+  }
+}
+
+
+void colvar::inertia_z::apply_force (colvarvalue const &force)
+{
+  if (!atoms.noforce)
+    atoms.apply_colvar_force (force.real_value);
+}
+
+
+
 colvar::rmsd::rmsd (std::string const &conf)
   : orientation (conf)
 {
@@ -1111,3 +1203,6 @@ void colvar::eigenvector::calc_Jacobian_derivative()
 
   jd.real_value = sum * std::sqrt (eigenvec_invnorm2); 
 }
+
+
+

lib/colvars/colvargrid.h

@@ -45,6 +45,9 @@ protected:
   /// Colvars collected in this grid
   std::vector<colvar *> cv;
 
+  /// Do we request actual value (for extended-system colvars)?
+  std::vector<bool> actual_value;
+
   /// Get the low-level index corresponding to an index
   inline size_t address (std::vector<int> const &ix) const
   {
@@ -70,6 +73,12 @@ public:
   /// Whether some colvars are periodic
   std::vector<bool>        periodic;
 
+  /// Whether some colvars have hard lower boundaries
+  std::vector<bool>        hard_lower_boundaries;
+
+  /// Whether some colvars have hard upper boundaries
+  std::vector<bool>        hard_upper_boundaries;
+
   /// Widths of the colvars in this grid
   std::vector<cvm::real>   widths;
 
@@ -170,8 +179,11 @@ public:
                                           cv (g.cv),
                                           lower_boundaries (g.lower_boundaries),
                                           upper_boundaries (g.upper_boundaries),
+                                          hard_lower_boundaries (g.hard_lower_boundaries),
+                                          hard_upper_boundaries (g.hard_upper_boundaries),
                                           periodic (g.periodic),
                                           widths (g.widths),
+                                          actual_value (g.actual_value),
                                           data()
   {
     save_delimiters = false;
@@ -224,8 +236,19 @@ public:
       }
 
       widths.push_back (cv[i]->width);
+      hard_lower_boundaries.push_back (cv[i]->hard_lower_boundary);
+      hard_upper_boundaries.push_back (cv[i]->hard_upper_boundary);
       periodic.push_back (cv[i]->periodic_boundaries());
 
+      // By default, get reported colvar value (for extended Lagrangian colvars)
+      actual_value.push_back (false);
+
+      // except if a colvar is specified twice in a row
+      // then the first instance is the actual value
+      if (i > 0 && cv[i-1] == cv[i]) {
+        actual_value[i-1] = true;
+      }
+
       if (margin) {
         if (periodic[i]) {
           // Shift the grid by half the bin width (values at edges instead of center of bins)
@@ -288,7 +311,7 @@ public:
   /// \brief Report the bin corresponding to the current value of variable i
   inline int current_bin_scalar(int const i) const
   {
-    return value_to_bin_scalar (cv[i]->value(), i);
+    return value_to_bin_scalar (actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
   }
 
   /// \brief Use the lower boundary and the width to report which bin
@@ -381,7 +404,8 @@ public:
   /// \brief Get the minimal distance (in number of bins) from the
   /// boundaries; a negative number is returned if the given point is
   /// off-grid
-  inline cvm::real bin_distance_from_boundaries (std::vector<colvarvalue> const &values)
+  inline cvm::real bin_distance_from_boundaries (std::vector<colvarvalue> const &values,
+                                                 bool skip_hard_boundaries = false)
   {
     cvm::real minimum = 1.0E+16;
     for (size_t i = 0; i < nd; i++) {
@@ -396,9 +420,9 @@ public:
       if (values[i].real_value > upper_boundaries[i])
         du *= -1.0;
 
-      if (dl < minimum) 
+      if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum))
         minimum = dl;
-      if (du < minimum) 
+      if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum))
         minimum = du;
     }
 

lib/colvars/colvarmodule.cpp

@@ -88,16 +88,18 @@ colvarmodule::colvarmodule (char const  *config_filename,
   cvm::log ("The trajectory file will be \""+
             cv_traj_name+"\".\n");
 
-  // open trajectory file
+  if (cv_traj_freq) {
-  if (cv_traj_append) {
+    // open trajectory file
-    cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+
+    if (cv_traj_append) {
-              "\".\n");
+      cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+
-    cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
+                "\".\n");
-  } else {
+      cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
-    proxy->backup_file (cv_traj_name.c_str());
+    } else {
-    cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
+      proxy->backup_file (cv_traj_name.c_str());
+      cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
+    }
+    cv_traj_os.setf (std::ios::scientific, std::ios::floatfield);
   }
-  cv_traj_os.setf (std::ios::scientific, std::ios::floatfield);
 
   // parse the options for collective variables
   init_colvars (conf);

lib/colvars/colvarmodule.h

@@ -2,7 +2,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2012-03-23"
+#define COLVARS_VERSION "2012-06-20"
 #endif
 
 #ifndef COLVARS_DEBUG

lib/colvars/colvartypes.h

@@ -103,7 +103,8 @@ public:
 
   inline cvm::rvector unit() const
   {
-    return cvm::rvector (x, y, z)/this->norm();
+    const cvm::real n = this->norm();
+    return (n > 0. ? cvm::rvector (x, y, z)/n : cvm::rvector (1., 0., 0.));
   }
 
   static inline size_t output_width (size_t const &real_width)

lib/colvars/colvarvalue.h

@@ -636,17 +636,13 @@ inline cvm::real colvarvalue::dist2 (colvarvalue const &x2) const
     colvarvalue::check_types (*this, x2);
 
   switch (this->value_type) {
-  case colvarvalue::type_scalar: {
+  case colvarvalue::type_scalar:
-    const cvm::real tmp = this->real_value - x2.real_value;
+    return (this->real_value - x2.real_value)*(this->real_value - x2.real_value);
-    return tmp*tmp;
-  }
   case colvarvalue::type_vector:
     return (this->rvector_value - x2.rvector_value).norm2();
-  case colvarvalue::type_unitvector: {
+  case colvarvalue::type_unitvector:
     // angle between (*this) and x2 is the distance
-    const cvm::real tmp = std::acos (this->rvector_value * x2.rvector_value);
+    return std::acos (this->rvector_value * x2.rvector_value) * std::acos (this->rvector_value * x2.rvector_value);
-    return tmp*tmp;
-  }
   case colvarvalue::type_quaternion:
     // angle between (*this) and x2 is the distance, the quaternion
     // object has it implemented internally