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

lib/colvars/Makefile.g++

@@ -1,4 +1,4 @@
-# library build makefile for colvars module
+# library build -*- makefile -*- for colvars module
 
 # which file will be copied to Makefile.lammps
 
@@ -14,9 +14,9 @@ SHELL =		/bin/sh
 
 # ------ DEFINITIONS ------
 
-SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp	\
+SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \
- colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp	\
+ colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp   \
- colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp	\
+ colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp   \
  colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp
 
 LIB = libcolvars.a
@@ -25,11 +25,13 @@ EXE = #colvars_standalone
 
 # ------ MAKE PROCEDURE ------
 
-default: $(LIB) $(EXE)
+default: $(LIB) $(EXE) Makefile.lammps
+
+Makefile.lammps:
+	@cp $(EXTRAMAKE) Makefile.lammps
 
 $(LIB):	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
-	@cp $(EXTRAMAKE) Makefile.lammps
 
 colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
 	$(CXX) -o $@ $(CXXFLAGS) $^

lib/colvars/Makefile.mingw32-cross

@@ -1,4 +1,4 @@
-# library build makefile for colvars module
+# library build -*- makefile -*- for colvars module
 
 # which file will be copied to Makefile.lammps
 
@@ -6,35 +6,42 @@ EXTRAMAKE = Makefile.lammps.empty
 
 # ------ SETTINGS ------
 
-CXX =		i686-pc-mingw32-g++
+CXX =		i686-w64-mingw32-g++
 CXXFLAGS =	-O2 -march=i686 -mtune=generic -mfpmath=387 -mpc64	\
 		-fno-rtti -fno-exceptions -finline-functions		\
                 -ffast-math -funroll-loops -fstrict-aliasing		\
                 -Wall -W -Wno-uninitialized
-ARCHIVE =	i686-pc-mingw32-ar
+ARCHIVE =	i686-w64-mingw32-ar
 ARCHFLAG =	-rscv
 SHELL =		/bin/sh
 
 # ------ DEFINITIONS ------
 
-SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp	\
+SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \
- colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp	\
+ colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp   \
- colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp	\
+ colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp   \
  colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp
 
-LIB = libcolvars.a
+DIR = Obj_mingw32/
-OBJ = $(SRC:.cpp=.o)
+LIB = $(DIR)libcolvars.a
+OBJ = $(SRC:%.cpp=$(DIR)%.o)
 EXE = #colvars_standalone
 
 # ------ MAKE PROCEDURE ------
 
-default: $(LIB) $(EXE)
+default: $(DIR) $(LIB) $(EXE) Makefile.lammps
 
-$(LIB):	$(OBJ)
+$(DIR):
+	mkdir $(DIR)
+
+Makefile.lammps:
+	@cp $(EXTRAMAKE) Makefile.lammps
+
+$(LIB): $(DIR)	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
 	@cp $(EXTRAMAKE) Makefile.lammps
 
-colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
+$(DIR)colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
 	$(CXX) -o $@ $(CXXFLAGS) $^
 
 # ------ MAKE FLAGS ------
@@ -46,58 +53,59 @@ colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB)
 
 # ------ COMPILE RULES ------
 
-.cpp.o:
+$(DIR)%.o: %.cpp
-	$(CXX) $(CXXFLAGS) -c $<
+	$(CXX) $(CXXFLAGS) -c $< -o $@
 
 # ------ DEPENDENCIES ------
 #
-colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarproxy_standalone.h colvaratoms.h colvarparse.h colvarvalue.h
-colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \
+$(DIR)colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \
  colvartypes.h colvarproxy.h colvaratoms.h colvarparse.h colvarvalue.h \
  colvarproxy_standalone.h
-colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarparse.h colvarvalue.h colvaratoms.h
-colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
+$(DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \
  colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \
  colvarbias.h colvargrid.h
-colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarbias.h colvar.h colvarparse.h
-colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
+$(DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h \
  colvarbias_meta.h colvarbias.h colvargrid.h
-colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h
-colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
+$(DIR)colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \
  colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarcomp.h \
  colvaratoms.h
-colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h
-colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
+$(DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \
  colvartypes.h colvarproxy.h colvarparse.h colvarvalue.h colvaratoms.h \
  colvar.h colvarcomp.h
-colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
+$(DIR)colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \
  colvarcomp.h colvaratoms.h
-colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
+$(DIR)colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \
  colvarproxy.h colvarvalue.h colvarparse.h colvar.h colvarcomp.h \
  colvaratoms.h
-colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
+$(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \
  colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \
  colvarcomp.h colvaratoms.h
-colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \
  colvargrid.h
-colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarparse.h colvarvalue.h colvar.h colvarbias.h colvarbias_meta.h \
  colvargrid.h colvarbias_abf.h
-colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h colvarparse.h
-colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
+$(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \
  colvarvalue.h
 
 # ------ CLEAN ------
 
 clean:
-	-rm *.o *~ $(LIB)
+	-rm $(DIR)*.o *~ $(LIB)
+	-rmdir $(DIR)
 

lib/colvars/README

@@ -4,9 +4,24 @@ that allows enhanced sampling in molecular dynamics simulations.
 The module is written to maximize performance, portability,
 flexibility of usage for the user, and extensibility for the developer.
 
-The following publication describes the principles of
+The development of the colvars library is now hosted on github at:
+http://colvars.github.io/
+You can use this site to get access to the latest development sources
+and the up-to-date documentation.
+
+Copy of the specific documentation is also in
+ doc/PDF/colvars-refman-lammps.pdf
+
+Please report bugs and request new features at:
+https://github.com/colvars/colvars/issues
+
+The following publications describe the principles of
 the implementation of this library:
 
+  Using collective variables to drive molecular dynamics simulations,
+  Giacomo Fiorin , Michael L. Klein & Jérôme Hénin (2013):
+  Molecular Physics DOI:10.1080/00268976.2013.813594
+
   Exploring Multidimensional Free Energy Landscapes Using 
   Time-Dependent Biases on Collective Variables, 
    J. Hénin, G. Fiorin, C. Chipot, and M. L. Klein, 

lib/colvars/colvar.cpp

@@ -14,7 +14,7 @@
 colvar::colvar (std::string const &conf)
 {
   cvm::log ("Initializing a new collective variable.\n");
-  
+
   get_keyval (conf, "name", this->name,
               (std::string ("colvar")+cvm::to_str (cvm::colvars.size()+1)));
 
@@ -26,7 +26,7 @@ colvar::colvar (std::string const &conf)
                         "\", as another colvar.\n");
   }
 
-  // all tasks disabled by default 
+  // all tasks disabled by default
   for (size_t i = 0; i < task_ntot; i++) {
     tasks[i] = false;
   }
@@ -152,7 +152,7 @@ colvar::colvar (std::string const &conf)
 
   // Decide whether the colvar is periodic
   // Used to wrap extended DOF if extendedLagrangian is on
-  if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1 
+  if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1
                                                 && (cvcs[0])->sup_coeff == 1.0 ) {
     this->b_periodic = true;
     this->period = (cvcs[0])->period;
@@ -167,6 +167,9 @@ colvar::colvar (std::string const &conf)
   // check the available features of each cvc
   for (size_t i = 0; i < cvcs.size(); i++) {
 
+    if ((cvcs[i])->b_debug_gradients)
+      enable (task_gradients);
+
     if ((cvcs[i])->sup_np != 1) {
       if (cvm::debug() && b_linear)
         cvm::log ("Warning: You are using a non-linear polynomial "
@@ -307,7 +310,7 @@ colvar::colvar (std::string const &conf)
 
       cvm::log ("Enabling the extended Lagrangian term for colvar \""+
                 this->name+"\".\n");
-    
+
       enable (task_extended_lagrangian);
 
       xr.type (this->type());
@@ -525,10 +528,10 @@ void colvar::enable (colvar::task const &t)
     if ( !tasks[task_extended_lagrangian] ) {
       enable (task_gradients);
 
-      if (!b_Jacobian_force) 
+      if (!b_Jacobian_force)
         cvm::fatal_error ("Error: colvar \""+this->name+
                           "\" does not have Jacobian forces implemented.\n");
-      if (!b_linear) 
+      if (!b_linear)
         cvm::fatal_error ("Error: colvar \""+this->name+
                           "\" must be defined as a linear combination "
                           "to calculate the Jacobian force.\n");
@@ -664,12 +667,24 @@ void colvar::disable (colvar::task const &t)
 }
 
 
+void colvar::setup() {
+  // loop over all components to reset masses of all groups
+  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.reset_mass(name,i,ig);
+    }
+  }
+}
+
+
 colvar::~colvar()
 {
   for (size_t i = 0; i < cvcs.size(); i++) {
     delete cvcs[i];
   }
-}  
+}
 
 
 
@@ -682,6 +697,8 @@ void colvar::calc()
     cvm::log ("Calculating colvar \""+this->name+"\".\n");
 
   // prepare atom groups for calculation
+  if (cvm::debug())
+    cvm::log ("Collecting data from atom groups.\n");
   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]);
@@ -697,19 +714,21 @@ void colvar::calc()
     for (size_t i = 0; i < cvcs.size(); i++) {
       for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
         cvcs[i]->atom_groups[ig]->read_velocities();
-      }   
+      }
     }
   }
   if (tasks[task_system_force]) {
     for (size_t i = 0; i < cvcs.size(); i++) {
       for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
         cvcs[i]->atom_groups[ig]->read_system_forces();
-      }   
+      }
     }
   }
 
   // calculate the value of the colvar
 
+  if (cvm::debug())
+    cvm::log ("Calculating colvar components.\n");
   x.reset();
   if (x.type() == colvarvalue::type_scalar) {
     // polynomial combination allowed
@@ -727,7 +746,7 @@ void colvar::calc()
         ( ((cvcs[i])->sup_np != 1) ?
           std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np) :
           (cvcs[i])->value().real_value );
-    } 
+    }
   } else {
     // only linear combination allowed
 
@@ -741,7 +760,7 @@ void colvar::calc()
                   cvm::to_str ((cvcs[i])->value(),
                                cvm::cv_width, cvm::cv_prec)+".\n");
       x += (cvcs[i])->sup_coeff * (cvcs[i])->value();
-    } 
+    }
   }
 
   if (cvm::debug())
@@ -762,9 +781,9 @@ void colvar::calc()
       // if requested, propagate (via chain rule) the gradients above
       // to the atoms used to define the roto-translation
       for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
-        if (cvcs[i]->atom_groups[ig]->b_fit_gradients) 
+        if (cvcs[i]->atom_groups[ig]->b_fit_gradients)
           cvcs[i]->atom_groups[ig]->calc_fit_gradients();
-      }   
+      }
 
       cvm::decrease_depth();
     }
@@ -801,7 +820,7 @@ void colvar::calc()
             for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) {
               int a = std::lower_bound (atom_ids.begin(), atom_ids.end(),
                   cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin();
-              atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad; 
+              atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad;
             }
           }
         }
@@ -925,7 +944,7 @@ cvm::real colvar::update()
   }
 
   if (tasks[task_Jacobian_force]) {
-    
+
     cvm::increase_depth();
     for (size_t i = 0; i < cvcs.size(); i++) {
       (cvcs[i])->calc_Jacobian_derivative();
@@ -941,8 +960,8 @@ cvm::real colvar::update()
     fj *= cvm::boltzmann() * cvm::temperature();
 
     // the instantaneous Jacobian force was not included in the reported system force;
-    // instead, it is subtracted from the applied force (silent Jacobian correction) 
+    // instead, it is subtracted from the applied force (silent Jacobian correction)
-    if (! tasks[task_report_Jacobian_force]) 
+    if (! tasks[task_report_Jacobian_force])
       f -= fj;
   }
 
@@ -966,9 +985,9 @@ cvm::real colvar::update()
     // leap to v_(i+1/2)
     if (tasks[task_langevin]) {
       vr -= dt * ext_gamma * vr.real_value;
-      vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass; 
+      vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass;
-    } 
+    }
-    vr  += (0.5 * dt) * fr / ext_mass; 
+    vr  += (0.5 * dt) * fr / ext_mass;
     xr  += dt * vr;
     xr.apply_constraints();
     if (this->b_periodic) this->wrap (xr);
@@ -989,13 +1008,13 @@ cvm::real colvar::update()
 void colvar::communicate_forces()
 {
   if (cvm::debug())
-    cvm::log ("Communicating forces from colvar \""+this->name+"\".\n"); 
+    cvm::log ("Communicating forces from colvar \""+this->name+"\".\n");
 
   if (x.type() == colvarvalue::type_scalar) {
 
     for (size_t i = 0; i < cvcs.size(); i++) {
       cvm::increase_depth();
-      (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff * 
+      (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff *
                               cvm::real ((cvcs[i])->sup_np) *
                               (std::pow ((cvcs[i])->value().real_value,
                                       (cvcs[i])->sup_np-1)) );
@@ -1012,7 +1031,7 @@ void colvar::communicate_forces()
   }
 
   if (cvm::debug())
-    cvm::log ("Done communicating forces from colvar \""+this->name+"\".\n"); 
+    cvm::log ("Done communicating forces from colvar \""+this->name+"\".\n");
 }
 
 
@@ -1200,7 +1219,7 @@ std::istream & colvar::read_traj (std::istream &is)
 
     is >> ft;
 
-    if (tasks[task_extended_lagrangian]) { 
+    if (tasks[task_extended_lagrangian]) {
       is >> fr;
       ft_reported = fr;
     } else {
@@ -1248,7 +1267,7 @@ std::ostream & colvar::write_restart (std::ostream &os) {
   os << "}\n\n";
 
   return os;
-}  
+}
 
 
 std::ostream & colvar::write_traj_label (std::ostream & os)
@@ -1418,7 +1437,7 @@ inline void history_add_value (size_t const           &history_length,
 inline void history_incr (std::list< std::list<colvarvalue> >           &history,
                           std::list< std::list<colvarvalue> >::iterator &history_p)
 {
-  if ((++history_p) == history.end()) 
+  if ((++history_p) == history.end())
     history_p = history.begin();
 }
 
@@ -1435,7 +1454,7 @@ void colvar::calc_acf()
 
     // first-step operations
 
-    colvar *cfcv = (acf_colvar_name.size() ? 
+    colvar *cfcv = (acf_colvar_name.size() ?
                     cvm::colvar_p (acf_colvar_name) :
                     this);
     if (cfcv->type() != this->type())
@@ -1474,7 +1493,7 @@ void colvar::calc_acf()
 
   } else {
 
-    colvar *cfcv = (acf_colvar_name.size() ? 
+    colvar *cfcv = (acf_colvar_name.size() ?
                     cvm::colvar_p (acf_colvar_name) :
                     this);
 
@@ -1539,7 +1558,7 @@ void colvar::calc_vel_acf (std::list<colvarvalue> &v_list,
 
     // inner products of previous velocities with current (acf_i and
     // vs_i are updated)
-    colvarvalue::inner_opt (v, vs_i, v_list.end(), acf_i); 
+    colvarvalue::inner_opt (v, vs_i, v_list.end(), acf_i);
 
     acf_nframes++;
   }
@@ -1559,7 +1578,7 @@ void colvar::calc_coor_acf (std::list<colvarvalue> &x_list,
 
     *(acf_i++) += x.norm2();
 
-    colvarvalue::inner_opt (x, xs_i, x_list.end(), acf_i); 
+    colvarvalue::inner_opt (x, xs_i, x_list.end(), acf_i);
 
     acf_nframes++;
   }
@@ -1581,7 +1600,7 @@ void colvar::calc_p2coor_acf (std::list<colvarvalue> &x_list,
     // value of P2(0) = 1
     *(acf_i++) += 1.0;
 
-    colvarvalue::p2leg_opt (x, xs_i, x_list.end(), acf_i); 
+    colvarvalue::p2leg_opt (x, xs_i, x_list.end(), acf_i);
 
     acf_nframes++;
   }

lib/colvars/colvar.h

@@ -215,7 +215,7 @@ protected:
   cvm::real ext_gamma;
   /// Amplitude of Gaussian white noise for Langevin extended dynamics
   cvm::real ext_sigma;
-  
+
   /// \brief Harmonic restraint force
   colvarvalue fr;
 
@@ -288,6 +288,9 @@ public:
   /// Disable the specified task
   void disable (colvar::task const &t);
 
+  /// Get ready for a run and possibly re-initialize internal data
+  void setup();
+
   /// Destructor
   ~colvar();
 
@@ -380,7 +383,7 @@ protected:
   colvarvalue            x_old;
 
   /// Time series of values and velocities used in correlation
-  /// functions 
+  /// functions
   std::list< std::list<colvarvalue> > acf_x_history, acf_v_history;
   /// Time series of values and velocities used in correlation
   /// functions (pointers)x

lib/colvars/colvaratoms.cpp

@@ -26,8 +26,8 @@ cvm::atom_group::atom_group (std::string const &conf,
 
 
 cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
-  : b_dummy (false), b_center (false), b_rotate (false), 
+  : b_dummy (false), b_center (false), b_rotate (false),
-    ref_pos_group (NULL), b_fit_gradients (false), 
+    ref_pos_group (NULL), b_fit_gradients (false),
     noforce (false)
 {
   this->reserve (atoms.size());
@@ -43,8 +43,8 @@ cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
 
 
 cvm::atom_group::atom_group()
-  : b_dummy (false), b_center (false), b_rotate (false), 
+  : b_dummy (false), b_center (false), b_rotate (false),
-    ref_pos_group (NULL), b_fit_gradients (false), 
+    ref_pos_group (NULL), b_fit_gradients (false),
     noforce (false)
 {
   total_mass = 0.0;
@@ -71,6 +71,18 @@ void cvm::atom_group::add_atom (cvm::atom const &a)
 }
 
 
+void cvm::atom_group::reset_mass(std::string &name, int i, int j)
+{
+  total_mass = 0.0;
+  for (cvm::atom_iter ai = this->begin();
+       ai != this->end(); ai++) {
+    total_mass += ai->mass;
+  }
+  cvm::log ("Re-initialized atom group "+name+":"+cvm::to_str (i)+"/"+
+            cvm::to_str (j)+". "+ cvm::to_str (this->size())+
+            " atoms: total mass = "+cvm::to_str (this->total_mass)+".\n");
+}
+
 void cvm::atom_group::parse (std::string const &conf,
                              char const        *key)
 {
@@ -129,6 +141,25 @@ void cvm::atom_group::parse (std::string const &conf,
 
       atom_indexes.clear();
     }
+
+    std::string index_group_name;
+    if (get_keyval (group_conf, "indexGroup", index_group_name)) {
+      // use an index group from the index file read globally
+      std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
+      std::list<std::vector<int> >::iterator index_groups_i = cvm::index_groups.begin();
+      for ( ; names_i != cvm::index_group_names.end() ; names_i++, index_groups_i++) {
+        if (*names_i == index_group_name)
+          break;
+      }
+      if (names_i == cvm::index_group_names.end()) {
+        cvm::fatal_error ("Error: could not find index group "+
+                          index_group_name+" among those provided by the index file.\n");
+      }
+      this->reserve (index_groups_i->size());
+      for (size_t i = 0; i < index_groups_i->size(); i++) {
+        this->push_back (cvm::atom ((*index_groups_i)[i]));
+      }
+    }
   }
 
   {
@@ -234,7 +265,7 @@ void cvm::atom_group::parse (std::string const &conf,
     }
   }
 
-  for (std::vector<cvm::atom>::iterator a1 = this->begin(); 
+  for (std::vector<cvm::atom>::iterator a1 = this->begin();
        a1 != this->end(); a1++) {
     std::vector<cvm::atom>::iterator a2 = a1;
     ++a2;
@@ -253,10 +284,10 @@ void cvm::atom_group::parse (std::string const &conf,
   if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
     b_dummy = true;
     this->total_mass = 1.0;
-  } else 
+  } else
     b_dummy = false;
 
-  if (b_dummy && (this->size())) 
+  if (b_dummy && (this->size()))
     cvm::fatal_error ("Error: cannot set up group \""+
                       std::string (key)+"\" as a dummy atom "
                       "and provide it with atom definitions.\n");
@@ -328,7 +359,7 @@ void cvm::atom_group::parse (std::string const &conf,
 
       std::string ref_pos_col;
       double ref_pos_col_value;
-      
+
       if (get_keyval (group_conf, "refPositionsCol", ref_pos_col, std::string (""), mode)) {
         // if provided, use PDB column to select coordinates
         bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode);
@@ -486,7 +517,7 @@ void cvm::atom_group::calc_apply_roto_translation()
   }
 }
 
-void cvm::atom_group::apply_translation (cvm::rvector const &t) 
+void cvm::atom_group::apply_translation (cvm::rvector const &t)
 {
   if (b_dummy) return;
 
@@ -496,7 +527,7 @@ void cvm::atom_group::apply_translation (cvm::rvector const &t)
   }
 }
 
-void cvm::atom_group::apply_rotation (cvm::rotation const &rot) 
+void cvm::atom_group::apply_rotation (cvm::rotation const &rot)
 {
   if (b_dummy) return;
 
@@ -617,11 +648,11 @@ void cvm::atom_group::calc_fit_gradients()
   }
 
   if (b_rotate) {
-    
+
     // add the rotation matrix contribution to the gradients
     cvm::rotation const rot_inv = rot.inverse();
     cvm::atom_pos const cog = this->center_of_geometry();
-    
+
     for (size_t i = 0; i < this->size(); i++) {
 
       cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));
@@ -785,7 +816,7 @@ void cvm::atom_group::apply_force (cvm::rvector const &force)
     for (cvm::atom_iter ai = this->begin();
          ai != this->end(); ai++) {
       ai->apply_force ((ai->mass/this->total_mass) * force);
-    } 
+    }
   }
 }
 

lib/colvars/colvaratoms.h

@@ -8,7 +8,7 @@
 
 /// \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
 /// data (id, mass, position and collective variable derivatives)
 /// altogether.  There may be multiple instances with identical
@@ -23,7 +23,7 @@ protected:
 
   /// \brief Index in the list of atoms involved by the colvars (\b
   /// NOT in the global topology!)
-  size_t        index;
+  int           index;
 
 public:
 
@@ -47,7 +47,7 @@ public:
 
   /// \brief Gradient of a scalar collective variable with respect
   /// to this atom
-  /// 
+  ///
   /// This can only handle a scalar collective variable (i.e. when
   /// the \link colvarvalue::real_value \endlink member is used
   /// from the \link colvarvalue \endlink class), which is also the
@@ -57,8 +57,8 @@ public:
   /// implementation
   cvm::rvector   grad;
 
-  /// \brief Default constructor, setting id to a non-valid value
+  /// \brief Default constructor, setting id and index to invalid numbers
-  inline atom() {}
+  atom() : id (-1), index (-1) { reset_data(); }
 
   /// \brief Initialize an atom for collective variable calculation
   /// and get its internal identifier \param atom_number Atom index in
@@ -121,7 +121,7 @@ class colvarmodule::atom_group
     public colvarparse
 {
 public:
-  // Note: all members here are kept public, to make possible to any
+  // Note: all members here are kept public, to allow any
   // object accessing and manipulating them
 
 
@@ -138,7 +138,7 @@ public:
 
   /// Allocates and populates the sorted list of atom ids
   void create_sorted_ids (void);
-  
+
 
   /// \brief When updating atomic coordinates, translate them to align with the
   /// center of mass of the reference coordinates
@@ -175,7 +175,7 @@ public:
   /// \brief If b_center or b_rotate is true, use this group to
   /// define the transformation (default: this group itself)
   atom_group                *ref_pos_group;
-  
+
   /// Total mass of the atom group
   cvm::real total_mass;
 
@@ -199,9 +199,14 @@ public:
   /// \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  
+  /// \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();
 
@@ -216,7 +221,7 @@ public:
   void calc_apply_roto_translation();
 
   /// \brief Save center of geometry fo ref positions,
-  /// then subtract it 
+  /// then subtract it
   void center_ref_pos();
 
   /// \brief Move all positions

lib/colvars/colvarbias.cpp

@@ -41,10 +41,11 @@ colvarbias::colvarbias (std::string const &conf, char const *key)
       add_colvar (colvars_str[i]);
     }
   }
-
   if (!colvars.size()) {
     cvm::fatal_error ("Error: no collective variables specified.\n");
   }
+
+  get_keyval (conf, "outputEnergy", b_output_energy, false);
 }
 
 
@@ -59,7 +60,7 @@ void colvarbias::add_colvar (std::string const &cv_name)
     cvp->enable (colvar::task_gradients);
     if (cvm::debug())
       cvm::log ("Applying this bias to collective variable \""+
-                cvp->name+"\".\n"); 
+                cvp->name+"\".\n");
     colvars.push_back (cvp);
     colvar_forces.push_back (colvarvalue (cvp->type()));
   } else {
@@ -79,13 +80,47 @@ void colvarbias::communicate_forces()
     }
     colvars[i]->add_bias_force (colvar_forces[i]);
   }
-}    
+}
+
+
+void colvarbias::change_configuration(std::string const &conf)
+{
+  cvm::fatal_error ("Error: change_configuration() not implemented.\n");
+}
+
+
+cvm::real colvarbias::energy_difference(std::string const &conf)
+{
+  cvm::fatal_error ("Error: energy_difference() not implemented.\n");
+  return 0.;
+}
+
+
+std::ostream & colvarbias::write_traj_label (std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy)
+    os << " E_"
+       << cvm::wrap_string (this->name, cvm::en_width-2);
+  return os;
+}
+
+
+std::ostream & colvarbias::write_traj (std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy)
+    os << " "
+       << bias_energy;
+  return os;
+}
+
 
 
 
 colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
                                           char const *key)
-  : colvarbias (conf, key), 
+  : colvarbias (conf, key),
     target_nsteps (0),
     target_nstages (0)
 {
@@ -171,25 +206,22 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf,
     }
   }
 
+  get_keyval (conf, "outputCenters", b_output_centers, false);
+  get_keyval (conf, "outputAccumulatedWork", b_output_acc_work, false);
+  acc_work = 0.0;
+
   if (cvm::debug())
     cvm::log ("Done initializing a new harmonic restraint bias.\n");
 }
 
-
+colvarbias_harmonic::~colvarbias_harmonic ()
-void colvarbias::change_configuration(std::string const &conf)
 {
-  cvm::fatal_error ("Error: change_configuration() not implemented.\n");
+  if (cvm::n_harm_biases > 0)
+    cvm::n_harm_biases -= 1;
 }
 
 
-cvm::real colvarbias::energy_difference(std::string const &conf)
+void colvarbias_harmonic::change_configuration (std::string const &conf)
-{
-  cvm::fatal_error ("Error: energy_difference() not implemented.\n");
-  return 0.;
-}
-
-
-void colvarbias_harmonic::change_configuration(std::string const &conf)
 {
   get_keyval (conf, "forceConstant", force_k, force_k);
   if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) {
@@ -317,7 +349,7 @@ cvm::real colvarbias_harmonic::update()
 
       if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) {
         // Start averaging after equilibration period, if requested
-        
+
         // Square distance normalized by square colvar width
         cvm::real dist_sq = 0.0;
         for (size_t i = 0; i < colvars.size(); i++) {
@@ -335,7 +367,7 @@ cvm::real colvarbias_harmonic::update()
 
           cvm::log ("Lambda= " + cvm::to_str (lambda) + " dA/dLambda= "
               + cvm::to_str (restraint_FE / cvm::real(target_nsteps - target_equil_steps)));
-      
+
         //  ...and move on to the next one
         if (stage < target_nstages) {
 
@@ -363,7 +395,7 @@ 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] =
@@ -379,6 +411,15 @@ cvm::real colvarbias_harmonic::update()
                                colvar_centers[i]))+"\n");
   }
 
+  if (b_output_acc_work) {
+    if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) {
+      for (size_t i = 0; i < colvars.size(); i++) {
+        // project forces on the calculated increments at this step
+        acc_work += colvar_forces[i] * centers_incr[i];
+      }
+    }
+  }
+
   if (cvm::debug())
     cvm::log ("Current forces for the harmonic bias \""+
               this->name+"\": "+cvm::to_str (colvar_forces)+".\n");
@@ -408,7 +449,7 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is)
     return is;
   }
 
-//   int id = -1; 
+//   int id = -1;
   std::string name = "";
 //   if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) &&
 //          (id != this->id) ) ||
@@ -442,6 +483,11 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is)
       cvm::fatal_error ("Error: current stage is missing from the restart.\n");
   }
 
+  if (b_output_acc_work) {
+    if (!get_keyval (conf, "accumulatedWork", acc_work))
+      cvm::fatal_error ("Error: accumulatedWork is missing from the restart.\n");
+  }
+
   is >> brace;
   if (brace != "}") {
     cvm::fatal_error ("Error: corrupt restart information for harmonic bias \""+
@@ -484,9 +530,61 @@ std::ostream & colvarbias_harmonic::write_restart (std::ostream &os)
        << stage << "\n";
   }
 
+  if (b_output_acc_work) {
+    os << "    accumulatedWork " << acc_work << "\n";
+  }
+
   os << "  }\n"
      << "}\n\n";
 
   return os;
 }
 
+
+std::ostream & colvarbias_harmonic::write_traj_label (std::ostream &os)
+{
+  os << " ";
+
+  if (b_output_energy)
+    os << " E_"
+       << cvm::wrap_string (this->name, cvm::en_width-2);
+
+  if (b_output_centers)
+    for (size_t i = 0; i < colvars.size(); i++) {
+      size_t const this_cv_width = (colvars[i]->value()).output_width (cvm::cv_width);
+      os << " x0_"
+         << cvm::wrap_string (colvars[i]->name, this_cv_width-3);
+    }
+
+  if (b_output_acc_work)
+    os << " W_"
+       << cvm::wrap_string (this->name, cvm::en_width-2);
+
+  return os;
+}
+
+
+std::ostream & colvarbias_harmonic::write_traj (std::ostream &os)
+{
+  os << " ";
+
+  if (b_output_energy)
+    os << " "
+       << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width)
+       << bias_energy;
+
+  if (b_output_centers)
+    for (size_t i = 0; i < colvars.size(); i++) {
+      os << " "
+         << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width)
+         << colvar_centers[i];
+    }
+
+  if (b_output_acc_work)
+    os << " "
+       << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width)
+       << acc_work;
+
+  return os;
+}
+

lib/colvars/colvarbias.h

@@ -14,7 +14,7 @@ public:
 
   /// Name of this bias
   std::string    name;
-  
+
   /// Add a new collective variable to this bias
   void add_colvar (std::string const &cv_name);
 
@@ -35,7 +35,7 @@ public:
   void communicate_forces();
 
   /// \brief Constructor
-  /// 
+  ///
   /// The constructor of the colvarbias base class is protected, so
   /// that it can only be called from inherited classes
   colvarbias (std::string const &conf, char const *key);
@@ -52,6 +52,13 @@ public:
   /// Write the bias configuration to a restart file
   virtual std::ostream & write_restart (std::ostream &os) = 0;
 
+  /// Write a label to the trajectory file (comment line)
+  virtual std::ostream & write_traj_label (std::ostream &os);
+
+  /// Output quantities such as the bias energy to the trajectory file
+  virtual std::ostream & write_traj (std::ostream &os);
+
+
 protected:
 
   /// \brief Pointers to collective variables to which the bias is
@@ -62,10 +69,12 @@ protected:
   /// \brief Current forces from this bias to the colvars
   std::vector<colvarvalue> colvar_forces;
 
-  /// \brief Current energy of this bias (colvar_forces should be
+  /// \brief Current energy of this bias (colvar_forces should be obtained by deriving this)
-  /// obtained by deriving this)
   cvm::real                bias_energy;
 
+  /// Whether to write the current bias energy from this bias to the trajectory file
+  bool                     b_output_energy;
+
   /// \brief Whether this bias has already accumulated information
   /// (when relevant)
   bool                     has_data;
@@ -94,11 +103,17 @@ public:
   /// Write the bias configuration to a restart file
   virtual std::ostream & write_restart (std::ostream &os);
 
+  /// Write a label to the trajectory file (comment line)
+  virtual std::ostream & write_traj_label (std::ostream &os);
+
+  /// Output quantities such as the bias energy to the trajectory file
+  virtual std::ostream & write_traj (std::ostream &os);
+
   /// \brief Constructor
   colvarbias_harmonic (std::string const &conf, char const *key);
 
   /// Destructor
-  virtual inline ~colvarbias_harmonic() {}
+  virtual ~colvarbias_harmonic();
 
 
 protected:
@@ -109,27 +124,9 @@ protected:
   /// \brief Restraint centers without wrapping or constraints applied
   std::vector<colvarvalue> colvar_centers_raw;
 
-  /// \brief Restraint force constant
-  cvm::real force_k;
-
   /// \brief Moving target?
   bool b_chg_centers;
 
-  /// \brief Changing force constant?
-  bool b_chg_force_k;
-
-  /// \brief Restraint force constant (target value)
-  cvm::real target_force_k;
-
-  /// \brief Equilibration steps for restraint FE calculation through TI
-  cvm::real target_equil_steps;
-
-  /// \brief Restraint force constant (starting value)
-  cvm::real starting_force_k;
-
-  /// \brief Lambda-schedule for custom varying force constant
-  std::vector<cvm::real> lambda_schedule;
-
   /// \brief New restraint centers
   std::vector<colvarvalue> target_centers;
 
@@ -137,12 +134,41 @@ protected:
   /// (or stage) towards the new values (calculated from target_nsteps)
   std::vector<colvarvalue> centers_incr;
 
+  /// Whether to write the current restraint centers to the trajectory file
+  bool b_output_centers;
+
+  /// Whether to write the current accumulated work to the trajectory file
+  bool b_output_acc_work;
+
+  /// \brief Accumulated work
+  cvm::real acc_work;
+
+
+  /// \brief Restraint force constant
+  cvm::real force_k;
+
+  /// \brief Changing force constant?
+  bool b_chg_force_k;
+
+  /// \brief Restraint force constant (target value)
+  cvm::real target_force_k;
+
+  /// \brief Restraint force constant (starting value)
+  cvm::real starting_force_k;
+
+  /// \brief Lambda-schedule for custom varying force constant
+  std::vector<cvm::real> lambda_schedule;
+
   /// \brief Exponent for varying the force constant
   cvm::real force_k_exp;
 
-  /// \brief Number of steps required to reach the target force constant
+  /// \brief Intermediate quantity to compute the restraint free energy
-  /// or restraint centers
+  /// (in TI, would be the accumulating FE derivative)
-  size_t target_nsteps;
+  cvm::real restraint_FE;
+
+
+  /// \brief Equilibration steps for restraint FE calculation through TI
+  cvm::real target_equil_steps;
 
   /// \brief Number of stages over which to perform the change
   /// If zero, perform a continuous change
@@ -150,10 +176,10 @@ protected:
 
   /// \brief Number of current stage of the perturbation
   int stage;
-  
+
-  /// \brief Intermediate quantity to compute the restraint free energy
+  /// \brief Number of steps required to reach the target force constant
-  /// (in TI, would be the accumulating FE derivative)
+  /// or restraint centers
-  cvm::real restraint_FE;
+  size_t target_nsteps;
 };
 
 

lib/colvars/colvarbias_abf.cpp

@@ -15,7 +15,7 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
 {
   if (cvm::temperature() == 0.0)
     cvm::log ("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n");
-  
+
   // ************* parsing general ABF options ***********************
 
   get_keyval (conf, "applyBias",  apply_bias, true);
@@ -32,7 +32,7 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
   }
 
   get_keyval (conf, "fullSamples", full_samples, 200);
-  if ( full_samples <= 1 ) full_samples = 1; 
+  if ( full_samples <= 1 ) full_samples = 1;
   min_samples = full_samples / 2;
   // full_samples - min_samples >= 1 is guaranteed
 
@@ -79,7 +79,21 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
     }
 
     // Here we could check for orthogonality of the Cartesian coordinates
-    // and make it just a warning if some parameter is set? 
+    // and make it just a warning if some parameter is set?
+  }
+
+  if (get_keyval (conf, "maxForce", max_force)) {
+    if (max_force.size() != colvars.size()) {
+      cvm::fatal_error ("Error: Number of parameters to maxForce does not match number of colvars.");
+    }
+    for (size_t i=0; i<colvars.size(); i++) {
+      if (max_force[i] < 0.0) {
+        cvm::fatal_error ("Error: maxForce should be non-negative.");
+      }
+    }
+    cap_force = true;
+  } else {
+    cap_force = false;
   }
 
   bin.assign (colvars.size(), 0);
@@ -96,7 +110,7 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
   if ( input_prefix.size() > 0 ) {
     read_gradients_samples ();
   }
-  
+
   cvm::log ("Finished ABF setup.\n");
 }
 
@@ -114,6 +128,9 @@ colvarbias_abf::~colvarbias_abf()
   }
 
   delete [] force;
+
+  if (cvm::n_abf_biases > 0)
+    cvm::n_abf_biases -= 1;
 }
 
 
@@ -125,7 +142,7 @@ cvm::real colvarbias_abf::update()
   if (cvm::debug()) cvm::log ("Updating ABF bias " + this->name);
 
   if (cvm::step_relative() == 0) {
-	
+
     // At first timestep, do only:
     // initialization stuff (file operations relying on n_abf_biases
     // compute current value of colvars
@@ -158,7 +175,7 @@ cvm::real colvarbias_abf::update()
   }
 
   // save bin for next timestep
-  force_bin = bin;  
+  force_bin = bin;
 
   // Reset biasing forces from previous timestep
   for (size_t i=0; i<colvars.size(); i++) {
@@ -176,19 +193,25 @@ cvm::real colvarbias_abf::update()
       fact = ( count < min_samples) ? 0.0 :
         (cvm::real (count - min_samples)) / (cvm::real (full_samples - min_samples));
     }
-	
+
     const cvm::real * grad  = &(gradients->value (bin));
 
     if ( fact != 0.0 ) {
-
       if ( (colvars.size() == 1) && colvars[0]->periodic_boundaries() ) {
         // Enforce a zero-mean bias on periodic, 1D coordinates
-        colvar_forces[0].real_value += fact * (grad[0] / cvm::real (count) - gradients->average ());
+        // in other words: boundary condition is that the biasing potential is periodic
+        colvar_forces[0].real_value = fact * (grad[0] / cvm::real (count) - gradients->average ());
       } else {
         for (size_t i=0; i<colvars.size(); i++) {
           // subtracting the mean force (opposite of the FE gradient) means adding the gradient
-          colvar_forces[i].real_value += fact * grad[i] / cvm::real (count);
+          colvar_forces[i].real_value = fact * grad[i] / cvm::real (count);
-          // without .real_value, the above would do (cheap) runtime type checking
+        }
+      }
+      if (cap_force) {
+        for (size_t i=0; i<colvars.size(); i++) {
+          if ( colvar_forces[i].real_value * colvar_forces[i].real_value > max_force[i] * max_force[i] ) {
+            colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]);
+          }
         }
       }
     }
@@ -203,7 +226,7 @@ cvm::real colvarbias_abf::update()
     // otherwise, backup and replace
     write_gradients_samples (output_prefix + ".hist", (cvm::step_absolute() > 0));
   }
-  return 0.0; // TODO compute bias energy whenever possible (i.e. 1D with updateBias off)
+  return 0.0;
 }
 
 
@@ -250,7 +273,7 @@ void colvarbias_abf::read_gradients_samples ()
     samples_in_name = input_prefix[i] + ".count";
     gradients_in_name = input_prefix[i] + ".grad";
     // For user-provided files, the per-bias naming scheme may not apply
-    
+
     std::ifstream is;
 
     cvm::log ("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);
@@ -272,7 +295,7 @@ void colvarbias_abf::read_gradients_samples ()
 std::ostream & colvarbias_abf::write_restart (std::ostream& os)
 {
 
-  std::ios::fmtflags flags (os.flags ()); 
+  std::ios::fmtflags flags (os.flags ());
   os.setf(std::ios::fmtflags (0), std::ios::floatfield); // default floating-point format
 
   os << "abf {\n"
@@ -325,7 +348,7 @@ std::istream & colvarbias_abf::read_restart (std::istream& is)
   if ( name == "" ) {
     cvm::fatal_error ("Error: \"abf\" block in the restart file has no name.\n");
   }
-  
+
   if ( !(is >> key)   || !(key == "samples")) {
     cvm::log ("Error: in reading restart configuration for ABF bias \""+
               this->name+"\" at position "+
@@ -381,7 +404,7 @@ colvarbias_histogram::colvarbias_histogram (std::string const &conf, char const
   bin.assign (colvars.size(), 0);
 
   out_name = cvm::output_prefix + "." + this->name + ".dat";
-  cvm::log ("Histogram will be written to file " + out_name); 
+  cvm::log ("Histogram will be written to file " + out_name);
 
   cvm::log ("Finished histogram setup.\n");
 }
@@ -395,6 +418,9 @@ colvarbias_histogram::~colvarbias_histogram()
     delete grid;
     grid = NULL;
   }
+
+  if (cvm::n_histo_biases > 0)
+    cvm::n_histo_biases -= 1;
 }
 
 /// Update the grid
@@ -452,7 +478,7 @@ std::istream & colvarbias_histogram::read_restart (std::istream& is)
     cvm::fatal_error ("Error: \"histogram\" block in the restart file "
                       "has no name.\n");
   }
-  
+
   if ( !(is >> key)   || !(key == "grid")) {
     cvm::log ("Error: in reading restart configuration for histogram \""+
               this->name+"\" at position "+

lib/colvars/colvarbias_abf.h

@@ -46,6 +46,10 @@ private:
   bool      b_history_files;
   size_t    history_freq;
 
+  /// Cap applied biasing force?
+  bool                    cap_force;
+  std::vector<cvm::real>  max_force;
+
   // Internal data and methods
 
   std::vector<int>  bin, force_bin;
@@ -84,7 +88,7 @@ private:
   std::vector<int>  bin;
   std::string	  out_name;
 
-  int		  output_freq;  
+  int		  output_freq;
   void		  write_grid ();
   std::ofstream	  grid_os;  /// Stream for writing grid to disk
 

lib/colvars/colvarbias_meta.cpp

@@ -54,7 +54,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
     get_keyval (conf, "multipleReplicas", b_replicas, false);
     if (b_replicas)
       comm = multiple_replicas;
-    else 
+    else
       comm = single_replica;
   }
 
@@ -76,7 +76,8 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
     }
 
     get_keyval (conf, "keepHills", keep_hills, false);
-    get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true);
+    if (! get_keyval (conf, "writeFreeEnergyFile", dump_fes, true))
+      get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
     get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false);
 
     for (size_t i = 0; i < colvars.size(); i++) {
@@ -113,7 +114,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
                         "when using more than one replica.\n");
 
     get_keyval (conf, "replicasRegistry",
-                replicas_registry_file, 
+                replicas_registry_file,
                 (this->name+".replicas.registry.txt"));
 
     get_keyval (conf, "replicaUpdateFrequency",
@@ -171,7 +172,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
     // if this replica was not included yet, we should generate a
     // new record for it: but first, we write this replica's files,
     // for the others to read
-      
+
     // open the "hills" buffer file
     replica_hills_os.open (replica_hills_file.c_str());
     if (!replica_hills_os.good())
@@ -188,7 +189,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
 
     // if we're running without grids, use a growing list of "hills" files
     // otherwise, just one state file and one "hills" file as buffer
-    std::ofstream list_os (replica_list_file.c_str(), 
+    std::ofstream list_os (replica_list_file.c_str(),
                            (use_grids ? std::ios::trunc : std::ios::app));
     if (! list_os.good())
       cvm::fatal_error ("Error: in opening file \""+
@@ -232,7 +233,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
     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");
@@ -258,6 +259,9 @@ colvarbias_meta::~colvarbias_meta()
 
   if (hills_traj_os.good())
     hills_traj_os.close();
+
+  if (cvm::n_meta_biases > 0)
+    cvm::n_meta_biases -= 1;
 }
 
 
@@ -266,7 +270,7 @@ colvarbias_meta::~colvarbias_meta()
 // Hill management member functions
 // **********************************************************************
 
-std::list<colvarbias_meta::hill>::const_iterator 
+std::list<colvarbias_meta::hill>::const_iterator
 colvarbias_meta::create_hill (colvarbias_meta::hill const &h)
 {
   hill_iter const hills_end = hills.end();
@@ -324,7 +328,7 @@ colvarbias_meta::delete_hill (hill_iter &h)
   }
 
   if (hills_traj_os.good()) {
-    // output to the trajectory 
+    // output to the trajectory
     hills_traj_os << "# DELETED this hill: "
                   << (hills.back()).output_traj()
                   << "\n";
@@ -459,9 +463,9 @@ cvm::real colvarbias_meta::update()
         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 { 
+      } else {
         create_hill (hill (hill_weight, colvars, hill_width));
-      } 
+      }
       break;
 
     case multiple_replicas:
@@ -470,9 +474,9 @@ cvm::real colvarbias_meta::update()
         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 { 
+      } else {
         create_hill (hill (hill_weight, colvars, hill_width, replica_id));
-      } 
+      }
       if (replica_hills_os.good()) {
         replica_hills_os << hills.back();
       } else {
@@ -591,11 +595,11 @@ cvm::real colvarbias_meta::update()
     calc_hills_force (ic, new_hills_begin, hills.end(), colvar_forces);
   }
 
-  if (cvm::debug()) 
+  if (cvm::debug())
     cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
               ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
 
-  if (cvm::debug()) 
+  if (cvm::debug())
     cvm::log ("Metadynamics bias \""+this->name+"\""+
               ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
               ": adding the forces from the other replicas.\n");
@@ -611,7 +615,7 @@ cvm::real colvarbias_meta::update()
                           replicas[ir]->hills.end(),
                           colvar_forces);
       }
-      if (cvm::debug()) 
+      if (cvm::debug())
         cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
                   ", hills forces = "+cvm::to_str (colvar_forces)+".\n");
     }
@@ -641,7 +645,7 @@ void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter      h_first,
   }
 
   for (hill_iter h = h_first; h != h_last; h++) {
-      
+
     // compute the gaussian exponent
     cvm::real cv_sqdev = 0.0;
     for (size_t i = 0; i < colvars.size(); i++) {
@@ -684,8 +688,8 @@ void colvarbias_meta::calc_hills_force (size_t const &i,
       if (h->value() == 0.0) continue;
       colvarvalue const &center = h->centers[i];
       cvm::real const    half_width = 0.5 * h->widths[i];
-      forces[i].real_value += 
+      forces[i].real_value +=
-        ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * 
+        ( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
           (colvars[i]->dist2_lgrad (x, center)).real_value );
     }
     break;
@@ -993,7 +997,7 @@ void colvarbias_meta::read_replica_files()
     }
 
     // (re)read the state file if necessary
-    if ( (! (replicas[ir])->has_data) || 
+    if ( (! (replicas[ir])->has_data) ||
          (! (replicas[ir])->replica_state_file_in_sync) ) {
 
       cvm::log ("Metadynamics bias \""+this->name+"\""+
@@ -1013,11 +1017,11 @@ void colvarbias_meta::read_replica_files()
       }
       is.close();
     }
-    
+
     // now read the hills added after writing the state file
     if ((replicas[ir])->replica_hills_file.size()) {
 
-      if (cvm::debug()) 
+      if (cvm::debug())
         cvm::log ("Metadynamics bias \""+this->name+"\""+
                   ": checking for new hills from replica \""+
                   (replicas[ir])->replica_id+"\" in the file \""+
@@ -1064,7 +1068,7 @@ void colvarbias_meta::read_replica_files()
                       "\" at position "+
                       cvm::to_str ((replicas[ir])->replica_hills_file_pos)+".\n");
 
-          // test whether this is the end of the file 
+          // test whether this is the end of the file
           is.seekg (0, std::ios::end);
           if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) {
             (replicas[ir])->update_status++;
@@ -1108,7 +1112,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
   size_t const start_pos = is.tellg();
 
   if (comm == single_replica) {
-    // if using a multiple replicas scheme, output messages 
+    // if using a multiple replicas scheme, output messages
     // are printed before and after calling this function
     cvm::log ("Restarting metadynamics bias \""+this->name+"\""+
               ".\n");
@@ -1118,7 +1122,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
        !(is >> brace) || !(brace == "{") ||
        !(is >> colvarparse::read_block ("configuration", conf)) ) {
 
-    if (comm == single_replica) 
+    if (comm == single_replica)
       cvm::log ("Error: in reading restart configuration for metadynamics bias \""+
                 this->name+"\""+
                 ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
@@ -1239,7 +1243,7 @@ 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 (!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+"\""+
@@ -1283,7 +1287,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
 
   // read the hills explicitly written (if there are any)
   while (read_hill (is)) {
-    if (cvm::debug()) 
+    if (cvm::debug())
       cvm::log ("Read a previously saved hill under the "
                 "metadynamics bias \""+
                 this->name+"\", created at step "+
@@ -1364,7 +1368,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
     hills.erase (hills.begin(), old_hills_end);
     hills_off_grid.erase (hills_off_grid.begin(), old_hills_off_grid_end);
   }
-  
+
   has_data = true;
 
   if (comm != single_replica) {
@@ -1372,7 +1376,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is)
   }
 
   return is;
-}  
+}
 
 
 std::istream & colvarbias_meta::read_hill (std::istream &is)
@@ -1404,7 +1408,7 @@ std::istream & colvarbias_meta::read_hill (std::istream &is)
 
   std::vector<colvarvalue> h_centers (colvars.size());
   for (size_t i = 0; i < colvars.size(); i++) {
-    h_centers[i].type ((colvars[i]->value()).type()); 
+    h_centers[i].type ((colvars[i]->value()).type());
   }
   {
     // it is safer to read colvarvalue objects one at a time;
@@ -1421,7 +1425,7 @@ std::istream & colvarbias_meta::read_hill (std::istream &is)
   get_keyval (data, "widths", h_widths,
               std::vector<cvm::real> (colvars.size(), (std::sqrt (2.0 * PI) / 2.0)),
               parse_silent);
-  
+
   std::string h_replica = "";
   if (comm != single_replica) {
     get_keyval (data, "replicaID", h_replica, replica_id, parse_silent);
@@ -1517,7 +1521,7 @@ std::ostream & colvarbias_meta::write_restart (std::ostream& os)
   }
 
   return os;
-}  
+}
 
 
 void colvarbias_meta::write_pmf()
@@ -1527,7 +1531,7 @@ void colvarbias_meta::write_pmf()
   pmf->create();
 
   std::string fes_file_name_prefix (cvm::output_prefix);
-  
+
   if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) {
     // if this is not the only free energy integrator, append
     // this bias's name, to distinguish it from the output of the other
@@ -1598,7 +1602,7 @@ void colvarbias_meta::write_replica_state_file()
     cvm::fatal_error ("Error: in opening file \""+
                       replica_state_file+"\" for writing.\n");
 
-  rep_state_os.setf (std::ios::scientific, std::ios::floatfield); 
+  rep_state_os.setf (std::ios::scientific, std::ios::floatfield);
   rep_state_os << "\n"
                << "metadynamics {\n"
                << "  configuration {\n"
@@ -1672,7 +1676,7 @@ std::string colvarbias_meta::hill::output_traj()
      << std::setw (cvm::en_width) << W << "\n";
 
   return os.str();
-}  
+}
 
 
 std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)

lib/colvars/colvarbias_meta.h

@@ -33,7 +33,7 @@ public:
 
   /// Destructor
   virtual ~colvarbias_meta();
-  
+
   virtual cvm::real update();
 
   virtual std::istream & read_restart (std::istream &is);
@@ -84,7 +84,7 @@ protected:
   std::istream & read_hill (std::istream &is);
 
   /// \brief step present in a state file
-  /// 
+  ///
   /// When using grids and reading state files containing them
   /// (multiple replicas), this is used to check whether a hill is
   /// newer or older than the grids
@@ -147,10 +147,10 @@ protected:
   /// time steps, appending the step number to each file
   bool       dump_fes_save;
 
-  /// \brief Whether to use well-tempered metadynamics 
+  /// \brief Whether to use well-tempered metadynamics
-  bool       well_tempered; 
+  bool       well_tempered;
 
-  /// \brief Biasing temperature in well-tempered metadynamics 
+  /// \brief Biasing temperature in well-tempered metadynamics
   cvm::real  bias_temperature;
 
   /// \brief Try to read the restart information by allocating new
@@ -282,7 +282,7 @@ public:
       centers[i] = cv[i]->value();
       widths[i] = cv[i]->width * hill_width;
     }
-    if (cvm::debug()) 
+    if (cvm::debug())
       cvm::log ("New hill, applied to "+cvm::to_str (cv.size())+
                 " collective variables, with centers "+
                 cvm::to_str (centers)+", widths "+
@@ -322,7 +322,7 @@ public:
   /// Destructor
   inline ~hill()
   {}
-  
+
   /// Get the energy
   inline cvm::real energy()
   {

lib/colvars/colvarcomp.cpp

@@ -39,7 +39,7 @@ colvar::cvc::cvc (std::string const &conf)
 }
 
 
-void colvar::cvc::parse_group (std::string const &conf, 
+void colvar::cvc::parse_group (std::string const &conf,
                                char const *group_key,
                                cvm::atom_group &group,
                                bool optional)
@@ -120,7 +120,7 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group)
     for (size_t id = 0; id < 3; id++) {
       // (re)read original positions
       group.read_positions();
-      // change one coordinate 
+      // change one coordinate
       group[ia].pos[id] += cvm::debug_gradients_step_size;
       // (re)do the fit (if defined)
       if (group.b_center || group.b_rotate) {

lib/colvars/colvarcomp.h

@@ -1,6 +1,16 @@
 #ifndef COLVARCOMP_H
 #define COLVARCOMP_H
 
+// Declaration of colvar::cvc base class and derived ones.
+//
+// Future cvc's could be declared on additional header files.
+// After the declaration of a new derived class, its metric
+// functions must be reimplemented as well.
+// If the new cvc has no symmetry or periodicity,
+// this can be done straightforwardly by using the macro:
+// simple_scalar_dist_functions (derived_class)
+
+
 #include <fstream>
 #include <cmath>
 
@@ -40,7 +50,7 @@
 /// 2. add a call to the parser in colvar.C, within the function colvar::colvar() \par
 /// 3. declare the class in colvarcomp.h \par
 /// 4. implement the class in one of the files colvarcomp_*.C
-/// 
+///
 /// </b>
 /// The cvm::atom and cvm::atom_group classes are available to
 /// transparently communicate with the simulation program.  However,
@@ -58,7 +68,7 @@ public:
   std::string name;
 
   /// \brief Description of the type of collective variable
-  /// 
+  ///
   /// Normally this string is set by the parent \link colvar \endlink
   /// object within its constructor, when all \link cvc \endlink
   /// objects are initialized; therefore the main "config string"
@@ -94,7 +104,7 @@ public:
 
   /// \brief Within the constructor, make a group parse its own
   /// options from the provided configuration string
-  void parse_group (std::string const &conf, 
+  void parse_group (std::string const &conf,
                     char const *group_key,
                     cvm::atom_group &group,
                     bool optional = false);
@@ -164,7 +174,7 @@ public:
   /// \brief Square distance between x1 and x2 (can be redefined to
   /// transparently implement constraints, symmetries and
   /// periodicities)
-  /// 
+  ///
   /// colvar::cvc::dist2() and the related functions are
   /// declared as "const" functions, but not "static", because
   /// additional parameters defining the metrics (e.g. the
@@ -182,7 +192,7 @@ public:
   /// to provide a gradient which is compatible with the constraint,
   /// i.e. already deprived of its component normal to the constraint
   /// hypersurface.
-  /// 
+  ///
   /// Finally, another useful application, if you are performing very
   /// many operations with these functions, could be to override the
   /// \link colvarvalue \endlink member functions and access directly
@@ -591,7 +601,7 @@ protected:
   std::vector<cvm::atom_pos>  ref_pos;
 
   /// Geometric center of the reference coordinates
-  cvm::rvector                ref_pos_center;
+  cvm::atom_pos               ref_pos_center;
 
   /// Eigenvector (of a normal or essential mode): will always have zero center
   std::vector<cvm::rvector>   eigenvec;
@@ -691,7 +701,7 @@ protected:
   /// \brief Compute system force on first site only to avoid unwanted
   /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
   bool b_1site_force;
-  
+
 public:
 
   /// Initialize by parsing the configuration
@@ -759,7 +769,7 @@ public:
   static cvm::real switching_function (cvm::real const &r0,
                                        int const &exp_num, int const &exp_den,
                                        cvm::atom &A1, cvm::atom &A2);
-  
+
   template<bool b_gradients>
   /// \brief Calculate a coordination number through the function
   /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
@@ -808,7 +818,7 @@ public:
   static cvm::real switching_function (cvm::real const &r0,
                                        int const &exp_num, int const &exp_den,
                                        cvm::atom &A1, cvm::atom &A2);
-  
+
   virtual cvm::real dist2 (colvarvalue const &x1,
                            colvarvalue const &x2) const;
   virtual colvarvalue dist2_lgrad (colvarvalue const &x1,
@@ -822,7 +832,7 @@ public:
 /// \brief Colvar component: hydrogen bond, defined as the product of
 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
 /// (colvarvalue::type_scalar type, range [0:1])
-class colvar::h_bond 
+class colvar::h_bond
   : public colvar::cvc
 {
 protected:
@@ -887,7 +897,7 @@ public:
 //   alpha_dihedrals();
 //   virtual inline ~alpha_dihedrals() {}
 //   virtual void calc_value();
-//   virtual void calc_gradients(); 
+//   virtual void calc_gradients();
 //   virtual void apply_force (colvarvalue const &force);
 //   virtual cvm::real dist2 (colvarvalue const &x1,
 //                            colvarvalue const &x2) const;
@@ -921,7 +931,7 @@ protected:
   /// List of hydrogen bonds
   std::vector<h_bond *>   hb;
 
-  /// Contribution of the hb terms 
+  /// Contribution of the hb terms
   cvm::real hb_coeff;
 
 public:
@@ -1140,6 +1150,62 @@ public:
 };
 
 
+
+// metrics functions for cvc implementations
+
+// simple definitions of the distance functions; these are useful only
+// for optimization (the type check performed in the default
+// colvarcomp functions is skipped)
+
+// definitions assuming the scalar type
+
+#define simple_scalar_dist_functions(TYPE)                              \
+                                                                        \
+  inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1,          \
+                                        colvarvalue const &x2) const    \
+  {                                                                     \
+    return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
+  }                                                                     \
+                                                                        \
+  inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1,  \
+                                                colvarvalue const &x2) const \
+  {                                                                     \
+    return 2.0 * (x1.real_value - x2.real_value);                       \
+  }                                                                     \
+                                                                        \
+  inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1,  \
+                                                colvarvalue const &x2) const \
+  {                                                                     \
+    return this->dist2_lgrad (x2, x1);                                  \
+  }                                                                     \
+                                                                        \
+  inline cvm::real colvar::TYPE::compare (colvarvalue const &x1,        \
+                                          colvarvalue const &x2) const  \
+  {                                                                     \
+    return this->dist2_lgrad (x1, x2);                                  \
+  }                                                                     \
+                                                                        \
+
+  simple_scalar_dist_functions (distance)
+  // NOTE: distance_z has explicit functions, see below
+  simple_scalar_dist_functions (distance_xy)
+  simple_scalar_dist_functions (distance_inv)
+  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 (orientation_angle)
+  simple_scalar_dist_functions (tilt)
+  simple_scalar_dist_functions (eigenvector)
+  //  simple_scalar_dist_functions (alpha_dihedrals)
+  simple_scalar_dist_functions (alpha_angles)
+  simple_scalar_dist_functions (dihedPC)
+
+
 // metrics functions for cvc implementations with a periodicity
 
 inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1,
@@ -1177,12 +1243,12 @@ inline void colvar::dihedral::wrap (colvarvalue &x) const
   if ((x.real_value - wrap_center) >= 180.0) {
     x.real_value -= 360.0;
     return;
-  } 
+  }
 
   if ((x.real_value - wrap_center) < -180.0) {
     x.real_value += 360.0;
     return;
-  } 
+  }
 
   return;
 }
@@ -1222,68 +1288,16 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const
   if ((x.real_value - wrap_center) >= 180.0) {
     x.real_value -= 360.0;
     return;
-  } 
+  }
 
   if ((x.real_value - wrap_center) < -180.0) {
     x.real_value += 360.0;
     return;
-  } 
+  }
 
   return;
 }
 
-// simple definitions of the distance functions; these are useful only
-// for optimization (the type check performed in the default
-// colvarcomp functions is skipped)
-
-// definitions assuming the scalar type
-
-#define simple_scalar_dist_functions(TYPE)                              \
-                                                                        \
-  inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1,          \
-                                        colvarvalue const &x2) const    \
-  {                                                                     \
-    return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \
-  }                                                                     \
-                                                                        \
-  inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1,  \
-                                                colvarvalue const &x2) const \
-  {                                                                     \
-    return 2.0 * (x1.real_value - x2.real_value);                       \
-  }                                                                     \
-                                                                        \
-  inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1,  \
-                                                colvarvalue const &x2) const \
-  {                                                                     \
-    return this->dist2_lgrad (x2, x1);                                  \
-  }                                                                     \
-                                                                        \
-  inline cvm::real colvar::TYPE::compare (colvarvalue const &x1,        \
-                                          colvarvalue const &x2) const  \
-  {                                                                     \
-    return this->dist2_lgrad (x1, x2);                                  \
-  }                                                                     \
-                                                                        \
-
-  simple_scalar_dist_functions (distance)
-  // NOTE: distance_z has explicit functions, see below 
-  simple_scalar_dist_functions (distance_xy)
-  simple_scalar_dist_functions (distance_inv)
-  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 (orientation_angle)
-  simple_scalar_dist_functions (tilt)
-  simple_scalar_dist_functions (eigenvector)
-  //  simple_scalar_dist_functions (alpha_dihedrals)
-  simple_scalar_dist_functions (alpha_angles)
-  simple_scalar_dist_functions (dihedPC)
-
 
 // Projected distance
 // Differences should always be wrapped around 0 (ignoring wrap_center)
@@ -1422,9 +1436,6 @@ inline cvm::real colvar::orientation::compare (colvarvalue const &x1,
 }
 
 
-
-
-
 #endif
 
 

lib/colvars/colvarcomp_angles.cpp

@@ -77,7 +77,7 @@ void colvar::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 );
 
@@ -246,7 +246,7 @@ void colvar::dihedral::calc_gradients()
   cvm::real   rA = A.norm();
   cvm::rvector B = cvm::rvector::outer (r23, r34);
   cvm::real   rB = B.norm();
-  cvm::rvector C = cvm::rvector::outer (r23, A); 
+  cvm::rvector C = cvm::rvector::outer (r23, A);
   cvm::real   rC = C.norm();
 
   cvm::real const cos_phi = (A*B)/(rA*rB);

lib/colvars/colvarcomp_coordnums.cpp

@@ -72,11 +72,14 @@ cvm::real colvar::coordnum::switching_function (cvm::rvector const &r0_vec,
 
 colvar::coordnum::coordnum (std::string const &conf)
   : distance (conf), b_anisotropic (false), b_group2_center_only (false)
-{ 
+{
   function_type = "coordnum";
   x.type (colvarvalue::type_scalar);
 
   // group1 and group2 are already initialized by distance()
+  if (group1.b_dummy)
+    cvm::fatal_error ("Error: only group2 is allowed to be a dummy atom\n");
+
 
   // need to specify this explicitly because the distance() constructor
   // has set it to true
@@ -105,7 +108,7 @@ colvar::coordnum::coordnum (std::string const &conf)
     cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n");
   }
 
-  get_keyval (conf, "group2CenterOnly", b_group2_center_only, false);
+  get_keyval (conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy);
 }
 
 
@@ -121,18 +124,10 @@ void colvar::coordnum::calc_value()
 {
   x.real_value = 0.0;
 
-  // these are necessary: for each atom, gradients are summed together
-  // by multiple calls to switching_function()
-  group1.reset_atoms_data();
-  group2.reset_atoms_data();
-
-  group1.read_positions();
-  group2.read_positions();
-
   if (b_group2_center_only) {
 
     // create a fake atom to hold the group2 com coordinates
-    cvm::atom group2_com_atom (group2[0]);
+    cvm::atom group2_com_atom;
     group2_com_atom.pos = group2.center_of_mass();
 
     if (b_anisotropic) {
@@ -165,9 +160,10 @@ void colvar::coordnum::calc_gradients()
   if (b_group2_center_only) {
 
     // create a fake atom to hold the group2 com coordinates
-    cvm::atom group2_com_atom (group2[0]);
+    cvm::atom group2_com_atom;
     group2_com_atom.pos = group2.center_of_mass();
 
+
     if (b_anisotropic) {
       for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
         switching_function<true> (r0_vec, en, ed, *ai1, group2_com_atom);
@@ -288,14 +284,6 @@ colvar::h_bond::~h_bond()
 
 void colvar::h_bond::calc_value()
 {
-  // this is necessary, because switching_function() will sum the new
-  // gradient to the current one
-  acceptor.reset_data();
-  donor.reset_data();
-
-  acceptor.read_position();
-  donor.read_position();
-
   x.real_value = colvar::coordnum::switching_function<false> (r0, en, ed, acceptor, donor);
 }
 
@@ -315,40 +303,11 @@ void colvar::h_bond::apply_force (colvarvalue const &force)
 }
 
 
-// Self-coordination number for a group
-
-template<bool calculate_gradients>
-cvm::real colvar::selfcoordnum::switching_function (cvm::real const &r0,
-                                                int const &en,
-                                                int const &ed,
-                                                cvm::atom &A1,
-                                                cvm::atom &A2)
-{
-  cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos);
-  cvm::real const l2 = diff.norm2()/(r0*r0);
-
-  // Assume en and ed are even integers, and avoid sqrt in the following
-  int const en2 = en/2;
-  int const ed2 = ed/2;
-
-  cvm::real const xn = std::pow (l2, en2);
-  cvm::real const xd = std::pow (l2, ed2);
-  cvm::real const func = (1.0-xn)/(1.0-xd);
-
-  if (calculate_gradients) {
-    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
-    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
-    A1.grad += (-1.0)*dFdl2*dl2dx;
-    A2.grad +=        dFdl2*dl2dx;
-  }
-
-  return func;
-}
 
 
 colvar::selfcoordnum::selfcoordnum (std::string const &conf)
  : distance (conf, false)
-{ 
+{
   function_type = "selfcoordnum";
   x.type (colvarvalue::type_scalar);
 
@@ -379,13 +338,9 @@ void colvar::selfcoordnum::calc_value()
 {
   x.real_value = 0.0;
 
-  // for each atom, gradients are summed by multiple calls to switching_function()
-  group1.reset_atoms_data();
-  group1.read_positions();
-
   for (size_t i = 0; i < group1.size() - 1; i++)
     for (size_t j = i + 1; j < group1.size(); j++)
-      x.real_value += switching_function<false> (r0, en, ed, group1[i], group1[j]);
+      x.real_value += colvar::coordnum::switching_function<false> (r0, en, ed, group1[i], group1[j]);
 }
 
 
@@ -393,7 +348,7 @@ void colvar::selfcoordnum::calc_gradients()
 {
   for (size_t i = 0; i < group1.size() - 1; i++)
     for (size_t j = i + 1; j < group1.size(); j++)
-      switching_function<true> (r0, en, ed, group1[i], group1[j]);
+      colvar::coordnum::switching_function<true> (r0, en, ed, group1[i], group1[j]);
 }
 
 void colvar::selfcoordnum::apply_force (colvarvalue const &force)

lib/colvars/colvarcomp_distances.cpp

@@ -113,7 +113,7 @@ void colvar::distance_vec::calc_value()
 }
 
 void colvar::distance_vec::calc_gradients()
-{ 
+{
   // gradients are not stored: a 3x3 matrix for each atom would be
   // needed to store just the identity matrix
 }
@@ -140,7 +140,7 @@ colvar::distance_z::distance_z (std::string const &conf)
   // TODO detect PBC from MD engine (in simple cases)
   // and then update period in real time
   if (period != 0.0)
-    b_periodic = true;    
+    b_periodic = true;
 
   if ((wrap_center != 0.0) && (period == 0.0)) {
     cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component,"
@@ -153,7 +153,7 @@ colvar::distance_z::distance_z (std::string const &conf)
   atom_groups.push_back (&ref1);
   // this group is optional
   parse_group (conf, "ref2", ref2, true);
- 
+
   if (ref2.size()) {
     atom_groups.push_back (&ref2);
     cvm::log ("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
@@ -522,7 +522,7 @@ colvar::gyration::gyration (std::string const &conf)
     cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
   } else {
     atoms.b_center = true;
-    atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0));
+    atoms.ref_pos.assign (1, cvm::atom_pos (0.0, 0.0, 0.0));
   }
 
   x.type (colvarvalue::type_scalar);
@@ -700,7 +700,7 @@ colvar::rmsd::rmsd (std::string const &conf)
   parse_group (conf, "atoms", atoms);
   atom_groups.push_back (&atoms);
 
-  if (atoms.b_dummy) 
+  if (atoms.b_dummy)
     cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom.");
 
   if (atoms.ref_pos_group != NULL) {
@@ -732,7 +732,7 @@ colvar::rmsd::rmsd (std::string const &conf)
 
       std::string ref_pos_col;
       double ref_pos_col_value;
-      
+
       if (get_keyval (conf, "refPositionsCol", ref_pos_col, std::string (""))) {
         // if provided, use PDB column to select coordinates
         bool found = get_keyval (conf, "refPositionsColValue", ref_pos_col_value, 0.0);
@@ -771,7 +771,7 @@ colvar::rmsd::rmsd (std::string const &conf)
   }
 }
 
-  
+
 void colvar::rmsd::calc_value()
 {
   // rotational-translational fit is handled by the atom group
@@ -813,9 +813,9 @@ void colvar::rmsd::calc_force_invgrads()
 {
   atoms.read_system_forces();
   ft.real_value = 0.0;
-    
+
   // Note: gradient square norm is 1/N_atoms
-          
+
   for (size_t ia = 0; ia < atoms.size(); ia++) {
     ft.real_value += atoms[ia].grad * atoms[ia].system_force;
   }
@@ -832,7 +832,7 @@ void colvar::rmsd::calc_Jacobian_derivative()
 
     // gradient of the rotation matrix
     cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
-    // gradients of products of 2 quaternion components 
+    // gradients of products of 2 quaternion components
     cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
     for (size_t ia = 0; ia < atoms.size(); ia++) {
 
@@ -850,15 +850,15 @@ void colvar::rmsd::calc_Jacobian_derivative()
       g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2];
 
       // Gradient of the rotation matrix wrt current Cartesian position
-      grad_rot_mat[0][0] = -2.0 * (g22 + g33); 
+      grad_rot_mat[0][0] = -2.0 * (g22 + g33);
-      grad_rot_mat[1][0] =  2.0 * (g12 + g03); 
+      grad_rot_mat[1][0] =  2.0 * (g12 + g03);
-      grad_rot_mat[2][0] =  2.0 * (g13 - g02); 
+      grad_rot_mat[2][0] =  2.0 * (g13 - g02);
-      grad_rot_mat[0][1] =  2.0 * (g12 - g03); 
+      grad_rot_mat[0][1] =  2.0 * (g12 - g03);
-      grad_rot_mat[1][1] = -2.0 * (g11 + g33); 
+      grad_rot_mat[1][1] = -2.0 * (g11 + g33);
-      grad_rot_mat[2][1] =  2.0 * (g01 + g23); 
+      grad_rot_mat[2][1] =  2.0 * (g01 + g23);
-      grad_rot_mat[0][2] =  2.0 * (g02 + g13); 
+      grad_rot_mat[0][2] =  2.0 * (g02 + g13);
-      grad_rot_mat[1][2] =  2.0 * (g23 - g01); 
+      grad_rot_mat[1][2] =  2.0 * (g23 - g01);
-      grad_rot_mat[2][2] = -2.0 * (g11 + g22); 
+      grad_rot_mat[2][2] = -2.0 * (g11 + g22);
 
       cvm::atom_pos &y = ref_pos[ia];
 
@@ -991,7 +991,7 @@ colvar::eigenvector::eigenvector (std::string const &conf)
                       "and eigenvector must be defined.\n");
   }
 
-  cvm::rvector eig_center (0.0, 0.0, 0.0);
+  cvm::atom_pos eig_center (0.0, 0.0, 0.0);
   for (size_t i = 0; i < atoms.size(); i++) {
     eig_center += eigenvec[i];
   }
@@ -1053,7 +1053,7 @@ colvar::eigenvector::eigenvector (std::string const &conf)
   }
 }
 
-  
+
 void colvar::eigenvector::calc_value()
 {
   x.real_value = 0.0;
@@ -1087,7 +1087,7 @@ void colvar::eigenvector::calc_force_invgrads()
 {
   atoms.read_system_forces();
   ft.real_value = 0.0;
-    
+
   for (size_t ia = 0; ia < atoms.size(); ia++) {
     ft.real_value += eigenvec_invnorm2 * atoms[ia].grad *
       atoms[ia].system_force;
@@ -1101,10 +1101,10 @@ void colvar::eigenvector::calc_Jacobian_derivative()
   cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
   cvm::quaternion &quat0 = atoms.rot.q;
 
-  // gradients of products of 2 quaternion components 
+  // gradients of products of 2 quaternion components
   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
 
-  cvm::atom_pos x_relative; 
+  cvm::atom_pos x_relative;
   cvm::real sum = 0.0;
 
   for (size_t ia = 0; ia < atoms.size(); ia++) {
@@ -1126,15 +1126,15 @@ void colvar::eigenvector::calc_Jacobian_derivative()
 
     // Gradient of the inverse rotation matrix wrt current Cartesian position
     // (transpose of the gradient of the direct rotation)
-    grad_rot_mat[0][0] = -2.0 * (g22 + g33); 
+    grad_rot_mat[0][0] = -2.0 * (g22 + g33);
-    grad_rot_mat[0][1] =  2.0 * (g12 + g03); 
+    grad_rot_mat[0][1] =  2.0 * (g12 + g03);
-    grad_rot_mat[0][2] =  2.0 * (g13 - g02); 
+    grad_rot_mat[0][2] =  2.0 * (g13 - g02);
-    grad_rot_mat[1][0] =  2.0 * (g12 - g03); 
+    grad_rot_mat[1][0] =  2.0 * (g12 - g03);
-    grad_rot_mat[1][1] = -2.0 * (g11 + g33); 
+    grad_rot_mat[1][1] = -2.0 * (g11 + g33);
-    grad_rot_mat[1][2] =  2.0 * (g01 + g23); 
+    grad_rot_mat[1][2] =  2.0 * (g01 + g23);
-    grad_rot_mat[2][0] =  2.0 * (g02 + g13); 
+    grad_rot_mat[2][0] =  2.0 * (g02 + g13);
-    grad_rot_mat[2][1] =  2.0 * (g23 - g01); 
+    grad_rot_mat[2][1] =  2.0 * (g23 - g01);
-    grad_rot_mat[2][2] = -2.0 * (g11 + g22); 
+    grad_rot_mat[2][2] = -2.0 * (g11 + g22);
 
     for (size_t i = 0; i < 3; i++) {
       for (size_t j = 0; j < 3; j++) {
@@ -1143,7 +1143,7 @@ void colvar::eigenvector::calc_Jacobian_derivative()
     }
   }
 
-  jd.real_value = sum * std::sqrt (eigenvec_invnorm2); 
+  jd.real_value = sum * std::sqrt (eigenvec_invnorm2);
 }
 
 

lib/colvars/colvarcomp_protein.cpp

@@ -122,7 +122,7 @@ void colvar::alpha_angles::calc_value()
 
   if (theta.size()) {
 
-    cvm::real const theta_norm = 
+    cvm::real const theta_norm =
       (1.0-hb_coeff) / cvm::real (theta.size());
 
     for (size_t i = 0; i < theta.size(); i++) {
@@ -162,7 +162,7 @@ void colvar::alpha_angles::calc_value()
 
 void colvar::alpha_angles::calc_gradients()
 {
-  for (size_t i = 0; i < theta.size(); i++) 
+  for (size_t i = 0; i < theta.size(); i++)
     (theta[i])->calc_gradients();
 
   for (size_t i = 0; i < hb.size(); i++)
@@ -175,9 +175,9 @@ void colvar::alpha_angles::apply_force (colvarvalue const &force)
 
   if (theta.size()) {
 
-    cvm::real const theta_norm = 
+    cvm::real const theta_norm =
       (1.0-hb_coeff) / cvm::real (theta.size());
-    
+
     for (size_t i = 0; i < theta.size(); i++) {
 
       cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol;
@@ -185,10 +185,10 @@ void colvar::alpha_angles::apply_force (colvarvalue const &force)
                             (1.0 - std::pow (t, (int) 4)) );
 
       cvm::real const dfdt =
-        1.0/(1.0 - std::pow (t, (int) 4)) * 
+        1.0/(1.0 - std::pow (t, (int) 4)) *
         ( (-2.0 * t) + (-1.0*f)*(-4.0 * std::pow (t, (int) 3)) );
 
-      (theta[i])->apply_force (theta_norm * 
+      (theta[i])->apply_force (theta_norm *
                                dfdt * (1.0/theta_tol) *
                                force.real_value );
     }
@@ -216,7 +216,7 @@ void colvar::alpha_angles::apply_force (colvarvalue const &force)
     // of this cvc (alpha)
     // This is true of all cvcs with sub-cvcs, and those
     // that do not calculate explicit gradients
-    // SO: we need a flag giving the availability of 
+    // SO: we need a flag giving the availability of
     // atomic gradients
 colvar::dihedPC::dihedPC (std::string const &conf)
   : cvc (conf)
@@ -260,7 +260,7 @@ colvar::dihedPC::dihedPC (std::string const &conf)
   std::string vecFileName;
   int         vecNumber;
   if (get_keyval (conf, "vectorFile", vecFileName, vecFileName)) {
-    get_keyval (conf, "vectorNumber", vecNumber, 0); 
+    get_keyval (conf, "vectorNumber", vecNumber, 0);
     if (vecNumber < 1)
       cvm::fatal_error ("A positive value of vectorNumber is required.");
 

lib/colvars/colvarcomp_rotations.cpp

@@ -74,7 +74,7 @@ colvar::orientation::orientation (std::string const &conf)
 
 }
 
-  
+
 colvar::orientation::orientation()
   : cvc ()
 {
@@ -160,7 +160,7 @@ void colvar::orientation_angle::calc_value()
 void colvar::orientation_angle::calc_gradients()
 {
   cvm::real const dxdq0 =
-    ( ((rot.q).q0 * (rot.q).q0 < 1.0) ? 
+    ( ((rot.q).q0 * (rot.q).q0 < 1.0) ?
       ((180.0 / PI) * (-2.0) / std::sqrt (1.0 - ((rot.q).q0 * (rot.q).q0))) :
       0.0 );
 

lib/colvars/colvargrid.cpp

@@ -141,7 +141,7 @@ void colvar_grid_gradient::write_1D_integral (std::ostream &os)
 {
   cvm::real bin, min, integral;
   std::vector<cvm::real> int_vals;
-  
+
   os << "#       xi            A(xi)\n";
 
   if ( cv.size() != 1 ) {
@@ -162,7 +162,7 @@ void colvar_grid_gradient::write_1D_integral (std::ostream &os)
   }
 
   for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix), bin += 1.0 ) {
-	  
+
     if (samples) {
       size_t const samples_here = samples->value (ix);
       if (samples_here)

lib/colvars/colvargrid.h

@@ -266,11 +266,11 @@ public:
       }
 
 
-      { 
+      {
         cvm::real nbins = ( upper_boundaries[i].real_value -
                             lower_boundaries[i].real_value ) / widths[i];
         int nbins_round = (int)(nbins+0.5);
-        
+
         if (std::fabs (nbins - cvm::real (nbins_round)) > 1.0E-10) {
           cvm::log ("Warning: grid interval ("+
                     cvm::to_str (lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
@@ -284,7 +284,7 @@ public:
         if (cvm::debug())
           cvm::log ("Number of points is "+cvm::to_str ((int) nbins_round)+
                     " for the colvar no. "+cvm::to_str (i+1)+".\n");
-        
+
         nx_i.push_back (nbins_round);
       }
 
@@ -336,7 +336,7 @@ public:
   {
     return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
   }
-  
+
   /// \brief Same as the standard version, but uses different parameters
   inline colvarvalue bin_to_value_scalar (int const &i_bin,
                                           colvarvalue const &new_offset,
@@ -349,7 +349,7 @@ public:
   inline void set_value (std::vector<int> const &ix,
                          T const &t,
                          size_t const &imult = 0)
-  { 
+  {
     data[this->address (ix)+imult] = t;
     has_data = true;
   }
@@ -367,7 +367,7 @@ public:
   /// \brief Add a constant to all elements (fast loop)
   inline void add_constant (T const &t)
   {
-    for (size_t i = 0; i < nt; i++) 
+    for (size_t i = 0; i < nt; i++)
       data[i] += t;
     has_data = true;
   }
@@ -375,11 +375,11 @@ public:
   /// \brief Multiply all elements by a scalar constant (fast loop)
   inline void multiply_constant (cvm::real const &a)
   {
-    for (size_t i = 0; i < nt; i++) 
+    for (size_t i = 0; i < nt; i++)
       data[i] *= a;
   }
 
- 
+
   /// \brief Get the bin indices corresponding to the provided values of
   /// the colvars
   inline std::vector<int> const get_colvars_index (std::vector<colvarvalue> const &values) const
@@ -432,7 +432,7 @@ public:
 
 
   /// \brief Add data from another grid of the same type
-  /// 
+  ///
   /// Note: this function maps other_grid inside this one regardless
   /// of whether it fits or not.
   void map_grid (colvar_grid<T> const &other_grid)
@@ -482,11 +482,11 @@ public:
     if (other_grid.multiplicity() != this->multiplicity())
       cvm::fatal_error ("Error: trying to sum togetehr two grids with values of "
                         "different multiplicity.\n");
-    if (scale_factor != 1.0) 
+    if (scale_factor != 1.0)
       for (size_t i = 0; i < data.size(); i++) {
         data[i] += scale_factor * other_grid.data[i];
       }
-    else 
+    else
       // skip multiplication if possible
       for (size_t i = 0; i < data.size(); i++) {
         data[i] += other_grid.data[i];
@@ -571,30 +571,30 @@ public:
   std::ostream & write_params (std::ostream &os)
   {
     os << "grid_parameters {\n  n_colvars " << nd << "\n";
-      
+
     os << "  lower_boundaries ";
-    for (size_t i = 0; i < nd; i++) 
+    for (size_t i = 0; i < nd; i++)
       os << " " << lower_boundaries[i];
     os << "\n";
 
     os << "  upper_boundaries ";
-    for (size_t i = 0; i < nd; i++) 
+    for (size_t i = 0; i < nd; i++)
       os << " " << upper_boundaries[i];
     os << "\n";
 
     os << "  widths ";
-    for (size_t i = 0; i < nd; i++) 
+    for (size_t i = 0; i < nd; i++)
       os << " " << widths[i];
     os << "\n";
 
     os << "  sizes ";
-    for (size_t i = 0; i < nd; i++) 
+    for (size_t i = 0; i < nd; i++)
       os << " " << nx[i];
     os << "\n";
 
     os << "}\n";
     return os;
-  } 
+  }
 
 
   bool parse_params (std::string const &conf)
@@ -648,9 +648,9 @@ public:
   {
     for (size_t i = 0; i < nd; i++) {
       if ( (std::sqrt (cv[i]->dist2 (cv[i]->lower_boundary,
-                                     lower_boundaries[i])) > 1.0E-10) || 
+                                     lower_boundaries[i])) > 1.0E-10) ||
            (std::sqrt (cv[i]->dist2 (cv[i]->upper_boundary,
-                                     upper_boundaries[i])) > 1.0E-10) || 
+                                     upper_boundaries[i])) > 1.0E-10) ||
            (std::sqrt (cv[i]->dist2 (cv[i]->width,
                                      widths[i])) > 1.0E-10) ) {
         cvm::fatal_error ("Error: restart information for a grid is "
@@ -669,11 +669,11 @@ public:
       // matter: boundaries should be EXACTLY the same (otherwise,
       // map_grid() should be used)
       if ( (std::fabs (other_grid.lower_boundaries[i] -
-                       lower_boundaries[i]) > 1.0E-10) || 
+                       lower_boundaries[i]) > 1.0E-10) ||
            (std::fabs (other_grid.upper_boundaries[i] -
-                       upper_boundaries[i]) > 1.0E-10) || 
+                       upper_boundaries[i]) > 1.0E-10) ||
            (std::fabs (other_grid.widths[i] -
-                       widths[i]) > 1.0E-10) || 
+                       widths[i]) > 1.0E-10) ||
            (data.size() != other_grid.data.size()) ) {
       cvm::fatal_error ("Error: inconsistency between "
                         "two grids that are supposed to be equal, "
@@ -681,7 +681,7 @@ public:
     }
   }
 }
-  
+
 
 /// \brief Write the grid data without labels, as they are
 /// represented in memory
@@ -714,7 +714,7 @@ public:
 std::istream & read_raw (std::istream &is)
 {
   size_t const start_pos = is.tellg();
-    
+
   for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix)) {
     for (size_t imult = 0; imult < mult; imult++) {
       T new_value;
@@ -845,7 +845,7 @@ std::istream & read_multicol (std::istream &is, bool add = false)
         if ( !(is >> x) ) end_of_file = true;
         bin[i] = value_to_bin_scalar (x, i);
       }
-      if (end_of_file) break;	
+      if (end_of_file) break;
 
       for (size_t imult = 0; imult < mult; imult++) {
         is >> new_value[imult];
@@ -904,7 +904,7 @@ public:
     ++(data[this->address (ix)]);
   }
 
-  /// \brief Get the binned count indexed by ix from the newly read data 
+  /// \brief Get the binned count indexed by ix from the newly read data
   inline size_t const & new_count (std::vector<int> const &ix,
                                    size_t const &imult = 0)
   {
@@ -993,7 +993,7 @@ public:
       A0 += data[address (ix)];
       grad[n] = 0.5 * (A1 - A0) / widths[n];
     }
-    return grad; 
+    return grad;
   }
 
   /// \brief Return the value of the function at ix divided by its

lib/colvars/colvarmodule.cpp

@@ -9,14 +9,14 @@
 
 colvarmodule::colvarmodule (char const  *config_filename,
                             colvarproxy *proxy_in)
-{ 
+{
   // pointer to the proxy object
   if (proxy == NULL) {
     proxy = proxy_in;
     parse = new colvarparse();
   } else {
-    cvm::fatal_error ("Error: trying to allocate twice the collective "
+    cvm::fatal_error ("Error: trying to allocate the collective "
-                      "variable module.\n");
+                      "variable module twice.\n");
   }
 
   cvm::log (cvm::line_marker);
@@ -44,6 +44,11 @@ colvarmodule::colvarmodule (char const  *config_filename,
     config_s.close();
   }
 
+  std::string index_file_name;
+  if (parse->get_keyval (conf, "indexFile", index_file_name)) {
+    read_index_file (index_file_name.c_str());
+  }
+
   parse->get_keyval (conf, "analysis", b_analysis, false);
 
   parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07,
@@ -80,12 +85,10 @@ colvarmodule::colvarmodule (char const  *config_filename,
              std::string (output_prefix+".colvars.state") :
              std::string ("colvars.state"))+"\".\n");
 
-  cv_traj_name = 
+  cv_traj_name =
     (output_prefix.size() ?
      std::string (output_prefix+".colvars.traj") :
      std::string ("colvars.traj"));
-  cvm::log ("The trajectory file will be \""+
-            cv_traj_name+"\".\n");
 
   if (cv_traj_freq) {
     // open trajectory file
@@ -94,6 +97,8 @@ colvarmodule::colvarmodule (char const  *config_filename,
                 "\".\n");
       cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
     } else {
+      cvm::log ("Writing to colvar trajectory file \""+cv_traj_name+
+                "\".\n");
       proxy->backup_file (cv_traj_name.c_str());
       cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
     }
@@ -135,7 +140,7 @@ colvarmodule::colvarmodule (char const  *config_filename,
 }
 
 
-std::istream & colvarmodule::read_restart (std::istream &is) 
+std::istream & colvarmodule::read_restart (std::istream &is)
 {
   {
     // read global restart information
@@ -192,7 +197,7 @@ void colvarmodule::init_colvars (std::string const &conf)
       (colvars.back())->check_keywords (colvar_conf, "colvar");
       cvm::decrease_depth();
     } else {
-      cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n"); 
+      cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n");
     }
     colvar_conf = "";
   }
@@ -347,13 +352,13 @@ void colvarmodule::calc() {
   }
 
   // calculate collective variables and their gradients
-  if (cvm::debug()) 
+  if (cvm::debug())
     cvm::log ("Calculating collective variables.\n");
   cvm::increase_depth();
   for (std::vector<colvar *>::iterator cvi = colvars.begin();
        cvi != colvars.end();
        cvi++) {
-    (*cvi)->calc(); 
+    (*cvi)->calc();
   }
   cvm::decrease_depth();
 
@@ -365,7 +370,7 @@ void colvarmodule::calc() {
   for (std::vector<colvarbias *>::iterator bi = biases.begin();
        bi != biases.end();
        bi++) {
-    total_bias_energy += (*bi)->update(); 
+    total_bias_energy += (*bi)->update();
   }
   cvm::decrease_depth();
 
@@ -376,12 +381,12 @@ void colvarmodule::calc() {
   for (std::vector<colvar *>::iterator cvi = colvars.begin();
        cvi != colvars.end();
        cvi++) {
-    (*cvi)->reset_bias_force(); 
+    (*cvi)->reset_bias_force();
   }
   for (std::vector<colvarbias *>::iterator bi = biases.begin();
        bi != biases.end();
        bi++) {
-    (*bi)->communicate_forces(); 
+    (*bi)->communicate_forces();
   }
   cvm::decrease_depth();
 
@@ -393,12 +398,12 @@ void colvarmodule::calc() {
     for (std::vector<colvar *>::iterator cvi = colvars.begin();
          cvi != colvars.end();
          cvi++) {
-      (*cvi)->analyse(); 
+      (*cvi)->analyse();
     }
     for (std::vector<colvarbias *>::iterator bi = biases.begin();
          bi != biases.end();
          bi++) {
-      (*bi)->analyse(); 
+      (*bi)->analyse();
     }
     cvm::decrease_depth();
   }
@@ -426,7 +431,7 @@ void colvarmodule::calc() {
        cvi != colvars.end();
        cvi++) {
     if ((*cvi)->tasks[colvar::task_gradients])
-      (*cvi)->communicate_forces(); 
+      (*cvi)->communicate_forces();
   }
   cvm::decrease_depth();
 
@@ -442,7 +447,7 @@ void colvarmodule::calc() {
       if (!write_restart (restart_out_os))
         cvm::fatal_error ("Error: in writing restart file.\n");
       restart_out_os.close();
-    }    
+    }
   }
 
   // write trajectory file, if needed
@@ -471,15 +476,20 @@ void colvarmodule::calc() {
     if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
       cv_traj_os << "# " << cvm::wrap_string ("step", cvm::it_width-2)
                  << " ";
-      if (cvm::debug()) 
+      if (cvm::debug())
         cv_traj_os.flush();
       for (std::vector<colvar *>::iterator cvi = colvars.begin();
 	   cvi != colvars.end();
 	   cvi++) {
 	(*cvi)->write_traj_label (cv_traj_os);
       }
+      for (std::vector<colvarbias *>::iterator bi = biases.begin();
+	   bi != biases.end();
+	   bi++) {
+	(*bi)->write_traj_label (cv_traj_os);
+      }
       cv_traj_os << "\n";
-      if (cvm::debug()) 
+      if (cvm::debug())
         cv_traj_os.flush();
     }
     cvm::decrease_depth();
@@ -494,13 +504,18 @@ void colvarmodule::calc() {
 	   cvi++) {
         (*cvi)->write_traj (cv_traj_os);
       }
+      for (std::vector<colvarbias *>::iterator bi = biases.begin();
+	   bi != biases.end();
+	   bi++) {
+	(*bi)->write_traj (cv_traj_os);
+      }
       cv_traj_os << "\n";
       if (cvm::debug())
         cv_traj_os.flush();
     }
     cvm::decrease_depth();
 
-    if (restart_out_freq) { 
+    if (restart_out_freq) {
       // flush the trajectory file if we are at the restart frequency
       if ( (cvm::step_relative() > 0) &&
            ((cvm::step_absolute() % restart_out_freq) == 0) ) {
@@ -527,7 +542,7 @@ void colvarmodule::analyze()
        cvi != colvars.end();
        cvi++) {
     cvm::increase_depth();
-    (*cvi)->analyse(); 
+    (*cvi)->analyse();
     cvm::decrease_depth();
   }
 
@@ -536,12 +551,20 @@ void colvarmodule::analyze()
        bi != biases.end();
        bi++) {
     cvm::increase_depth();
-    (*bi)->analyse(); 
+    (*bi)->analyse();
     cvm::decrease_depth();
   }
 
 }
 
+void colvarmodule::setup()
+{
+  // loop over all components of all colvars to reset masses of all groups
+  for (std::vector<colvar *>::iterator cvi = colvars.begin();
+       cvi != colvars.end();  cvi++) {
+    (*cvi)->setup();
+  }
+}
 
 colvarmodule::~colvarmodule()
 {
@@ -565,7 +588,7 @@ colvarmodule::~colvarmodule()
 
   delete parse;
   proxy = NULL;
-}  
+}
 
 
 void colvarmodule::write_output_files()
@@ -590,7 +613,7 @@ void colvarmodule::write_output_files()
     (*cvi)->write_output_files();
   }
   cvm::decrease_depth();
-  
+
   // do not close to avoid problems with multiple NAMD runs
   cv_traj_os.flush();
 }
@@ -631,7 +654,7 @@ bool colvarmodule::read_traj (char const *traj_filename,
 
         continue;
 
-      } else { 
+      } else {
 
         if ((it % 1000) == 0)
           std::cerr << "Reading from trajectory, step = " << it
@@ -721,6 +744,55 @@ void cvm::exit (std::string const &message)
 }
 
 
+void cvm::read_index_file (char const *filename)
+{
+  std::ifstream is (filename);
+  if (!is.good())
+    fatal_error ("Error: in opening index file \""+
+                 std::string (filename)+"\".\n");
+  // std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
+  // std::list<std::vector<int> >::iterator lists_i = cvm::index_groups.begin();
+  while (is.good()) {
+    char open, close;
+    std::string group_name;
+    if ( (is >> open) && (open == '[') &&
+         (is >> group_name) &&
+         (is >> close) && (close == ']') ) {
+      cvm::index_group_names.push_back (group_name);
+      cvm::index_groups.push_back (std::vector<int> ());
+    } else {
+      cvm::fatal_error ("Error: in parsing index file \""+
+                 std::string (filename)+"\".\n");
+    }
+
+    int atom_number = 1;
+    size_t pos = is.tellg();
+    while ( (is >> atom_number) && (atom_number > 0) ) {
+      (cvm::index_groups.back()).push_back (atom_number);
+      pos = is.tellg();
+    }
+    is.clear();
+    is.seekg (pos, std::ios::beg);
+    std::string delim;
+    if ( (is >> delim) && (delim == "[") ) {
+      // new group
+      is.clear();
+      is.seekg (pos, std::ios::beg);
+    } else {
+      break;
+    }
+  }
+
+  cvm::log ("The following index groups were read from the index file \""+
+            std::string (filename)+"\":\n");
+  std::list<std::string>::iterator names_i = cvm::index_group_names.begin();
+  std::list<std::vector<int> >::iterator lists_i = cvm::index_groups.begin();
+  for ( ; names_i != cvm::index_group_names.end() ; names_i++, lists_i++) {
+    cvm::log ("  "+(*names_i)+" ("+cvm::to_str (lists_i->size())+" atoms).\n");
+  }
+
+}
+
 
 // static pointers
 std::vector<colvar *>     colvarmodule::colvars;
@@ -741,6 +813,8 @@ size_t    colvarmodule::cv_traj_freq = 0;
 size_t    colvarmodule::depth = 0;
 bool      colvarmodule::b_analysis = false;
 cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-04;
+std::list<std::string> colvarmodule::index_group_names;
+std::list<std::vector<int> > colvarmodule::index_groups;
 
 
 // file name prefixes
@@ -827,7 +901,7 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
                          std::string ("euler")) ) {
 
     // parse the Euler angles
-    
+
     char sep;
     cvm::real phi, theta, psi;
     if ( !(is >> sep)   || !(sep == '(') ||
@@ -846,7 +920,7 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
 
     // parse the quaternion components
 
-    is.seekg (start_pos, std::ios::beg);      
+    is.seekg (start_pos, std::ios::beg);
     char sep;
     if ( !(is >> sep)  || !(sep == '(') ||
          !(is >> q.q0) || !(is >> sep)  || !(sep == ',') ||
@@ -932,15 +1006,13 @@ cvm::quaternion::position_derivative_inner (cvm::rvector const &pos,
 
 
 
+
+
 // Calculate the optimal rotation between two groups, and implement it
 // as a quaternion.  The method is the one documented in: Coutsias EA,
 // Seok C, Dill KA.  Using quaternions to calculate RMSD.  J Comput
 // Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254
 
-
-
-
-
 void colvarmodule::rotation::build_matrix (std::vector<cvm::atom_pos> const &pos1,
                                            std::vector<cvm::atom_pos> const &pos2,
                                            matrix2d<cvm::real, 4, 4>        &S)
@@ -959,7 +1031,7 @@ void colvarmodule::rotation::build_matrix (std::vector<cvm::atom_pos> const &pos
     C.zx() += pos1[i].z * pos2[i].x;
     C.zy() += pos1[i].z * pos2[i].y;
     C.zz() += pos1[i].z * pos2[i].z;
-  } 
+  }
 
   // build the "overlap" matrix, whose eigenvectors are stationary
   // points of the RMSD in the space of rotations
@@ -1068,11 +1140,12 @@ void colvarmodule::rotation::calc_optimal_rotation
   if (q_old.norm2() > 0.0) {
     q.match (q_old);
     if (q_old.inner (q) < (1.0 - crossing_threshold)) {
-      cvm::log ("Warning: discontinuous rotation!\n");
+      cvm::log ("Warning: one molecular orientation has changed by more than "+
+                cvm::to_str (crossing_threshold)+": discontinuous rotation ?\n");
     }
   }
   q_old = q;
-                     
+
   if (cvm::debug()) {
     if (b_debug_gradients) {
       cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+
@@ -1143,8 +1216,8 @@ void colvarmodule::rotation::calc_optimal_rotation
       for (size_t i = 0; i < 4; i++) {
         for (size_t j = 0; j < 4; j++) {
           dq0_1[p] +=
-            (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] + 
+            (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
-            (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] + 
+            (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
             (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p];
         }
       }
@@ -1193,15 +1266,15 @@ void colvarmodule::rotation::calc_optimal_rotation
       for (size_t i = 0; i < 4; i++) {
         for (size_t j = 0; j < 4; j++) {
           dq0_2[p] +=
-            (Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] + 
+            (Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
-            (Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] + 
+            (Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
             (Q3[i] * ds_2[i][j] * Q0[j]) / (L0-L3) * Q3[p];
         }
       }
     }
 
     if (cvm::debug()) {
-      
+
       if (b_debug_gradients) {
 
       matrix2d<cvm::real, 4, 4> S_new;
@@ -1216,7 +1289,7 @@ void colvarmodule::rotation::calc_optimal_rotation
           // diagonalize the new overlap matrix
           for (size_t i = 0; i < 4; i++) {
             for (size_t j = 0; j < 4; j++) {
-              S_new[i][j] += 
+              S_new[i][j] +=
                 colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
             }
           }

lib/colvars/colvarmodule.h

@@ -2,16 +2,16 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2013-04-17"
+#define COLVARS_VERSION "2013-06-19"
 #endif
 
 #ifndef COLVARS_DEBUG
 #define COLVARS_DEBUG false
 #endif
 
-/// \file colvarmodule.h 
+/// \file colvarmodule.h
 /// \brief Collective variables main module
-/// 
+///
 /// This file declares the main class for defining and manipulating
 /// collective variables: there should be only one instance of this
 /// class, because several variables are made static (i.e. they are
@@ -40,7 +40,7 @@ class colvarproxy;
 /// Class to control the collective variables calculation.  An object
 /// (usually one) of this class is spawned from the MD program,
 /// containing all i/o routines and general interface.
-/// 
+///
 /// At initialization, the colvarmodule object creates a proxy object
 /// to provide a transparent interface between the MD program and the
 /// child objects
@@ -54,7 +54,7 @@ private:
 public:
 
   friend class colvarproxy;
-  
+
   /// Defining an abstract real number allows to switch precision
   typedef  double    real;
   /// Residue identifier
@@ -147,7 +147,7 @@ public:
   /// \brief Number of histograms initialized (no limit on the
   /// number)
   static size_t n_histo_biases;
-   
+
   /// \brief Whether debug output should be enabled (compile-time option)
   static inline bool debug()
   {
@@ -157,8 +157,9 @@ public:
 
   /// \brief Constructor \param config_name Configuration file name
   /// \param restart_name (optional) Restart file name
-  colvarmodule (char const *config_name, 
+  colvarmodule (char const *config_name,
                 colvarproxy *proxy_in);
+
   /// Destructor
   ~colvarmodule();
 
@@ -168,6 +169,11 @@ public:
   /// Initialize collective variable biases
   void init_biases (std::string const &conf);
 
+  /// Re-initialize data at the beginning of a run. For use with
+  /// MD codes that can change system parameters like atom masses
+  /// between run commands.
+  void setup();
+
   /// Load new configuration for the given bias -
   /// currently works for harmonic (force constant and/or centers)
   void change_configuration (std::string const &bias_name, std::string const &conf);
@@ -194,16 +200,16 @@ public:
   bool read_traj (char const *traj_filename,
                   size_t      traj_read_begin,
                   size_t      traj_read_end);
- 
+
   /// Get the pointer of a colvar from its name (returns NULL if not found)
   static colvar * colvar_p (std::string const &name);
 
   /// Quick conversion of an object to a string
-  template<typename T> static std::string to_str (T const &x, 
+  template<typename T> static std::string to_str (T const &x,
                                                   size_t const &width = 0,
                                                   size_t const &prec = 0);
   /// Quick conversion of a vector of objects to a string
-  template<typename T> static std::string to_str (std::vector<T> const &x, 
+  template<typename T> static std::string to_str (std::vector<T> const &x,
                                                   size_t const &width = 0,
                                                   size_t const &prec = 0);
 
@@ -247,10 +253,10 @@ public:
 
   /// \brief Time step of MD integrator (fs)
   static real dt();
-  
+
   /// Request calculation of system force from MD engine
   static void request_system_force();
-  
+
   /// Print a message to the main log
   static void log (std::string const &message);
 
@@ -278,7 +284,7 @@ public:
 
   /// \brief Get the closest periodic image to a reference position
   /// \param pos The position to look for the closest periodic image
-  /// \param ref_pos (optional) The reference position 
+  /// \param ref_pos (optional) The reference position
   static void select_closest_image (atom_pos &pos,
                                     atom_pos const &ref_pos);
 
@@ -290,6 +296,16 @@ public:
                                      atom_pos const &ref_pos);
 
 
+  /// \brief Names of groups from a Gromacs .ndx file to be read at startup
+  static std::list<std::string> index_group_names;
+
+  /// \brief Groups from a Gromacs .ndx file read at startup
+  static std::list<std::vector<int> > index_groups;
+
+  /// \brief Read a Gromacs .ndx file
+  static void read_index_file (char const *filename);
+
+
   /// \brief Create atoms from a file \param filename name of the file
   /// (usually a PDB) \param atoms array of the atoms to be allocated
   /// \param pdb_field (optiona) if "filename" is a PDB file, use this
@@ -371,7 +387,7 @@ std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
 std::istream & operator >> (std::istream &is, cvm::rvector &v);
 
 
-template<typename T> std::string cvm::to_str (T const &x, 
+template<typename T> std::string cvm::to_str (T const &x,
                                               size_t const &width,
                                               size_t const &prec) {
   std::ostringstream os;
@@ -384,7 +400,7 @@ template<typename T> std::string cvm::to_str (T const &x,
   return os.str();
 }
 
-template<typename T> std::string cvm::to_str (std::vector<T> const &x, 
+template<typename T> std::string cvm::to_str (std::vector<T> const &x,
                                               size_t const &width,
                                               size_t const &prec) {
   if (!x.size()) return std::string ("");
@@ -434,7 +450,7 @@ inline void cvm::request_system_force()
 {
   proxy->request_system_force (true);
 }
-  
+
 inline void cvm::select_closest_image (atom_pos &pos,
                                        atom_pos const &ref_pos)
 {

lib/colvars/colvarparse.cpp

@@ -241,11 +241,11 @@ _get_keyval_vector_ (colvarvalue);
 
 
 
-bool colvarparse::get_keyval (std::string const &conf,        
+bool colvarparse::get_keyval (std::string const &conf,
-                              char const *key,                
+                              char const *key,
-                              bool &value,                    
+                              bool &value,
-                              bool const &def_value,          
+                              bool const &def_value,
-                              Parse_Mode const parse_mode)    
+                              Parse_Mode const parse_mode)
 {
   std::string data;
   bool b_found = false, b_found_any = false;
@@ -277,7 +277,7 @@ 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::fatal_error ("Error: boolean values only are allowed "
                         "for \""+std::string (key)+"\".\n");
     if (parse_mode != parse_silent) {
       cvm::log ("# "+std::string (key)+" = "+
@@ -371,8 +371,8 @@ void colvarparse::check_keywords (std::string &conf, char const *key)
     std::string uk;
     std::istringstream line_is (line);
     line_is >> uk;
-    if (cvm::debug())
+    // if (cvm::debug())
-      cvm::log ("Checking the validity of \""+uk+"\" from line:\n" + line);
+    //   cvm::log ("Checking the validity of \""+uk+"\" from line:\n" + line);
     uk = to_lower_cppstr (uk);
 
     bool found_keyword = false;
@@ -426,8 +426,8 @@ bool colvarparse::key_lookup (std::string const &conf,
   colvarparse::dummy_pos = 0;
 
   // start from the first occurrence of key
-  size_t pos = conf_lower.find (key, save_pos); 
+  size_t pos = conf_lower.find (key, save_pos);
-  
+
   // iterate over all instances until it finds the isolated keyword
   while (true) {
 
@@ -442,7 +442,7 @@ bool colvarparse::key_lookup (std::string const &conf,
       if ( std::string ("\n"+white_space+
                         "}").find (conf[pos-1]) ==
            std::string::npos ) {
-        // none of the valid delimiting characters is on the left of key 
+        // none of the valid delimiting characters is on the left of key
         b_isolated_left = false;
       }
     }
@@ -451,7 +451,7 @@ bool colvarparse::key_lookup (std::string const &conf,
       if ( std::string ("\n"+white_space+
                         "{").find (conf[pos+key.size()]) ==
            std::string::npos ) {
-        // none of the valid delimiting characters is on the right of key 
+        // none of the valid delimiting characters is on the right of key
         b_isolated_right = false;
       }
     }
@@ -461,7 +461,7 @@ bool colvarparse::key_lookup (std::string const &conf,
 
     bool const b_isolated = (b_isolated_left && b_isolated_right &&
                              b_not_within_block);
-    
+
     if (b_isolated) {
       // found it
       break;
@@ -490,7 +490,7 @@ bool colvarparse::key_lookup (std::string const &conf,
   size_t line_begin = (pl == std::string::npos) ? 0 : pos;
   size_t nl = conf.find  ("\n", pos);
   size_t line_end = (nl == std::string::npos) ? conf.size() : nl;
-  std::string line (conf, line_begin, (line_end-line_begin)); 
+  std::string line (conf, line_begin, (line_end-line_begin));
 
   size_t data_begin = (to_lower_cppstr (line)).find (key) + key.size();
   data_begin = line.find_first_not_of (white_space, data_begin+1);
@@ -527,9 +527,9 @@ bool colvarparse::key_lookup (std::string const &conf,
 
           line_begin = line_end;
           nl = conf.find ('\n', line_begin+1);
-          if (nl == std::string::npos) 
+          if (nl == std::string::npos)
             line_end = conf.size();
-          else 
+          else
             line_end = nl;
           line.append (conf, line_begin, (line_end-line_begin));
 
@@ -565,7 +565,7 @@ bool colvarparse::key_lookup (std::string const &conf,
       //                 << "\n";
     }
   }
-  
+
   save_pos = line_end;
 
   return true;

lib/colvars/colvarparse.h

@@ -53,15 +53,15 @@ public:
   {}
 
   /// How a keyword is parsed in a string
-  enum Parse_Mode { 
+  enum Parse_Mode {
     /// \brief (default) Read the first instance of a keyword (if
     /// any), report its value, and print a warning when there is more
     /// than one
-    parse_normal, 
+    parse_normal,
     /// \brief Like parse_normal, but don't send any message to the log
     /// (useful e.g. in restart files when such messages are very
     /// numerous and redundant)
-    parse_silent 
+    parse_silent
   };
 
   /// \fn get_keyval bool const get_keyval (std::string const &conf,
@@ -112,7 +112,7 @@ public:
                    std::vector<_type_> const &def_values =              \
                    std::vector<_type_> (0, static_cast<_type_>(_def_value_)),                \
                    Parse_Mode const parse_mode = parse_normal)
-  
+
     _get_keyval_vector_proto_ (int, 0);
     _get_keyval_vector_proto_ (size_t, 0);
     _get_keyval_vector_proto_ (std::string, std::string (""));
@@ -169,7 +169,7 @@ public:
   /// data (optional) holds the string provided after "key", if any
   /// \param save_pos (optional) stores the position of the keyword
   /// within "conf", useful when doing multiple calls \param
-  /// save_delimiters (optional) 
+  /// save_delimiters (optional)
   bool key_lookup (std::string const &conf,
                    char const *key,
                    std::string &data = dummy_string,

lib/colvars/colvarproxy.h

@@ -16,6 +16,12 @@ public:
   /// Pointer to the instance of colvarmodule
   colvarmodule *colvars;
 
+  /// Default destructor
+  virtual inline ~colvarproxy() {}
+
+
+  // **************** SYSTEM-WIDE PHYSICAL QUANTITIES ****************
+
   /// \brief Value of the unit for atomic coordinates with respect to
   /// angstroms (used by some variables for hard-coded default values)
   virtual cvm::real unit_angstrom() = 0;
@@ -29,21 +35,11 @@ public:
   /// \brief Time step of the simulation (fs)
   virtual cvm::real dt() = 0;
 
-  /// Pass restraint energy value for current timestep to MD engine
+  /// \brief Pseudo-random number with Gaussian distribution
-  virtual void add_energy (cvm::real energy) = 0;
+  virtual cvm::real rand_gaussian (void) = 0;
 
-  /// Tell the proxy whether system forces are needed
-  virtual void request_system_force (bool yesno) = 0;
-
-  /// Print a message to the main log
-  virtual void log (std::string const &message) = 0;
-
-  /// Print a message to the main log and exit with error code
-  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;
 
+  // **************** SIMULATION PARAMETERS ****************
 
   /// \brief Prefix to be used for input files (restarts, not
   /// configuration)
@@ -59,25 +55,32 @@ public:
   /// \brief Restarts will be fritten 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 simple distance vector between two positions
+  /// \brief Get the PBC-aware distance vector between two positions
-  /// (with periodic boundary conditions handled transparently)
   virtual cvm::rvector position_distance (cvm::atom_pos const &pos1,
                                           cvm::atom_pos const &pos2) = 0;
 
-  /// \brief Get the square distance between two positions (with
+  /// \brief Get the PBC-aware square distance between two positions;
-  /// periodic boundary conditions handled transparently)
+  /// may be implemented independently from position_distance() for optimization purposes
-  ///
-  /// Note: in the case of periodic boundary conditions, this provides
-  /// an analytical square distance (while taking the square of
-  /// position_distance() would produce leads to a cusp)
   virtual cvm::real position_dist2 (cvm::atom_pos const &pos1,
-                                    cvm::atom_pos const &pos2) = 0;
+                                    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 
+  /// \param ref_pos The reference position
   virtual void select_closest_image (cvm::atom_pos &pos,
                                      cvm::atom_pos const &ref_pos) = 0;
 
@@ -89,6 +92,18 @@ public:
                               cvm::atom_pos const &ref_pos);
 
 
+
+  // **************** INPUT/OUTPUT ****************
+
+  /// Print a message to the main log
+  virtual void log (std::string const &message) = 0;
+
+  /// Print a message to the main log and exit with error code
+  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;
+
   /// \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
@@ -97,7 +112,7 @@ public:
   virtual void load_atoms (char const *filename,
                            std::vector<cvm::atom> &atoms,
                            std::string const pdb_field,
-                           double const pdb_field_value = 0.0) = 0;
+                           double const pdb_field_value = 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
@@ -108,14 +123,9 @@ public:
                             std::string const pdb_field,
                             double const pdb_field_value = 0.0) = 0;
 
-  /// \brief Rename the given file, under the convention provided by
+  /// \brief Rename the given file, before overwriting it
-  /// the MD program
+  virtual void backup_file (char const *filename) {}
-  virtual void backup_file (char const *filename) = 0;
 
-  /// \brief Pseudo-random number with Gaussian distribution
-  virtual cvm::real rand_gaussian (void) = 0;
-
-  virtual inline ~colvarproxy() {}
 };
 
 
@@ -128,6 +138,12 @@ 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)
+{
+  return (position_distance (pos1, pos2)).norm2();
+}
+
 #endif
 
 

lib/colvars/colvartypes.h

@@ -19,7 +19,7 @@ class colvarmodule::rvector {
 public:
 
   cvm::real x, y, z;
-     
+
   inline rvector()
     : x (0.0), y (0.0), z (0.0)
   {}
@@ -55,7 +55,7 @@ public:
   }
 
 
-  inline cvm::rvector & operator = (cvm::real const &v) 
+  inline cvm::rvector & operator = (cvm::real const &v)
   {
     x = v;
     y = v;
@@ -63,28 +63,28 @@ public:
     return *this;
   }
 
-  inline void operator += (cvm::rvector const &v) 
+  inline void operator += (cvm::rvector const &v)
   {
     x += v.x;
     y += v.y;
     z += v.z;
   }
 
-  inline void operator -= (cvm::rvector const &v) 
+  inline void operator -= (cvm::rvector const &v)
   {
     x -= v.x;
     y -= v.y;
     z -= v.z;
   }
 
-  inline void operator *= (cvm::real const &v) 
+  inline void operator *= (cvm::real const &v)
   {
     x *= v;
     y *= v;
     z *= v;
   }
 
-  inline void operator /= (cvm::real const& v) 
+  inline void operator /= (cvm::real const& v)
   {
     x /= v;
     y /= v;
@@ -113,57 +113,57 @@ public:
   }
 
 
-  static inline cvm::rvector outer (cvm::rvector const &v1, cvm::rvector const &v2) 
+  static inline cvm::rvector outer (cvm::rvector const &v1, cvm::rvector const &v2)
   {
     return cvm::rvector ( v1.y*v2.z - v2.y*v1.z,
                          -v1.x*v2.z + v2.x*v1.z,
                           v1.x*v2.y - v2.x*v1.y);
   }
 
-  friend inline cvm::rvector operator - (cvm::rvector const &v) 
+  friend inline cvm::rvector operator - (cvm::rvector const &v)
   {
     return cvm::rvector (-v.x, -v.y, -v.z);
   }
 
-  friend inline int operator == (cvm::rvector const &v1, cvm::rvector const &v2) 
+  friend inline int operator == (cvm::rvector const &v1, cvm::rvector const &v2)
   {
     return (v1.x == v2.x) && (v1.y == v2.y) && (v1.z == v2.z);
   }
 
-  friend inline int operator != (cvm::rvector const &v1, cvm::rvector const &v2) 
+  friend inline int operator != (cvm::rvector const &v1, cvm::rvector const &v2)
   {
     return (v1.x != v2.x) || (v1.y != v2.y) || (v1.z != v2.z);
   }
 
-  friend inline cvm::rvector operator + (cvm::rvector const &v1, cvm::rvector const &v2) 
+  friend inline cvm::rvector operator + (cvm::rvector const &v1, cvm::rvector const &v2)
   {
     return cvm::rvector (v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
   }
-  friend inline cvm::rvector operator - (cvm::rvector const &v1, cvm::rvector const &v2) 
+  friend inline cvm::rvector operator - (cvm::rvector const &v1, cvm::rvector const &v2)
   {
     return cvm::rvector (v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
   }
 
-  friend inline cvm::real operator * (cvm::rvector const &v1, cvm::rvector const &v2) 
+  friend inline cvm::real operator * (cvm::rvector const &v1, cvm::rvector const &v2)
   {
     return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
   }
 
-  friend inline cvm::rvector operator * (cvm::real const &a, cvm::rvector const &v) 
+  friend inline cvm::rvector operator * (cvm::real const &a, cvm::rvector const &v)
   {
     return cvm::rvector (a*v.x, a*v.y, a*v.z);
   }
 
-  friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real const &a) 
+  friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real const &a)
   {
     return cvm::rvector (a*v.x, a*v.y, a*v.z);
   }
 
-  friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real const &a) 
+  friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real const &a)
   {
     return cvm::rvector (v.x/a, v.y/a, v.z/a);
   }
-     
+
 
 };
 
@@ -186,7 +186,7 @@ public:
   {
     return length;
   }
- 
+
   /// Default constructor
   inline vector1d (T const &t = T())
   {
@@ -251,7 +251,7 @@ public:
     return prod;
   }
 
-  /// Formatted output 
+  /// Formatted output
   friend std::ostream & operator << (std::ostream &os,
                                      vector1d<T, length> const &v)
   {
@@ -310,7 +310,7 @@ public:
     this->alloc();
     reset();
   }
- 
+
   /// Constructor from a 2-d C array
   inline matrix2d (T const **m)
   {
@@ -380,7 +380,7 @@ public:
 //     }
 //   }
 
-  /// Formatted output 
+  /// Formatted output
   friend std::ostream & operator << (std::ostream &os,
                                      matrix2d<T, outer_length, inner_length> const &m)
   {
@@ -454,16 +454,16 @@ public:
   inline cvm::real zz() const { return array[2][2]; }
 
   /// Constructor from a 2-d C array
-  inline rmatrix (cvm::real const **m) 
+  inline rmatrix (cvm::real const **m)
-    : cvm::matrix2d<cvm::real, 3, 3> (m) 
+    : cvm::matrix2d<cvm::real, 3, 3> (m)
   {}
 
   /// Default constructor
-  inline rmatrix() 
+  inline rmatrix()
     : cvm::matrix2d<cvm::real, 3, 3>()
   {}
 
-  /// Constructor component by component 
+  /// Constructor component by component
   inline rmatrix (cvm::real const &xxi,
                   cvm::real const &xyi,
                   cvm::real const &xzi,
@@ -472,7 +472,7 @@ public:
                   cvm::real const &yzi,
                   cvm::real const &zxi,
                   cvm::real const &zyi,
-                  cvm::real const &zzi) 
+                  cvm::real const &zzi)
     : cvm::matrix2d<cvm::real, 3, 3>()
   {
     this->xx() = xxi;
@@ -488,7 +488,7 @@ public:
 
   /// Destructor
   inline ~rmatrix()
-  {}    
+  {}
 
   /// Return the determinant
   inline cvm::real determinant() const
@@ -528,7 +528,7 @@ public:
 
 };
 
-                    
+
 inline cvm::rvector operator * (cvm::rmatrix const &m,
                                 cvm::rvector const &r)
 {
@@ -680,7 +680,7 @@ public:
     return std::sqrt (this->norm2());
   }
 
-  /// Return the conjugate quaternion 
+  /// Return the conjugate quaternion
   inline cvm::quaternion conjugate() const
   {
     return cvm::quaternion (q0, -q1, -q2, -q3);
@@ -720,7 +720,7 @@ public:
     return cvm::quaternion (0.0, v.x, v.y, v.z);
   }
   /// Return the vector component
-  inline cvm::rvector get_vector() const 
+  inline cvm::rvector get_vector() const
   {
     return cvm::rvector (q1, q2, q3);
   }
@@ -1002,11 +1002,11 @@ public:
     if (q.q0 != 0.0) {
 
       // cvm::real const x = iprod/q.q0;
-      
+
       cvm::real const dspindx = (180.0/PI) * 2.0 * (1.0 / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
 
-      return 
+      return
-        cvm::quaternion ( dspindx * (iprod * (-1.0) / (q.q0*q.q0)), 
+        cvm::quaternion ( dspindx * (iprod * (-1.0) / (q.q0*q.q0)),
                           dspindx * ((1.0/q.q0) * axis.x),
                           dspindx * ((1.0/q.q0) * axis.y),
                           dspindx * ((1.0/q.q0) * axis.z));
@@ -1023,11 +1023,11 @@ public:
   inline cvm::real cos_theta (cvm::rvector const &axis) const
   {
     cvm::rvector const q_vec = q.get_vector();
-    cvm::real const alpha = 
+    cvm::real const alpha =
       (180.0/PI) * 2.0 * std::atan2 (axis * q_vec, q.q0);
 
     cvm::real const cos_spin_2 = std::cos (alpha * (PI/180.0) * 0.5);
-    cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ? 
+    cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ?
                                     (q.q0 / cos_spin_2) :
                                     (0.0) );
     // cos(2t) = 2*cos(t)^2 - 1
@@ -1044,7 +1044,7 @@ public:
 
     if (q.q0 != 0.0)  {
 
-      cvm::real const d_cos_theta_dq0 = 
+      cvm::real const d_cos_theta_dq0 =
         (4.0 * q.q0 / (cos_spin_2*cos_spin_2)) *
         (1.0 - (iprod*iprod)/(q.q0*q.q0) / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
 
@@ -1052,7 +1052,7 @@ public:
         (4.0 * q.q0 / (cos_spin_2*cos_spin_2) *
          (iprod/q.q0) / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
 
-      return cvm::quaternion (d_cos_theta_dq0, 
+      return cvm::quaternion (d_cos_theta_dq0,
                               d_cos_theta_dqn * axis.x,
                               d_cos_theta_dqn * axis.y,
                               d_cos_theta_dqn * axis.z);

lib/colvars/colvarvalue.cpp

@@ -28,7 +28,7 @@ void colvarvalue::error_rside
   cvm::fatal_error ("Trying to assign a colvar value with type \""+
                     type_desc[this->value_type]+"\" to one with type \""+
                     type_desc[vt]+"\".\n");
-}    
+}
 
 void colvarvalue::error_lside
 (colvarvalue::Type const &vt) const
@@ -36,7 +36,7 @@ void colvarvalue::error_lside
   cvm::fatal_error ("Trying to use a colvar value with type \""+
                     type_desc[vt]+"\" as one of type \""+
                     type_desc[this->value_type]+"\".\n");
-}    
+}
 
 
 
@@ -125,7 +125,7 @@ void colvarvalue::p2leg_opt (colvarvalue const                        &x,
     break;
   case colvarvalue::type_vector:
     while (xvi != xv_end) {
-      cvm::real const cosine = 
+      cvm::real const cosine =
         ((xvi)->rvector_value * x.rvector_value) /
         ((xvi)->rvector_value.norm() * x.rvector_value.norm());
       xvi++;
@@ -167,7 +167,7 @@ void colvarvalue::p2leg_opt (colvarvalue const                        &x,
     break;
   case colvarvalue::type_vector:
     while (xvi != xv_end) {
-      cvm::real const cosine = 
+      cvm::real const cosine =
         ((xvi)->rvector_value * x.rvector_value) /
         ((xvi)->rvector_value.norm() * x.rvector_value.norm());
       xvi++;
@@ -225,11 +225,11 @@ std::istream & operator >> (std::istream &is, colvarvalue &x)
 
   switch (x.type()) {
   case colvarvalue::type_scalar:
-    is >> x.real_value; 
+    is >> x.real_value;
     break;
   case colvarvalue::type_vector:
   case colvarvalue::type_unitvector:
-    is >> x.rvector_value; 
+    is >> x.rvector_value;
     break;
   case colvarvalue::type_quaternion:
     is >> x.quaternion_value;
@@ -242,7 +242,7 @@ std::istream & operator >> (std::istream &is, colvarvalue &x)
 }
 
 
-size_t colvarvalue::output_width (size_t const &real_width)
+size_t colvarvalue::output_width (size_t const &real_width) const
 {
   switch (this->value_type) {
   case colvarvalue::type_scalar:

lib/colvars/colvarvalue.h

@@ -37,7 +37,7 @@
 /// the problem, because \link colvarvalue \endlink objects are first
 /// initialized in the configuration, and the stream operation will be
 /// performed only when reading restart files.
-/// 
+///
 /// No problem of course with the output streams: \code os << x;
 /// \endcode will print a different output according to the value of
 /// colvarvalue::value_type, and the width of such output is returned
@@ -92,7 +92,7 @@ public:
   cvm::quaternion quaternion_value;
 
   /// Current type of this colvarvalue object
-  Type  value_type; 
+  Type  value_type;
 
   static inline bool type_checking()
   {
@@ -146,7 +146,7 @@ public:
     case type_scalar:
       real_value = x.real_value;
       break;
-    case type_vector: 
+    case type_vector:
     case type_unitvector:
       rvector_value = x.rvector_value;
       break;
@@ -254,7 +254,7 @@ public:
   /// Ensure that the two types are the same within a binary operator
   void static check_types (colvarvalue const &x1, colvarvalue const &x2);
 
-  /// Undefined operation 
+  /// Undefined operation
   void undef_op() const;
 
   /// Trying to assign this \link colvarvalue \endlink object to
@@ -265,10 +265,10 @@ public:
   /// with a different type to this object
   void error_rside (Type const &vt) const;
 
-  ///<EFBFBD>Give the number of characters required to output this
+  /// Give the number of characters required to output this
   /// colvarvalue, given the current type assigned and the number of
   /// characters for a real number
-  size_t output_width (size_t const &real_width);
+  size_t output_width (size_t const &real_width) const;
 
 
   // optimized routines for operations with an array; xv and inner as
@@ -324,7 +324,7 @@ inline void colvarvalue::reset()
   case colvarvalue::type_notset:
   default:
     break;
-  }  
+  }
 }
 
 
@@ -344,7 +344,7 @@ inline void colvarvalue::apply_constraints()
   case colvarvalue::type_notset:
   default:
     break;
-  }  
+  }
 }
 
 
@@ -362,7 +362,7 @@ inline cvm::real colvarvalue::norm2() const
   case colvarvalue::type_notset:
   default:
     return 0.0;
-  }  
+  }
 }
 
 
@@ -383,7 +383,7 @@ inline colvarvalue colvarvalue::inverse() const
   case colvarvalue::type_notset:
   default:
     undef_op();
-  }  
+  }
   return colvarvalue();
 }
 
@@ -391,11 +391,11 @@ inline colvarvalue colvarvalue::inverse() const
 inline colvarvalue & colvarvalue::operator = (colvarvalue const &x)
 {
   if (this->value_type != type_notset)
-    if (this->value_type != x.value_type) 
+    if (this->value_type != x.value_type)
       error_lside (x.value_type);
 
   this->value_type = x.value_type;
-  
+
   switch (this->value_type) {
   case colvarvalue::type_scalar:
     this->real_value = x.real_value;