Update Colvars to version 2018-03-14 49/3849/1
authorGiacomo Fiorin <giacomo.fiorin@gmail.com>
Wed, 14 Mar 2018 18:04:45 +0000 (14:04 -0400)
committerGiacomo Fiorin <giacomo.fiorin@gmail.com>
Wed, 14 Mar 2018 18:04:45 +0000 (14:04 -0400)
This update fixes an issue reported by Marcelo Melo with the order of atoms in
Colvars groups being given manually by number, but out of order.  The intended
behavior is that Colvars atom groups can be defined out-of-order, but
coordinates are read from files in-order.  This fix introduces a back-mapping
to the original order, retaining backward compatibility as much as possible.

As stated already to various UIUC members in private, given the increased
frequent use of reinitatoms (e.g. in constant pH) I'd really like to phase out
most selection commands internal to Colvars.  These "selections" should be
managed centrally by NAMD, kept updated between topology changes, and made
available to the objects as Tcl procs.  These may be VMD-like selections (a
long standing issue) but they don't have to.

Below are the relevant commits from https://github.com/Colvars/colvars/

9b85d5f 2018-03-13 Fix issue with out-of-order atom selections; clarify format for ref positions [Giacomo Fiorin]

Change-Id: Ia84f3241b0b41227ae05ae784c9de8eed83cd9e0

18 files changed:
colvars/src/colvar.cpp
colvars/src/colvaratoms.cpp
colvars/src/colvaratoms.h
colvars/src/colvarbias_abf.cpp
colvars/src/colvarcomp_distances.cpp
colvars/src/colvarcomp_rotations.cpp
colvars/src/colvargrid.cpp
colvars/src/colvarmodule.cpp
colvars/src/colvarmodule.h
colvars/src/colvarproxy.cpp
colvars/src/colvarproxy.h
colvars/src/colvars_version.h
colvars/src/colvartypes.h
lib/abf_integrate/abf_data.h
src/colvarproxy_namd.C
src/colvarproxy_namd.h
src/colvarproxy_namd_version.h
ug/ug_colvars.tex

index 5a4e8b6..142517f 100644 (file)
@@ -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"
 #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");
index 1be6f42..c8dcfc0 100644 (file)
@@ -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);
-  }
-  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;
+  // 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]);
   }
-  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++;
+  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);
   }
-  return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
+
+  // 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 COLVARS_OK;
 }
 
 
index 0dda6ab..93ae594 100644 (file)
@@ -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;
+  }
+
+  std::string const print_atom_ids() const;
+
+  /// Allocates and populates sorted_ids and sorted_ids_map
+  int create_sorted_ids();
 
-  /// Allocates and populates the sorted list of atom ids
-  int create_sorted_ids(void);
+  /// 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
index b3b5b3e..771955a 100644 (file)
@@ -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.
index 9911f4c..081ec4a 100644 (file)
@@ -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();
 }
 
 
@@ -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);
     }
   }
 
index 2650a9f..498ef7c 100644 (file)
@@ -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);
     }
   }
 
index 1ac4aae..407b336 100644 (file)
@@ -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];
     }
index ded32e3..82bc72c 100644 (file)
@@ -432,13 +432,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 +452,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 +472,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();
@@ -1694,40 +1694,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::vector<cvm::rvector> *pos,
+                     cvm::atom_group *atoms,
                      std::string const &pdb_field,
-                     double const pdb_field_value)
+                     double pdb_field_value)
 {
-  // Differentiate between PDB and XYZ files
-  // for XYZ files, use CVM internal parser
-  // otherwise call proxy function for PDB
+  int error_code = COLVARS_OK;
 
-  std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name);
+  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;
@@ -1742,12 +1761,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);
@@ -1757,7 +1776,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];
     }
@@ -1798,17 +1817,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)
 {
index 8afdabc..f3d99a2 100644 (file)
@@ -311,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)
@@ -484,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();
@@ -500,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;
 
@@ -518,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);
-
-  /// \brief Load the coordinates for a group of atoms from a file
-  /// (PDB or XYZ)
+                        double pdb_field_value = 0.0);
+
+  /// \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;
index a5f50b3..079c8a3 100644 (file)
 
 
 
-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 std::floor(x+0.5);
+}
+
+
+void colvarproxy_system::update_pbc_lattice()
 {
-  return (position_distance(pos1, pos2)).norm2();
+  // 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",
index 5fbc66f..845f93c 100644 (file)
@@ -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);
-
-  /// \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"
+                         double pdb_field_value = 0.0);
+
+  /// \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();
index 7ef89ef..8873b80 100644 (file)
@@ -1,5 +1,5 @@
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2018-02-24"
+#define COLVARS_VERSION "2018-03-14"
 // 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
index 97257d1..4ef9557 100644 (file)
@@ -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);
     }
   }
 
index 7ff1cc5..e082979 100644 (file)
@@ -57,7 +57,7 @@ class ABFdata {
     int *PBC;
 
     /// Constructor: reads from a file
-     ABFdata(const char *gradFileName);
+    explicit ABFdata(const char *gradFileName);
     ~ABFdata();
 
     /// \brief Returns an offset for scalar fields based on a n-index.
index b475cef..6b323b9 100644 (file)
@@ -273,6 +273,29 @@ void colvarproxy_namd::calculate()
 
   previous_NAMD_step = step;
 
+  {
+    Vector const a = lattice->a();
+    Vector const b = lattice->b();
+    Vector const c = lattice->c();
+    unit_cell_x.set(a.x, a.y, a.z);
+    unit_cell_y.set(b.x, b.y, c.z);
+    unit_cell_z.set(c.x, c.y, c.z);
+  }
+
+  if (!lattice->a_p() && !lattice->b_p() && !lattice->c_p()) {
+    boundaries_type = boundaries_non_periodic;
+    reset_pbc_lattice();
+  } else if (lattice->a_p() && lattice->b_p() && lattice->c_p()) {
+    if (lattice->orthogonal()) {
+      boundaries_type = boundaries_pbc_ortho;
+    } else {
+      boundaries_type = boundaries_pbc_triclinic;
+    }
+    colvarproxy_system::update_pbc_lattice();
+  } else {
+    boundaries_type = boundaries_unsupported;
+  }
+
   if (cvm::debug()) {
     log(std::string(cvm::line_marker)+
         "colvarproxy_namd, step no. "+cvm::to_str(colvars->it)+"\n"+
@@ -654,6 +677,19 @@ void colvarproxy_namd::clear_atom(int index)
 }
 
 
+cvm::rvector colvarproxy_namd::position_distance(cvm::atom_pos const &pos1,
+                                                 cvm::atom_pos const &pos2)
+  const
+{
+  Position const p1(pos1.x, pos1.y, pos1.z);
+  Position const p2(pos2.x, pos2.y, pos2.z);
+  // return p2 - p1
+  Vector const d = this->lattice->delta(p2, p1);
+  return cvm::rvector(d.x, d.y, d.z);
+}
+
+
+
 enum e_pdb_field {
   e_pdb_none,
   e_pdb_occ,
index 0fb11f9..15279ee 100644 (file)
@@ -75,7 +75,8 @@ public:
   int reset();
 
   // synchronize the local arrays with requested or forced atoms
-  int update_atoms_map(AtomIDList::const_iterator begin, AtomIDList::const_iterator end);
+  int update_atoms_map(AtomIDList::const_iterator begin,
+                       AtomIDList::const_iterator end);
 
   void calculate();
 
@@ -253,13 +254,7 @@ public:
   }
 
   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 load_atoms(char const *filename,
                  cvm::atom_group &atoms,
@@ -291,40 +286,4 @@ public:
 };
 
 
-inline cvm::rvector colvarproxy_namd::position_distance(cvm::atom_pos const &pos1,
-                                                        cvm::atom_pos const &pos2)
-{
-  Position const p1(pos1.x, pos1.y, pos1.z);
-  Position const p2(pos2.x, pos2.y, pos2.z);
-  // return p2 - p1
-  Vector const d = this->lattice->delta(p2, p1);
-  return cvm::rvector(d.x, d.y, d.z);
-}
-
-
-inline cvm::real colvarproxy_namd::position_dist2(cvm::atom_pos const &pos1,
-                                                  cvm::atom_pos const &pos2)
-{
-  Lattice const *l = this->lattice;
-  Vector const p1(pos1.x, pos1.y, pos1.z);
-  Vector const p2(pos2.x, pos2.y, pos2.z);
-  Vector const d = l->delta(p1, p2);
-  return cvm::real(d.x*d.x + d.y*d.y + d.z*d.z);
-}
-
-
-inline void colvarproxy_namd::select_closest_image(cvm::atom_pos &pos,
-                                                   cvm::atom_pos const &ref_pos)
-{
-  Position const p(pos.x, pos.y, pos.z);
-  Position const rp(ref_pos.x, ref_pos.y, ref_pos.z);
-  ScaledPosition const srp = this->lattice->scale(rp);
-  Position const np = this->lattice->nearest(p, srp);
-  pos.x = np.x;
-  pos.y = np.y;
-  pos.z = np.z;
-}
-
-
-
 #endif
index fd8860f..70d6e04 100644 (file)
@@ -1,5 +1,5 @@
 #ifndef COLVARPROXY_VERSION
-#define COLVARPROXY_VERSION "2018-02-24"
+#define COLVARPROXY_VERSION "2018-03-14"
 // 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
index 64ebe53..50f32d9 100644 (file)
@@ -1017,7 +1017,7 @@ The complete list of selection keywords available in \MDENGINE{} is:
     This option sets the PSF segment identifier for
     \texttt{atomNameResidueRange}.  Multiple values may be provided,
     which correspond to multiple instances of
-    \texttt{atomNameResidueRange}, in the order of their occurrence.
+    \texttt{atomNameResidueRange}, in order of their occurrence.
     This option is only necessary if a PSF topology file is used.}
 
 \item %
@@ -1098,38 +1098,42 @@ At this point, only the RMSD component will correctly adjust the Jacobian deriva
     boolean}{%
     \texttt{off}}{%
     If this option is \texttt{on}, the coordinates of this group will be optimally superimposed to the reference positions provided by \cvnamebasedonly{either} \texttt{refPositions} or \texttt{refPositionsFile}.
-    The rotation will be performed around the center of geometry if \texttt{centerReference} is \texttt{on}, around the origin otherwise.
+    The rotation will be performed around the center of geometry if \texttt{centerReference} is \texttt{on}, or around the origin otherwise.
     The algorithm used is the same employed by the \texttt{orientation} colvar component~\cite{Coutsias2004}.
     Forces applied to the atoms of this group will also be implicitly rotated back to the original frame.
     \textbf{Note}: unless otherwise specified, \texttt{rmsd} and \texttt{eigenvector} set this option to \texttt{on} \emph{by default}.
 }
 
 \item %
+  \labelkey{atom-group|refPositions}
   \key
     {refPositions}{%
     atom group}{%
     Reference positions for fitting (\AA)}{%
     space-separated list of \texttt{(x, y, z)} triplets}{%
     \label{key:colvars:atom_group:refPositions}
-    This option provides a list of reference coordinates for \texttt{centerReference} or \texttt{rotateReference}.
-    If only \texttt{centerReference} is \texttt{on}, the list may contain a single (x, y, z) triplet; if also \texttt{rotateReference} is \texttt{on}, the list should be as long as the atom group.
+    This option provides a list of reference coordinates for \texttt{centerReference} and/or \texttt{rotateReference}, and is mutually exclusive with \texttt{refPositionsFile}.
+    If only \texttt{centerReference} is \texttt{on}, the list may contain a single (x, y, z) triplet; if also \texttt{rotateReference} is \texttt{on}, the list should be as long as the atom group, and \emph{its order must match the order in which atoms were defined}.
 }
 
 \item %
+  \labelkey{atom-group|refPositionsFile}
   \key
     {refPositionsFile}{%
     atom group}{%
     File containing the reference positions for fitting}{%
     UNIX filename}{%
     \label{key:colvars:atom_group:refPositionsFile}
-    This keyword provides the reference coordinates for fitting from the given file, and is mutually exclusive with \texttt{refPositions}.
-    The acceptable file format is XYZ\cvnamebasedonly{or PDB}.
-    Atomic positions are read differently depending on the three following scenarios:
-    \emph{i)} the file contains exactly as many records as the atoms in the group: all positions are read in sequence;
-    \emph{ii)} the file contains coordinates for the entire system: only the positions corresponding to the numeric indices of the atom group are read\cvnamebasedonly{;
-    \emph{iii)} if the file is a PDB file and \texttt{refPositionsCol} is specified, positions are read according to the value of the column \texttt{refPositionsCol} (which may be the same as \texttt{atomsCol})}.
+    This option provides a list of reference coordinates for \texttt{centerReference} and/or \texttt{rotateReference}, and is mutually exclusive with \texttt{refPositions}.
+    The acceptable file format is XYZ, which is read in double precision\cvnamebasedonly{, or PDB; \emph{the latter is discouraged if the precision of the reference coordinates is a concern}}.
+    Atomic positions are read differently depending on the following scenarios:
+    \textbf{(i)} the file contains exactly as many records as the atoms in the group: all positions are read in sequence;
+    \textbf{(ii)} (most common case) the file contains coordinates for the entire system: only the positions corresponding to the numeric indices of the atom group are read\cvnamebasedonly{;
+      \textbf{(iii)} if the file is a PDB file and \texttt{refPositionsCol} is specified, positions are read according to the value of the column \texttt{refPositionsCol} (which may be the same as \texttt{atomsCol})}.
+    In each case, atoms are read from the file \emph{in order of increasing number}.
 }
 
+
 \cvnamebasedonly{
 \item %
   \key
@@ -1211,8 +1215,8 @@ The following two options have default values appropriate for the vast majority
     Apply forces from this colvar to this group}{%
     boolean}{%
     \texttt{on}}{%
-    If this option is \texttt{off}, no forces are applied from this
-    colvar to this group. Other forces are not affected (i.e. those
+    If this option is \texttt{off}, no forces are applied the atoms in the group.
+    Other forces are not affected (i.e. those
     from the MD engine, from other colvars, and other external forces).
     For dummy atoms, this option is \texttt{off} by default.
  }
@@ -1239,8 +1243,8 @@ The following two options have default values appropriate for the vast majority
 }
 \cvlammpsonly{
 In simulations with periodic boundary conditions, many of the implemented colvar components rely on the fact that each position within a group of atoms is at the nearest periodic image from the center of geometry of the group itself.
-However, due to the internal wrapping of individual atomic positions done by LAMMPS, this assumption is inaccurate if groups lies at one of the unit cell's boundaries.
-For this reason, within the Colvars module coordinates are unwrapped by default to avoid discontinuities due to coordinate wrapping (see \texttt{unwrap} keyword in \ref{sec:colvars_mdengine_parameters}).
+However, due to the internal wrapping of individual atomic positions done by LAMMPS, this assumption is broken if the group straddles one of the unit cell's boundaries.
+For this reason, within the Colvars module all coordinates are unwrapped by default to avoid discontinuities (see \texttt{unwrap} keyword in \ref{sec:colvars_mdengine_parameters}).
 
 The user should determine whether maintaining the default value of \texttt{unwrap}, depending on the specifics of each system.
 In general, internal coordinate wrapping by LAMMPS does not affect the calculation of colvars if each atom group satisfies one or more of the following:
@@ -1962,11 +1966,8 @@ with respect to $\{\mathbf{x}_{i}(t)\}$.
     \texttt{rmsd}}{%
     Reference coordinates}{%
     space-separated list of \texttt{(x, y, z)} triplets}{%
-    This option (mutually exclusive with \texttt{refPositionsFile}) sets the reference coordinates.  If only \texttt{centerReference} is \texttt{on}, the list can be a single (x, y, z) triplet; if also \texttt{rotateReference} is \texttt{on}, the list should be as long as the atom group.  This option 
-    is independent from that with the same keyword within the
-    \texttt{atoms~\{...\}} block (see \ref{key:colvars:atom_group:refPositions}).  The latter (and related fitting
-    options for the atom group) are normally not needed,
-    and should be omitted altogether except for advanced usage cases.
+    This option (mutually exclusive with \texttt{refPositionsFile}) sets the reference coordinates for RMSD calculation, and uses these to compute the roto-translational fit.  
+    It is functionally equivalent to the option \refkey{refPositions}{atom-group|refPositions} in the atom group definition, which also supports more advanced fitting options.
     }
 
 \item %
@@ -1976,7 +1977,8 @@ with respect to $\{\mathbf{x}_{i}(t)\}$.
     \texttt{rmsd}}{%
     Reference coordinates file}{%
     UNIX filename}{%
-    This option (mutually exclusive with \texttt{refPositions}) sets the file name for the reference coordinates to be compared with.  The format is the same as that provided by \texttt{refPositionsFile} within an atom group's definition (see \ref{key:colvars:atom_group:refPositionsFile}).
+    This option (mutually exclusive with \texttt{refPositions}) sets the reference coordinates for RMSD calculation, and uses these to compute the roto-translational fit.  
+    It is functionally equivalent to the option \refkey{refPositionsFile}{atom-group|refPositionsFile} in the atom group definition, which also supports more advanced fitting options.
     }
 
 \cvnamebasedonly{
@@ -2086,23 +2088,20 @@ automatically when reading them from the configuration.
     \texttt{eigenvector}}{%
     Vector components}{%
     space-separated list of \texttt{(x, y, z)} triplets}{%
-    This option \cvnamebasedonly{(mutually exclusive with \texttt{vectorFile})} sets the values of the vector components.}
+    This option (mutually exclusive with \texttt{vectorFile}) sets the values of the vector components.}
 
-\cvnamebasedonly{
 \item %
   \key
     {vectorFile}{%
     \texttt{eigenvector}}{%
-    PDB file containing vector components}{%
+    file containing vector components}{%
     UNIX filename}{%
-    This option (mutually exclusive with \texttt{vector}) sets the
-    name of a PDB file where the vector components will be read from the
-    X, Y, and Z fields.
-    \textbf{Note:} \emph{The PDB file has limited precision and fixed
-      point numbers: in some cases, the vector may not be
-      accurately represented, and }\texttt{vector}\emph{ should be
-      used instead.}}
+    This option (mutually exclusive with \texttt{vector}) sets the name of a coordinate file containing the vector components; the file is read according to the same format used for \texttt{refPositionsFile}.
+    \cvnamebasedonly{For a PDB file specifically, the components are read from the X, Y and Z fields.
+      \textbf{Note:} \emph{The PDB file has limited precision and fixed-point numbers: in some cases, the vector components may not be accurately represented; a XYZ file should be used instead, containing floating-point numbers.}}
+  }
 
+\cvnamebasedonly{
 \item %
   \key
     {vectorCol}{%