diff --git a/doc/PDF/colvars-refman-lammps.pdf b/doc/PDF/colvars-refman-lammps.pdf
index a9a1b5624b..7a5817f6d7 100644
Binary files a/doc/PDF/colvars-refman-lammps.pdf and b/doc/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index e046a0bb12..132793ab2e 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -161,11 +161,33 @@ colvar::colvar(std::string const &conf)
     enable(task_scripted);
     cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
 
-    // Only accept scalar scripted colvars
-    // might accept other types when the infrastructure is in place
-    // for derivatives of vectors wrt vectors
-    x.type(colvarvalue::type_scalar);
-    x_reported.type(x);
+    std::string type_str;
+    get_keyval(conf, "scriptedFunctionType", type_str, "scalar");
+
+    x.type(colvarvalue::type_notset);
+    int t;
+    for (t = 0; t < colvarvalue::type_all; t++) {
+      if (type_str == colvarvalue::type_keyword(colvarvalue::Type(t))) {
+        x.type(colvarvalue::Type(t));
+        break;
+      }
+    }
+    if (x.type() == colvarvalue::type_notset) {
+      cvm::error("Could not parse scripted colvar type.");
+      return;
+    }
+    x_reported.type (x.type());
+    cvm::log(std::string("Expecting colvar value of type ")
+      + colvarvalue::type_desc(x.type()));
+
+    if (x.type() == colvarvalue::type_vector) {
+      int size;
+      if (!get_keyval(conf, "scriptedFunctionVectorSize", size)) {
+        cvm::error("Error: no size specified for vector scripted function.");
+        return;
+      }
+      x.vector1d_value.resize(size);
+    }
 
     // Sort array of cvcs based on values of componentExp
     std::vector<cvc *> temp_vec;
@@ -190,15 +212,6 @@ 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;
@@ -1176,25 +1189,28 @@ void colvar::communicate_forces()
     cvm::log("Communicating forces from colvar \""+this->name+"\".\n");
 
   if (tasks[task_scripted]) {
-    std::vector<colvarvalue> func_grads(cvcs.size());
+    std::vector<cvm::matrix2d<cvm::real> > func_grads;
+    for (i = 0; i < cvcs.size(); i++) {
+      func_grads.push_back(cvm::matrix2d<cvm::real> (x.size(),
+                                                     cvcs[i]->value().size()));
+    }
     int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads);
 
-    if (res == COLVARS_NOT_IMPLEMENTED) {
-      cvm::error("Colvar gradient scripts are not implemented.");
-      return;
-    }
     if (res != COLVARS_OK) {
-      cvm::error("Error running colvar gradient script");
+      if (res == COLVARS_NOT_IMPLEMENTED) {
+        cvm::error("Colvar gradient scripts are not implemented.");
+      } else {
+        cvm::error("Error running colvar gradient script");
+      }
       return;
     }
 
     for (i = 0; i < cvcs.size(); i++) {
       cvm::increase_depth();
-      // Force is scalar times colvarvalue (scalar or vector)
-      // Note: this can only handle scalar colvars (scalar values of f)
-      // A non-scalar colvar would need the gradient to be expressed
-      // as an order-2 tensor
-      (cvcs[i])->apply_force(f.real_value * func_grads[i]);
+      // cvc force is colvar force times colvar/cvc Jacobian
+      // (vector-matrix product)
+      (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[i],
+                             cvcs[i]->value().type()));
       cvm::decrease_depth();
     }
   } else if (x.type() == colvarvalue::type_scalar) {
diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp
index 309a089258..78d7341af9 100644
--- a/lib/colvars/colvarbias_alb.cpp
+++ b/lib/colvars/colvarbias_alb.cpp
@@ -1,8 +1,20 @@
+// -*- c++ -*-
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+
 #include "colvarmodule.h"
 #include "colvarbias_alb.h"
 #include "colvarbias.h"
-#include <stdio.h>
-#include <stdlib.h>
+
+#ifdef _MSC_VER
+#if _MSC_VER <= 1700
+#define copysign(A,B) _copysign(A,B)
+double fmax(double A, double B) { return ( A > B ? A : B ); }
+double fmin(double A, double B) { return ( A < B ? A : B ); }
+#endif
+#endif
 
 /* Note about nomenclature. Force constant is called a coupling
  * constant here to emphasize its changing in the code. Outwards,
diff --git a/lib/colvars/colvarbias_alb.h b/lib/colvars/colvarbias_alb.h
index 0b32cb8545..d8c46b499f 100644
--- a/lib/colvars/colvarbias_alb.h
+++ b/lib/colvars/colvarbias_alb.h
@@ -1,3 +1,5 @@
+// -*- c++ -*-
+
 #ifndef COLVARBIAS_ALB_H
 #define COLVARBIAS_ALB_H
 
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 55f49cd21e..2528088926 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -1,6 +1,8 @@
 /// -*- c++ -*-
 
 #include <sstream>
+#include <string.h>
+
 #include "colvarmodule.h"
 #include "colvarparse.h"
 #include "colvarproxy.h"
@@ -255,6 +257,7 @@ int colvarmodule::parse_biases_type(std::string const &conf,
     }
     bias_conf = "";
   }
+  return COLVARS_OK;
 }
 
 
@@ -573,7 +576,7 @@ int colvarmodule::calc() {
   // write trajectory file, if needed
   if (cv_traj_freq && cv_traj_name.size()) {
 
-    if (!cv_traj_os.good()) {
+    if (!cv_traj_os.is_open()) {
       open_traj_file(cv_traj_name);
     }
 
@@ -734,21 +737,10 @@ int colvarmodule::setup_output()
      std::string(""));
 
   if (cv_traj_freq && cv_traj_name.size()) {
-    // open trajectory file
-    if (cv_traj_append) {
-      cvm::log("Appending to colvar trajectory file \""+cv_traj_name+
-                "\".\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);
-    }
-    cv_traj_os.setf(std::ios::scientific, std::ios::floatfield);
+    open_traj_file(cv_traj_name);
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 
+  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
 
 
@@ -928,21 +920,27 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
 
 int colvarmodule::open_traj_file(std::string const &file_name)
 {
+  if (cv_traj_os.is_open()) {
+    return COLVARS_OK;
+  }
+
   // (re)open trajectory file
   if (cv_traj_append) {
     cvm::log("Appending to colvar trajectory file \""+file_name+
               "\".\n");
     cv_traj_os.open(file_name.c_str(), std::ios::app);
   } else {
-    cvm::log("Overwriting colvar trajectory file \""+file_name+
+    cvm::log("Writing to colvar trajectory file \""+file_name+
               "\".\n");
     proxy->backup_file(file_name.c_str());
     cv_traj_os.open(file_name.c_str(), std::ios::out);
   }
-  if (!cv_traj_os.good()) {
+
+  if (!cv_traj_os.is_open()) {
     cvm::error("Error: cannot write to file \""+file_name+"\".\n",
                 FILE_ERROR);
   }
+
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
 
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index edf62f1f88..bcfc6a1d9e 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2014-11-06"
+#define COLVARS_VERSION "2014-11-13"
 #endif
 
 #ifndef COLVARS_DEBUG
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index 9d3c73e109..e9388f68fb 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -170,7 +170,7 @@ public:
 
   virtual int run_colvar_gradient_callback(std::string const &name,
                                          std::vector<const colvarvalue *> const &cvcs,
-                                         std::vector<colvarvalue> &gradient)
+                                         std::vector<cvm::matrix2d<cvm::real> > &gradient)
   { return COLVARS_NOT_IMPLEMENTED; }
 
   // **************** INPUT/OUTPUT ****************
diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp
index a550a09915..0c3a90efbf 100644
--- a/lib/colvars/colvarscript.cpp
+++ b/lib/colvars/colvarscript.cpp
@@ -1,6 +1,9 @@
 // -*- c++ -*-
 
 #include <cstdlib>
+#include <stdlib.h>
+#include <string.h>
+
 #include "colvarscript.h"
 
 
@@ -290,10 +293,6 @@ int colvarscript::proc_bias(int argc, char const *argv[]) {
   // Subcommands for MW ABF
   if (subcmd == "bin") {
     int r = b->current_bin();
-    if (r < 0) {
-      result = "Error: calling current_bin() for bias " + b->name;
-      return COLVARSCRIPT_ERROR;
-    }
     result = cvm::to_str(r);
     return COLVARSCRIPT_OK;
   }
diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp
index 120342e71f..2ec448bd14 100644
--- a/lib/colvars/colvartypes.cpp
+++ b/lib/colvars/colvartypes.cpp
@@ -1,4 +1,7 @@
-/// -*- c++ -*-
+// -*- c++ -*-
+
+#include <stdlib.h>
+#include <string.h>
 
 #include "colvarmodule.h"
 #include "colvartypes.h"
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index 9bed239962..5680fb93e8 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -29,12 +29,9 @@ protected:
 
 public:
 
-  size_t length;
-
   /// Default constructor
   inline vector1d(size_t const n = 0)
   {
-    length = n;
     data.resize(n);
     reset();
   }
@@ -42,7 +39,6 @@ public:
   /// Constructor from C array
   inline vector1d(size_t const n, T const *t)
   {
-    length = n;
     data.resize(n);
     reset();
     size_t i;
@@ -405,6 +401,11 @@ public:
     data.assign(data.size(), T(0.0));
   }
 
+  inline size_t size() const
+  {
+    return data.size();
+  }
+
   /// Default constructor
   inline matrix2d()
     : outer_length(0), inner_length(0)
@@ -557,26 +558,63 @@ public:
   }
 
   /// Matrix multiplication
-  inline friend matrix2d<T> const & operator * (matrix2d<T> const &m1, matrix2d<T> const &m2)
+//   inline friend matrix2d<T> const & operator * (matrix2d<T> const &m1, matrix2d<T> const &m2)
+//   {
+//     matrix2d<T> result(m1.outer_length, m2.inner_length);
+//     if (m1.inner_length != m2.outer_length) {
+//       cvm::error("Error: trying to multiply two matrices of incompatible sizes, "+
+//                  cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
+//                  cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
+//     } else {
+//       size_t i, j, k;
+//       for (i = 0; i < m1.outer_length; i++) {
+//         for (j = 0; j < m2.inner_length; j++) {
+//           for (k = 0; k < m1.inner_length; k++) {
+//             result[i][j] += m1[i][k] * m2[k][j];
+//           }
+//         }
+//       }
+//     }
+//     return result;
+//   }
+
+  /// vector-matrix multiplication
+  inline friend vector1d<T> operator * (vector1d<T> const &v, matrix2d<T> const &m)
   {
-    matrix2d<T> result(m1.outer_length, m2.inner_length);
-    if (m1.inner_length != m2.outer_length) {
-      cvm::error("Error: trying to multiply two matrices of incompatible sizes, "+
-                 cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
-                 cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
+    vector1d<T> result(m.inner_length);
+    if (m.outer_length != v.size()) {
+      cvm::error("Error: trying to multiply a vector and a matrix of incompatible sizes, "+
+                  cvm::to_str(v.size()) + " and " + cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length) + ".\n");
     } else {
-      size_t i, j, k;
-      for (i = 0; i < m1.outer_length; i++) {
-        for (j = 0; j < m2.inner_length; j++) {
-          for (k = 0; k < m1.inner_length; k++) {
-            result[i][j] += m1[i][k] * m2[k][j];
-          }
+      size_t i, k;
+      for (i = 0; i < m.inner_length; i++) {
+        for (k = 0; k < m.outer_length; k++) {
+          result[i] += m[k][i] * v[k];
         }
       }
     }
     return result;
   }
 
+//   /// matrix-vector multiplication (unused for now)
+//   inline friend vector1d<T> const & operator * (matrix2d<T> const &m, vector1d<T> const &v)
+//   {
+//     vector1d<T> result(m.outer_length);
+//     if (m.inner_length != v.size()) {
+//       cvm::error("Error: trying to multiply a matrix and a vector of incompatible sizes, "+
+//                   cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length)
+//                   + " and " + cvm::to_str(v.length) + ".\n");
+//     } else {
+//       size_t i, k;
+//       for (i = 0; i < m.outer_length; i++) {
+//         for (k = 0; k < m.inner_length; k++) {
+//           result[i] += m[i][k] * v[k];
+//         }
+//       }
+//     }
+//     return result;
+//   }
+
   /// Formatted output
   friend std::ostream & operator << (std::ostream &os,
                                      matrix2d<T> const &m)
diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h
index 16fc5737ca..3f74fedb47 100644
--- a/lib/colvars/colvarvalue.h
+++ b/lib/colvars/colvarvalue.h
@@ -16,20 +16,21 @@
 /// treat different data types.  By default, a \link colvarvalue
 /// \endlink variable is a scalar number.  To use it as
 /// another type, declare and initialize it as
-/// \code colvarvalue x(colvarvalue::type_xxx), use \link x.type
+/// \code colvarvalue x(colvarvalue::type_xxx)\endcode, use \link x.type
 /// (colvarvalue::type_xxx) \endlink at a later stage, or if unset,
-//  assign the type with \code x = y; \endcode, provided y is correctly set.
+///  assign the type with \code x = y; \endcode, provided y is correctly set.
 ///
 /// All operators (either unary or binary) on a \link
 /// colvarvalue \endlink object performs one or more checks on the
 /// \link Type \endlink, except when reading from a stream, when there is no way to
 /// detect the \link Type \endlink.  To use  \code is >> x; \endcode x \b MUST
 /// already have a type correcly set up for properly parsing the
-/// stream. No problem of course with the output streams: \code os << x;
+/// stream. No problem of course with the output streams: \code os << x; \endcode
 ///
 /// \em Note \em on \em performance: to avoid type checks in a long array of \link
 /// colvarvalue \endlink objects, use one of the existing "_opt" functions or implement a new one
 
+
 class colvarvalue {
 
 public:
@@ -55,7 +56,9 @@ public:
     /// 4-dimensional vector that is a derivative of a quaternion
     type_quaternionderiv,
     /// vector (arbitrary dimension)
-    type_vector
+    type_vector,
+    /// Needed to iterate through enum
+    type_all
   };
 
   /// Current type of this colvarvalue object
@@ -501,9 +504,6 @@ inline size_t colvarvalue::size() const
 inline colvarvalue::colvarvalue(colvarvalue const &x)
   : value_type(x.type())
 {
-  // reset the value based on the previous type
-  reset();
-
   switch (x.type()) {
   case type_scalar:
     real_value = x.real_value;
@@ -530,9 +530,7 @@ inline colvarvalue::colvarvalue(colvarvalue const &x)
 
 inline colvarvalue::colvarvalue(cvm::vector1d<cvm::real> const &v, Type const &vti)
 {
-  // reset the value based on the previous type
-  reset();
-  if ((vti != type_vector) && (v.size() != num_dimensions(value_type))) {
+  if ((vti != type_vector) && (v.size() != num_dimensions(vti))) {
     cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+
                "\" using a vector of size "+cvm::to_str(v.size())+
                ".\n");
@@ -770,7 +768,7 @@ inline cvm::vector1d<cvm::real> const colvarvalue::as_vector() const
   case colvarvalue::type_scalar:
     {
       cvm::vector1d<cvm::real> v(1);
-      v = real_value;
+      v[0] = real_value;
       return v;
     }
   case colvarvalue::type_3vector: