diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 339ba95f7d..20654966b2 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -1104,7 +1104,6 @@ void colvar::eigenvector::calc_Jacobian_derivative()
   // gradients of products of 2 quaternion components
   cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
 
-  cvm::atom_pos x_relative;
   cvm::real sum = 0.0;
 
   for (size_t ia = 0; ia < atoms.size(); ia++) {
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 423894bd22..b6f57c16c2 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -1,3 +1,4 @@
+#include <sstream>
 #include "colvarmodule.h"
 #include "colvarparse.h"
 #include "colvarproxy.h"
@@ -320,6 +321,27 @@ void colvarmodule::change_configuration (std::string const &bias_name,
   cvm::decrease_depth();
 }
 
+
+std::string colvarmodule::read_colvar(std::string const &name)
+{
+  cvm::increase_depth();
+  int found = 0;
+  std::stringstream ss;
+  for (std::vector<colvar *>::iterator cvi = colvars.begin();
+       cvi != colvars.end();
+       cvi++) {
+    if ( (*cvi)->name == name ) {
+      ++found;
+      ss << (*cvi)->value();
+    }
+  }
+  if ( found < 1 ) cvm::fatal_error ("Error: colvar not found");
+  if ( found > 1 ) cvm::fatal_error ("Error: duplicate colvar");
+  cvm::decrease_depth();
+  return ss.str();
+}
+
+
 cvm::real colvarmodule::energy_difference (std::string const &bias_name,
                                            std::string const &conf)
 {
@@ -793,6 +815,44 @@ void cvm::read_index_file (char const *filename)
 
 }
 
+void cvm::load_coords_xyz (char const *filename,
+                           std::vector<atom_pos> &pos,
+                           const std::vector<int> &indices)
+{
+  std::ifstream xyz_is (filename);
+  int natoms;
+  char symbol[256];
+  std::string line;
+
+  if ( ! (xyz_is >> natoms) ) {
+    cvm::fatal_error ("Error: cannot parse XYZ file "
+                      + std::string (filename) + ".\n");
+  }
+  // skip comment line
+  std::getline (xyz_is, line);
+  std::getline (xyz_is, line);
+  xyz_is.width (255);
+  std::vector<atom_pos>::iterator pos_i = pos.begin();
+
+  if (pos.size() != natoms) { // Use specified indices
+    int next = 0; // indices are zero-based
+    std::vector<int>::const_iterator index = indices.begin();
+    for ( ; pos_i != pos.end() ; pos_i++, index++) {
+
+      while ( next < *index ) {
+        std::getline (xyz_is, line);
+        next++;
+      }
+      xyz_is >> symbol;
+      xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
+    }
+  } else {          // Use all positions
+    for ( ; pos_i != pos.end() ; pos_i++) {
+      xyz_is >> symbol;
+      xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
+    }
+  }
+}
 
 // static pointers
 std::vector<colvar *>     colvarmodule::colvars;
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 2587f95d88..2146b2cb94 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -2,7 +2,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2013-06-19"
+#define COLVARS_VERSION "2013-10-22"
 #endif
 
 #ifndef COLVARS_DEBUG
@@ -22,6 +22,7 @@
 #include <iostream>
 #include <iomanip>
 #include <string>
+#include <cstring>
 #include <sstream>
 #include <fstream>
 #include <cmath>
@@ -178,6 +179,9 @@ public:
   /// currently works for harmonic (force constant and/or centers)
   void change_configuration (std::string const &bias_name, std::string const &conf);
 
+  /// Read a colvar value
+  std::string read_colvar(std::string const &name);
+
   /// Calculate change in energy from using alt. config. for the given bias -
   /// currently works for harmonic (force constant and/or centers)
   real energy_difference (std::string const &bias_name, std::string const &conf);
@@ -316,14 +320,18 @@ public:
                           double const pdb_field_value = 0.0);
 
   /// \brief Load the coordinates for a group of atoms from a file
-  /// (usually a PDB); the number of atoms in "filename" must match
-  /// the number of elements in "pos"
+  /// (PDB or XYZ)
   static void load_coords (char const *filename,
                            std::vector<atom_pos> &pos,
                            const std::vector<int> &indices,
                            std::string const &pdb_field,
                            double const pdb_field_value = 0.0);
 
+  /// \brief Load the coordinates for a group of atoms from an
+  /// XYZ file
+  static void load_coords_xyz (char const *filename,
+                              std::vector<atom_pos> &pos,
+                              const std::vector<int> &indices);
 
   /// Frequency for collective variables trajectory output
   static size_t cv_traj_freq;
@@ -489,7 +497,19 @@ inline void cvm::load_coords (char const *file_name,
                               std::string const &pdb_field,
                               double const pdb_field_value)
 {
-  proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value);
+  // Differentiate between PDB and XYZ files
+  // for XYZ files, use CVM internal parser
+  // otherwise call proxy function for PDB
+
+  char const *ext = strlen(file_name) > 4 ? file_name + (strlen(file_name) - 4) : file_name;
+  if ( !strncmp(ext, ".xyz", 4) || !strncmp(ext, ".XYZ", 4) ) {
+    if ( pdb_field.size() > 0 ) {
+      cvm::fatal_error ("Error: PDB column may not be specified for XYZ coordinate file.\n");
+    }
+    cvm::load_coords_xyz (file_name, pos, indices);
+  } else {
+    proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value);
+  }
 }
 
 inline void cvm::backup_file (char const *filename)