another day, another colvars update

lib/colvars/colvar.cpp

@@ -188,6 +188,15 @@ colvar::colvar(std::string const &conf)
       sorted_cvc_values.push_back(&(cvcs[j]->value()));
     }
 
+    // // these two are the vector value and the Jacobian matrix of the scripted function, respectively
+    // x_cvc.type(type_vector);
+    // dx_cvc.type(type_matrix); // TODO: not implemented yet
+    // for (j = 0; j < cvcs.size(); j++) {
+    //   x_cvc.add_elem(cvcs[j]->value());
+    //   dx_cvc.add_elem(cvcs[j]->value());
+    // }
+
+
     b_homogeneous = false;
     // Scripted functions are deemed non-periodic
     b_periodic = false;

lib/colvars/colvar.h

@@ -190,6 +190,13 @@ protected:
   /// Value of the colvar
   colvarvalue x;
 
+  // TODO: implement functionality to treat these
+  // /// Vector of individual values from CVCs
+  // colvarvalue x_cvc;
+
+  // /// Jacobian matrix of individual values from CVCs
+  // colvarvalue dx_cvc;
+
   /// Cached reported value (x may be manipulated)
   colvarvalue x_reported;
 

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2014-10-27"
+#define COLVARS_VERSION "2014-10-28"
 #endif
 
 #ifndef COLVARS_DEBUG

lib/colvars/colvartypes.cpp

@@ -4,6 +4,7 @@
 #include "colvartypes.h"
 #include "colvarparse.h"
 
+
 std::string cvm::rvector::to_simple_string() const
 {
   std::ostringstream os;
@@ -13,6 +14,7 @@ std::string cvm::rvector::to_simple_string() const
   return os.str();
 }
 
+
 int cvm::rvector::from_simple_string(std::string const &s)
 {
   std::stringstream stream(s);
@@ -24,6 +26,7 @@ int cvm::rvector::from_simple_string(std::string const &s)
   return COLVARS_OK;
 }
 
+
 std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v)
 {
   std::streamsize const w = os.width();
@@ -277,11 +280,14 @@ void colvarmodule::rotation::diagonalize_matrix(cvm::matrix2d<cvm::real> &S,
 {
   // diagonalize
   int jac_nrot = 0;
-  jacobi(S, S_eigval, S_eigvec, &jac_nrot);
+  jacobi(S, S_eigval.c_array(), S_eigvec, &jac_nrot);
-  eigsrt(S_eigval, S_eigvec);
+  eigsrt(S_eigval.c_array(), S_eigvec);
   // jacobi saves eigenvectors by columns
   transpose(S_eigvec);
 
+  S_eigval.update_from_c_array();
+  S_eigval.delete_c_array();
+
   // normalize eigenvectors
   for (size_t ie = 0; ie < 4; ie++) {
     cvm::real norm2 = 0.0;
@@ -301,7 +307,7 @@ void colvarmodule::rotation::calc_optimal_rotation
 {
   cvm::matrix2d<cvm::real> S(4, 4, 0.0);
   cvm::matrix2d<cvm::real> S_backup(4, 4, 0.0);
-  cvm::vector1d<cvm::real> S_eigval(4, 0.0);
+  cvm::vector1d<cvm::real> S_eigval(4);
   cvm::matrix2d<cvm::real> S_eigvec(4, 4, 0.0);
 
   build_matrix(pos1, pos2, S);
@@ -470,7 +476,7 @@ void colvarmodule::rotation::calc_optimal_rotation
       if (b_debug_gradients) {
 
         cvm::matrix2d<cvm::real> S_new(4, 4, 0.0);
-        cvm::vector1d<cvm::real> S_new_eigval(4, 0.0);
+        cvm::vector1d<cvm::real> S_new_eigval(4);
         cvm::matrix2d<cvm::real> S_new_eigvec(4, 4, 0.0);
 
         // make an infitesimal move along each cartesian coordinate of
@@ -519,6 +525,7 @@ void colvarmodule::rotation::calc_optimal_rotation
   h=a[k][l];                \
   a[i][j]=g-s*(h+g*tau);    \
   a[k][l]=h+s*(g-h*tau);
+
 #define n 4
 
 void jacobi(cvm::real **a, cvm::real d[], cvm::real **v, int *nrot)

lib/colvars/colvartypes.h

@@ -18,145 +18,158 @@
 /// \brief Arbitrary size array (one dimensions) suitable for linear
 /// algebra operations (i.e. for floating point numbers it can be used
 /// with library functions)
-template <class T> class colvarmodule::vector1d
+template <class T> class colvarmodule::vector1d : public std::vector<T>
 {
 protected:
 
-  /// Underlying C-array
+  // member used to exchange with functions expecting a C-style array
   T *array;
   size_t length;
 
 public:
 
-  /// Allocation routine, used by all constructors
-  inline void alloc() {
-    if (length > 0) {
-      array = new T [length];
-    } else {
-      array = NULL;
-    }
-  }
-
-  /// Deallocation routine
-  inline void dealloc() {
-    if (array != NULL) {
-      delete [] array;
-    }
-  }
-
-  /// Length of the array
-  inline size_t size() const
-  {
-    return length;
-  }
-
   /// Default constructor
-  inline vector1d(size_t const n = 0, T const &t = T())
+  inline vector1d(size_t const n = 0) : array (NULL)
   {
     length = n;
-    this->alloc();
+    this->resize(n);
     reset();
   }
 
-  /// Constructor from a 1-d C array
+  /// Create a copy (unless it exists already), and return it
-  inline vector1d(size_t const n, T const *v)
+  inline T * c_array()
   {
-    length = n;
+    if (array == NULL) {
-    this->alloc();
+      array = new T[this->size()];
-    for (size_t i = 0; i < length; i++) {
+      for (size_t i = 0; i < this->size(); i++) {
-      array[i] = v[i];
+        array[i] = (*this)[i];
+      }
+    }
+    return array;
+  }
+
+  inline void update_from_c_array()
+  {
+    for (size_t i = 0; i < this->size(); i++) {
+      (*this)[i] = array[i];
     }
   }
 
-  /// Copy constructor
+  inline void delete_c_array()
-  inline vector1d(vector1d<T> const &v)
   {
-    length = v.length;
+    delete [] array;
-    this->alloc();
+    array = NULL;
-    for (size_t i = 0; i < length; i++) {
+  }
-      array[i] = v.array[i];
+
-    }
+  inline ~vector1d()
+  {
+    this->delete_c_array();
   }
 
   /// Set all elements to zero
   inline void reset()
   {
-    for (size_t i = 0; i < length; i++) {
+    this->assign(this->size(), T(0.0));
-      array[i] = T(0.0);
-    }
   }
 
-  /// Assignment
+  inline static void check_sizes(size_t const n1, size_t const n2)
-  inline vector1d<T> & operator = (vector1d<T> const &v)
   {
-    if (length != v.length) {
+    if (n1 != n2) {
-      this->dealloc();
+      cvm::error("Error: trying to perform an operation between vectors of different sizes, "+
-      length = v.length;
+                 cvm::to_str(n1)+" and "+cvm::to_str(n2)+".\n");
-      this->alloc();
     }
-    for (size_t i = 0; i < length; i++) {
+  }
-      this->array[i] = v.array[i];
+
+  inline void operator += (vector1d<T> const &v)
+  {
+    check_sizes(this->size(), v.size());
+    for (size_t i = 0; i < this->size(); i++) {
+      (*this)[i] += v[i];
     }
-    return *this;
   }
 
-  /// Destructor
+  inline void operator -= (vector1d<T> const &v)
-  inline ~vector1d() {
+  {
-    this->dealloc();
+    check_sizes(this->size(), v.size());
+    for (size_t i = 0; i < this->size(); i++) {
+      (*this)[i] -= v[i];
+    }
   }
 
-  /// Return the 1-d C array
+  inline void operator *= (cvm::real const &a)
-  inline operator T *() { return array; }
+  {
-
+    for (size_t i = 0; i < this->size(); i++) {
-  inline T & operator [] (int const &i) {
+      (*this)[i] *= a;
-    return array[i];
+    }
   }
 
-  inline T const operator [] (int const &i) const {
+  inline void operator /= (cvm::real const &a)
-    return array[i];
+  {
+    for (size_t i = 0; i < this->size(); i++) {
+      (*this)[i] /= a;
+    }
+  }
+
+  inline friend vector1d<T> operator + (vector1d<T> const &v1, vector1d<T> const &v2)
+  {
+    check_sizes(v1.size(), v2.size());
+    vector1d<T> result(v1.size());
+    for (size_t i = 0; i < v1.size(); i++) {
+      result[i] = v1[i] + v2[i];
+    }
+    return result;
+  }
+
+  inline friend vector1d<T> operator - (vector1d<T> const &v1, vector1d<T> const &v2)
+  {
+    check_sizes(v1.size(), v2.size());
+    vector1d<T> result(v1.size());
+    for (size_t i = 0; i < v1.size(); i++) {
+      result[i] = v1[i] - v2[i];
+    }
+    return result;
+  }
+
+  inline friend vector1d<T> operator * (vector1d<T> const &v, cvm::real const &a)
+  {
+    vector1d<T> result(v.size());
+    for (size_t i = 0; i < v.size(); i++) {
+      result[i] = v[i] * a;
+    }
+    return result;
+  }
+
+  inline friend vector1d<T> operator * (cvm::real const &a, vector1d<T> const &v)
+  {
+    return v * a;
+  }
+
+  inline friend vector1d<T> operator / (vector1d<T> const &v, cvm::real const &a)
+  {
+    vector1d<T> result(v.size());
+    for (size_t i = 0; i < v.size(); i++) {
+      result[i] = v[i] / a;
+    }
+    return result;
   }
 
   /// Inner product
-  inline friend T operator * (vector1d<T> const &v1,
+  inline friend T operator * (vector1d<T> const &v1, vector1d<T> const &v2)
-                              vector1d<T> const &v2)
   {
+    check_sizes(v1.size(), v2.size());
     T prod(0.0);
-    size_t const min_length = (v1.length > v2.length) ? v1.length : v2.length;
+    for (size_t i = 0; i < v1.size(); i++) {
-    for (size_t i = 0; i < min_length; i++) {
+      prod += v1[i] * v2[i];
-      prod += v1.array[i] * v2.array[i];
     }
     return prod;
   }
 
-  /// Slicing
-  inline vector1d<T> const slice(size_t const i1, size_t const i2) const
-  {
-    if ((i2 < i1) || (i1 < 0) || (i2 >= length)) {
-      cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
-    }
-    vector1d<T> result(i2 - i1, T());
-    for (size_t i = 0; i < (i2 - i1); i++) {
-      result[i] = (*this)[i1+i];
-    }
-  }
-
-  /// Assign to a slice
-  inline void sliceassign(size_t const i1, size_t const i2, vector1d<T> const &v)
-  {
-    if ((i2 < i1) || (i1 < 0) || (i2 >= length)) {
-      cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
-    }
-    for (size_t i = 0; i < (i2 - i1); i++) {
-      (*this)[i1+i] = v[i];
-    }
-  }
-
   /// Squared norm
   inline cvm::real norm2() const
   {
     cvm::real result = 0.0;
-    for (size_t i = 0; i < length; i++) {
+    for (size_t i = 0; i < this->size(); i++) {
-      result += array[i] * array[i];
+      result += (*this)[i] * (*this)[i];
     }
     return result;
   }
@@ -166,23 +179,75 @@ public:
     return std::sqrt(this->norm2());
   }
 
+  /// Slicing
+  inline vector1d<T> const slice(size_t const i1, size_t const i2) const
+  {
+    if ((i2 < i1) || (i1 < 0) || (i2 >= this->size())) {
+      cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
+    }
+    vector1d<T> result(i2 - i1);
+    for (size_t i = 0; i < (i2 - i1); i++) {
+      result[i] = (*this)[i1+i];
+    }
+    return result;
+  }
+
+  /// Assign a vector to a slice of this vector
+  inline void sliceassign(size_t const i1, size_t const i2, vector1d<T> const &v)
+  {
+    if ((i2 < i1) || (i1 < 0) || (i2 >= this->size())) {
+      cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
+    }
+    for (size_t i = 0; i < (i2 - i1); i++) {
+      (*this)[i1+i] = v[i];
+    }
+  }
+
   /// Formatted output
-  friend std::ostream & operator << (std::ostream &os,
+
-                                     vector1d<T> const &v)
+  inline friend std::ostream & operator << (std::ostream &os, cvm::vector1d<T> const &v)
   {
     std::streamsize const w = os.width();
     std::streamsize const p = os.precision();
 
     os << "( ";
-    for (size_t i = 0; i < v.length-1; i++) {
+    for (size_t i = 0; i < v.size()-1; i++) {
       os.width(w); os.precision(p);
-      os << v.array[i] << " , ";
+      os << v[i] << " , ";
     }
     os.width(w); os.precision(p);
-    os << v.array[v.length-1] << " )";
+    os << v[v.size()-1] << " )";
     return os;
   }
 
+
+  inline std::string to_simple_string() const
+  {
+    if (this->size() == 0) return std::string("");
+    std::ostringstream os;
+    os.setf(std::ios::scientific, std::ios::floatfield);
+    os.precision(cvm::cv_prec);
+    os << (*this)[0];
+    for (size_t i = 1; i < this->size(); i++) {
+      os << " " << (*this)[i];
+    }
+    return os.str();
+  }
+
+
+  inline int from_simple_string(std::string const &s)
+  {
+    std::stringstream stream(s);
+    size_t i = 0;
+    while ((stream >> (*this)[i]) && (i < this->size())) {
+      i++;
+    }
+    if (i < this->size()) {
+      return COLVARS_ERROR;
+    }
+    return COLVARS_OK;
+  }
+
 };
 
 
@@ -242,10 +307,11 @@ public:
 
   inline cvm::vector1d<cvm::real> const as_vector() const
   {
-    cvm::vector1d<cvm::real> result(3, 0.0);
+    cvm::vector1d<cvm::real> result(3);
     result[0] = this->x;
     result[1] = this->y;
     result[2] = this->z;
+    return result;
   }
 
   inline cvm::rvector & operator = (cvm::real const &v)
@@ -762,11 +828,12 @@ public:
 
   inline cvm::vector1d<cvm::real> const as_vector() const
   {
-    cvm::vector1d<cvm::real> result(4, 0.0);
+    cvm::vector1d<cvm::real> result(4);
     result[0] = q0;
     result[1] = q1;
     result[2] = q2;
     result[3] = q3;
+    return result;
   }
 
   /// Square norm of the quaternion
@@ -1014,7 +1081,7 @@ public:
   {
     dS_1.resize(n, cvm::matrix2d<cvm::rvector>(4, 4, cvm::rvector(0.0, 0.0, 0.0)));
     dL0_1.resize(n, cvm::rvector(0.0, 0.0, 0.0));
-    dQ0_1.resize(n, cvm::vector1d<cvm::rvector>(4, cvm::rvector(0.0, 0.0, 0.0)));
+    dQ0_1.resize(n, cvm::vector1d<cvm::rvector>(4));
   }
 
   /// Allocate space for the derivatives of the rotation
@@ -1022,7 +1089,7 @@ public:
   {
     dS_2.resize(n, cvm::matrix2d<cvm::rvector>(4, 4, cvm::rvector(0.0, 0.0, 0.0)));
     dL0_2.resize(n, cvm::rvector(0.0, 0.0, 0.0));
-    dQ0_2.resize(n, cvm::vector1d<cvm::rvector>(4, cvm::rvector(0.0, 0.0, 0.0)));
+    dQ0_2.resize(n, cvm::vector1d<cvm::rvector>(4));
   }
 
   /// \brief Calculate the optimal rotation and store the

lib/colvars/colvarvalue.cpp

@@ -1,9 +1,606 @@
 /// -*- c++ -*-
 
 #include <vector>
+#include <sstream>
+#include <iostream>
 
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 
 
 
+void colvarvalue::add_elem(colvarvalue const &x)
+{
+  if (this->value_type != type_vector) {
+    cvm::error("Error: trying to set an element for a variable that is not set to be a vector.\n");
+    return;
+  }
+  size_t const n = vector1d_value.size();
+  size_t const nd = num_dimensions(x.value_type);
+  elem_types.push_back(x.value_type);
+  elem_indices.push_back(n);
+  elem_sizes.push_back(nd);
+  vector1d_value.resize(n + nd);
+  set_elem(n, x);
+}
+
+
+colvarvalue const colvarvalue::get_elem(int const i_begin, int const i_end, Type const vt) const
+{
+  if (vector1d_value.size() > 0) {
+    cvm::vector1d<cvm::real> const v(vector1d_value.slice(i_begin, i_end));
+    return colvarvalue(v, vt);
+  } else {
+    cvm::error("Error: trying to get an element from a variable that is not a vector.\n");
+    return colvarvalue(type_notset);
+  }
+}
+
+
+void colvarvalue::set_elem(int const i_begin, int const i_end, colvarvalue const &x)
+{
+  if (vector1d_value.size() > 0) {
+    vector1d_value.sliceassign(i_begin, i_end, x.as_vector());
+  } else {
+    cvm::error("Error: trying to set an element for a variable that is not a vector.\n");
+  }
+}
+
+
+colvarvalue const colvarvalue::get_elem(int const icv) const
+{
+  if (elem_types.size() > 0) {
+    return get_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv],
+                    elem_types[icv]);
+  } else {
+    cvm::error("Error: trying to get a colvarvalue element from a vector colvarvalue that was initialized as a plain array.\n");
+    return colvarvalue(type_notset);
+  }
+}
+
+
+void colvarvalue::set_elem(int const icv, colvarvalue const &x)
+{
+  if (elem_types.size() > 0) {
+    check_types_assign(elem_types[icv], x.value_type);
+    set_elem(elem_indices[icv], elem_indices[icv] + elem_sizes[icv], x);
+  } else {
+    cvm::error("Error: trying to set a colvarvalue element for a colvarvalue that was initialized as a plain array.\n");
+  }
+}
+
+
+colvarvalue colvarvalue::inverse() const
+{
+  switch (value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue(1.0/real_value);
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return colvarvalue(cvm::rvector(1.0/rvector_value.x,
+                                    1.0/rvector_value.y,
+                                    1.0/rvector_value.z));
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0,
+                                       1.0/quaternion_value.q1,
+                                       1.0/quaternion_value.q2,
+                                       1.0/quaternion_value.q3));
+    break;
+  case colvarvalue::type_vector:
+    {
+      cvm::vector1d<cvm::real> result(vector1d_value);
+      if (elem_types.size() > 0) {
+        // if we have information about non-scalar types, use it
+        for (size_t i = 0; i < elem_types.size(); i++) {
+          result.sliceassign(elem_indices[i], elem_indices[i]+elem_sizes[i],
+                             cvm::vector1d<cvm::real>((this->get_elem(i)).inverse()));
+        }
+      } else {
+        for (size_t i = 0; i < result.size(); i++) {
+          if (result[i] != 0.0) {
+            result = 1.0/result[i];
+          }
+        }
+      }
+      return colvarvalue(result, type_vector);
+    }
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+    break;
+  }
+  return colvarvalue();
+}
+
+
+// binary operations between two colvarvalues
+
+colvarvalue operator + (colvarvalue const &x1,
+                        colvarvalue const &x2)
+{
+  colvarvalue::check_types(x1, x2);
+
+  switch (x1.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue(x1.real_value + x2.real_value);
+  case colvarvalue::type_3vector:
+    return colvarvalue(x1.rvector_value + x2.rvector_value);
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return colvarvalue(x1.rvector_value + x2.rvector_value,
+                       colvarvalue::type_unit3vector);
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return colvarvalue(x1.quaternion_value + x2.quaternion_value);
+  case colvarvalue::type_vector:
+    return colvarvalue(x1.vector1d_value + x2.vector1d_value, colvarvalue::type_vector);
+  case colvarvalue::type_notset:
+  default:
+    x1.undef_op();
+    return colvarvalue(colvarvalue::type_notset);
+  };
+}
+
+
+colvarvalue operator - (colvarvalue const &x1,
+                        colvarvalue const &x2)
+{
+  colvarvalue::check_types(x1, x2);
+
+  switch (x1.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue(x1.real_value - x2.real_value);
+  case colvarvalue::type_3vector:
+    return colvarvalue(x1.rvector_value - x2.rvector_value);
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return colvarvalue(x1.rvector_value - x2.rvector_value,
+                       colvarvalue::type_unit3vector);
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return colvarvalue(x1.quaternion_value - x2.quaternion_value);
+  case colvarvalue::type_vector:
+    return colvarvalue(x1.vector1d_value - x2.vector1d_value, colvarvalue::type_vector);
+  case colvarvalue::type_notset:
+  default:
+    x1.undef_op();
+    return colvarvalue(colvarvalue::type_notset);
+  };
+}
+
+
+// binary operations with real numbers
+
+colvarvalue operator * (cvm::real const &a,
+                        colvarvalue const &x)
+{
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue(a * x.real_value);
+  case colvarvalue::type_3vector:
+    return colvarvalue(a * x.rvector_value);
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return colvarvalue(a * x.rvector_value,
+                       colvarvalue::type_unit3vector);
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return colvarvalue(a * x.quaternion_value);
+  case colvarvalue::type_vector:
+    return colvarvalue(x.vector1d_value * a, colvarvalue::type_vector);
+  case colvarvalue::type_notset:
+  default:
+    x.undef_op();
+    return colvarvalue(colvarvalue::type_notset);
+  }
+}
+
+
+colvarvalue operator * (colvarvalue const &x,
+                        cvm::real const &a)
+{
+  return a * x;
+}
+
+
+colvarvalue operator / (colvarvalue const &x,
+                        cvm::real const &a)
+{
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    return colvarvalue(x.real_value / a);
+  case colvarvalue::type_3vector:
+    return colvarvalue(x.rvector_value / a);
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return colvarvalue(x.rvector_value / a,
+                       colvarvalue::type_unit3vector);
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return colvarvalue(x.quaternion_value / a);
+  case colvarvalue::type_vector:
+    return colvarvalue(x.vector1d_value / a, colvarvalue::type_vector);
+  case colvarvalue::type_notset:
+  default:
+    x.undef_op();
+    return colvarvalue(colvarvalue::type_notset);
+  }
+}
+
+
+// inner product between two colvarvalues
+
+cvm::real operator * (colvarvalue const &x1,
+                      colvarvalue const &x2)
+{
+  colvarvalue::check_types(x1, x2);
+
+  switch (x1.value_type) {
+  case colvarvalue::type_scalar:
+    return (x1.real_value * x2.real_value);
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return (x1.rvector_value * x2.rvector_value);
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    // the "*" product is the quaternion product, here the inner
+    // member function is used instead
+    return (x1.quaternion_value.inner(x2.quaternion_value));
+  case colvarvalue::type_vector:
+    return (x1.vector1d_value * x2.vector1d_value);
+  case colvarvalue::type_notset:
+  default:
+    x1.undef_op();
+    return 0.0;
+  };
+}
+
+
+colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const
+{
+  colvarvalue::check_types(*this, x2);
+
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    return 2.0 * (this->real_value - x2.real_value);
+  case colvarvalue::type_3vector:
+    return 2.0 * (this->rvector_value - x2.rvector_value);
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    {
+      cvm::rvector const &v1 = this->rvector_value;
+      cvm::rvector const &v2 = x2.rvector_value;
+      cvm::real const cos_t = v1 * v2;
+      cvm::real const sin_t = std::sqrt(1.0 - cos_t*cos_t);
+      return colvarvalue( 2.0 * sin_t *
+                          cvm::rvector((-1.0) * sin_t * v2.x +
+                                       cos_t/sin_t * (v1.x - cos_t*v2.x),
+                                       (-1.0) * sin_t * v2.y +
+                                       cos_t/sin_t * (v1.y - cos_t*v2.y),
+                                       (-1.0) * sin_t * v2.z +
+                                       cos_t/sin_t * (v1.z - cos_t*v2.z)
+                                       ),
+                          colvarvalue::type_unit3vector );
+    }
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return this->quaternion_value.dist2_grad(x2.quaternion_value);
+  case colvarvalue::type_vector:
+    return colvarvalue(2.0 * (this->vector1d_value - x2.vector1d_value), colvarvalue::type_vector);
+    break;
+  case colvarvalue::type_notset:
+  default:
+    this->undef_op();
+    return colvarvalue(colvarvalue::type_notset);
+  };
+}
+
+
+std::string colvarvalue::to_simple_string() const
+{
+  switch (type()) {
+  case colvarvalue::type_scalar:
+    return cvm::to_str(real_value, 0, cvm::cv_prec);
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return rvector_value.to_simple_string();
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return quaternion_value.to_simple_string();
+    break;
+  case colvarvalue::type_vector:
+    return vector1d_value.to_simple_string();
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+    break;
+  }
+  return std::string();
+}
+
+
+int colvarvalue::from_simple_string(std::string const &s)
+{
+  switch (type()) {
+  case colvarvalue::type_scalar:
+    return ((std::istringstream(s) >> real_value)
+            ? COLVARS_OK : COLVARS_ERROR);
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return rvector_value.from_simple_string(s);
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return quaternion_value.from_simple_string(s);
+    break;
+  case colvarvalue::type_vector:
+    return vector1d_value.from_simple_string(s);
+    break;
+  case colvarvalue::type_notset:
+  default:
+    undef_op();
+    break;
+  }
+  return COLVARS_ERROR;
+}
+
+std::ostream & operator << (std::ostream &os, colvarvalue const &x)
+{
+  switch (x.type()) {
+  case colvarvalue::type_scalar:
+    os << x.real_value;
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    os << x.rvector_value;
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    os << x.quaternion_value;
+    break;
+  case colvarvalue::type_notset:
+    os << "not set"; break;
+  }
+  return os;
+}
+
+
+std::ostream & operator << (std::ostream &os, std::vector<colvarvalue> const &v)
+{
+  for (size_t i = 0; i < v.size(); i++) {
+    os << v[i];
+  }
+  return os;
+}
+
+
+std::istream & operator >> (std::istream &is, colvarvalue &x)
+{
+  if (x.type() == colvarvalue::type_notset) {
+    cvm::error("Trying to read from a stream a colvarvalue, "
+               "which has not yet been assigned a data type.\n");
+    return is;
+  }
+
+  switch (x.type()) {
+  case colvarvalue::type_scalar:
+    is >> x.real_value;
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    is >> x.rvector_value;
+    break;
+  case colvarvalue::type_unit3vector:
+    is >> x.rvector_value;
+    x.apply_constraints();
+    break;
+  case colvarvalue::type_quaternion:
+    is >> x.quaternion_value;
+    x.apply_constraints();
+    break;
+  case colvarvalue::type_quaternionderiv:
+    is >> x.quaternion_value;
+    break;
+  default:
+    x.undef_op();
+  }
+  return is;
+}
+
+
+size_t colvarvalue::output_width(size_t const &real_width) const
+{
+  switch (this->value_type) {
+  case colvarvalue::type_scalar:
+    return real_width;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    return cvm::rvector::output_width(real_width);
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    return cvm::quaternion::output_width(real_width);
+  case colvarvalue::type_notset:
+  default:
+    return 0;
+  }
+}
+
+
+void colvarvalue::inner_opt(colvarvalue                        const &x,
+                            std::vector<colvarvalue>::iterator       &xv,
+                            std::vector<colvarvalue>::iterator const &xv_end,
+                            std::vector<cvm::real>::iterator         &result)
+{
+  // doing type check only once, here
+  colvarvalue::check_types(x, *xv);
+
+  std::vector<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator    &ii = result;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->real_value * x.real_value;
+    }
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    while (xvi != xv_end) {
+      *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value);
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+
+void colvarvalue::inner_opt(colvarvalue const                      &x,
+                            std::list<colvarvalue>::iterator       &xv,
+                            std::list<colvarvalue>::iterator const &xv_end,
+                            std::vector<cvm::real>::iterator       &result)
+{
+  // doing type check only once, here
+  colvarvalue::check_types(x, *xv);
+
+  std::list<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator  &ii = result;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->real_value * x.real_value;
+    }
+    break;
+  case colvarvalue::type_3vector:
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    while (xvi != xv_end) {
+      *(ii++) += (xvi++)->rvector_value * x.rvector_value;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    while (xvi != xv_end) {
+      *(ii++) += ((xvi++)->quaternion_value).cosine(x.quaternion_value);
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+
+void colvarvalue::p2leg_opt(colvarvalue const                        &x,
+                            std::vector<colvarvalue>::iterator       &xv,
+                            std::vector<colvarvalue>::iterator const &xv_end,
+                            std::vector<cvm::real>::iterator         &result)
+{
+  // doing type check only once, here
+  colvarvalue::check_types(x, *xv);
+
+  std::vector<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator    &ii = result;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    cvm::error("Error: cannot calculate Legendre polynomials "
+               "for scalar variables.\n");
+    return;
+    break;
+  case colvarvalue::type_3vector:
+    while (xvi != xv_end) {
+      cvm::real const cosine =
+        ((xvi)->rvector_value * x.rvector_value) /
+        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
+      xvi++;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value);
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+
+void colvarvalue::p2leg_opt(colvarvalue const                        &x,
+                            std::list<colvarvalue>::iterator         &xv,
+                            std::list<colvarvalue>::iterator const   &xv_end,
+                            std::vector<cvm::real>::iterator         &result)
+{
+  // doing type check only once, here
+  colvarvalue::check_types(x, *xv);
+
+  std::list<colvarvalue>::iterator &xvi = xv;
+  std::vector<cvm::real>::iterator  &ii = result;
+
+  switch (x.value_type) {
+  case colvarvalue::type_scalar:
+    cvm::error("Error: cannot calculate Legendre polynomials "
+               "for scalar variables.\n");
+    break;
+  case colvarvalue::type_3vector:
+    while (xvi != xv_end) {
+      cvm::real const cosine =
+        ((xvi)->rvector_value * x.rvector_value) /
+        ((xvi)->rvector_value.norm() * x.rvector_value.norm());
+      xvi++;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_unit3vector:
+  case colvarvalue::type_unit3vectorderiv:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value;
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  case colvarvalue::type_quaternion:
+  case colvarvalue::type_quaternionderiv:
+    while (xvi != xv_end) {
+      cvm::real const cosine = (xvi++)->quaternion_value.cosine(x.quaternion_value);
+      *(ii++) += 1.5*cosine*cosine - 0.5;
+    }
+    break;
+  default:
+    x.undef_op();
+  };
+}
+
+
+