diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf
index 2d28758819..7f52367c82 100644
Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ
diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp
index 5a4e8b617a..c5b6d906c1 100644
--- a/lib/colvars/colvar.cpp
+++ b/lib/colvars/colvar.cpp
@@ -7,6 +7,10 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
+#include <list>
+#include <vector>
+#include <algorithm>
+
 #include "colvarmodule.h"
 #include "colvarvalue.h"
 #include "colvarparse.h"
@@ -14,15 +18,6 @@
 #include "colvarcomp.h"
 #include "colvarscript.h"
 
-// used in build_atom_list()
-#include <algorithm>
-
-
-/// Compare two cvcs using their names
-/// Used to sort CVC array in scripted coordinates
-bool compare(colvar::cvc *i, colvar::cvc *j) {
-  return i->name < j->name;
-}
 
 
 colvar::colvar()
@@ -34,6 +29,15 @@ colvar::colvar()
 }
 
 
+namespace {
+  /// Compare two cvcs using their names
+  /// Used to sort CVC array in scripted coordinates
+  bool compare(colvar::cvc *i, colvar::cvc *j)
+  {
+    return i->name < j->name;
+  }
+}
+
 int colvar::init(std::string const &conf)
 {
   cvm::log("Initializing a new collective variable.\n");
@@ -143,6 +147,9 @@ int colvar::init(std::string const &conf)
     x.type(cvc_value);
     x_reported.type(cvc_value);
   }
+
+  set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
+
   // If using scripted biases, any colvar may receive bias forces
   // and will need its gradient
   if (cvm::scripted_forces()) {
@@ -198,6 +205,7 @@ int colvar::init(std::string const &conf)
   if (is_enabled(f_cv_homogeneous) && cvcs[0]->b_periodic) { // TODO make this a CVC feature
     bool b_periodic = true;
     period = cvcs[0]->period;
+    wrap_center = cvcs[0]->wrap_center;
     for (i = 1; i < cvcs.size(); i++) {
       if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
         b_periodic = false;
@@ -211,6 +219,14 @@ int colvar::init(std::string const &conf)
     set_enabled(f_cv_periodic, b_periodic);
   }
 
+  // Allow scripted/custom functions to be defined as periodic
+  if ( (is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && is_enabled(f_cv_scalar) ) {
+    if (get_keyval(conf, "period", period, 0.)) {
+      set_enabled(f_cv_periodic, true);
+      get_keyval(conf, "wrapAround", wrap_center, 0.);
+    }
+  }
+
   // check that cvcs are compatible
 
   for (i = 0; i < cvcs.size(); i++) {
@@ -443,8 +459,6 @@ int colvar::init_grid_parameters(std::string const &conf)
   upper_boundary.type(value());
   upper_wall.type(value());
 
-  set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
-
   if (is_enabled(f_cv_scalar)) {
 
     if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
@@ -1503,7 +1517,7 @@ cvm::real colvar::update_forces_energy()
     vr  += (0.5 * dt) * f_ext / ext_mass;
     xr  += dt * vr;
     xr.apply_constraints();
-    if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
+    this->wrap(xr);
   }
 
   // Now adding the force on the actual colvar (for those biases that
@@ -1714,9 +1728,18 @@ colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
 
 void colvar::wrap(colvarvalue &x) const
 {
-  if (is_enabled(f_cv_homogeneous)) {
-    (cvcs[0])->wrap(x);
+  if ( !is_enabled(f_cv_periodic) ) {
+    return;
   }
+
+  if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
+    // Scripted functions do their own wrapping, as cvcs might not be periodic
+    cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
+    x.real_value -= shift * period;
+  } else {
+    cvcs[0]->wrap(x);
+  }
+
   return;
 }
 
@@ -2244,7 +2267,7 @@ void colvar::calc_runave()
     runave.type(value().type());
     runave.reset();
 
-    // first-step operations
+    // first-step operationsf
 
     if (cvm::debug())
       cvm::log("Colvar \""+this->name+
diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h
index 32c329460d..271231fb2a 100644
--- a/lib/colvars/colvar.h
+++ b/lib/colvars/colvar.h
@@ -216,6 +216,7 @@ public:
 
   /// Period, if this variable is periodic
   cvm::real period;
+  cvm::real wrap_center;
 
 
   /// \brief Expand the boundaries of multiples of width, to keep the
diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 1be6f42e59..c8dcfc0766 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -7,6 +7,10 @@
 // If you wish to distribute your changes, please submit them to the
 // Colvars repository at GitHub.
 
+#include <list>
+#include <vector>
+#include <algorithm>
+
 #include "colvarmodule.h"
 #include "colvarproxy.h"
 #include "colvarparse.h"
@@ -436,7 +440,7 @@ int cvm::atom_group::parse(std::string const &group_conf)
   }
 
   bool b_print_atom_ids = false;
-  get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false, colvarparse::parse_silent);
+  get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false);
 
   // Calculate all required properties (such as total mass)
   setup();
@@ -715,13 +719,10 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
                      "if provided, must be non-zero.\n", INPUT_ERROR);
           return COLVARS_ERROR;
         }
-      } else {
-        // if not, rely on existing atom indices for the group
-        group_for_fit->create_sorted_ids();
-        ref_pos.resize(group_for_fit->size());
       }
 
-      cvm::load_coords(ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids,
+      ref_pos.resize(group_for_fit->size());
+      cvm::load_coords(ref_pos_file.c_str(), &ref_pos, group_for_fit,
                        ref_pos_col, ref_pos_col_value);
     }
 
@@ -789,33 +790,39 @@ void cvm::atom_group::do_feature_side_effects(int id)
 }
 
 
-int cvm::atom_group::create_sorted_ids(void)
+int cvm::atom_group::create_sorted_ids()
 {
   // Only do the work if the vector is not yet populated
-  if (sorted_ids.size())
+  if (sorted_atoms_ids.size())
     return COLVARS_OK;
 
-  std::list<int> temp_id_list;
-  for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
-    temp_id_list.push_back(ai->id);
+  // Sort the internal IDs
+  std::list<int> sorted_atoms_ids_list;
+  for (size_t i = 0; i < this->size(); i++) {
+    sorted_atoms_ids_list.push_back(atoms_ids[i]);
   }
-  temp_id_list.sort();
-  temp_id_list.unique();
-  if (temp_id_list.size() != this->size()) {
-    cvm::error("Error: duplicate atom IDs in atom group? (found " +
-               cvm::to_str(temp_id_list.size()) +
-               " unique atom IDs instead of" +
-               cvm::to_str(this->size()) + ").\n");
-    return COLVARS_ERROR;
+  sorted_atoms_ids_list.sort();
+  sorted_atoms_ids_list.unique();
+  if (sorted_atoms_ids_list.size() != this->size()) {
+    return cvm::error("Error: duplicate atom IDs in atom group? (found " +
+                      cvm::to_str(sorted_atoms_ids_list.size()) +
+                      " unique atom IDs instead of " +
+                      cvm::to_str(this->size()) + ").\n", BUG_ERROR);
   }
-  sorted_ids = std::vector<int> (temp_id_list.size());
-  unsigned int id_i = 0;
-  std::list<int>::iterator li;
-  for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
-    sorted_ids[id_i] = *li;
-    id_i++;
+
+  // Compute map between sorted and unsorted elements
+  sorted_atoms_ids.resize(this->size());
+  sorted_atoms_ids_map.resize(this->size());
+  std::list<int>::iterator lsii = sorted_atoms_ids_list.begin();
+  size_t ii = 0;
+  for ( ; ii < this->size(); lsii++, ii++) {
+    sorted_atoms_ids[ii] = *lsii;
+    size_t const pos = std::find(atoms_ids.begin(), atoms_ids.end(), *lsii) -
+      atoms_ids.begin();
+    sorted_atoms_ids_map[ii] = pos;
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+
+  return COLVARS_OK;
 }
 
 
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 0dda6ab792..93ae594aae 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -227,9 +227,16 @@ protected:
   /// \brief Array of atom objects
   std::vector<cvm::atom> atoms;
 
-  /// \brief Array of atom identifiers for the MD program (0-based)
+  /// \brief Internal atom IDs for host code
   std::vector<int> atoms_ids;
 
+  /// Sorted list of internal atom IDs (populated on-demand by
+  /// create_sorted_ids); used to read coordinate files
+  std::vector<int> sorted_atoms_ids;
+
+  /// Map entries of sorted_atoms_ids onto the original positions in the group
+  std::vector<int> sorted_atoms_ids_map;
+
   /// \brief Dummy atom position
   cvm::atom_pos dummy_atom_pos;
 
@@ -273,19 +280,34 @@ public:
     return atoms.size();
   }
 
-  std::string const print_atom_ids() const;
-
   /// \brief If this option is on, this group merely acts as a wrapper
   /// for a fixed position; any calls to atoms within or to
   /// functions that return disaggregated data will fail
   bool b_dummy;
 
-  /// Sorted list of zero-based (internal) atom ids
-  /// (populated on-demand by create_sorted_ids)
-  std::vector<int> sorted_ids;
+  /// Internal atom IDs (populated during initialization)
+  inline std::vector<int> const &ids() const
+  {
+    return atoms_ids;
+  }
 
-  /// Allocates and populates the sorted list of atom ids
-  int create_sorted_ids(void);
+  std::string const print_atom_ids() const;
+
+  /// Allocates and populates sorted_ids and sorted_ids_map
+  int create_sorted_ids();
+
+  /// Sorted internal atom IDs (populated on-demand by create_sorted_ids);
+  /// used to read coordinate files
+  inline std::vector<int> const &sorted_ids() const
+  {
+    return sorted_atoms_ids;
+  }
+
+  /// Map entries of sorted_atoms_ids onto the original positions in the group
+  inline std::vector<int> const &sorted_ids_map() const
+  {
+    return sorted_atoms_ids_map;
+  }
 
   /// Detect whether two groups share atoms
   /// If yes, returns 1-based number of a common atom; else, returns 0
diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp
index b3b5b3eb16..771955a4fa 100644
--- a/lib/colvars/colvarbias_abf.cpp
+++ b/lib/colvars/colvarbias_abf.cpp
@@ -17,17 +17,17 @@ colvarbias_abf::colvarbias_abf(char const *key)
   : colvarbias(key),
     b_UI_estimator(false),
     b_CZAR_estimator(false),
+    pabf_freq(0),
     system_force(NULL),
     gradients(NULL),
-    pmf(NULL),
     samples(NULL),
+    pmf(NULL),
     z_gradients(NULL),
     z_samples(NULL),
     czar_gradients(NULL),
     czar_pmf(NULL),
     last_gradients(NULL),
-    last_samples(NULL),
-    pabf_freq(0)
+    last_samples(NULL)
 {
 }
 
@@ -315,8 +315,6 @@ colvarbias_abf::~colvarbias_abf()
 
 int colvarbias_abf::update()
 {
-  int       iter;
-
   if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
 
   size_t i;
@@ -368,7 +366,12 @@ int colvarbias_abf::update()
     if ( b_integrate ) {
       if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
         cvm::real err;
-        iter = pmf->integrate(integrate_steps, integrate_tol, err);
+        int iter = pmf->integrate(integrate_steps, integrate_tol, err);
+        if ( iter == integrate_steps ) {
+          cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol)
+            + " in " + cvm::to_str(integrate_steps)
+            + " steps. Residual error: " +  cvm::to_str(err));
+        }
         pmf->set_zero_minimum(); // TODO: do this only when necessary
       }
     }
@@ -485,7 +488,6 @@ int colvarbias_abf::update()
 
 
 int colvarbias_abf::replica_share() {
-  int p;
 
   if ( !cvm::replica_enabled() ) {
     cvm::error("Error: shared ABF: No replicas.\n");
@@ -507,6 +509,7 @@ int colvarbias_abf::replica_share() {
   char* msg_data = new char[msg_total];
 
   if (cvm::replica_index() == 0) {
+    int p;
     // Replica 0 collects the delta gradient and count from the others.
     for (p = 1; p < cvm::replica_num(); p++) {
       // Receive the deltas.
diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp
index 589de1d32a..c717e8f4c6 100644
--- a/lib/colvars/colvarcomp.cpp
+++ b/lib/colvars/colvarcomp.cpp
@@ -21,6 +21,9 @@ colvar::cvc::cvc()
     b_try_scalable(true)
 {
   init_cvc_requires();
+  sup_coeff = 1.0;
+  period = 0.0;
+  wrap_center = 0.0;
 }
 
 
@@ -30,40 +33,47 @@ colvar::cvc::cvc(std::string const &conf)
     b_periodic(false),
     b_try_scalable(true)
 {
+  init_cvc_requires();
+  sup_coeff = 1.0;
+  period = 0.0;
+  wrap_center = 0.0;
+  init(conf);
+}
+
+
+int colvar::cvc::init(std::string const &conf)
+{
+  int error_code = COLVARS_OK;
   if (cvm::debug())
     cvm::log("Initializing cvc base object.\n");
 
-  init_cvc_requires();
-
-  if (get_keyval(conf, "name", this->name, std::string(""), parse_silent)) {
+  get_keyval(conf, "name", this->name, this->name);
+  if (name.size() > 0) {
     // Temporary description until child object is initialized
     description = "cvc " + name;
   } else {
     description = "uninitialized cvc";
   }
 
-  get_keyval(conf, "componentCoeff", sup_coeff, 1.0);
-  get_keyval(conf, "componentExp", sup_np, 1);
+  get_keyval(conf, "componentCoeff", sup_coeff, sup_coeff);
+  get_keyval(conf, "componentExp", sup_np, sup_np);
 
-  get_keyval(conf, "period", period, 0.0);
-  get_keyval(conf, "wrapAround", wrap_center, 0.0);
+  get_keyval(conf, "period", period, period);
+  get_keyval(conf, "wrapAround", wrap_center, wrap_center);
 
-  get_keyval_feature((colvarparse *)this, conf, "debugGradients",
+  get_keyval_feature(dynamic_cast<colvarparse *>(this), conf, "debugGradients",
                      f_cvc_debug_gradient, false, parse_silent);
 
-  {
-    bool b_no_PBC = false;
-    get_keyval(conf, "forceNoPBC", b_no_PBC, false);
-    if (b_no_PBC) {
-      disable(f_cvc_pbc_minimum_image);
-    } else {
-      enable(f_cvc_pbc_minimum_image);
-    }
-    // this does not use get_keyval_feature() only for backward compatibility
+  bool b_no_PBC = !is_enabled(f_cvc_pbc_minimum_image); // Enabled by default
+  get_keyval(conf, "forceNoPBC", b_no_PBC, b_no_PBC);
+  if (b_no_PBC) {
+    disable(f_cvc_pbc_minimum_image);
+  } else {
+    enable(f_cvc_pbc_minimum_image);
   }
 
   // Attempt scalable calculations when in parallel? (By default yes, if available)
-  get_keyval(conf, "scalable", b_try_scalable, true);
+  get_keyval(conf, "scalable", b_try_scalable, b_try_scalable);
 
   if (cvm::debug())
     cvm::log("Done initializing cvc base object.\n");
diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h
index 52078a3a30..b5c3b16098 100644
--- a/lib/colvars/colvarcomp.h
+++ b/lib/colvars/colvarcomp.h
@@ -98,12 +98,14 @@ public:
 
   /// \brief Constructor
   ///
-  /// At least one constructor which reads a string should be defined
-  /// for every class inheriting from cvc \param conf Contents
-  /// of the configuration file pertaining to this \link cvc
-  /// \endlink
+  /// Calls the init() function of the class
   cvc(std::string const &conf);
 
+  /// An init function should be defined for every class inheriting from cvc
+  /// \param conf Contents of the configuration file pertaining to this \link
+  /// cvc \endlink
+  virtual int init(std::string const &conf);
+
   /// \brief Within the constructor, make a group parse its own
   /// options from the provided configuration string
   /// Returns reference to new group
@@ -231,7 +233,7 @@ public:
   virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
                                   colvarvalue const &x2) const;
 
-  /// \brief Wrapp value (for periodic/symmetric cvcs)
+  /// \brief Wrap value (for periodic/symmetric cvcs)
   virtual void wrap(colvarvalue &x) const;
 
   /// \brief Pointers to all atom groups, to let colvars collect info
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 9911f4c87e..57b2a9a625 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -148,7 +148,7 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
 cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
                                       colvarvalue const &x2) const
 {
-  return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
+  return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2();
 }
 
 
@@ -192,7 +192,7 @@ colvar::distance_z::distance_z(std::string const &conf)
   // this group is optional
   ref2 = parse_group(conf, "ref2", true);
 
-  if (ref2 && ref2->size()) {
+  if ( ref2 ) {
     cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
     fixed_axis = false;
     if (key_lookup(conf, "axis"))
@@ -306,7 +306,7 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
   if (!ref1->noforce)
     ref1->apply_colvar_force(force.real_value);
 
-  if (ref2 && ref2->size() && !ref2->noforce)
+  if (ref2 && !ref2->noforce)
     ref2->apply_colvar_force(force.real_value);
 
   if (!main->noforce)
@@ -464,7 +464,7 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
   if (!ref1->noforce)
     ref1->apply_colvar_force(force.real_value);
 
-  if (ref2 && ref2->size() && !ref2->noforce)
+  if (ref2 && !ref2->noforce)
     ref2->apply_colvar_force(force.real_value);
 
   if (!main->noforce)
@@ -979,14 +979,12 @@ colvar::rmsd::rmsd(std::string const &conf)
                      "if provided, must be non-zero.\n");
           return;
         }
-      } else {
-        // if not, rely on existing atom indices for the group
-        atoms->create_sorted_ids();
-        ref_pos.resize(atoms->size());
       }
 
-      cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms->sorted_ids,
-                        ref_pos_col, ref_pos_col_value);
+      ref_pos.resize(atoms->size());
+
+      cvm::load_coords(ref_pos_file.c_str(), &ref_pos, atoms,
+                       ref_pos_col, ref_pos_col_value);
     }
   }
 
@@ -1172,13 +1170,11 @@ colvar::eigenvector::eigenvector(std::string const &conf)
                             "if provided, must be non-zero.\n");
           return;
         }
-      } else {
-        // if not, use atom indices
-        atoms->create_sorted_ids();
       }
 
       ref_pos.resize(atoms->size());
-      cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value);
+      cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
+                       file_col, file_col_value);
     }
   }
 
@@ -1249,13 +1245,11 @@ colvar::eigenvector::eigenvector(std::string const &conf)
           cvm::error("Error: vectorColValue, if provided, must be non-zero.\n");
           return;
         }
-      } else {
-        // if not, use atom indices
-        atoms->create_sorted_ids();
       }
 
       eigenvec.resize(atoms->size());
-      cvm::load_coords(file_name.c_str(), eigenvec, atoms->sorted_ids, file_col, file_col_value);
+      cvm::load_coords(file_name.c_str(), &eigenvec, atoms,
+                       file_col, file_col_value);
     }
   }
 
diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp
index 2650a9fe18..498ef7c2f5 100644
--- a/lib/colvars/colvarcomp_rotations.cpp
+++ b/lib/colvars/colvarcomp_rotations.cpp
@@ -50,12 +50,11 @@ colvar::orientation::orientation(std::string const &conf)
                             "if provided, must be non-zero.\n");
           return;
         }
-      } else {
-        // if not, use atom indices
-        atoms->create_sorted_ids();
       }
+
       ref_pos.resize(atoms->size());
-      cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value);
+      cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
+                       file_col, file_col_value);
     }
   }
 
diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp
index a058ad55c2..80bd6670d3 100644
--- a/lib/colvars/colvardeps.cpp
+++ b/lib/colvars/colvardeps.cpp
@@ -549,7 +549,7 @@ void colvardeps::init_cv_requires() {
     f_req_exclude(f_cv_custom_function, f_cv_scripted);
 
     init_feature(f_cv_periodic, "periodic", f_type_static);
-    f_req_self(f_cv_periodic, f_cv_homogeneous);
+    f_req_self(f_cv_periodic, f_cv_scalar);
     init_feature(f_cv_scalar, "scalar", f_type_static);
     init_feature(f_cv_linear, "linear", f_type_static);
     init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp
index 1ac4aae133..407b336afd 100644
--- a/lib/colvars/colvargrid.cpp
+++ b/lib/colvars/colvargrid.cpp
@@ -329,7 +329,6 @@ void integrate_potential::update_div_local(const std::vector<int> &ix0)
   const int linear_index = address(ix0);
   int i, j, k;
   std::vector<int> ix = ix0;
-  const cvm::real * g;
 
   if (nd == 2) {
     // gradients at grid points surrounding the current scalar grid point
@@ -783,7 +782,7 @@ void integrate_potential::atimes(const std::vector<cvm::real> &A, std::vector<cv
 /// Inversion of preconditioner matrix (e.g. diagonal of the Laplacian)
 void integrate_potential::asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x)
 {
-  for (size_t i=0; i<nt; i++) {
+  for (size_t i=0; i<int(nt); i++) {
     x[i] = b[i] * inv_lap_diag[i]; // Jacobi preconditioner - little benefit in tests so far
   }
   return;
@@ -803,7 +802,7 @@ void integrate_potential::nr_linbcg_sym(const std::vector<cvm::real> &b, std::ve
 
   iter=0;
   atimes(x,r);
-  for (j=0;j<nt;j++) {
+  for (j=0;j<int(nt);j++) {
     r[j]=b[j]-r[j];
   }
   bnrm=l2norm(b);
@@ -814,26 +813,26 @@ void integrate_potential::nr_linbcg_sym(const std::vector<cvm::real> &b, std::ve
   bkden = 1.0;
   while (iter < itmax) {
     ++iter;
-    for (bknum=0.0,j=0;j<nt;j++) {
+    for (bknum=0.0,j=0;j<int(nt);j++) {
       bknum += r[j]*r[j];  // precon: z[j]*r[j]
     }
     if (iter == 1) {
-      for (j=0;j<nt;j++) {
+      for (j=0;j<int(nt);j++) {
         p[j] = r[j];  // precon: p[j] = z[j]
       }
     } else {
       bk=bknum/bkden;
-      for (j=0;j<nt;j++) {
+      for (j=0;j<int(nt);j++) {
         p[j] = bk*p[j] + r[j];  // precon:  bk*p[j] + z[j]
       }
     }
     bkden = bknum;
     atimes(p,z);
-    for (akden=0.0,j=0;j<nt;j++) {
+    for (akden=0.0,j=0;j<int(nt);j++) {
       akden += z[j]*p[j];
     }
     ak = bknum/akden;
-    for (j=0;j<nt;j++) {
+    for (j=0;j<int(nt);j++) {
       x[j] += ak*p[j];
       r[j] -= ak*z[j];
     }
diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp
index 9898e2d5e8..c9e2d14116 100644
--- a/lib/colvars/colvarmodule.cpp
+++ b/lib/colvars/colvarmodule.cpp
@@ -30,10 +30,10 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   depth_s = 0;
   cv_traj_os = NULL;
 
-  // pointer to the proxy object
   if (proxy == NULL) {
-    proxy = proxy_in;
-    parse = new colvarparse();
+    proxy = proxy_in; // Pointer to the proxy object
+    parse = new colvarparse(); // Parsing object for global options
+    version_int = proxy->get_version_from_string(COLVARS_VERSION);
   } else {
     // TODO relax this error to handle multiple molecules in VMD
     // once the module is not static anymore
@@ -58,22 +58,24 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
   // "it_restart" will be set by the input state file, if any;
   // "it" should be updated by the proxy
   colvarmodule::it = colvarmodule::it_restart = 0;
-  colvarmodule::it_restart_from_state_file = true;
 
-  colvarmodule::use_scripted_forces = false;
+  use_scripted_forces = false;
+  scripting_after_biases = false;
 
-  colvarmodule::b_analysis = false;
+  b_analysis = false;
 
   colvarmodule::debug_gradients_step_size = 1.0e-07;
 
   colvarmodule::rotation::monitor_crossings = false;
   colvarmodule::rotation::crossing_threshold = 1.0e-02;
 
-  colvarmodule::cv_traj_freq = 100;
-  colvarmodule::restart_out_freq = proxy->restart_frequency();
+  cv_traj_freq = 100;
+  restart_out_freq = proxy->restart_frequency();
 
   // by default overwrite the existing trajectory file
-  colvarmodule::cv_traj_append = false;
+  cv_traj_append = false;
+
+  cv_traj_write_labels = true;
 }
 
 
@@ -189,26 +191,27 @@ int colvarmodule::parse_config(std::string &conf)
 {
   extra_conf.clear();
 
-  // parse global options
+  // Parse global options
   if (catch_input_errors(parse_global_params(conf))) {
     return get_error();
   }
 
-  // parse the options for collective variables
+  // Parse the options for collective variables
   if (catch_input_errors(parse_colvars(conf))) {
     return get_error();
   }
 
-  // parse the options for biases
+  // Parse the options for biases
   if (catch_input_errors(parse_biases(conf))) {
     return get_error();
   }
 
-  // done parsing known keywords, check that all keywords found were valid ones
+  // Done parsing known keywords, check that all keywords found were valid ones
   if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) {
     return get_error();
   }
 
+  // Parse auto-generated configuration (e.g. for back-compatibility)
   if (extra_conf.size()) {
     catch_input_errors(parse_global_params(extra_conf));
     catch_input_errors(parse_colvars(extra_conf));
@@ -222,13 +225,11 @@ int colvarmodule::parse_config(std::string &conf)
   cvm::log("Collective variables module (re)initialized.\n");
   cvm::log(cvm::line_marker);
 
-  // update any necessary proxy data
+  // Update any necessary proxy data
   proxy->setup();
 
-  if (cv_traj_os != NULL) {
-    // configuration might have changed, better redo the labels
-    write_traj_label(*cv_traj_os);
-  }
+  // configuration might have changed, better redo the labels
+  cv_traj_write_labels = true;
 
   return get_error();
 }
@@ -279,15 +280,18 @@ int colvarmodule::parse_global_params(std::string const &conf)
   parse->get_keyval(conf, "colvarsTrajAppend",
                     cv_traj_append, cv_traj_append, colvarparse::parse_silent);
 
-  parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false);
+  parse->get_keyval(conf, "scriptedColvarForces",
+                    use_scripted_forces, use_scripted_forces);
 
-  parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true);
+  parse->get_keyval(conf, "scriptingAfterBiases",
+                    scripting_after_biases, scripting_after_biases);
 
   if (use_scripted_forces && !proxy->force_script_defined) {
-    cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
+    return cvm::error("User script for scripted colvar forces not found.",
+                      INPUT_ERROR);
   }
 
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+  return cvm::get_error();
 }
 
 
@@ -432,13 +436,13 @@ int colvarmodule::parse_biases(std::string const &conf)
 }
 
 
-int colvarmodule::num_variables() const
+size_t colvarmodule::num_variables() const
 {
   return colvars.size();
 }
 
 
-int colvarmodule::num_variables_feature(int feature_id) const
+size_t colvarmodule::num_variables_feature(int feature_id) const
 {
   size_t n = 0;
   for (std::vector<colvar *>::const_iterator cvi = colvars.begin();
@@ -452,13 +456,13 @@ int colvarmodule::num_variables_feature(int feature_id) const
 }
 
 
-int colvarmodule::num_biases() const
+size_t colvarmodule::num_biases() const
 {
   return biases.size();
 }
 
 
-int colvarmodule::num_biases_feature(int feature_id) const
+size_t colvarmodule::num_biases_feature(int feature_id) const
 {
   size_t n = 0;
   for (std::vector<colvarbias *>::const_iterator bi = biases.begin();
@@ -472,7 +476,7 @@ int colvarmodule::num_biases_feature(int feature_id) const
 }
 
 
-int colvarmodule::num_biases_type(std::string const &type) const
+size_t colvarmodule::num_biases_type(std::string const &type) const
 {
   size_t n = 0;
   for (std::vector<colvarbias *>::const_iterator bi = biases.begin();
@@ -971,13 +975,20 @@ int colvarmodule::write_restart_file(std::string const &out_name)
 int colvarmodule::write_traj_files()
 {
   if (cv_traj_os == NULL) {
-    open_traj_file(cv_traj_name);
+    if (open_traj_file(cv_traj_name) != COLVARS_OK) {
+      return cvm::get_error();
+    } else {
+      cv_traj_write_labels = true;
+    }
   }
 
   // write labels in the traj file every 1000 lines and at first timestep
-  if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
+  if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 ||
+      cvm::step_relative() == 0 ||
+      cv_traj_write_labels) {
     write_traj_label(*cv_traj_os);
   }
+  cv_traj_write_labels = false;
 
   if ((cvm::step_absolute() % cv_traj_freq) == 0) {
     write_traj(*cv_traj_os);
@@ -1064,7 +1075,7 @@ int colvarmodule::reset()
 {
   parse->init();
 
-  cvm::log("Resetting the Collective Variables Module.\n");
+  cvm::log("Resetting the Collective Variables module.\n");
 
   // Iterate backwards because we are deleting the elements as we go
   for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
@@ -1073,6 +1084,7 @@ int colvarmodule::reset()
     delete *bi; // the bias destructor updates the biases array
   }
   biases.clear();
+  biases_active_.clear();
 
   // Iterate backwards because we are deleting the elements as we go
   for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
@@ -1088,7 +1100,7 @@ int colvarmodule::reset()
   proxy->reset();
 
   if (cv_traj_os != NULL) {
-    // Do not close file here, as we might not be done with it yet.
+    // Do not close traj file here, as we might not be done with it yet.
     proxy->flush_output_stream(cv_traj_os);
   }
 
@@ -1188,12 +1200,10 @@ std::istream & colvarmodule::read_restart(std::istream &is)
     // read global restart information
     std::string restart_conf;
     if (is >> colvarparse::read_block("configuration", restart_conf)) {
-      if (it_restart_from_state_file) {
-        parse->get_keyval(restart_conf, "step",
-                          it_restart, (size_t) 0,
-                          colvarparse::parse_silent);
+      parse->get_keyval(restart_conf, "step",
+                        it_restart, (size_t) 0,
+                        colvarparse::parse_silent);
         it = it_restart;
-      }
       std::string restart_version;
       parse->get_keyval(restart_conf, "version",
                         restart_version, std::string(""),
@@ -1688,40 +1698,59 @@ int cvm::read_index_file(char const *filename)
   return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
 }
 
+
 int cvm::load_atoms(char const *file_name,
                     cvm::atom_group &atoms,
                     std::string const &pdb_field,
-                    double const pdb_field_value)
+                    double pdb_field_value)
 {
   return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value);
 }
 
-int cvm::load_coords(char const *file_name,
-                     std::vector<cvm::atom_pos> &pos,
-                     const std::vector<int> &indices,
-                     std::string const &pdb_field,
-                     double const pdb_field_value)
-{
-  // Differentiate between PDB and XYZ files
-  // for XYZ files, use CVM internal parser
-  // otherwise call proxy function for PDB
 
-  std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name);
+int cvm::load_coords(char const *file_name,
+                     std::vector<cvm::rvector> *pos,
+                     cvm::atom_group *atoms,
+                     std::string const &pdb_field,
+                     double pdb_field_value)
+{
+  int error_code = COLVARS_OK;
+
+  std::string const ext(strlen(file_name) > 4 ?
+                        (file_name + (strlen(file_name) - 4)) :
+                        file_name);
+
+  atoms->create_sorted_ids();
+
+  std::vector<cvm::rvector> sorted_pos(atoms->size(), cvm::rvector(0.0));
+
+  // Differentiate between PDB and XYZ files
   if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) {
-    if ( pdb_field.size() > 0 ) {
-      cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR);
-      return COLVARS_ERROR;
+    if (pdb_field.size() > 0) {
+      return cvm::error("Error: PDB column may not be specified "
+                        "for XYZ coordinate files.\n", INPUT_ERROR);
     }
-    return cvm::load_coords_xyz(file_name, pos, indices);
+    // For XYZ files, use internal parser
+    error_code |= cvm::load_coords_xyz(file_name, &sorted_pos, atoms);
   } else {
-    return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value);
+    // Otherwise, call proxy function for PDB
+    error_code |= proxy->load_coords(file_name,
+                                     sorted_pos, atoms->sorted_ids(),
+                                     pdb_field, pdb_field_value);
   }
+
+  std::vector<int> const &map = atoms->sorted_ids_map();
+  for (size_t i = 0; i < atoms->size(); i++) {
+    (*pos)[map[i]] = sorted_pos[i];
+  }
+
+  return error_code;
 }
 
 
 int cvm::load_coords_xyz(char const *filename,
-                         std::vector<atom_pos> &pos,
-                         const std::vector<int> &indices)
+                         std::vector<rvector> *pos,
+                         cvm::atom_group *atoms)
 {
   std::ifstream xyz_is(filename);
   unsigned int natoms;
@@ -1736,12 +1765,12 @@ int cvm::load_coords_xyz(char const *filename,
   cvm::getline(xyz_is, line);
   cvm::getline(xyz_is, line);
   xyz_is.width(255);
-  std::vector<atom_pos>::iterator pos_i = pos.begin();
+  std::vector<atom_pos>::iterator pos_i = pos->begin();
 
-  if (pos.size() != natoms) { // Use specified indices
+  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++) {
+    std::vector<int>::const_iterator index = atoms->sorted_ids().begin();
+    for ( ; pos_i != pos->end() ; pos_i++, index++) {
 
       while ( next < *index ) {
         cvm::getline(xyz_is, line);
@@ -1751,7 +1780,7 @@ int cvm::load_coords_xyz(char const *filename,
       xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
     }
   } else {          // Use all positions
-    for ( ; pos_i != pos.end() ; pos_i++) {
+    for ( ; pos_i != pos->end() ; pos_i++) {
       xyz_is >> symbol;
       xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
     }
@@ -1792,17 +1821,13 @@ void cvm::request_total_force()
   proxy->request_total_force(true);
 }
 
-cvm::rvector cvm::position_distance(atom_pos const &pos1,
-                                            atom_pos const &pos2)
+
+cvm::rvector cvm::position_distance(cvm::atom_pos const &pos1,
+                                    cvm::atom_pos const &pos2)
 {
   return proxy->position_distance(pos1, pos2);
 }
 
-cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
-                                      cvm::atom_pos const &pos2)
-{
-  return proxy->position_dist2(pos1, pos2);
-}
 
 cvm::real cvm::rand_gaussian(void)
 {
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 64a98b77a0..f3d99a2b23 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -73,8 +73,17 @@ private:
   /// Impossible to initialize the main object without arguments
   colvarmodule();
 
+  /// Integer representing the version string (allows comparisons)
+  int version_int;
+
 public:
 
+  /// Get the version number (higher = more recent)
+  int version_number() const
+  {
+    return version_int;
+  }
+
   friend class colvarproxy;
   // TODO colvarscript should be unaware of colvarmodule's internals
   friend class colvarscript;
@@ -158,10 +167,6 @@ public:
     return it;
   }
 
-  /// If true, get it_restart from the state file; if set to false,
-  /// the MD program is providing it
-  bool it_restart_from_state_file;
-
   /// \brief Finite difference step size (if there is no dynamics, or
   /// if gradients need to be tested independently from the size of
   /// dt)
@@ -306,19 +311,19 @@ private:
 public:
 
   /// Return how many variables are defined
-  int num_variables() const;
+  size_t num_variables() const;
 
   /// Return how many variables have this feature enabled
-  int num_variables_feature(int feature_id) const;
+  size_t num_variables_feature(int feature_id) const;
 
   /// Return how many biases are defined
-  int num_biases() const;
+  size_t num_biases() const;
 
   /// Return how many biases have this feature enabled
-  int num_biases_feature(int feature_id) const;
+  size_t num_biases_feature(int feature_id) const;
 
   /// Return how many biases of this type are defined
-  int num_biases_type(std::string const &type) const;
+  size_t num_biases_type(std::string const &type) const;
 
   /// Return the names of time-dependent biases with forces enabled (ABF,
   /// metadynamics, etc)
@@ -479,9 +484,6 @@ public:
   /// Print a message to the main log and set global error code
   static int error(std::string const &message, int code = COLVARS_ERROR);
 
-  /// Print a message to the main log and exit normally
-  static void exit(std::string const &message);
-
   // Replica exchange commands.
   static bool replica_enabled();
   static int replica_index();
@@ -495,15 +497,6 @@ public:
   static rvector position_distance(atom_pos const &pos1,
                                    atom_pos const &pos2);
 
-  /// \brief Get the square distance between two positions (with
-  /// periodic boundary conditions handled transparently)
-  ///
-  /// 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)
-  static real position_dist2(atom_pos const &pos1,
-                             atom_pos const &pos2);
-
   /// \brief Names of groups from a Gromacs .ndx file to be read at startup
   std::list<std::string> index_group_names;
 
@@ -513,29 +506,36 @@ public:
   /// \brief Read a Gromacs .ndx file
   int 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
-  /// field to determine which are the atoms to be set
+  /// \brief Select atom IDs from a file (usually PDB) \param filename name of
+  /// the file \param atoms array into which atoms read from "filename" will be
+  /// appended \param pdb_field (optional) if the file is a PDB and this
+  /// string is non-empty, select atoms for which this field is non-zero
+  /// \param pdb_field_value (optional) if non-zero, select only atoms whose
+  /// pdb_field equals this
   static int load_atoms(char const *filename,
                         atom_group &atoms,
                         std::string const &pdb_field,
-                        double const pdb_field_value = 0.0);
+                        double pdb_field_value = 0.0);
 
-  /// \brief Load the coordinates for a group of atoms from a file
-  /// (PDB or XYZ)
+  /// \brief Load coordinates for a group of atoms from a file (PDB or XYZ);
+  /// if "pos" is already allocated, the number of its elements must match the
+  /// number of entries in "filename" \param filename name of the file \param
+  /// pos array of coordinates \param atoms group containing the atoms (used
+  /// to obtain internal IDs) \param pdb_field (optional) if the file is a PDB
+  /// and this string is non-empty, select atoms for which this field is
+  /// non-zero \param pdb_field_value (optional) if non-zero, select only
+  /// atoms whose pdb_field equals this
   static int load_coords(char const *filename,
-                         std::vector<atom_pos> &pos,
-                         const std::vector<int> &indices,
+                         std::vector<rvector> *pos,
+                         atom_group *atoms,
                          std::string const &pdb_field,
-                         double const pdb_field_value = 0.0);
+                         double pdb_field_value = 0.0);
 
   /// \brief Load the coordinates for a group of atoms from an
   /// XYZ file
   static int load_coords_xyz(char const *filename,
-                              std::vector<atom_pos> &pos,
-                              const std::vector<int> &indices);
+                             std::vector<rvector> *pos,
+                             atom_group *atoms);
 
   /// Frequency for collective variables trajectory output
   static size_t cv_traj_freq;
@@ -568,6 +568,9 @@ protected:
   /// Appending to the existing trajectory file?
   bool cv_traj_append;
 
+  /// Write labels at the next iteration
+  bool cv_traj_write_labels;
+
 private:
 
   /// Counter for the current depth in the object hierarchy (useg e.g. in output)
diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp
index 8767d5f459..c2992c309c 100644
--- a/lib/colvars/colvarproxy.cpp
+++ b/lib/colvars/colvarproxy.cpp
@@ -26,7 +26,10 @@
 
 
 
-colvarproxy_system::colvarproxy_system() {}
+colvarproxy_system::colvarproxy_system()
+{
+  reset_pbc_lattice();
+}
 
 
 colvarproxy_system::~colvarproxy_system() {}
@@ -55,10 +58,73 @@ bool colvarproxy_system::total_forces_same_step() const
 }
 
 
-cvm::real colvarproxy_system::position_dist2(cvm::atom_pos const &pos1,
-                                             cvm::atom_pos const &pos2)
+inline int round_to_integer(cvm::real x)
 {
-  return (position_distance(pos1, pos2)).norm2();
+  return std::floor(x+0.5);
+}
+
+
+void colvarproxy_system::update_pbc_lattice()
+{
+  // Periodicity is assumed in all directions
+
+  if (boundaries_type == boundaries_unsupported ||
+      boundaries_type == boundaries_non_periodic) {
+    cvm::error("Error: setting PBC lattice with unsupported boundaries.\n",
+               BUG_ERROR);
+    return;
+  }
+
+  {
+    cvm::rvector const v = cvm::rvector::outer(unit_cell_y, unit_cell_z);
+    reciprocal_cell_x = v/(v*unit_cell_x);
+  }
+  {
+    cvm::rvector const v = cvm::rvector::outer(unit_cell_z, unit_cell_x);
+    reciprocal_cell_y = v/(v*unit_cell_y);
+  }
+  {
+    cvm::rvector const v = cvm::rvector::outer(unit_cell_x, unit_cell_y);
+    reciprocal_cell_z = v/(v*unit_cell_z);
+  }
+}
+
+
+void colvarproxy_system::reset_pbc_lattice()
+{
+  unit_cell_x.reset();
+  unit_cell_y.reset();
+  unit_cell_z.reset();
+  reciprocal_cell_x.reset();
+  reciprocal_cell_y.reset();
+  reciprocal_cell_z.reset();
+}
+
+
+cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1,
+                                                   cvm::atom_pos const &pos2)
+  const
+{
+  if (boundaries_type == boundaries_unsupported) {
+    cvm::error("Error: unsupported boundary conditions.\n", INPUT_ERROR);
+  }
+
+  cvm::rvector diff = (pos2 - pos1);
+
+  if (boundaries_type == boundaries_non_periodic) return diff;
+
+  cvm::real const x_shift = round_to_integer(reciprocal_cell_x*diff);
+  cvm::real const y_shift = round_to_integer(reciprocal_cell_y*diff);
+  cvm::real const z_shift = round_to_integer(reciprocal_cell_z*diff);
+
+  diff.x -= x_shift*unit_cell_x.x + y_shift*unit_cell_y.x +
+    z_shift*unit_cell_z.x;
+  diff.y -= x_shift*unit_cell_x.y + y_shift*unit_cell_y.y +
+    z_shift*unit_cell_z.y;
+  diff.z -= x_shift*unit_cell_x.z + y_shift*unit_cell_y.z +
+    z_shift*unit_cell_z.z;
+
+  return diff;
 }
 
 
@@ -132,7 +198,7 @@ void colvarproxy_atoms::clear_atom(int index)
 int colvarproxy_atoms::load_atoms(char const *filename,
                                   cvm::atom_group &atoms,
                                   std::string const &pdb_field,
-                                  double const)
+                                  double)
 {
   return cvm::error("Error: loading atom identifiers from a file "
                     "is currently not implemented.\n",
@@ -142,9 +208,9 @@ int colvarproxy_atoms::load_atoms(char const *filename,
 
 int colvarproxy_atoms::load_coords(char const *filename,
                                    std::vector<cvm::atom_pos> &pos,
-                                   const std::vector<int> &indices,
+                                   std::vector<int> const &sorted_ids,
                                    std::string const &pdb_field,
-                                   double const)
+                                   double)
 {
   return cvm::error("Error: loading atomic coordinates from a file "
                     "is currently not implemented.\n",
@@ -661,6 +727,7 @@ int colvarproxy_io::close_output_stream(std::string const &output_name)
   for ( ; osi != output_files.end(); osi++, osni++) {
     if (*osni == output_name) {
       ((std::ofstream *) (*osi))->close();
+      delete *osi;
       output_files.erase(osi);
       output_stream_names.erase(osni);
       return COLVARS_OK;
@@ -729,3 +796,13 @@ size_t colvarproxy::restart_frequency()
   return 0;
 }
 
+
+int colvarproxy::get_version_from_string(char const *version_string)
+{
+  std::string const v(version_string);
+  std::istringstream is(v.substr(0, 4) + v.substr(5, 2) + v.substr(8, 2));
+  int newint;
+  is >> newint;
+  return newint;
+}
+
diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h
index bf29048270..845f93c9bc 100644
--- a/lib/colvars/colvarproxy.h
+++ b/lib/colvars/colvarproxy.h
@@ -14,6 +14,7 @@
 #include <list>
 
 #include "colvarmodule.h"
+#include "colvartypes.h"
 #include "colvarvalue.h"
 
 
@@ -29,7 +30,7 @@
 /// To interface to a new MD engine, the simplest solution is to derive a new
 /// class from \link colvarproxy \endlink.  Currently implemented are: \link
 /// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link
-/// colvarproxy_vmd, \endlink.
+/// colvarproxy_vmd \endlink.
 
 
 // forward declarations
@@ -68,14 +69,16 @@ public:
 
   /// \brief Get the PBC-aware distance vector between two positions
   virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
-                                         cvm::atom_pos const &pos2) = 0;
+                                         cvm::atom_pos const &pos2) const;
 
-  /// \brief Get the PBC-aware square distance between two positions;
-  /// may need to be reimplemented independently from position_distance() for optimization purposes
-  virtual cvm::real position_dist2(cvm::atom_pos const &pos1,
-                                   cvm::atom_pos const &pos2);
+  /// Recompute PBC reciprocal lattice (assumes XYZ periodicity)
+  void update_pbc_lattice();
 
-  /// Tell the proxy whether total forces are needed (may not always be available)
+  /// Set the lattice vectors to zero
+  void reset_pbc_lattice();
+
+  /// \brief Tell the proxy whether total forces are needed (they may not
+  /// always be available)
   virtual void request_total_force(bool yesno);
 
   /// Are total forces being used?
@@ -83,6 +86,29 @@ public:
 
   /// Are total forces from the current step available?
   virtual bool total_forces_same_step() const;
+
+protected:
+
+  /// \brief Type of boundary conditions
+  ///
+  /// Orthogonal and triclinic cells are made available to objects.
+  /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS)
+  /// minimum-image distances are computed by the host engine regardless.
+  enum Boundaries_type {
+    boundaries_non_periodic,
+    boundaries_pbc_ortho,
+    boundaries_pbc_triclinic,
+    boundaries_unsupported
+  };
+
+  /// Type of boundary conditions
+  Boundaries_type boundaries_type;
+
+  /// Bravais lattice vectors
+  cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
+
+  /// Reciprocal lattice vectors
+  cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
 };
 
 
@@ -121,24 +147,30 @@ public:
   /// (costly) set the corresponding atoms_ncopies to zero
   virtual void clear_atom(int index);
 
-  /// \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
-  /// "filename" is a PDB file, use this field to determine which are
-  /// the atoms to be set
+  /// \brief Select atom IDs from a file (usually PDB) \param filename name of
+  /// the file \param atoms array to which atoms read from "filename" will be
+  /// appended \param pdb_field (optional) if the file is a PDB and this
+  /// string is non-empty, select atoms for which this field is non-zero
+  /// \param pdb_field_value (optional) if non-zero, select only atoms whose
+  /// pdb_field equals this
   virtual int load_atoms(char const *filename,
                          cvm::atom_group &atoms,
                          std::string const &pdb_field,
-                         double const pdb_field_value = 0.0);
+                         double 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
-  /// elements must match the number of atoms in "filename"
+  /// \brief Load a set of coordinates from a file (usually PDB); if "pos" is
+  /// already allocated, the number of its elements must match the number of
+  /// entries in "filename" \param filename name of the file \param pos array
+  /// of coordinates \param sorted_ids array of sorted internal IDs, used to
+  /// loop through the file only once \param pdb_field (optional) if the file
+  /// is a PDB and this string is non-empty, select atoms for which this field
+  /// is non-zero \param pdb_field_value (optional) if non-zero, select only
+  /// atoms whose pdb_field equals this
   virtual int load_coords(char const *filename,
                           std::vector<cvm::atom_pos> &pos,
-                          const std::vector<int> &indices,
+                          std::vector<int> const &sorted_ids,
                           std::string const &pdb_field,
-                          double const pdb_field_value = 0.0);
+                          double pdb_field_value = 0.0);
 
   /// Clear atomic data
   int reset();
@@ -636,6 +668,15 @@ public:
     return b_simulation_running;
   }
 
+  /// Convert a version string "YYYY-MM-DD" into an integer
+  int get_version_from_string(char const *version_string);
+
+  /// Get the version number (higher = more recent)
+  int version_number() const
+  {
+    return version_int;
+  }
+
 protected:
 
   /// Whether a simulation is running (warn against irrecovarable errors)
@@ -644,6 +685,9 @@ protected:
   /// Whether the entire module should be deallocated by the host engine
   bool b_delete_requested;
 
+  /// Integer representing the version string (allows comparisons)
+  int version_int;
+
 };
 
 
diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h
index dc4b8bd07e..68f5cd13ab 100644
--- a/lib/colvars/colvars_version.h
+++ b/lib/colvars/colvars_version.h
@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-01-17"
+#define COLVARS_VERSION "2018-04-29"
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars
diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h
index 97257d18ad..4ef9557dcb 100644
--- a/lib/colvars/colvartypes.h
+++ b/lib/colvars/colvartypes.h
@@ -130,7 +130,7 @@ public:
     }
   }
 
-  inline void operator *= (cvm::real const &a)
+  inline void operator *= (cvm::real a)
   {
     size_t i;
     for (i = 0; i < this->size(); i++) {
@@ -138,7 +138,7 @@ public:
     }
   }
 
-  inline void operator /= (cvm::real const &a)
+  inline void operator /= (cvm::real a)
   {
     size_t i;
     for (i = 0; i < this->size(); i++) {
@@ -146,7 +146,8 @@ public:
     }
   }
 
-  inline friend vector1d<T> operator + (vector1d<T> const &v1, vector1d<T> const &v2)
+  inline friend vector1d<T> operator + (vector1d<T> const &v1,
+                                        vector1d<T> const &v2)
   {
     check_sizes(v1.size(), v2.size());
     vector1d<T> result(v1.size());
@@ -157,7 +158,8 @@ public:
     return result;
   }
 
-  inline friend vector1d<T> operator - (vector1d<T> const &v1, vector1d<T> const &v2)
+  inline friend vector1d<T> operator - (vector1d<T> const &v1,
+                                        vector1d<T> const &v2)
   {
     check_sizes(v1.size(), v2.size());
     vector1d<T> result(v1.size());
@@ -168,7 +170,7 @@ public:
     return result;
   }
 
-  inline friend vector1d<T> operator * (vector1d<T> const &v, cvm::real const &a)
+  inline friend vector1d<T> operator * (vector1d<T> const &v, cvm::real a)
   {
     vector1d<T> result(v.size());
     size_t i;
@@ -178,12 +180,12 @@ public:
     return result;
   }
 
-  inline friend vector1d<T> operator * (cvm::real const &a, vector1d<T> const &v)
+  inline friend vector1d<T> operator * (cvm::real a, vector1d<T> const &v)
   {
     return v * a;
   }
 
-  inline friend vector1d<T> operator / (vector1d<T> const &v, cvm::real const &a)
+  inline friend vector1d<T> operator / (vector1d<T> const &v, cvm::real a)
   {
     vector1d<T> result(v.size());
     size_t i;
@@ -246,7 +248,8 @@ public:
   }
 
   /// Assign a vector to a slice of this vector
-  inline void sliceassign(size_t const i1, size_t const i2, vector1d<T> const &v)
+  inline void sliceassign(size_t const i1, size_t const i2,
+                          vector1d<T> const &v)
   {
     if ((i2 < i1) || (i2 >= this->size())) {
       cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
@@ -259,12 +262,13 @@ public:
 
   /// Formatted output
 
-  inline size_t output_width(size_t const &real_width) const
+  inline size_t output_width(size_t real_width) const
   {
     return real_width*(this->size()) + 3*(this->size()-1) + 4;
   }
 
-  inline friend std::istream & operator >> (std::istream &is, cvm::vector1d<T> &v)
+  inline friend std::istream & operator >> (std::istream &is,
+                                            cvm::vector1d<T> &v)
   {
     if (v.size() == 0) return is;
     size_t const start_pos = is.tellg();
@@ -288,7 +292,8 @@ public:
     return is;
   }
 
-  inline friend std::ostream & operator << (std::ostream &os, cvm::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();
@@ -377,6 +382,15 @@ protected:
     {
       return vector1d<T>(length, data);
     }
+    inline int set(cvm::vector1d<T> const &v) const
+    {
+      if (v.size() != length) {
+        return cvm::error("Error: setting a matrix row from a vector of "
+                          "incompatible size.\n", BUG_ERROR);
+      }
+      for (size_t i = 0; i < length; i++) data[i] = v[i];
+      return COLVARS_OK;
+    }
   };
 
   std::vector<T> data;
@@ -515,9 +529,12 @@ public:
   {
     if ((m1.outer_length != m2.outer_length) ||
         (m1.inner_length != m2.inner_length)) {
-      cvm::error("Error: trying to perform an operation between matrices of different 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");
+      cvm::error("Error: trying to perform an operation between "
+                 "matrices of different 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");
     }
   }
 
@@ -539,7 +556,7 @@ public:
     }
   }
 
-  inline void operator *= (cvm::real const &a)
+  inline void operator *= (cvm::real a)
   {
     size_t i;
     for (i = 0; i < data.size(); i++) {
@@ -547,7 +564,7 @@ public:
     }
   }
 
-  inline void operator /= (cvm::real const &a)
+  inline void operator /= (cvm::real a)
   {
     size_t i;
     for (i = 0; i < data.size(); i++) {
@@ -555,7 +572,8 @@ public:
     }
   }
 
-  inline friend matrix2d<T> operator + (matrix2d<T> const &m1, matrix2d<T> const &m2)
+  inline friend matrix2d<T> operator + (matrix2d<T> const &m1,
+                                        matrix2d<T> const &m2)
   {
     check_sizes(m1, m2);
     matrix2d<T> result(m1.outer_length, m1.inner_length);
@@ -566,7 +584,8 @@ public:
     return result;
   }
 
-  inline friend matrix2d<T> operator - (matrix2d<T> const &m1, matrix2d<T> const &m2)
+  inline friend matrix2d<T> operator - (matrix2d<T> const &m1,
+                                        matrix2d<T> const &m2)
   {
     check_sizes(m1, m2);
     matrix2d<T> result(m1.outer_length, m1.inner_length);
@@ -577,7 +596,7 @@ public:
     return result;
   }
 
-  inline friend matrix2d<T> operator * (matrix2d<T> const &m, cvm::real const &a)
+  inline friend matrix2d<T> operator * (matrix2d<T> const &m, cvm::real a)
   {
     matrix2d<T> result(m.outer_length, m.inner_length);
     size_t i;
@@ -587,12 +606,12 @@ public:
     return result;
   }
 
-  inline friend matrix2d<T> operator * (cvm::real const &a, matrix2d<T> const &m)
+  inline friend matrix2d<T> operator * (cvm::real a, matrix2d<T> const &m)
   {
     return m * a;
   }
 
-  inline friend matrix2d<T> operator / (matrix2d<T> const &m, cvm::real const &a)
+  inline friend matrix2d<T> operator / (matrix2d<T> const &m, cvm::real a)
   {
     matrix2d<T> result(m.outer_length, m.inner_length);
     size_t i;
@@ -602,34 +621,17 @@ public:
     return result;
   }
 
-  /// Matrix multiplication
-//   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)
+  inline friend vector1d<T> operator * (vector1d<T> const &v,
+                                        matrix2d<T> const &m)
   {
     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");
+      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, k;
       for (i = 0; i < m.inner_length; i++) {
@@ -641,25 +643,6 @@ public:
     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)
@@ -725,49 +708,52 @@ public:
   cvm::real x, y, z;
 
   inline rvector()
-    : x(0.0), y(0.0), z(0.0)
-  {}
+  {
+    reset();
+  }
 
-  inline rvector(cvm::real const &x_i,
-                 cvm::real const &y_i,
-                 cvm::real const &z_i)
-    : x(x_i), y(y_i), z(z_i)
-  {}
+  /// \brief Set all components to zero
+  inline void reset()
+  {
+    set(0.0);
+  }
+
+  inline rvector(cvm::real x_i, cvm::real y_i, cvm::real z_i)
+  {
+    set(x_i, y_i, z_i);
+  }
 
   inline rvector(cvm::vector1d<cvm::real> const &v)
-    : x(v[0]), y(v[1]), z(v[2])
-  {}
+  {
+    set(v[0], v[1], v[2]);
+  }
 
   inline rvector(cvm::real t)
-    : x(t), y(t), z(t)
-  {}
+  {
+    set(t);
+  }
 
-  /// \brief Set all components to a scalar value
-  inline void set(cvm::real const &value) {
+  /// \brief Set all components to a scalar
+  inline void set(cvm::real value)
+  {
     x = y = z = value;
   }
 
   /// \brief Assign all components
-  inline void set(cvm::real const &x_i,
-                  cvm::real const &y_i,
-                  cvm::real const &z_i) {
+  inline void set(cvm::real x_i, cvm::real y_i, cvm::real z_i)
+  {
     x = x_i;
     y = y_i;
     z = z_i;
   }
 
-  /// \brief Set all components to zero
-  inline void reset() {
-    x = y = z = 0.0;
-  }
-
   /// \brief Access cartesian components by index
-  inline cvm::real & operator [] (int const &i) {
+  inline cvm::real & operator [] (int i) {
     return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : x;
   }
 
   /// \brief Access cartesian components by index
-  inline cvm::real const & operator [] (int const &i) const {
+  inline cvm::real  operator [] (int i) const {
     return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : x;
   }
 
@@ -780,14 +766,6 @@ public:
     return result;
   }
 
-  inline cvm::rvector & operator = (cvm::real const &v)
-  {
-    x = v;
-    y = v;
-    z = v;
-    return *this;
-  }
-
   inline void operator += (cvm::rvector const &v)
   {
     x += v.x;
@@ -802,7 +780,7 @@ public:
     z -= v.z;
   }
 
-  inline void operator *= (cvm::real const &v)
+  inline void operator *= (cvm::real v)
   {
     x *= v;
     y *= v;
@@ -832,13 +810,14 @@ public:
     return (n > 0. ? cvm::rvector(x, y, z)/n : cvm::rvector(1., 0., 0.));
   }
 
-  static inline size_t output_width(size_t const &real_width)
+  static inline size_t output_width(size_t real_width)
   {
     return 3*real_width + 10;
   }
 
 
-  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,
@@ -850,41 +829,35 @@ public:
     return cvm::rvector(-v.x, -v.y, -v.z);
   }
 
-  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)
-  {
-    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)
+  /// Inner (dot) product
+  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 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 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 a)
   {
     return cvm::rvector(v.x/a, v.y/a, v.z/a);
   }
@@ -946,15 +919,9 @@ public:
   {}
 
   /// Constructor component by component
-  inline rmatrix(cvm::real const &xxi,
-                 cvm::real const &xyi,
-                 cvm::real const &xzi,
-                 cvm::real const &yxi,
-                 cvm::real const &yyi,
-                 cvm::real const &yzi,
-                 cvm::real const &zxi,
-                 cvm::real const &zyi,
-                 cvm::real const &zzi)
+  inline rmatrix(cvm::real xxi, cvm::real xyi, cvm::real xzi,
+                 cvm::real yxi, cvm::real yyi, cvm::real yzi,
+                 cvm::real zxi, cvm::real zyi, cvm::real zzi)
     : cvm::matrix2d<cvm::real>(3, 3)
   {
     this->xx() = xxi;
@@ -983,31 +950,13 @@ public:
 
   inline cvm::rmatrix transpose() const
   {
-    return cvm::rmatrix(this->xx(),
-                        this->yx(),
-                        this->zx(),
-                        this->xy(),
-                        this->yy(),
-                        this->zy(),
-                        this->xz(),
-                        this->yz(),
-                        this->zz());
+    return cvm::rmatrix(this->xx(), this->yx(), this->zx(),
+                        this->xy(), this->yy(), this->zy(),
+                        this->xz(), this->yz(), this->zz());
   }
 
   friend cvm::rvector operator * (cvm::rmatrix const &m, cvm::rvector const &r);
 
-  //   friend inline cvm::rmatrix const operator * (cvm::rmatrix const &m1, cvm::rmatrix const &m2) {
-  //     return cvm::rmatrix (m1.xx()*m2.xx() + m1.xy()*m2.yx() + m1.xz()*m2.yz(),
-  //                     m1.xx()*m2.xy() + m1.xy()*m2.yy() + m1.xz()*m2.zy(),
-  //                     m1.xx()*m2.xz() + m1.xy()*m2.yz() + m1.xz()*m2.zz(),
-  //                     m1.yx()*m2.xx() + m1.yy()*m2.yx() + m1.yz()*m2.yz(),
-  //                     m1.yx()*m2.xy() + m1.yy()*m2.yy() + m1.yz()*m2.yy(),
-  //                     m1.yx()*m2.xz() + m1.yy()*m2.yz() + m1.yz()*m2.yz(),
-  //                     m1.zx()*m2.xx() + m1.zy()*m2.yx() + m1.zz()*m2.yz(),
-  //                     m1.zx()*m2.xy() + m1.zy()*m2.yy() + m1.zz()*m2.yy(),
-  //                     m1.zx()*m2.xz() + m1.zy()*m2.yz() + m1.zz()*m2.yz());
-  //   }
-
 };
 
 
@@ -1031,7 +980,7 @@ public:
   cvm::real q0, q1, q2, q3;
 
   /// Constructor from a 3-d vector
-  inline quaternion(cvm::real const &x, cvm::real const &y, cvm::real const &z)
+  inline quaternion(cvm::real x, cvm::real y, cvm::real z)
     : q0(0.0), q1(x), q2(y), q3(z)
   {}
 
@@ -1041,10 +990,10 @@ public:
   {}
 
   /// Constructor component by component
-  inline quaternion(cvm::real const &q0i,
-                    cvm::real const &q1i,
-                    cvm::real const &q2i,
-                    cvm::real const &q3i)
+  inline quaternion(cvm::real q0i,
+                    cvm::real q1i,
+                    cvm::real q2i,
+                    cvm::real q3i)
     : q0(q0i), q1(q1i), q2(q2i), q3(q3i)
   {}
 
@@ -1055,9 +1004,9 @@ public:
   /// "Constructor" after Euler angles (in radians)
   ///
   /// http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
-  inline void set_from_euler_angles(cvm::real const &phi_in,
-                                    cvm::real const &theta_in,
-                                    cvm::real const &psi_in)
+  inline void set_from_euler_angles(cvm::real phi_in,
+                                    cvm::real theta_in,
+                                    cvm::real psi_in)
   {
     q0 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) +
            (std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
@@ -1079,7 +1028,7 @@ public:
   }
 
   /// \brief Set all components to a scalar
-  inline void set(cvm::real const &value = 0.0)
+  inline void set(cvm::real value)
   {
     q0 = q1 = q2 = q3 = value;
   }
@@ -1087,7 +1036,7 @@ public:
   /// \brief Set all components to zero (null quaternion)
   inline void reset()
   {
-    q0 = q1 = q2 = q3 = 0.0;
+    set(0.0);
   }
 
   /// \brief Set the q0 component to 1 and the others to 0 (quaternion
@@ -1099,7 +1048,7 @@ public:
   }
 
   /// Tell the number of characters required to print a quaternion, given that of a real number
-  static inline size_t output_width(size_t const &real_width)
+  static inline size_t output_width(size_t real_width)
   {
     return 4*real_width + 13;
   }
@@ -1113,7 +1062,7 @@ public:
   friend std::istream & operator >> (std::istream &is, cvm::quaternion &q);
 
   /// Access the quaternion as a 4-d array (return a reference)
-  inline cvm::real & operator [] (int const &i) {
+  inline cvm::real & operator [] (int i) {
     switch (i) {
     case 0:
       return this->q0;
@@ -1130,7 +1079,7 @@ public:
   }
 
   /// Access the quaternion as a 4-d array (return a value)
-  inline cvm::real operator [] (int const &i) const {
+  inline cvm::real operator [] (int i) const {
     switch (i) {
     case 0:
       return this->q0;
@@ -1175,12 +1124,12 @@ public:
     return cvm::quaternion(q0, -q1, -q2, -q3);
   }
 
-  inline void operator *= (cvm::real const &a)
+  inline void operator *= (cvm::real a)
   {
     q0 *= a; q1 *= a; q2 *= a; q3 *= a;
   }
 
-  inline void operator /= (cvm::real const &a)
+  inline void operator /= (cvm::real a)
   {
     q0 /= a; q1 /= a; q2 /= a; q3 /= a;
   }
@@ -1215,19 +1164,22 @@ public:
   }
 
 
-  friend inline cvm::quaternion operator + (cvm::quaternion const &h, cvm::quaternion const &q)
+  friend inline cvm::quaternion operator + (cvm::quaternion const &h,
+                                            cvm::quaternion const &q)
   {
     return cvm::quaternion(h.q0+q.q0, h.q1+q.q1, h.q2+q.q2, h.q3+q.q3);
   }
 
-  friend inline cvm::quaternion operator - (cvm::quaternion const &h, cvm::quaternion const &q)
+  friend inline cvm::quaternion operator - (cvm::quaternion const &h,
+                                            cvm::quaternion const &q)
   {
     return cvm::quaternion(h.q0-q.q0, h.q1-q.q1, h.q2-q.q2, h.q3-q.q3);
   }
 
   /// \brief Provides the quaternion product.  \b NOTE: for the inner
   /// product use: \code h.inner (q); \endcode
-  friend inline cvm::quaternion operator * (cvm::quaternion const &h, cvm::quaternion const &q)
+  friend inline cvm::quaternion operator * (cvm::quaternion const &h,
+                                            cvm::quaternion const &q)
   {
     return cvm::quaternion(h.q0*q.q0 - h.q1*q.q1 - h.q2*q.q2 - h.q3*q.q3,
                            h.q0*q.q1 + h.q1*q.q0 + h.q2*q.q3 - h.q3*q.q2,
@@ -1235,18 +1187,18 @@ public:
                            h.q0*q.q3 + h.q3*q.q0 + h.q1*q.q2 - h.q2*q.q1);
   }
 
-  friend inline cvm::quaternion operator * (cvm::real const &c,
+  friend inline cvm::quaternion operator * (cvm::real c,
                                             cvm::quaternion const &q)
   {
     return cvm::quaternion(c*q.q0, c*q.q1, c*q.q2, c*q.q3);
   }
   friend inline cvm::quaternion operator * (cvm::quaternion const &q,
-                                            cvm::real const &c)
+                                            cvm::real c)
   {
     return cvm::quaternion(q.q0*c, q.q1*c, q.q2*c, q.q3*c);
   }
   friend inline cvm::quaternion operator / (cvm::quaternion const &q,
-                                            cvm::real const &c)
+                                            cvm::real c)
   {
     return cvm::quaternion(q.q0/c, q.q1/c, q.q2/c, q.q3/c);
   }
@@ -1407,7 +1359,7 @@ public:
   std::vector< cvm::vector1d<cvm::rvector> > dQ0_1, dQ0_2;
 
   /// Allocate space for the derivatives of the rotation
-  inline void request_group1_gradients(size_t const &n)
+  inline void request_group1_gradients(size_t n)
   {
     dS_1.resize(n, cvm::matrix2d<cvm::rvector>(4, 4));
     dL0_1.resize(n, cvm::rvector(0.0, 0.0, 0.0));
@@ -1415,7 +1367,7 @@ public:
   }
 
   /// Allocate space for the derivatives of the rotation
-  inline void request_group2_gradients(size_t const &n)
+  inline void request_group2_gradients(size_t n)
   {
     dS_2.resize(n, cvm::matrix2d<cvm::rvector>(4, 4));
     dL0_2.resize(n, cvm::rvector(0.0, 0.0, 0.0));
@@ -1448,7 +1400,7 @@ public:
   }
 
   /// Constructor after an axis of rotation and an angle (in radians)
-  inline rotation(cvm::real const &angle, cvm::rvector const &axis)
+  inline rotation(cvm::real angle, cvm::rvector const &axis)
     : b_debug_gradients(false)
   {
     cvm::rvector const axis_n = axis.unit();
@@ -1500,20 +1452,18 @@ 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)));
 
-      cvm::real const dspindx = (180.0/PI) * 2.0 * (1.0 / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
-
-      return
-        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));
+      return 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));
     } else {
       // (1/(1+x^2)) ~ (1/x)^2
-      return
-        cvm::quaternion((180.0/PI) * 2.0 * ((-1.0)/iprod), 0.0, 0.0, 0.0);
-      // XX TODO: What if iprod == 0? XX
+      // The documentation of spinAngle discourages its use when q_vec and
+      // axis are not close
+      return cvm::quaternion((180.0/PI) * 2.0 * ((-1.0)/iprod), 0.0, 0.0, 0.0);
     }
   }
 
diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp
index c5b9e5a60c..ff74602b41 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.cpp
+++ b/src/USER-COLVARS/colvarproxy_lammps.cpp
@@ -133,6 +133,8 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
 
 void colvarproxy_lammps::init(const char *conf_file)
 {
+  version_int = get_version_from_string(COLVARPROXY_VERSION);
+
   // create the colvarmodule instance
   colvars = new colvarmodule(this);
 
@@ -203,6 +205,25 @@ double colvarproxy_lammps::compute()
   }
   previous_step = _lmp->update->ntimestep;
 
+  unit_cell_x.set(_lmp->domain->xprd, 0.0, 0.0);
+  unit_cell_y.set(0.0, _lmp->domain->yprd, 0.0);
+  unit_cell_z.set(0.0, 0.0, _lmp->domain->zprd);
+
+  if (_lmp->domain->xperiodic == 0 && _lmp->domain->yperiodic == 0 &&
+      _lmp->domain->zperiodic == 0) {
+    boundaries_type = boundaries_non_periodic;
+    reset_pbc_lattice();
+  } else if ((_lmp->domain->nonperiodic == 0) &&
+             (_lmp->domain->dimension == 3) &&
+             (_lmp->domain->triclinic == 0)) {
+    // Orthogonal unit cell
+    boundaries_type = boundaries_pbc_ortho;
+    colvarproxy_system::update_pbc_lattice();
+    // It is safer to let LAMMPS deal with high-tilt triclinic boxes
+  } else {
+    boundaries_type = boundaries_unsupported;
+  }
+
   if (cvm::debug()) {
     cvm::log(std::string(cvm::line_marker)+
              "colvarproxy_lammps, step no. "+cvm::to_str(colvars->it)+"\n"+
@@ -263,8 +284,10 @@ bool colvarproxy_lammps::deserialize_status(std::string &rst)
   }
 }
 
+
 cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1,
                                                    cvm::atom_pos const &pos2)
+  const
 {
   double xtmp = pos2.x - pos1.x;
   double ytmp = pos2.y - pos1.y;
@@ -273,28 +296,6 @@ cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1,
   return cvm::rvector(xtmp, ytmp, ztmp);
 }
 
-cvm::real colvarproxy_lammps::position_dist2(cvm::atom_pos const &pos1,
-                                             cvm::atom_pos const &pos2)
-{
-  double xtmp = pos2.x - pos1.x;
-  double ytmp = pos2.y - pos1.y;
-  double ztmp = pos2.z - pos1.z;
-  _lmp->domain->minimum_image(xtmp,ytmp,ztmp);
-  return cvm::real(xtmp*xtmp + ytmp*ytmp + ztmp*ztmp);
-}
-
-
-void colvarproxy_lammps::select_closest_image(cvm::atom_pos &pos,
-                                              cvm::atom_pos const &ref)
-{
-  double xtmp = pos.x - ref.x;
-  double ytmp = pos.y - ref.y;
-  double ztmp = pos.z - ref.z;
-  _lmp->domain->minimum_image(xtmp,ytmp,ztmp);
-  pos.x = ref.x + xtmp;
-  pos.y = ref.y + ytmp;
-  pos.z = ref.z + ztmp;
-}
 
 void colvarproxy_lammps::log(std::string const &message)
 {
@@ -308,12 +309,14 @@ void colvarproxy_lammps::log(std::string const &message)
   }
 }
 
+
 void colvarproxy_lammps::error(std::string const &message)
 {
   // In LAMMPS, all errors are fatal
   fatal_error(message);
 }
 
+
 void colvarproxy_lammps::fatal_error(std::string const &message)
 {
   log(message);
@@ -321,13 +324,6 @@ void colvarproxy_lammps::fatal_error(std::string const &message)
                    "Fatal error in the collective variables module.\n");
 }
 
-void colvarproxy_lammps::exit(std::string const &message)
-{
-  log(message);
-  log("Request to exit the simulation made.\n");
-  do_exit=true;
-}
-
 
 int colvarproxy_lammps::backup_file(char const *filename)
 {
@@ -342,10 +338,12 @@ int colvarproxy_lammps::backup_file(char const *filename)
 
 // multi-replica support
 
-void colvarproxy_lammps::replica_comm_barrier() {
+void colvarproxy_lammps::replica_comm_barrier()
+{
   MPI_Barrier(inter_comm);
 }
 
+
 int colvarproxy_lammps::replica_comm_recv(char* msg_data,
                                           int buf_len, int src_rep)
 {
@@ -359,6 +357,7 @@ int colvarproxy_lammps::replica_comm_recv(char* msg_data,
   return retval;
 }
 
+
 int colvarproxy_lammps::replica_comm_send(char* msg_data,
                                           int msg_len, int dest_rep)
 {
diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h
index af2aa04dfc..cdd86cbd16 100644
--- a/src/USER-COLVARS/colvarproxy_lammps.h
+++ b/src/USER-COLVARS/colvarproxy_lammps.h
@@ -120,14 +120,9 @@ class colvarproxy_lammps : public colvarproxy {
   void log(std::string const &message);
   void error(std::string const &message);
   void fatal_error(std::string const &message);
-  void exit(std::string const &message);
 
   cvm::rvector position_distance(cvm::atom_pos const &pos1,
-                                 cvm::atom_pos const &pos2);
-  cvm::real position_dist2(cvm::atom_pos const &pos1,
-                           cvm::atom_pos const &pos2);
-  void select_closest_image(cvm::atom_pos &pos,
-                            cvm::atom_pos const &ref_pos);
+                                 cvm::atom_pos const &pos2) const;
 
   int backup_file(char const *filename);
 
diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h
index 79f77dad4e..0d451212b2 100644
--- a/src/USER-COLVARS/colvarproxy_lammps_version.h
+++ b/src/USER-COLVARS/colvarproxy_lammps_version.h
@@ -1,5 +1,5 @@
 #ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2017-12-01"
+#define COLVARPROXY_VERSION "2018-04-29"
 // This file is part of the Collective Variables module (Colvars).
 // The original version of Colvars and its updates are located at:
 // https://github.com/colvars/colvars