Merge pull request #187 from akohlmey/colvars-update-2016-09-30

update colvars library to version 2016-09-30

doc/src/fix_colvars.txt

@@ -39,8 +39,8 @@ metadynamics, Steered Molecular Dynamics (SMD) and Umbrella Sampling
 (US) via a flexible harmonic restraint bias. The colvars library is
 hosted at "http://colvars.github.io/"_http://colvars.github.io/
 
-This documentation describes only the fix colvars command itself and 
-LAMMPS specific parts of the code.  The full documentation of the 
+This documentation describes only the fix colvars command itself and
+LAMMPS specific parts of the code.  The full documentation of the
 colvars library is available as "this supplementary PDF document"_PDF/colvars-refman-lammps.pdf
 
 A detailed discussion of the implementation of the portable collective
@@ -122,7 +122,7 @@ not a limitation of functionality.
 
 [Default:]
 
-The default options are input = NULL, output = out, seed = 1966, unwrap yes, 
+The default options are input = NULL, output = out, seed = 1966, unwrap yes,
 and tstat = NULL.
 
 :line

examples/cmap/charmm22.cmap

@@ -1,4 +1,5 @@
-# Title: charmm correction map
+# DATE: 2016-09-26 CONTRIBUTOR: Robert Latour, latourr@clemson.edu  CITATION: TBA
+# Title: charmm22/charmm27 dihedral correction map
 
 #  alanine map, type 1
 

examples/tad/README

@@ -1,5 +1,5 @@
 Run this example as:
 
-mpirun -np 3 lmp_linux -partition 3x1 -in in.tad
+mpirun -np 3 lmp_g++ -partition 3x1 -in in.tad
 
 You should be able to use any number of replicas >= 3.

lib/colvars/colvar.cpp

@@ -150,7 +150,7 @@ colvar::colvar(std::string const &conf)
     feature_states[f_cv_linear]->enabled = lin;
   }
 
-  // Colvar is homogeneous iff:
+  // Colvar is homogeneous if:
   // - it is linear (hence not scripted)
   // - all cvcs have coefficient 1 or -1
   // i.e. sum or difference of cvcs
@@ -375,28 +375,16 @@ colvar::colvar(std::string const &conf)
 
   {
     bool temp;
-    if (get_keyval(conf, "outputSystemForce", temp, false)) {
-      cvm::error("Colvar option outputSystemForce is deprecated.\n"
-        "Please use outputTotalForce, or outputSystemForce within an ABF bias.");
+    if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) {
+      cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n"
+                 "The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR);
       return;
     }
   }
 
-  {
-    bool b_output_total_force;
-    get_keyval(conf, "outputTotalForce", b_output_total_force, false);
-    if (b_output_total_force) {
-      enable(f_cv_output_total_force);
-    }
-  }
-
-  {
-    bool b_output_applied_force;
-    get_keyval(conf, "outputAppliedForce", b_output_applied_force, false);
-    if (b_output_applied_force) {
-      enable(f_cv_output_applied_force);
-    }
-  }
+  get_keyval_feature(this, conf, "outputTotalForce", f_cv_output_total_force, false);
+  get_keyval_feature(this, conf, "outputAppliedForce", f_cv_output_applied_force, false);
+  get_keyval_feature(this, conf, "subtractAppliedForce", f_cv_subtract_applied_force, false);
 
   // Start in active state by default
   enable(f_cv_active);
@@ -409,6 +397,8 @@ colvar::colvar(std::string const &conf)
   fj.type(value());
   ft.type(value());
   ft_reported.type(value());
+  f_old.type(value());
+  f_old.reset();
 
   if (cvm::b_analysis)
     parse_analysis(conf);
@@ -519,6 +509,8 @@ int colvar::init_components(std::string const &conf)
     "number", "coordNum");
   error_code |= init_components_type<selfcoordnum>(conf, "self-coordination "
     "number", "selfCoordNum");
+  error_code |= init_components_type<groupcoordnum>(conf, "group-coordination "
+    "number", "groupCoord");
   error_code |= init_components_type<angle>(conf, "angle", "angle");
   error_code |= init_components_type<dipole_angle>(conf, "dipole angle", "dipoleAngle");
   error_code |= init_components_type<dihedral>(conf, "dihedral", "dihedral");
@@ -1104,6 +1096,14 @@ int colvar::calc_colvar_properties()
 
   } else {
 
+    if (is_enabled(f_cv_subtract_applied_force)) {
+      // correct the total force only if it has been measured
+      // TODO add a specific test instead of relying on sq norm
+      if (ft.norm2() > 0.0) {
+        ft -= f_old;
+      }
+    }
+
     x_reported = x;
     ft_reported = ft;
   }
@@ -1210,6 +1210,10 @@ cvm::real colvar::update_forces_energy()
     x_old = x;
   }
 
+  if (is_enabled(f_cv_subtract_applied_force)) {
+    f_old = f;
+  }
+
   if (cvm::debug())
     cvm::log("Done updating colvar \""+this->name+"\".\n");
   return (potential_energy + kinetic_energy);
@@ -1640,15 +1644,9 @@ std::ostream & colvar::write_traj(std::ostream &os)
   }
 
   if (is_enabled(f_cv_output_applied_force)) {
-    if (is_enabled(f_cv_extended_Lagrangian)) {
-      os << " "
-         << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-         << fr;
-    } else {
-      os << " "
-         << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
-         << f;
-    }
+    os << " "
+       << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
+       << applied_force();
   }
 
   return os;

lib/colvars/colvar.h

@@ -175,6 +175,9 @@ public:
   /// (if defined) contribute to it
   colvarvalue f;
 
+  /// Applied force at the previous step (to be subtracted from total force if needed)
+  colvarvalue f_old;
+
   /// \brief Total force, as derived from the atomic trajectory;
   /// should equal the system force plus \link f \endlink
   colvarvalue ft;
@@ -272,10 +275,13 @@ public:
   /// \brief Calculate the quantities associated to the colvar (but not to the CVCs)
   int calc_colvar_properties();
 
-  /// Get the current biasing force
-  inline colvarvalue bias_force() const
+  /// Get the current applied force
+  inline colvarvalue const applied_force() const
   {
-    return fb;
+    if (is_enabled(f_cv_extended_Lagrangian)) {
+      return fr;
+    }
+    return f;
   }
 
   /// Set the total biasing force to zero
@@ -482,6 +488,7 @@ public:
   class dihedral;
   class coordnum;
   class selfcoordnum;
+  class groupcoordnum;
   class h_bond;
   class rmsd;
   class orientation_angle;

lib/colvars/colvarbias_histogram.cpp

@@ -67,7 +67,7 @@ int colvarbias_histogram::init(std::string const &conf)
 
   if (colvar_array_size > 0) {
     weights.assign(colvar_array_size, 1.0);
-    get_keyval(conf, "weights", weights, weights, colvarparse::parse_silent);
+    get_keyval(conf, "weights", weights, weights);
   }
 
   for (i = 0; i < colvars.size(); i++) {
@@ -79,7 +79,7 @@ int colvarbias_histogram::init(std::string const &conf)
 
   {
     std::string grid_conf;
-    if (key_lookup(conf, "grid", grid_conf)) {
+    if (key_lookup(conf, "histogramGrid", grid_conf)) {
       grid->parse_params(grid_conf);
     }
   }

lib/colvars/colvarbias_meta.cpp

@@ -92,7 +92,11 @@ int colvarbias_meta::init(std::string const &conf)
     get_keyval(conf, "keepHills", keep_hills, false);
     if (! get_keyval(conf, "writeFreeEnergyFile", dump_fes, true))
       get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
-    get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false);
+    if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) {
+      cvm::log("Option \"saveFreeEnergyFile\" is deprecated, "
+               "please use \"keepFreeEnergyFiles\" instead.");
+    }
+    get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save);
 
     hills_energy           = new colvar_grid_scalar(colvars);
     hills_energy_gradients = new colvar_grid_gradient(colvars);

lib/colvars/colvarbias_restraint.cpp

@@ -612,3 +612,250 @@ cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k,
 
 
 
+colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key)
+  : colvarbias(key)
+{
+  lower_boundary = 0.0;
+  upper_boundary = 0.0;
+  width = 0.0;
+  gaussian_width = 0.0;
+}
+
+
+int colvarbias_restraint_histogram::init(std::string const &conf)
+{
+  colvarbias::init(conf);
+
+  get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary);
+  get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary);
+  get_keyval(conf, "width", width, width);
+
+  if (width <= 0.0) {
+    cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR);
+  }
+
+  get_keyval(conf, "gaussianWidth", gaussian_width, 2.0 * width, colvarparse::parse_silent);
+  get_keyval(conf, "gaussianSigma", gaussian_width, 2.0 * width);
+
+  if (lower_boundary >= upper_boundary) {
+    cvm::error("Error: the upper boundary, "+
+               cvm::to_str(upper_boundary)+
+               ", is not higher than the lower boundary, "+
+               cvm::to_str(lower_boundary)+".\n",
+               INPUT_ERROR);
+  }
+
+  cvm::real const nbins = (upper_boundary - lower_boundary) / width;
+  int const nbins_round = (int)(nbins);
+
+  if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
+    cvm::log("Warning: grid interval ("+
+             cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+
+             cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+
+             ") is not commensurate to its bin width ("+
+             cvm::to_str(width, cvm::cv_width, cvm::cv_prec)+").\n");
+  }
+
+  p.resize(nbins_round);
+  ref_p.resize(nbins_round);
+  p_diff.resize(nbins_round);
+
+  bool const inline_ref_p =
+    get_keyval(conf, "refHistogram", ref_p.data_array(), ref_p.data_array());
+  std::string ref_p_file;
+  get_keyval(conf, "refHistogramFile", ref_p_file, std::string(""));
+  if (ref_p_file.size()) {
+    if (inline_ref_p) {
+      cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n",
+                 INPUT_ERROR);
+    } else {
+      std::ifstream is(ref_p_file.c_str());
+      std::string data_s = "";
+      std::string line;
+      while (getline_nocomments(is, line)) {
+        data_s.append(line+"\n");
+      }
+      if (data_s.size() == 0) {
+        cvm::error("Error: file \""+ref_p_file+"\" empty or unreadable.\n", FILE_ERROR);
+      }
+      is.close();
+      cvm::vector1d<cvm::real> data;
+      if (data.from_simple_string(data_s) != 0) {
+        cvm::error("Error: could not read histogram from file \""+ref_p_file+"\".\n");
+      }
+      if (data.size() == 2*ref_p.size()) {
+        // file contains both x and p(x)
+        size_t i;
+        for (i = 0; i < ref_p.size(); i++) {
+          ref_p[i] = data[2*i+1];
+        }
+      } else if (data.size() == ref_p.size()) {
+        ref_p = data;
+      } else {
+        cvm::error("Error: file \""+ref_p_file+"\" contains a histogram of different length.\n",
+                   INPUT_ERROR);
+      }
+    }
+  }
+  cvm::real const ref_integral = ref_p.sum() * width;
+  if (std::fabs(ref_integral - 1.0) > 1.0e-03) {
+    cvm::log("Reference distribution not normalized, normalizing to unity.\n");
+    ref_p /= ref_integral;
+  }
+
+  get_keyval(conf, "writeHistogram", b_write_histogram, false);
+  get_keyval(conf, "forceConstant", force_k, 1.0);
+
+  return COLVARS_OK;
+}
+
+
+colvarbias_restraint_histogram::~colvarbias_restraint_histogram()
+{
+  p.resize(0);
+  ref_p.resize(0);
+  p_diff.resize(0);
+}
+
+
+int colvarbias_restraint_histogram::update()
+{
+  if (cvm::debug())
+    cvm::log("Updating the histogram restraint bias \""+this->name+"\".\n");
+
+  size_t vector_size = 0;
+  size_t icv;
+  for (icv = 0; icv < colvars.size(); icv++) {
+    vector_size += colvars[icv]->value().size();
+  }
+
+  cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size);
+
+  // calculate the histogram
+  p.reset();
+  for (icv = 0; icv < colvars.size(); icv++) {
+    colvarvalue const &cv = colvars[icv]->value();
+    if (cv.type() == colvarvalue::type_scalar) {
+      cvm::real const cv_value = cv.real_value;
+      size_t igrid;
+      for (igrid = 0; igrid < p.size(); igrid++) {
+        cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+        p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                                    (2.0 * gaussian_width * gaussian_width));
+      }
+    } else if (cv.type() == colvarvalue::type_vector) {
+      size_t idim;
+      for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
+        cvm::real const cv_value = cv.vector1d_value[idim];
+        size_t igrid;
+        for (igrid = 0; igrid < p.size(); igrid++) {
+          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+          p[igrid] += norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                                      (2.0 * gaussian_width * gaussian_width));
+        }
+      }
+    } else {
+      // TODO
+    }
+  }
+
+  cvm::real const force_k_cv = force_k * vector_size;
+
+  // calculate the difference between current and reference
+  p_diff = p - ref_p;
+  bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
+
+  // calculate the forces
+  for (icv = 0; icv < colvars.size(); icv++) {
+    colvarvalue const &cv = colvars[icv]->value();
+    colvarvalue &cv_force = colvar_forces[icv];
+    cv_force.type(cv);
+    cv_force.reset();
+
+    if (cv.type() == colvarvalue::type_scalar) {
+      cvm::real const cv_value = cv.real_value;
+      cvm::real &force = cv_force.real_value;
+      size_t igrid;
+      for (igrid = 0; igrid < p.size(); igrid++) {
+        cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+        force += force_k_cv * p_diff[igrid] *
+          norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                          (2.0 * gaussian_width * gaussian_width)) *
+          (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
+      }
+    } else if (cv.type() == colvarvalue::type_vector) {
+      size_t idim;
+      for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
+        cvm::real const cv_value = cv.vector1d_value[idim];
+        cvm::real &force = cv_force.vector1d_value[idim];
+        size_t igrid;
+        for (igrid = 0; igrid < p.size(); igrid++) {
+          cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
+          force += force_k_cv * p_diff[igrid] *
+            norm * std::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
+                            (2.0 * gaussian_width * gaussian_width)) *
+            (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
+        }
+      }
+    } else {
+      // TODO
+    }
+  }
+
+  return COLVARS_OK;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os)
+{
+  if (b_write_histogram) {
+    std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat");
+    std::ofstream os(file_name.c_str());
+    os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width)
+       << "  " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3)
+       << ")\n";
+    size_t igrid;
+    for (igrid = 0; igrid < p.size(); igrid++) {
+      cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
+      os << "  "
+         << std::setprecision(cvm::cv_prec)
+         << std::setw(cvm::cv_width)
+         << x_grid
+         << "  "
+         << std::setprecision(cvm::cv_prec)
+         << std::setw(cvm::cv_width)
+         << p[igrid] << "\n";
+    }
+    os.close();
+  }
+  return os;
+}
+
+
+std::istream & colvarbias_restraint_histogram::read_restart(std::istream &is)
+{
+  return is;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_traj_label(std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy) {
+    os << " E_"
+       << cvm::wrap_string(this->name, cvm::en_width-2);
+  }
+  return os;
+}
+
+
+std::ostream & colvarbias_restraint_histogram::write_traj(std::ostream &os)
+{
+  os << " ";
+  if (b_output_energy) {
+    os << " "
+       << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
+       << bias_energy;
+  }
+  return os;
+}

lib/colvars/colvarbias_restraint.h

@@ -168,4 +168,52 @@ protected:
 };
 
 
+/// Restrain the 1D histogram of a set of variables (or of a multidimensional one)
+// TODO this could be reimplemented more cleanly as a derived class of both restraint and histogram
+class colvarbias_restraint_histogram : public colvarbias {
+
+public:
+
+  colvarbias_restraint_histogram(char const *key);
+  int init(std::string const &conf);
+  ~colvarbias_restraint_histogram();
+
+  virtual int update();
+
+  virtual std::istream & read_restart(std::istream &is);
+  virtual std::ostream & write_restart(std::ostream &os);
+  virtual std::ostream & write_traj_label(std::ostream &os);
+  virtual std::ostream & write_traj(std::ostream &os);
+
+protected:
+
+  /// Probability density
+  cvm::vector1d<cvm::real> p;
+
+  /// Reference probability density
+  cvm::vector1d<cvm::real> ref_p;
+
+  /// Difference between probability density and reference
+  cvm::vector1d<cvm::real> p_diff;
+
+  /// Lower boundary of the grid
+  cvm::real lower_boundary;
+
+  /// Upper boundary of the grid
+  cvm::real upper_boundary;
+
+  /// Resolution of the grid
+  cvm::real width;
+
+  /// Width of the Gaussians
+  cvm::real gaussian_width;
+
+  /// Restraint force constant
+  cvm::real force_k;
+
+  /// Write the histogram to a file
+  bool b_write_histogram;
+};
+
+
 #endif

lib/colvars/colvarcomp.cpp

@@ -54,6 +54,21 @@ colvar::cvc::cvc(std::string const &conf)
 }
 
 
+int colvar::cvc::init_total_force_params(std::string const &conf)
+{
+  if (get_keyval_feature(this, conf, "oneSiteSystemForce",
+                         f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
+    cvm::log("Warning: keyword \"oneSiteSystemForce\" is deprecated: "
+             "please use \"oneSiteTotalForce\" instead.\n");
+  }
+  if (get_keyval_feature(this, conf, "oneSiteTotalForce",
+                         f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) {
+    cvm::log("Computing total force on group 1 only");
+  }
+  return COLVARS_OK;
+}
+
+
 cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
                                           char const *group_key,
                                           bool optional)
@@ -77,8 +92,6 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
 
       if (is_enabled(f_cvc_scalable)) {
         cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n");
-      } else {
-        cvm::log("Scalable calculation is not available for group \""+group->key+"\" with the current configuration.\n");
       }
     }
 

lib/colvars/colvarcomp.h

@@ -108,6 +108,9 @@ public:
                    char const *group_key,
                    bool optional = false);
 
+  /// \brief Parse options pertaining to total force calculation
+  virtual int init_total_force_params(std::string const &conf);
+
   /// \brief After construction, set data related to dependency handling
   int setup();
 
@@ -306,9 +309,6 @@ protected:
   cvm::rvector     dist_v;
   /// Use absolute positions, ignoring PBCs when present
   bool b_no_PBC;
-  /// Compute total force on first site only to avoid unwanted
-  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
-  bool b_1site_force;
 public:
   distance(std::string const &conf);
   distance();
@@ -388,9 +388,6 @@ protected:
   cvm::atom_group  *ref2;
   /// Use absolute positions, ignoring PBCs when present
   bool b_no_PBC;
-  /// Compute total force on one site only to avoid unwanted
-  /// coupling to other colvars (see e.g. Ciccotti et al., 2005)
-  bool b_1site_force;
   /// Vector on which the distance vector is projected
   cvm::rvector axis;
   /// Norm of the axis
@@ -854,6 +851,62 @@ public:
                                   colvarvalue const &x2) const;
 };
 
+
+/// \brief Colvar component: coordination number between two groups
+/// (colvarvalue::type_scalar type, range [0:N1*N2])
+class colvar::groupcoordnum
+  : public colvar::distance
+{
+protected:
+  /// \brief "Cutoff" for isotropic calculation (default)
+  cvm::real     r0;
+  /// \brief "Cutoff vector" for anisotropic calculation
+  cvm::rvector  r0_vec;
+  /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
+  /// used
+  bool b_anisotropic;
+  /// Integer exponent of the function numerator
+  int en;
+  /// Integer exponent of the function denominator
+  int ed;
+public:
+  /// Constructor
+  groupcoordnum(std::string const &conf);
+  groupcoordnum();
+  virtual inline ~groupcoordnum() {}
+  virtual void calc_value();
+  virtual void calc_gradients();
+  virtual void apply_force(colvarvalue const &force);
+  template<bool b_gradients>
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
+  /// coordination number \param exp_num \i n exponent \param exp_den
+  /// \i m exponent \param A1 atom \param A2 atom
+  static cvm::real switching_function(cvm::real const &r0,
+                                      int const &exp_num, int const &exp_den,
+                                      cvm::atom &A1, cvm::atom &A2);
+
+  /*
+  template<bool b_gradients>
+  /// \brief Calculate a coordination number through the function
+  /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
+  /// vector of different cutoffs in the three directions \param
+  /// exp_num \i n exponent \param exp_den \i m exponent \param A1
+  /// atom \param A2 atom
+  static cvm::real switching_function(cvm::rvector const &r0_vec,
+                                      int const &exp_num, int const &exp_den,
+                                      cvm::atom &A1, cvm::atom &A2);
+
+  virtual cvm::real dist2(colvarvalue const &x1,
+                          colvarvalue const &x2) const;
+  virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
+                                  colvarvalue const &x2) const;
+  */
+};
+
+
 /// \brief Colvar component: hydrogen bond, defined as the product of
 /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol))
 /// (colvarvalue::type_scalar type, range [0:1])

lib/colvars/colvarcomp_angles.cpp

@@ -18,9 +18,9 @@ colvar::angle::angle(std::string const &conf)
   group1 = parse_group(conf, "group1");
   group2 = parse_group(conf, "group2");
   group3 = parse_group(conf, "group3");
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
+
+  init_total_force_params(conf);
+
   x.type(colvarvalue::type_scalar);
 }
 
@@ -33,7 +33,6 @@ colvar::angle::angle(cvm::atom const &a1,
   provide(f_cvc_inv_gradient);
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
-  b_1site_force = false;
 
   group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
   group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@@ -94,7 +93,7 @@ void colvar::angle::calc_force_invgrads()
   // centered on group2, which means group2 is kept fixed
   // when propagating changes in the angle)
 
-  if (b_1site_force) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
     group1->read_total_forces();
     cvm::real norm_fact = 1.0 / dxdr1.norm2();
     ft.real_value = norm_fact * dxdr1 * group1->total_force();
@@ -140,9 +139,8 @@ colvar::dipole_angle::dipole_angle(std::string const &conf)
   group2 = parse_group(conf, "group2");
   group3 = parse_group(conf, "group3");
 
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
+  init_total_force_params(conf);
+
   x.type(colvarvalue::type_scalar);
 }
 
@@ -152,7 +150,6 @@ colvar::dipole_angle::dipole_angle(cvm::atom const &a1,
                       cvm::atom const &a3)
 {
   function_type = "dipole_angle";
-  b_1site_force = false;
 
   group1 = new cvm::atom_group(std::vector<cvm::atom>(1, a1));
   group2 = new cvm::atom_group(std::vector<cvm::atom>(1, a2));
@@ -250,14 +247,13 @@ colvar::dihedral::dihedral(std::string const &conf)
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
 
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
   group1 = parse_group(conf, "group1");
   group2 = parse_group(conf, "group2");
   group3 = parse_group(conf, "group3");
   group4 = parse_group(conf, "group4");
 
+  init_total_force_params(conf);
+
   x.type(colvarvalue::type_scalar);
 }
 
@@ -422,7 +418,7 @@ void colvar::dihedral::calc_force_invgrads()
   cvm::real const fact4 = d34 * std::sqrt(1.0 - dot4 * dot4);
 
   group1->read_total_forces();
-  if ( b_1site_force ) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
     // This is only measuring the force on group 1
     ft.real_value = PI/180.0 * fact1 * (cross1 * group1->total_force());
   } else {

lib/colvars/colvarcomp_coordnums.cpp

@@ -338,3 +338,151 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
   }
 }
 
+
+// groupcoordnum member functions
+colvar::groupcoordnum::groupcoordnum(std::string const &conf)
+  : distance(conf), b_anisotropic(false)
+{
+  function_type = "groupcoordnum";
+  x.type(colvarvalue::type_scalar);
+
+  // group1 and group2 are already initialized by distance()
+  if (group1->b_dummy || group2->b_dummy)
+    cvm::fatal_error("Error: neither group can be a dummy atom\n");
+
+  bool const b_scale = get_keyval(conf, "cutoff", r0,
+                                  cvm::real(4.0 * cvm::unit_angstrom()));
+
+  if (get_keyval(conf, "cutoff3", r0_vec,
+                 cvm::rvector(4.0, 4.0, 4.0), parse_silent)) {
+
+    if (b_scale)
+      cvm::fatal_error("Error: cannot specify \"scale\" and "
+                       "\"scale3\" at the same time.\n");
+    b_anisotropic = true;
+    // remove meaningless negative signs
+    if (r0_vec.x < 0.0) r0_vec.x *= -1.0;
+    if (r0_vec.y < 0.0) r0_vec.y *= -1.0;
+    if (r0_vec.z < 0.0) r0_vec.z *= -1.0;
+  }
+
+  get_keyval(conf, "expNumer", en, int(6) );
+  get_keyval(conf, "expDenom", ed, int(12));
+
+  if ( (en%2) || (ed%2) ) {
+    cvm::fatal_error("Error: odd exponents provided, can only use even ones.\n");
+  }
+
+}
+
+
+colvar::groupcoordnum::groupcoordnum()
+  : b_anisotropic(false)
+{
+  function_type = "groupcoordnum";
+  x.type(colvarvalue::type_scalar);
+}
+
+
+template<bool calculate_gradients>
+cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
+                                                    int const &en,
+                                                    int const &ed,
+                                                    cvm::atom &A1,
+                                                    cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
+  cvm::real const l2 = diff.norm2()/(r0*r0);
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow(l2, en2);
+  cvm::real const xd = std::pow(l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx = (2.0/(r0*r0))*diff;
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+
+  return func;
+}
+
+
+#if 0  // AMG: I don't think there's any reason to support anisotropic,
+       //      and I don't have those flags below in calc_value, but
+       //      if I need them, I'll also need to uncomment this method
+template<bool calculate_gradients>
+cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
+                                                    int const &en,
+                                                    int const &ed,
+                                                    cvm::atom &A1,
+                                                    cvm::atom &A2)
+{
+  cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
+  cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
+  cvm::real const l2 = scal_diff.norm2();
+
+  // Assume en and ed are even integers, and avoid sqrt in the following
+  int const en2 = en/2;
+  int const ed2 = ed/2;
+
+  cvm::real const xn = std::pow(l2, en2);
+  cvm::real const xd = std::pow(l2, ed2);
+  cvm::real const func = (1.0-xn)/(1.0-xd);
+
+  if (calculate_gradients) {
+    cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
+    cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
+                             (2.0/(r0_vec.y*r0_vec.y))*diff.y,
+                             (2.0/(r0_vec.z*r0_vec.z))*diff.z);
+    A1.grad += (-1.0)*dFdl2*dl2dx;
+    A2.grad +=        dFdl2*dl2dx;
+  }
+  return func;
+}
+#endif
+
+
+
+void colvar::groupcoordnum::calc_value()
+{
+
+  // create fake atoms to hold the com coordinates
+  cvm::atom group1_com_atom;
+  cvm::atom group2_com_atom;
+  group1_com_atom.pos = group1->center_of_mass();
+  group2_com_atom.pos = group2->center_of_mass();
+
+  x.real_value = coordnum::switching_function<false>(r0, en, ed,
+                                                     group1_com_atom, group2_com_atom);
+
+}
+
+
+void colvar::groupcoordnum::calc_gradients()
+{
+  cvm::atom group1_com_atom;
+  cvm::atom group2_com_atom;
+  group1_com_atom.pos = group1->center_of_mass();
+  group2_com_atom.pos = group2->center_of_mass();
+
+  coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
+  group1->set_weighted_gradient(group1_com_atom.grad);
+  group2->set_weighted_gradient(group2_com_atom.grad);
+
+}
+
+
+void colvar::groupcoordnum::apply_force(colvarvalue const &force)
+{
+  if (!group1->noforce)
+    group1->apply_colvar_force(force.real_value);
+
+  if (!group2->noforce)
+    group2->apply_colvar_force(force.real_value);
+}

lib/colvars/colvarcomp_distances.cpp

@@ -18,14 +18,14 @@ colvar::distance::distance(std::string const &conf)
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
 
+  group1 = parse_group(conf, "group1");
+  group2 = parse_group(conf, "group2");
+
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group 1 only");
-  }
-  group1 = parse_group(conf, "group1");
-  group2 = parse_group(conf, "group2");
+
+  init_total_force_params(conf);
 
   x.type(colvarvalue::type_scalar);
 }
@@ -38,7 +38,6 @@ colvar::distance::distance()
   provide(f_cvc_inv_gradient);
   provide(f_cvc_Jacobian);
   provide(f_cvc_com_based);
-  b_1site_force = false;
   b_no_PBC = false;
   x.type(colvarvalue::type_scalar);
 }
@@ -67,7 +66,7 @@ void colvar::distance::calc_gradients()
 void colvar::distance::calc_force_invgrads()
 {
   group1->read_total_forces();
-  if ( b_1site_force ) {
+  if (is_enabled(f_cvc_one_site_total_force)) {
     ft.real_value = -1.0 * (group1->total_force() * dist_v.unit());
   } else {
     group2->read_total_forces();
@@ -97,6 +96,7 @@ colvar::distance_vec::distance_vec(std::string const &conf)
   : distance(conf)
 {
   function_type = "distance_vec";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_3vector);
 }
 
@@ -105,6 +105,7 @@ colvar::distance_vec::distance_vec()
   : distance()
 {
   function_type = "distance_vec";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_3vector);
 }
 
@@ -185,9 +186,9 @@ colvar::distance_z::distance_z(std::string const &conf)
   if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) {
     cvm::log("Computing distance using absolute positions (not minimal-image)");
   }
-  if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) {
-    cvm::log("Computing total force on group \"main\" only");
-  }
+
+  init_total_force_params(conf);
+
 }
 
 colvar::distance_z::distance_z()
@@ -251,7 +252,7 @@ void colvar::distance_z::calc_force_invgrads()
 {
   main->read_total_forces();
 
-  if (fixed_axis && !b_1site_force) {
+  if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
     ref1->read_total_forces();
     ft.real_value = 0.5 * ((main->total_force() - ref1->total_force()) * axis);
   } else {
@@ -351,7 +352,7 @@ void colvar::distance_xy::calc_force_invgrads()
 {
   main->read_total_forces();
 
-  if (fixed_axis && !b_1site_force) {
+  if (fixed_axis && !is_enabled(f_cvc_one_site_total_force)) {
     ref1->read_total_forces();
     ft.real_value = 0.5 / x.real_value * ((main->total_force() - ref1->total_force()) * dist_v_ortho);
   } else {
@@ -382,6 +383,7 @@ colvar::distance_dir::distance_dir(std::string const &conf)
   : distance(conf)
 {
   function_type = "distance_dir";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_unit3vector);
 }
 
@@ -390,6 +392,7 @@ colvar::distance_dir::distance_dir()
   : distance()
 {
   function_type = "distance_dir";
+  provide(f_cvc_com_based);
   x.type(colvarvalue::type_unit3vector);
 }
 
@@ -461,7 +464,6 @@ colvar::distance_inv::distance_inv()
 {
   function_type = "distance_inv";
   exponent = 6;
-  b_1site_force = false;
   x.type(colvarvalue::type_scalar);
 }
 

lib/colvars/colvardeps.cpp

@@ -293,6 +293,9 @@ void colvardeps::init_cv_requires() {
     f_description(f_cv_output_total_force, "output total force");
     f_req_self(f_cv_output_total_force, f_cv_total_force);
 
+    f_description(f_cv_subtract_applied_force, "subtract applied force from total force");
+    f_req_self(f_cv_subtract_applied_force, f_cv_total_force);
+
     f_description(f_cv_lower_boundary, "lower boundary");
     f_req_self(f_cv_lower_boundary, f_cv_scalar);
 
@@ -376,6 +379,11 @@ void colvardeps::init_cvc_requires() {
 
     f_description(f_cvc_com_based, "depends on group centers of mass");
 
+    // Compute total force on first site only to avoid unwanted
+    // coupling to other colvars (see e.g. Ciccotti et al., 2005)
+    f_description(f_cvc_one_site_total_force, "compute total collective force only from one group center");
+    f_req_self(f_cvc_one_site_total_force, f_cvc_com_based);
+
     f_description(f_cvc_scalable, "scalable calculation");
     f_req_self(f_cvc_scalable, f_cvc_scalable_com);
 

lib/colvars/colvardeps.h

@@ -176,6 +176,8 @@ public:
     f_cv_total_force,
     /// \brief Calculate total force from atomic forces
     f_cv_total_force_calc,
+    /// \brief Subtract the applied force from the total force
+    f_cv_subtract_applied_force,
     /// \brief Estimate Jacobian derivative
     f_cv_Jacobian,
     /// \brief Do not report the Jacobian force as part of the total force
@@ -236,6 +238,7 @@ public:
     /// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed
     f_cvc_debug_gradient,
     f_cvc_Jacobian,
+    f_cvc_one_site_total_force,
     f_cvc_com_based,
     f_cvc_scalable,
     f_cvc_scalable_com,

lib/colvars/colvargrid.h

@@ -382,8 +382,8 @@ public:
   inline int current_bin_scalar(int const i, int const iv) const
   {
     return value_to_bin_scalar(actual_value[i] ?
-        cv[i]->actual_value().vector1d_value[iv] :
-        cv[i]->value().vector1d_value[iv], i);
+                               cv[i]->actual_value().vector1d_value[iv] :
+                               cv[i]->value().vector1d_value[iv], i);
   }
 
   /// \brief Use the lower boundary and the width to report which bin
@@ -395,8 +395,8 @@ public:
 
   /// \brief Same as the standard version, but uses another grid definition
   inline int value_to_bin_scalar(colvarvalue const &value,
-                                  colvarvalue const &new_offset,
-                                  cvm::real   const &new_width) const
+                                 colvarvalue const &new_offset,
+                                 cvm::real   const &new_width) const
   {
     return (int) std::floor( (value.real_value - new_offset.real_value) / new_width );
   }
@@ -410,22 +410,22 @@ public:
 
   /// \brief Same as the standard version, but uses different parameters
   inline colvarvalue bin_to_value_scalar(int const &i_bin,
-                                          colvarvalue const &new_offset,
-                                          cvm::real const &new_width) const
+                                         colvarvalue const &new_offset,
+                                         cvm::real const &new_width) const
   {
     return new_offset.real_value + new_width * (0.5 + i_bin);
   }
 
   /// Set the value at the point with index ix
   inline void set_value(std::vector<int> const &ix,
-                         T const &t,
-                         size_t const &imult = 0)
+                        T const &t,
+                        size_t const &imult = 0)
   {
     data[this->address(ix)+imult] = t;
     has_data = true;
   }
 
- /// \brief Get the change from this to other_grid
+  /// \brief Get the change from this to other_grid
   /// and store the result in this.
   /// this_grid := other_grid - this_grid
   /// Grids must have the same dimensions.
@@ -434,13 +434,13 @@ public:
 
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to subtract two grids with "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
 
     if (other_grid.data.size() != this->data.size()) {
       cvm::error("Error: trying to subtract two grids with "
-                        "different size.\n");
+                 "different size.\n");
       return;
     }
 
@@ -457,13 +457,13 @@ public:
   {
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to copy two grids with "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
 
     if (other_grid.data.size() != this->data.size()) {
       cvm::error("Error: trying to copy two grids with "
-                        "different size.\n");
+                 "different size.\n");
       return;
     }
 
@@ -493,7 +493,7 @@ public:
   /// \brief Get the binned value indexed by ix, or the first of them
   /// if the multiplicity is larger than 1
   inline T const & value(std::vector<int> const &ix,
-                          size_t const &imult = 0) const
+                         size_t const &imult = 0) const
   {
     return data[this->address(ix) + imult];
   }
@@ -541,7 +541,7 @@ public:
   /// boundaries; a negative number is returned if the given point is
   /// off-grid
   inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
-                                                 bool skip_hard_boundaries = false)
+                                                bool skip_hard_boundaries = false)
   {
     cvm::real minimum = 1.0E+16;
     for (size_t i = 0; i < nd; i++) {
@@ -574,7 +574,7 @@ public:
   {
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to merge two grids with values of "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
 
@@ -593,8 +593,8 @@ public:
       for (size_t i = 0; i < nd; i++) {
         oix[i] =
           value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
-                               ogb[i],
-                               ogw[i]);
+                              ogb[i],
+                              ogw[i]);
       }
 
       if (! other_grid.index_ok(oix)) {
@@ -614,11 +614,11 @@ public:
   /// \brief Add data from another grid of the same type, AND
   /// identical definition (boundaries, widths)
   void add_grid(colvar_grid<T> const &other_grid,
-                 cvm::real scale_factor = 1.0)
+                cvm::real scale_factor = 1.0)
   {
     if (other_grid.multiplicity() != this->multiplicity()) {
       cvm::error("Error: trying to sum togetehr two grids with values of "
-                        "different multiplicity.\n");
+                 "different multiplicity.\n");
       return;
     }
     if (scale_factor != 1.0)
@@ -636,7 +636,7 @@ public:
   /// \brief Return the value suitable for output purposes (so that it
   /// may be rescaled or manipulated without changing it permanently)
   virtual inline T value_output(std::vector<int> const &ix,
-                                 size_t const &imult = 0)
+                                size_t const &imult = 0)
   {
     return value(ix, imult);
   }
@@ -645,9 +645,9 @@ public:
   /// into the internal representation (the two may be different,
   /// e.g. when using colvar_grid_count)
   virtual inline void value_input(std::vector<int> const &ix,
-                                   T const &t,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  T const &t,
+                                  size_t const &imult = 0,
+                                  bool add = false)
   {
     if ( add )
       data[address(ix) + imult] += t;
@@ -737,7 +737,8 @@ public:
   }
 
   /// Read a grid definition from a config string
-  int parse_params(std::string const &conf)
+  int parse_params(std::string const &conf,
+                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal)
   {
     if (cvm::debug()) cvm::log("Reading grid configuration from string.\n");
 
@@ -746,30 +747,33 @@ public:
 
     {
       size_t nd_in = 0;
+      // this is only used in state files
       colvarparse::get_keyval(conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
       if (nd_in != nd) {
         cvm::error("Error: trying to read data for a grid "
-                    "that contains a different number of colvars ("+
-                    cvm::to_str(nd_in)+") than the grid defined "
-                    "in the configuration file("+cvm::to_str(nd)+
-                    ").\n");
+                   "that contains a different number of colvars ("+
+                   cvm::to_str(nd_in)+") than the grid defined "
+                   "in the configuration file("+cvm::to_str(nd)+
+                   ").\n");
         return COLVARS_ERROR;
       }
     }
 
+    // underscore keywords are used in state file
     colvarparse::get_keyval(conf, "lower_boundaries",
-                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
+                            lower_boundaries, lower_boundaries, colvarparse::parse_silent);
     colvarparse::get_keyval(conf, "upper_boundaries",
-                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+                            upper_boundaries, upper_boundaries, colvarparse::parse_silent);
 
-    // support also camel case
+    // camel case keywords are used in config file
     colvarparse::get_keyval(conf, "lowerBoundaries",
-                             lower_boundaries, lower_boundaries, colvarparse::parse_silent);
+                            lower_boundaries, lower_boundaries, parse_mode);
     colvarparse::get_keyval(conf, "upperBoundaries",
-                             upper_boundaries, upper_boundaries, colvarparse::parse_silent);
+                            upper_boundaries, upper_boundaries, parse_mode);
 
-    colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent);
+    colvarparse::get_keyval(conf, "widths", widths, widths, parse_mode);
 
+    // only used in state file
     colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent);
 
     if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
@@ -808,13 +812,13 @@ public:
   {
     for (size_t i = 0; i < nd; i++) {
       if ( (std::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
-                                     lower_boundaries[i])) > 1.0E-10) ||
+                                   lower_boundaries[i])) > 1.0E-10) ||
            (std::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
-                                     upper_boundaries[i])) > 1.0E-10) ||
+                                   upper_boundaries[i])) > 1.0E-10) ||
            (std::sqrt(cv[i]->dist2(cv[i]->width,
-                                     widths[i])) > 1.0E-10) ) {
+                                   widths[i])) > 1.0E-10) ) {
         cvm::error("Error: restart information for a grid is "
-                    "inconsistent with that of its colvars.\n");
+                   "inconsistent with that of its colvars.\n");
         return;
       }
     }
@@ -830,19 +834,19 @@ public:
       // matter: boundaries should be EXACTLY the same (otherwise,
       // map_grid() should be used)
       if ( (std::fabs(other_grid.lower_boundaries[i] -
-                       lower_boundaries[i]) > 1.0E-10) ||
+                      lower_boundaries[i]) > 1.0E-10) ||
            (std::fabs(other_grid.upper_boundaries[i] -
-                       upper_boundaries[i]) > 1.0E-10) ||
+                      upper_boundaries[i]) > 1.0E-10) ||
            (std::fabs(other_grid.widths[i] -
-                       widths[i]) > 1.0E-10) ||
+                      widths[i]) > 1.0E-10) ||
            (data.size() != other_grid.data.size()) ) {
-      cvm::error("Error: inconsistency between "
-                  "two grids that are supposed to be equal, "
-                  "aside from the data stored.\n");
-      return;
+        cvm::error("Error: inconsistency between "
+                   "two grids that are supposed to be equal, "
+                   "aside from the data stored.\n");
+        return;
+      }
     }
   }
-}
 
 
   /// \brief Read grid entry in restart file
@@ -853,7 +857,7 @@ public:
     if ((is >> key) && (key == std::string("grid_parameters"))) {
       is.seekg(start_pos, std::ios::beg);
       is >> colvarparse::read_block("grid_parameters", conf);
-      parse_params(conf);
+      parse_params(conf, colvarparse::parse_silent);
     } else {
       cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n");
       is.seekg(start_pos, std::ios::beg);
@@ -871,11 +875,11 @@ public:
   }
 
 
-/// \brief Write the grid data without labels, as they are
-/// represented in memory
-/// \param buf_size Number of values per line
+  /// \brief Write the grid data without labels, as they are
+  /// represented in memory
+  /// \param buf_size Number of values per line
   std::ostream & write_raw(std::ostream &os,
-                            size_t const buf_size = 3)
+                           size_t const buf_size = 3)
   {
     std::streamsize const w = os.width();
     std::streamsize const p = os.precision();
@@ -935,10 +939,10 @@ public:
     os << std::setw(2) << "# " << nd << "\n";
     for (size_t i = 0; i < nd; i++) {
       os << "# "
-        << std::setw(10) << lower_boundaries[i]
-        << std::setw(10) << widths[i]
-        << std::setw(10) << nx[i] << "  "
-        << periodic[i] << "\n";
+         << std::setw(10) << lower_boundaries[i]
+         << std::setw(10) << widths[i]
+         << std::setw(10) << nx[i] << "  "
+         << periodic[i] << "\n";
     }
 
 
@@ -951,14 +955,14 @@ public:
 
       for (size_t i = 0; i < nd; i++) {
         os << " "
-          << std::setw(w) << std::setprecision(p)
-          << bin_to_value_scalar(ix[i], i);
+           << std::setw(w) << std::setprecision(p)
+           << bin_to_value_scalar(ix[i], i);
       }
       os << " ";
       for (size_t imult = 0; imult < mult; imult++) {
         os << " "
-          << std::setw(w) << std::setprecision(p)
-          << value_output(ix, imult);
+           << std::setw(w) << std::setprecision(p)
+           << value_output(ix, imult);
       }
       os << "\n";
     }
@@ -986,7 +990,7 @@ public:
 
     if ( !(is >> hash) || (hash != "#") ) {
       cvm::error("Error reading grid at position "+
-                  cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                 cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
       return is;
     }
 
@@ -1008,7 +1012,7 @@ public:
     for (size_t i = 0; i < nd; i++ ) {
       if ( !(is >> hash) || (hash != "#") ) {
         cvm::error("Error reading grid at position "+
-                    cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
+                   cvm::to_str(is.tellg())+" in stream(read \"" + hash + "\")\n");
         return is;
       }
 
@@ -1016,10 +1020,10 @@ public:
 
 
       if ( (std::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) ||
-          (std::fabs(width - widths[i] ) > 1.0e-10) ||
-          (nx_read[i] != nx[i]) ) {
+           (std::fabs(width - widths[i] ) > 1.0e-10) ||
+           (nx_read[i] != nx[i]) ) {
         cvm::log("Warning: reading from different grid definition (colvar "
-                  + cvm::to_str(i+1) + "); remapping data on new grid.\n");
+                 + cvm::to_str(i+1) + "); remapping data on new grid.\n");
         remap = true;
       }
     }
@@ -1063,7 +1067,6 @@ public:
 
   /// \brief Write the grid data without labels, as they are
   /// represented in memory
-  /// \param buf_size Number of values per line
   std::ostream & write_opendx(std::ostream &os)
   {
     // write the header
@@ -1122,11 +1125,11 @@ public:
 
   /// Constructor
   colvar_grid_count(std::vector<int> const &nx_i,
-                     size_t const           &def_count = 0);
+                    size_t const           &def_count = 0);
 
   /// Constructor from a vector of colvars
   colvar_grid_count(std::vector<colvar *>  &colvars,
-                     size_t const           &def_count = 0);
+                    size_t const           &def_count = 0);
 
   /// Increment the counter at given position
   inline void incr_count(std::vector<int> const &ix)
@@ -1136,7 +1139,7 @@ public:
 
   /// \brief Get the binned count indexed by ix from the newly read data
   inline size_t const & new_count(std::vector<int> const &ix,
-                                   size_t const &imult = 0)
+                                  size_t const &imult = 0)
   {
     return new_data[address(ix) + imult];
   }
@@ -1145,9 +1148,9 @@ public:
   /// into the internal representation (it may have been rescaled or
   /// manipulated)
   virtual inline void value_input(std::vector<int> const &ix,
-                                   size_t const &t,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  size_t const &t,
+                                  size_t const &imult = 0,
+                                  bool add = false)
   {
     if (add) {
       data[address(ix)] += t;
@@ -1164,7 +1167,7 @@ public:
   /// \brief Return the log-gradient from finite differences
   /// on the *same* grid for dimension n
   inline const cvm::real log_gradient_finite_diff( const std::vector<int> &ix0,
-    int n = 0)
+                                                   int n = 0)
   {
     cvm::real A0, A1;
     std::vector<int> ix;
@@ -1377,7 +1380,7 @@ public:
   /// \brief Return the value of the function at ix divided by its
   /// number of samples (if the count grid is defined)
   virtual inline cvm::real value_output(std::vector<int> const &ix,
-                                         size_t const &imult = 0)
+                                        size_t const &imult = 0)
   {
     if (samples)
       return (samples->value(ix) > 0) ?
@@ -1391,9 +1394,9 @@ public:
   /// into the internal representation (it may have been rescaled or
   /// manipulated)
   virtual inline void value_input(std::vector<int> const &ix,
-                                   cvm::real const &new_value,
-                                   size_t const &imult = 0,
-                                   bool add = false)
+                                  cvm::real const &new_value,
+                                  size_t const &imult = 0,
+                                  bool add = false)
   {
     if (add) {
       if (samples)

lib/colvars/colvarmodule.cpp

@@ -293,6 +293,9 @@ int colvarmodule::parse_biases(std::string const &conf)
   /// initialize histograms
   parse_biases_type<colvarbias_histogram>(conf, "histogram", n_histo_biases);
 
+  /// initialize histogram restraints
+  parse_biases_type<colvarbias_restraint_histogram>(conf, "histogramRestraint", n_rest_biases);
+
   /// initialize linear restraints
   parse_biases_type<colvarbias_restraint_linear>(conf, "linear", n_rest_biases);
 

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-09-14"
+#define COLVARS_VERSION "2016-09-30"
 #endif
 
 #ifndef COLVARS_DEBUG

lib/colvars/colvarscript.cpp

@@ -243,11 +243,17 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) {
   }
 
   if (subcmd == "getappliedforce") {
-    result = (cv->bias_force()).to_simple_string();
+    result = (cv->applied_force()).to_simple_string();
     return COLVARS_OK;
   }
 
   if (subcmd == "getsystemforce") {
+    // TODO warning here
+    result = (cv->total_force()).to_simple_string();
+    return COLVARS_OK;
+  }
+
+  if (subcmd == "gettotalforce") {
     result = (cv->total_force()).to_simple_string();
     return COLVARS_OK;
   }

lib/colvars/colvartypes.h

@@ -57,6 +57,12 @@ public:
     }
   }
 
+  /// Return a reference to the data
+  inline std::vector<T> &data_array()
+  {
+    return data;
+  }
+
   inline ~vector1d()
   {
     data.clear();
@@ -203,6 +209,16 @@ public:
     return std::sqrt(this->norm2());
   }
 
+  inline cvm::real sum() const
+  {
+    cvm::real result = 0.0;
+    size_t i;
+    for (i = 0; i < this->size(); i++) {
+      result += (*this)[i];
+    }
+    return result;
+  }
+
   /// Slicing
   inline vector1d<T> const slice(size_t const i1, size_t const i2) const
   {
@@ -295,11 +311,23 @@ public:
   {
     std::stringstream stream(s);
     size_t i = 0;
-    while ((stream >> (*this)[i]) && (i < this->size())) {
-      i++;
-    }
-    if (i < this->size()) {
-      return COLVARS_ERROR;
+    if (this->size()) {
+      while ((stream >> (*this)[i]) && (i < this->size())) {
+        i++;
+      }
+      if (i < this->size()) {
+        return COLVARS_ERROR;
+      }
+    } else {
+      T input;
+      while (stream >> input) {
+        if ((i % 100) == 0) {
+          data.reserve(data.size()+100);
+        }
+        data.resize(data.size()+1);
+        data[i] = input;
+        i++;
+      }
     }
     return COLVARS_OK;
   }
@@ -434,6 +462,12 @@ public:
     this->clear();
   }
 
+  /// Return a reference to the data
+  inline std::vector<T> &data_array()
+  {
+    return data;
+  }
+
   inline row & operator [] (size_t const i)
   {
     return rows[i];

src/.gitignore

@@ -303,6 +303,8 @@
 /fix_bond_create.h
 /fix_bond_swap.cpp
 /fix_bond_swap.h
+/fix_cmap.cpp
+/fix_cmap.h
 /fix_deposit.cpp
 /fix_deposit.h
 /fix_efield.cpp

src/KOKKOS/comm_tiled_kokkos.cpp

@@ -0,0 +1,266 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#include <string.h>
+#include "comm_tiled_kokkos.h"
+#include "comm_brick.h"
+#include "atom_kokkos.h"
+#include "atom_vec.h"
+#include "domain.h"
+#include "force.h"
+#include "pair.h"
+#include "neighbor.h"
+#include "modify.h"
+#include "fix.h"
+#include "compute.h"
+#include "output.h"
+#include "dump.h"
+#include "memory.h"
+#include "error.h"
+#include "atom_masks.h"
+
+using namespace LAMMPS_NS;
+
+#define BUFFACTOR 1.5
+#define BUFFACTOR 1.5
+#define BUFMIN 1000
+#define BUFEXTRA 1000
+#define EPSILON 1.0e-6
+
+#define DELTA_PROCS 16
+
+enum{SINGLE,MULTI};               // same as in Comm
+enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED};    // several files
+
+/* ---------------------------------------------------------------------- */
+
+CommTiledKokkos::CommTiledKokkos(LAMMPS *lmp) : CommTiled(lmp)
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+//IMPORTANT: we *MUST* pass "*oldcomm" to the Comm initializer here, as
+//           the code below *requires* that the (implicit) copy constructor
+//           for Comm is run and thus creating a shallow copy of "oldcomm".
+//           The call to Comm::copy_arrays() then converts the shallow copy
+//           into a deep copy of the class with the new layout.
+
+CommTiledKokkos::CommTiledKokkos(LAMMPS *lmp, Comm *oldcomm) : CommTiled(lmp,oldcomm)
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+CommTiledKokkos::~CommTiledKokkos()
+{
+
+}
+
+/* ---------------------------------------------------------------------- */
+
+
+/* ----------------------------------------------------------------------
+   forward communication of atom coords every timestep
+   other per-atom attributes may also be sent via pack/unpack routines
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm(int dummy)
+{
+  if (comm_x_only) {
+    atomKK->sync(Host,X_MASK);
+    atomKK->modified(Host,X_MASK);
+  } else if (ghost_velocity) {
+    atomKK->sync(Host,X_MASK | V_MASK);
+    atomKK->modified(Host,X_MASK | V_MASK);
+  } else {
+    atomKK->sync(Host,ALL_MASK);
+    atomKK->modified(Host,ALL_MASK);
+  }
+
+  CommTiled::forward_comm(dummy);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication of forces on atoms every timestep
+   other per-atom attributes may also be sent via pack/unpack routines
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm()
+{
+  if (comm_f_only)
+    atomKK->sync(Host,F_MASK);
+  else
+    atomKK->sync(Host,ALL_MASK);
+  CommTiled::reverse_comm();
+  if (comm_f_only)
+    atomKK->modified(Host,F_MASK);
+  else
+    atomKK->modified(Host,ALL_MASK);
+  atomKK->sync(Device,ALL_MASK);
+}
+
+/* ----------------------------------------------------------------------
+   exchange: move atoms to correct processors
+   atoms exchanged with procs that touch sub-box in each of 3 dims
+   send out atoms that have left my box, receive ones entering my box
+   atoms will be lost if not inside a touching proc's box
+     can happen if atom moves outside of non-periodic bounary
+     or if atom moves more than one proc away
+   this routine called before every reneighboring
+   for triclinic, atoms must be in lamda coords (0-1) before exchange is called
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::exchange()
+{
+  atomKK->sync(Host,ALL_MASK);
+  CommTiled::exchange();
+  atomKK->modified(Host,ALL_MASK);
+}
+
+/* ----------------------------------------------------------------------
+   borders: list nearby atoms to send to neighboring procs at every timestep
+   one list is created per swap/proc that will be made
+   as list is made, actually do communication
+   this does equivalent of a forward_comm(), so don't need to explicitly
+     call forward_comm() on reneighboring timestep
+   this routine is called before every reneighboring
+   for triclinic, atoms must be in lamda coords (0-1) before borders is called
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::borders()
+{
+  atomKK->sync(Host,ALL_MASK);
+  CommTiled::borders();
+  atomKK->modified(Host,ALL_MASK);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Pair
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_pair(Pair *pair)
+{
+  CommTiled::forward_comm_pair(pair);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Pair
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_pair(Pair *pair)
+{
+  CommTiled::reverse_comm_pair(pair);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Fix
+   size/nsize used only to set recv buffer limit
+   size = 0 (default) -> use comm_forward from Fix
+   size > 0 -> Fix passes max size per atom
+   the latter is only useful if Fix does several comm modes,
+     some are smaller than max stored in its comm_forward
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_fix(Fix *fix, int size)
+{
+  CommTiled::forward_comm_fix(fix,size);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Fix
+   size/nsize used only to set recv buffer limit
+   size = 0 (default) -> use comm_forward from Fix
+   size > 0 -> Fix passes max size per atom
+   the latter is only useful if Fix does several comm modes,
+     some are smaller than max stored in its comm_forward
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_fix(Fix *fix, int size)
+{
+  CommTiled::reverse_comm_fix(fix,size);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Fix with variable size data
+   query fix for all pack sizes to insure buf_send is big enough
+   handshake sizes before irregular comm to insure buf_recv is big enough
+   NOTE: how to setup one big buf recv with correct offsets ??
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_fix_variable(Fix *fix)
+{
+  CommTiled::reverse_comm_fix_variable(fix);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Compute
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_compute(Compute *compute)
+{
+  CommTiled::forward_comm_compute(compute);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Compute
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_compute(Compute *compute)
+{
+  CommTiled::reverse_comm_compute(compute);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication invoked by a Dump
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_dump(Dump *dump)
+{
+  CommTiled::forward_comm_dump(dump);
+}
+
+/* ----------------------------------------------------------------------
+   reverse communication invoked by a Dump
+   nsize used only to set recv buffer limit
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::reverse_comm_dump(Dump *dump)
+{
+  CommTiled::reverse_comm_dump(dump);
+}
+
+/* ----------------------------------------------------------------------
+   forward communication of Nsize values in per-atom array
+------------------------------------------------------------------------- */
+
+void CommTiledKokkos::forward_comm_array(int nsize, double **array)
+{
+  CommTiled::forward_comm_array(nsize,array);
+}
+
+/* ----------------------------------------------------------------------
+   exchange info provided with all 6 stencil neighbors
+   NOTE: this method is currently not used
+------------------------------------------------------------------------- */
+
+int CommTiledKokkos::exchange_variable(int n, double *inbuf, double *&outbuf)
+{
+  CommTiled::exchange_variable(n,inbuf,outbuf);
+}

src/KOKKOS/comm_tiled_kokkos.h

@@ -0,0 +1,59 @@
+/* -*- c++ -*- ----------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifndef LMP_COMM_TILED_KOKKOS_H
+#define LMP_COMM_TILED_KOKKOS_H
+
+#include "comm_tiled.h"
+#include "kokkos_type.h"
+
+namespace LAMMPS_NS {
+
+class CommTiledKokkos : public CommTiled {
+ public:
+  CommTiledKokkos(class LAMMPS *);
+  CommTiledKokkos(class LAMMPS *, class Comm *);
+  virtual ~CommTiledKokkos();
+
+  void forward_comm(int dummy = 0);    // forward comm of atom coords
+  void reverse_comm();                 // reverse comm of forces
+  void exchange();                     // move atoms to new procs
+  void borders();                      // setup list of atoms to comm
+
+  void forward_comm_pair(class Pair *);    // forward comm from a Pair
+  void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
+  void forward_comm_fix(class Fix *, int size=0);
+                                                   // forward comm from a Fix
+  void reverse_comm_fix(class Fix *, int size=0);
+                                                   // reverse comm from a Fix
+  void reverse_comm_fix_variable(class Fix *);
+                                     // variable size reverse comm from a Fix
+  void forward_comm_compute(class Compute *);  // forward from a Compute
+  void reverse_comm_compute(class Compute *);  // reverse from a Compute
+  void forward_comm_dump(class Dump *);    // forward comm from a Dump
+  void reverse_comm_dump(class Dump *);    // reverse comm from a Dump
+
+  void forward_comm_array(int, double **);          // forward comm of array
+  int exchange_variable(int, double *, double *&);  // exchange on neigh stencil
+
+ private:
+
+};
+
+}
+
+#endif
+
+/* ERROR/WARNING messages:
+
+*/

src/KOKKOS/neighbor_kokkos.cpp

@@ -280,15 +280,19 @@ void NeighborKokkos::choose_build(int index, NeighRequest *rq)
   if (rq->kokkos_host != 0) {
     PairPtrHost pb = NULL;
     if (rq->ghost) {
-      if (rq->full) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,1>;
-      else if (rq->half) &NeighborKokkos::full_bin_kokkos<LMPHostType,1,1>;
-      pair_build_host[index] = pb;
+      if (rq->full) {
+        if (rq->full_cluster) pb = &NeighborKokkos::full_bin_cluster_kokkos<LMPHostType>;
+        else pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,1>;
+      }
+      else if (rq->half) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,1,1>;
     } else {
-      if (rq->full) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,0>;
+      if (rq->full) {
+        if (rq->full_cluster) pb = &NeighborKokkos::full_bin_cluster_kokkos<LMPHostType>;
+        else pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,0,0>;
+      }
       else if (rq->half) pb = &NeighborKokkos::full_bin_kokkos<LMPHostType,1,0>;
-      pair_build_host[index] = pb;
     }
-    return;
+    pair_build_host[index] = pb;
   }
   if (rq->kokkos_device != 0) {
     PairPtrDevice pb = NULL;

src/MOLECULE/fix_cmap.cpp

@@ -13,11 +13,11 @@
 
 /* ----------------------------------------------------------------------
    Implementation of the CHARMM CMAP; adds an extra energy term for the
-   peptide backbone dihedrals.  The tools/ch2lmp.pl conversion script, which
-   generates an extra section in the LAMMPS data file, is needed in order to
-   generate the info used by this fix style.
+   peptide backbone dihedrals.  The tools/ch2lmp/charmm2lammps.pl
+   conversion script, which generates an extra section in the LAMMPS data
+   file, is needed in order to generate the info used by this fix style.
 
-   Contributing authors: 
+   Contributing authors:
    Xiaohu Hu, CMB/ORNL (hux2@ornl.gov)
    David Hyde-Volpe, Tigran Abramyan, and Robert A. Latour (Clemson University)
    Chris Lorenz (Kings College-London)
@@ -27,11 +27,11 @@
    - MacKerell et al., J. Comput. Chem. 25(2004):1400-1415.
  -------------------------------------------------------------------------*/
 
-#include "mpi.h"
-#include "math.h"
-#include "stdlib.h"
-#include "string.h"
-#include "stdio.h"
+#include <mpi.h>
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
 #include "fix_cmap.h"
 #include "atom.h"
 #include "atom_vec.h"
@@ -80,7 +80,7 @@ FixCMAP::FixCMAP(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg)
 
   MPI_Comm_rank(world,&me);
   MPI_Comm_size(world,&nprocs);
- 
+
   // allocate memory for CMAP data
 
   memory->create(g_axis,CMAPDIM,"cmap:g_axis");
@@ -123,8 +123,8 @@ FixCMAP::~FixCMAP()
 {
   // unregister callbacks to this fix from Atom class
 
-  atom->delete_callback(id,0); 
-  atom->delete_callback(id,1); 
+  atom->delete_callback(id,0);
+  atom->delete_callback(id,1);
 
   memory->destroy(g_axis);
   memory->destroy(cmapgrid);
@@ -173,7 +173,7 @@ void FixCMAP::init()
   }
 
   // pre-compute the derivatives of the maps
-  
+
   for (i = 0; i < 6; i++)
     set_map_derivatives(cmapgrid[i],d1cmapgrid[i],d2cmapgrid[i],d12cmapgrid[i]);
 
@@ -194,7 +194,7 @@ void FixCMAP::setup(int vflag)
     ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
     post_force_respa(vflag,nlevels_respa-1,0);
     ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
-  }  
+  }
 }
 
 /* --------------------------------------------------------------------- */
@@ -241,8 +241,8 @@ void FixCMAP::pre_neighbor()
       atom3 = atom->map(crossterm_atom3[i][m]);
       atom4 = atom->map(crossterm_atom4[i][m]);
       atom5 = atom->map(crossterm_atom5[i][m]);
-      
-      if (atom1 == -1 || atom2 == -1 || atom3 == -1 || 
+
+      if (atom1 == -1 || atom2 == -1 || atom3 == -1 ||
           atom4 == -1 || atom5 == -1) {
         char str[128];
         sprintf(str,"CMAP atoms "
@@ -260,7 +260,7 @@ void FixCMAP::pre_neighbor()
       atom4 = domain->closest_image(i,atom4);
       atom5 = domain->closest_image(i,atom5);
 
-      if (i <= atom1 && i <= atom2 && i <= atom3 && 
+      if (i <= atom1 && i <= atom2 && i <= atom3 &&
           i <= atom4 && i <= atom5) {
         if (ncrosstermlist == maxcrossterm) {
           maxcrossterm += LISTDELTA;
@@ -274,7 +274,7 @@ void FixCMAP::pre_neighbor()
         crosstermlist[ncrosstermlist][5] = crossterm_type[i][m];
         ncrosstermlist++;
       }
-    } 
+    }
   }
 }
 
@@ -304,7 +304,7 @@ void FixCMAP::post_force(int vflag)
   double dpr32r43,dpr45r43,r43,vb12x,vb12y,vb12z,vb54x,vb54y,vb54z;
   // cross-term dihedral angles
   double phi,psi,phi1,psi1;
-  double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6]; 
+  double f1[3],f2[3],f3[3],f4[3],f5[3],vcmap[6];
   double gs[4],d1gs[4],d2gs[4],d12gs[4];
   double engfraction;
   // vectors needed for the gradient/force calculation
@@ -316,12 +316,12 @@ void FixCMAP::post_force(int vflag)
   // Definition of cross-term dihedrals
 
   //         phi dihedral
-  //   |--------------------| 
+  //   |--------------------|
   //   a1-----a2-----a3-----a4-----a5    cross-term atoms
   //   C      N      CA     C      N     cross-term atom types
   //          |--------------------|
   //               psi dihedral
- 
+
   double **x = atom->x;
   double **f = atom->f;
   int nlocal = atom->nlocal;
@@ -342,7 +342,7 @@ void FixCMAP::post_force(int vflag)
     if (type == 0) continue;
 
     // calculate bond vectors for both dihedrals
-      
+
     // phi
     // vb21 = r2 - r1
 
@@ -369,50 +369,50 @@ void FixCMAP::post_force(int vflag)
       vb43x = -1.0*vb34x;
       vb43y = -1.0*vb34y;
       vb43z = -1.0*vb34z;
-    
-      vb45x = x[i4][0] - x[i5][0]; 
-      vb45y = x[i4][1] - x[i5][1]; 
-      vb45z = x[i4][2] - x[i5][2]; 
+
+      vb45x = x[i4][0] - x[i5][0];
+      vb45y = x[i4][1] - x[i5][1];
+      vb45z = x[i4][2] - x[i5][2];
       vb54x = -1.0*vb45x;
       vb54y = -1.0*vb45y;
       vb54z = -1.0*vb45z;
- 
+
       // calculate normal vectors for planes that define the dihedral angles
-      
+
       a1x = vb12y*vb23z - vb12z*vb23y;
       a1y = vb12z*vb23x - vb12x*vb23z;
       a1z = vb12x*vb23y - vb12y*vb23x;
-      
+
       b1x = vb43y*vb23z - vb43z*vb23y;
       b1y = vb43z*vb23x - vb43x*vb23z;
       b1z = vb43x*vb23y - vb43y*vb23x;
-      
+
       a2x = vb23y*vb34z - vb23z*vb34y;
       a2y = vb23z*vb34x - vb23x*vb34z;
       a2z = vb23x*vb34y - vb23y*vb34x;
-      
+
       b2x = vb45y*vb43z - vb45z*vb43y;
       b2y = vb45z*vb43x - vb45x*vb43z;
       b2z = vb45x*vb43y - vb45y*vb43x;
-      
+
       // calculate terms used later in calculations
-      
+
       r32 = sqrt(vb32x*vb32x + vb32y*vb32y + vb32z*vb32z);
       a1sq = a1x*a1x + a1y*a1y + a1z*a1z;
       b1sq = b1x*b1x + b1y*b1y + b1z*b1z;
-      
+
       r43 = sqrt(vb43x*vb43x + vb43y*vb43y + vb43z*vb43z);
       a2sq = a2x*a2x + a2y*a2y + a2z*a2z;
       b2sq = b2x*b2x + b2y*b2y + b2z*b2z;
-      //if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) 
+      //if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001)
       //  printf("a1sq b1sq a2sq b2sq: %f %f %f %f \n",a1sq,b1sq,a2sq,b2sq);
       if (a1sq<0.0001 || b1sq<0.0001 || a2sq<0.0001 || b2sq<0.0001) continue;
       dpr21r32 = vb21x*vb32x + vb21y*vb32y + vb21z*vb32z;
       dpr34r32 = vb34x*vb32x + vb34y*vb32y + vb34z*vb32z;
       dpr32r43 = vb32x*vb43x + vb32y*vb43y + vb32z*vb43z;
       dpr45r43 = vb45x*vb43x + vb45y*vb43y + vb45z*vb43z;
-      
-      // calculate the backbone dihedral angles as VMD and GROMACS 
+
+      // calculate the backbone dihedral angles as VMD and GROMACS
 
       phi = dihedral_angle_atan2(vb21x,vb21y,vb21z,a1x,a1y,a1z,b1x,b1y,b1z,r32);
       psi = dihedral_angle_atan2(vb32x,vb32y,vb32z,a2x,a2y,a2z,b2x,b2y,b2z,r43);
@@ -425,7 +425,7 @@ void FixCMAP::post_force(int vflag)
       psi1 = psi;
       if (psi1 < 0.0) psi1 += 360.0;
 
-      // find the neighbor grid point index 
+      // find the neighbor grid point index
 
       li1 = int(((phi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));
       li2 = int(((psi1+CMAPXMIN2)/CMAPDX)+((CMAPDIM*1.0)/2.0));
@@ -436,87 +436,87 @@ void FixCMAP::post_force(int vflag)
       mli4 = li4 % CMAPDIM;
       mli31 = (li3+1) % CMAPDIM;
       mli41 = (li4+1)  %CMAPDIM;
-      mli1 = li1 % CMAPDIM; 
+      mli1 = li1 % CMAPDIM;
       mli2 = li2 % CMAPDIM;
       mli11 = (li1+1) % CMAPDIM;
       mli21 = (li2+1)  %CMAPDIM;
       t1 = type-1;
       if (t1 < 0 || t1 > 5) error->all(FLERR,"Invalid CMAP crossterm_type");
-      
-      // determine the values and derivatives for the grid square points 
+
+      // determine the values and derivatives for the grid square points
 
       gs[0] = cmapgrid[t1][mli3][mli4];
       gs[1] = cmapgrid[t1][mli31][mli4];
       gs[2] = cmapgrid[t1][mli31][mli41];
       gs[3] = cmapgrid[t1][mli3][mli41];
-      d1gs[0] = d1cmapgrid[t1][mli1][mli2];  
+      d1gs[0] = d1cmapgrid[t1][mli1][mli2];
       d1gs[1] = d1cmapgrid[t1][mli11][mli2];
       d1gs[2] = d1cmapgrid[t1][mli11][mli21];
       d1gs[3] = d1cmapgrid[t1][mli1][mli21];
-      
-      d2gs[0] = d2cmapgrid[t1][mli1][mli2];  
+
+      d2gs[0] = d2cmapgrid[t1][mli1][mli2];
       d2gs[1] = d2cmapgrid[t1][mli11][mli2];
       d2gs[2] = d2cmapgrid[t1][mli11][mli21];
       d2gs[3] = d2cmapgrid[t1][mli1][mli21];
-      
-      d12gs[0] = d12cmapgrid[t1][mli1][mli2];  
+
+      d12gs[0] = d12cmapgrid[t1][mli1][mli2];
       d12gs[1] = d12cmapgrid[t1][mli11][mli2];
       d12gs[2] = d12cmapgrid[t1][mli11][mli21];
       d12gs[3] = d12cmapgrid[t1][mli1][mli21];
-      
-      // calculate the cmap energy and the gradient (dE/dphi,dE/dpsi) 
-      
-      bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs);
-      
-      // sum up cmap energy contributions
-      
-      engfraction = 0.2 * E;
-      if (i1 < nlocal) ecmap += engfraction; 
-      if (i2 < nlocal) ecmap += engfraction; 
-      if (i3 < nlocal) ecmap += engfraction; 
-      if (i4 < nlocal) ecmap += engfraction; 
-      if (i5 < nlocal) ecmap += engfraction; 
 
-      // calculate the derivatives dphi/dr_i 
-      
+      // calculate the cmap energy and the gradient (dE/dphi,dE/dpsi)
+
+      bc_interpol(phi,psi,li3,li4,gs,d1gs,d2gs,d12gs);
+
+      // sum up cmap energy contributions
+
+      engfraction = 0.2 * E;
+      if (i1 < nlocal) ecmap += engfraction;
+      if (i2 < nlocal) ecmap += engfraction;
+      if (i3 < nlocal) ecmap += engfraction;
+      if (i4 < nlocal) ecmap += engfraction;
+      if (i5 < nlocal) ecmap += engfraction;
+
+      // calculate the derivatives dphi/dr_i
+
       dphidr1x = 1.0*r32/a1sq*a1x;
       dphidr1y = 1.0*r32/a1sq*a1y;
       dphidr1z = 1.0*r32/a1sq*a1z;
-      
-      dphidr2x = -1.0*r32/a1sq*a1x - dpr21r32/a1sq/r32*a1x + 
+
+      dphidr2x = -1.0*r32/a1sq*a1x - dpr21r32/a1sq/r32*a1x +
         dpr34r32/b1sq/r32*b1x;
-      dphidr2y = -1.0*r32/a1sq*a1y - dpr21r32/a1sq/r32*a1y + 
+      dphidr2y = -1.0*r32/a1sq*a1y - dpr21r32/a1sq/r32*a1y +
         dpr34r32/b1sq/r32*b1y;
-      dphidr2z = -1.0*r32/a1sq*a1z - dpr21r32/a1sq/r32*a1z + 
+      dphidr2z = -1.0*r32/a1sq*a1z - dpr21r32/a1sq/r32*a1z +
         dpr34r32/b1sq/r32*b1z;
 
       dphidr3x = dpr34r32/b1sq/r32*b1x - dpr21r32/a1sq/r32*a1x - r32/b1sq*b1x;
       dphidr3y = dpr34r32/b1sq/r32*b1y - dpr21r32/a1sq/r32*a1y - r32/b1sq*b1y;
       dphidr3z = dpr34r32/b1sq/r32*b1z - dpr21r32/a1sq/r32*a1z - r32/b1sq*b1z;
-    
+
       dphidr4x = r32/b1sq*b1x;
       dphidr4y = r32/b1sq*b1y;
       dphidr4z = r32/b1sq*b1z;
 
-      // calculate the derivatives dpsi/dr_i 
-      
+      // calculate the derivatives dpsi/dr_i
+
       dpsidr1x = 1.0*r43/a2sq*a2x;
       dpsidr1y = 1.0*r43/a2sq*a2y;
       dpsidr1z = 1.0*r43/a2sq*a2z;
-      
+
       dpsidr2x = r43/a2sq*a2x + dpr32r43/a2sq/r43*a2x - dpr45r43/b2sq/r43*b2x;
       dpsidr2y = r43/a2sq*a2y + dpr32r43/a2sq/r43*a2y - dpr45r43/b2sq/r43*b2y;
       dpsidr2z = r43/a2sq*a2z + dpr32r43/a2sq/r43*a2z - dpr45r43/b2sq/r43*b2z;
-      
+
       dpsidr3x = dpr45r43/b2sq/r43*b2x - dpr32r43/a2sq/r43*a2x - r43/b2sq*b2x;
       dpsidr3y = dpr45r43/b2sq/r43*b2y - dpr32r43/a2sq/r43*a2y - r43/b2sq*b2y;
       dpsidr3z = dpr45r43/b2sq/r43*b2z - dpr32r43/a2sq/r43*a2z - r43/b2sq*b2z;
-      
+
       dpsidr4x = r43/b2sq*b2x;
       dpsidr4y = r43/b2sq*b2y;
       dpsidr4z = r43/b2sq*b2z;
 
-      // calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr) 
+      // calculate forces on cross-term atoms: F = -(dE/dPhi)*(dPhi/dr)
 
       f1[0] = dEdPhi*dphidr1x;
       f1[1] = dEdPhi*dphidr1y;
@@ -533,9 +533,9 @@ void FixCMAP::post_force(int vflag)
       f5[0] = -dEdPsi*dpsidr4x;
       f5[1] = -dEdPsi*dpsidr4y;
       f5[2] = -dEdPsi*dpsidr4z;
-      
+
       // apply force to each of the 5 atoms
-      
+
       if (i1 < nlocal) {
         f[i1][0] += f1[0];
         f[i1][1] += f1[1];
@@ -561,7 +561,7 @@ void FixCMAP::post_force(int vflag)
         f[i5][1] += f5[1];
         f[i5][2] += f5[2];
       }
-      
+
       // tally energy and/or virial
 
       if (evflag) {
@@ -622,15 +622,18 @@ double FixCMAP::compute_scalar()
 
 void FixCMAP::read_grid_map(char *cmapfile)
 {
-  char line[MAXLINE];
-  char *chunk;
+  char linebuf[MAXLINE];
+  char *chunk,*line;
   int i1, i2, i3, i4, i5, i6, j1, j2, j3, j4, j5, j6, counter;
-  
-  FILE *fp = fopen(cmapfile,"r");
-  if (fp == NULL) {
-    char str[128];
-    sprintf(str,"Cannot open fix cmap file %s",cmapfile);
-    error->one(FLERR,str);
+
+  FILE *fp = NULL;
+  if (comm->me == 0) {
+    fp = force->open_potential(cmapfile);
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open fix cmap file %s",cmapfile);
+      error->one(FLERR,str);
+    }
   }
 
   for (int ix1 = 0; ix1 < 6; ix1++)
@@ -642,8 +645,24 @@ void FixCMAP::read_grid_map(char *cmapfile)
   i1 = i2 = i3 = i4 = i5 = i6 = 0;
   j1 = j2 = j3 = j4 = j5 = j6 = 0;
 
-  while (fgets(line,MAXLINE,fp) != NULL) {
-    if (line == "" || line[0] == '#') { ;; }
+  int done = 0;
+
+  while (!done) {
+    // only read on rank 0 and broadcast to all other ranks
+    if (comm->me == 0)
+      done = (fgets(linebuf,MAXLINE,fp) == NULL);
+
+    MPI_Bcast(&done,1,MPI_INT,0,world);
+    if (done) continue;
+
+    MPI_Bcast(linebuf,MAXLINE,MPI_CHAR,0,world);
+
+    // remove leading whitespace
+    line = linebuf;
+    while (line && (*line == ' ' || *line == '\t' || *line == '\r')) ++line;
+
+    // skip if empty line or comment
+    if (!line || *line =='\n' || *line == '\0' || *line == '#') continue;
 
     // read in the cmap grid point values
     // NOTE: The order to read the 6 grid maps is HARD-CODED, thus errors
@@ -655,101 +674,99 @@ void FixCMAP::read_grid_map(char *cmapfile)
     // 3. Proline map
     // 4. Two adjacent prolines map
     // 5. Glycine map
-    // 6. Glycine before proline map	
-    
-    else {
-      chunk = strtok(line, " \r\n");
-      while (chunk != NULL) {
+    // 6. Glycine before proline map
 
-        // alanine map 
-        
-        if (counter < CMAPDIM*CMAPDIM) {
-          cmapgrid[0][i1][j1] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j1++;
-          if (j1 == CMAPDIM) {
-            j1 = 0;
-            i1++;
-          }
-          counter++;
+    chunk = strtok(line, " \r\n");
+    while (chunk != NULL) {
+
+      // alanine map
+
+      if (counter < CMAPDIM*CMAPDIM) {
+        cmapgrid[0][i1][j1] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j1++;
+        if (j1 == CMAPDIM) {
+          j1 = 0;
+          i1++;
         }
-        
-        // alanine-proline map 
-        
-        else if (counter >= CMAPDIM*CMAPDIM &&
-                 counter < 2*CMAPDIM*CMAPDIM) {
-          cmapgrid[1][i2][j2]= atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j2++;
-          if (j2 == CMAPDIM) {
-            j2 = 0;
-            i2++;
-          }
-          counter++;
-        }
-        
-        // proline map 
-        
-        else if (counter >= 2*CMAPDIM*CMAPDIM &&
-                 counter < 3*CMAPDIM*CMAPDIM) {
-          cmapgrid[2][i3][j3] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j3++;
-          if (j3 == CMAPDIM) {
-            j3 = 0;
-            i3++;
-          }
-          counter++;
-        }
-        
-        // 2 adjacent prolines map 
-        
-        else if (counter >= 3*CMAPDIM*CMAPDIM &&
-                 counter < 4*CMAPDIM*CMAPDIM) {
-          cmapgrid[3][i4][j4] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j4++;
-          if (j4 == CMAPDIM) {
-            j4 = 0;
-            i4++;
-          }
-          counter++;
-        }
-        
-        // glycine map 
-        
-        else if (counter >= 4*CMAPDIM*CMAPDIM &&
-                 counter < 5*CMAPDIM*CMAPDIM) {
-          cmapgrid[4][i5][j5] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j5++;
-          if (j5 == CMAPDIM) {
-            j5 = 0;
-            i5++;
-          }
-          counter++;
-        }
-        
-        // glycine-proline map 
-        
-        else if (counter >= 5*CMAPDIM*CMAPDIM &&
-                 counter < 6*CMAPDIM*CMAPDIM) {
-          cmapgrid[5][i6][j6] = atof(chunk);
-          chunk = strtok(NULL, " \r\n");
-          j6++;
-          if (j6 == CMAPDIM) {
-            j6 = 0;
-            i6++;
-          }
-          counter++;
-        }
-        
-        else break;
+        counter++;
       }
+
+      // alanine-proline map
+
+      else if (counter >= CMAPDIM*CMAPDIM &&
+               counter < 2*CMAPDIM*CMAPDIM) {
+        cmapgrid[1][i2][j2]= atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j2++;
+        if (j2 == CMAPDIM) {
+          j2 = 0;
+          i2++;
+        }
+        counter++;
+      }
+
+      // proline map
+
+      else if (counter >= 2*CMAPDIM*CMAPDIM &&
+               counter < 3*CMAPDIM*CMAPDIM) {
+        cmapgrid[2][i3][j3] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j3++;
+        if (j3 == CMAPDIM) {
+          j3 = 0;
+          i3++;
+        }
+        counter++;
+      }
+
+      // 2 adjacent prolines map
+
+      else if (counter >= 3*CMAPDIM*CMAPDIM &&
+               counter < 4*CMAPDIM*CMAPDIM) {
+        cmapgrid[3][i4][j4] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j4++;
+        if (j4 == CMAPDIM) {
+          j4 = 0;
+          i4++;
+        }
+        counter++;
+      }
+
+      // glycine map
+
+      else if (counter >= 4*CMAPDIM*CMAPDIM &&
+               counter < 5*CMAPDIM*CMAPDIM) {
+        cmapgrid[4][i5][j5] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j5++;
+        if (j5 == CMAPDIM) {
+          j5 = 0;
+          i5++;
+        }
+        counter++;
+      }
+
+      // glycine-proline map
+
+      else if (counter >= 5*CMAPDIM*CMAPDIM &&
+               counter < 6*CMAPDIM*CMAPDIM) {
+        cmapgrid[5][i6][j6] = atof(chunk);
+        chunk = strtok(NULL, " \r\n");
+        j6++;
+        if (j6 == CMAPDIM) {
+          j6 = 0;
+          i6++;
+        }
+        counter++;
+      }
+
+      else break;
     }
   }
-  
-  fclose(fp);
+
+  if (comm->me == 0) fclose(fp);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -783,18 +800,18 @@ void FixCMAP::spline(double *y, double *ddy, int n)
 /* ---------------------------------------------------------------------- */
 
 void FixCMAP::spl_interpolate(double x, double *y, double *ddy, double &yo,
-		double &dyo)
+                              double &dyo)
 {
   // perform a 1D cubic spline interpolation
 
   int ix;
   double a,b,a1,b1,a2,b2;
-  
+
   ix = int((x-CMAPXMIN)/CMAPDX-(1./2.));
 
   a = (CMAPXMIN+(ix*1.0)*CMAPDX-x)/CMAPDX;
   b = (x-CMAPXMIN-(((ix-1)*1.0)*CMAPDX))/CMAPDX;
-  
+
   a1 = a*a*a-a;
   b1 = b*b*b-b;
 
@@ -807,7 +824,7 @@ void FixCMAP::spl_interpolate(double x, double *y, double *ddy, double &yo,
 /* ---------------------------------------------------------------------- */
 
 void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
-                                   double **d12yo)
+                                  double **d12yo)
 {
   // precompute the gradient and cross-derivatives of the map grid points.
   // use the bicubic spline to calculate the derivatives
@@ -833,7 +850,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
   memory->create(tddmap,CMAPDIM*2,CMAPDIM*2,"cmap:tddmap");
 
   // periodically expand the original map
-  // use the expanded map for bicubic spline interpolation, 
+  // use the expanded map for bicubic spline interpolation,
   //   which is used to obtain the derivatives
   // actual interpolation is done with bicubic interpolation
 
@@ -844,7 +861,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
       tmap[i][j] = map[ii][jj];
     }
   }
-  
+
   for (i = 0; i < CMAPDIM*2; i++)
     spline(tmap[i], tddmap[i], CMAPDIM*2);
 
@@ -869,7 +886,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
         tmp_y[k] = tyyk;
         tmp_dy[k] = tdyk;
       }
-      
+
       spline(tmp_y,tmp_ddy,CMAPDIM+xm+xm);
       ix = int((phi-CMAPXMIN)/CMAPDX);
       a = (CMAPXMIN+((ix+1)*1.0)*CMAPDX-phi)/CMAPDX;
@@ -899,7 +916,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
       d12yo[i%p][j%p] = d12y;
     }
   }
-  
+
   memory->destroy(tmp_y);
   memory->destroy(tmp_dy);
   memory->destroy(tmp_ddy);
@@ -911,7 +928,7 @@ void FixCMAP::set_map_derivatives(double **map, double **d1yo, double **d2yo,
 
 double FixCMAP::dihedral_angle_atan2(double fx, double fy, double fz,
                                       double ax, double ay, double az,
-                                      double bx, double by, double bz, 
+                                      double bx, double by, double bz,
                                       double absg)
 {
   // calculate the dihedral angle
@@ -953,19 +970,19 @@ void FixCMAP::bc_coeff(double *gs, double *d1gs, double *d2gs, double *d12gs)
       2,-2, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
       0, 0, 0, 0, 0, 0, 0, 0, 2,-2, 0, 0, 1, 1, 0, 0,
       -6, 6,-6, 6,-3,-3, 3, 3,-4, 4, 2,-2,-2,-2,-1,-1,
-      4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1 
+      4,-4, 4,-4, 2, 2,-2,-2, 2,-2,-2, 2, 1, 1, 1, 1
     };
 
-  int i, j, k, l, in;
+  int i, j, k, in;
   double xx, x[16];
-  
+
   for (i = 0; i < 4; i++) {
     x[i] = gs[i];
     x[i+4] = d1gs[i]*CMAPDX;
     x[i+8] = d2gs[i]*CMAPDX;
     x[i+12] = d12gs[i]*CMAPDX*CMAPDX;
   }
-  
+
   in = 0;
   for (i = 0; i < 4; i++) {
     for (j = 0; j < 4; j++) {
@@ -1006,7 +1023,7 @@ void FixCMAP::bc_interpol(double x1, double x2, int low1, int low2, double *gs,
     dEdPhi = u*dEdPhi + (3.0*cij[3][i]*t+2.0*cij[2][i])*t+cij[1][i];
     dEdPsi = t*dEdPsi + (3.0*cij[i][3]*u+2.0*cij[i][2])*u+cij[i][1];
   }
-  
+
   dEdPhi *= (180.0/MY_PI/CMAPDX);
   dEdPsi *= (180.0/MY_PI/CMAPDX);
 }
@@ -1023,7 +1040,7 @@ void FixCMAP::read_data_header(char *line)
     sscanf(line,BIGINT_FORMAT,&ncmap);
   } else error->all(FLERR,"Invalid read data header line for fix cmap");
 
-  // didn't set in constructor b/c this fix could be defined 
+  // didn't set in constructor b/c this fix could be defined
   // before newton command
 
   newton_bond = force->newton_bond;
@@ -1117,7 +1134,7 @@ void FixCMAP::read_data_section(char *keyword, int n, char *buf,
       crossterm_atom5[m][num_crossterm[m]] = atom5;
       num_crossterm[m]++;
     }
-    
+
     if ((m = atom->map(atom5)) >= 0) {
       if (num_crossterm[m] == CMAPMAX)
         error->one(FLERR,"Too many CMAP crossterms for one atom");
@@ -1226,7 +1243,7 @@ void FixCMAP::write_data_section(int mth, FILE *fp,
                                   int n, double **buf, int index)
 {
   for (int i = 0; i < n; i++)
-    fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT 
+    fprintf(fp,"%d %d " TAGINT_FORMAT " " TAGINT_FORMAT
             " " TAGINT_FORMAT " " TAGINT_FORMAT " " TAGINT_FORMAT "\n",
             index+i,(int) ubuf(buf[i][0]).i,(tagint) ubuf(buf[i][1]).i,
             (tagint) ubuf(buf[i][2]).i,(tagint) ubuf(buf[i][3]).i,
@@ -1261,8 +1278,8 @@ void FixCMAP::restart(char *buf)
   ncmap = *((bigint *) buf);
 }
 
-/* ---------------------------------------------------------------------- 
-   pack values in local atom-based arrays for restart file 
+/* ----------------------------------------------------------------------
+   pack values in local atom-based arrays for restart file
 ------------------------------------------------------------------------- */
 
 int FixCMAP::pack_restart(int i, double *buf)
@@ -1275,14 +1292,14 @@ int FixCMAP::pack_restart(int i, double *buf)
     buf[n++] = ubuf(crossterm_atom3[i][m]).d;
     buf[n++] = ubuf(crossterm_atom4[i][m]).d;
     buf[n++] = ubuf(crossterm_atom5[i][m]).d;
-  }  
+  }
   buf[0] = n;
 
-  return n; 
+  return n;
 }
 
-/* ---------------------------------------------------------------------- 
-   unpack values from atom->extra array to restart the fix 
+/* ----------------------------------------------------------------------
+   unpack values from atom->extra array to restart the fix
 ------------------------------------------------------------------------- */
 
 void FixCMAP::unpack_restart(int nlocal, int nth)
@@ -1290,9 +1307,9 @@ void FixCMAP::unpack_restart(int nlocal, int nth)
   double **extra = atom->extra;
 
   // skip to Nth set of extra values
-   
+
    int n = 0;
-   for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]); 
+   for (int i = 0; i < nth; i++) n += static_cast<int> (extra[nlocal][n]);
 
    int count = static_cast<int> (extra[nlocal][n++]);
    num_crossterm[nlocal] = (count-1)/6;
@@ -1307,17 +1324,17 @@ void FixCMAP::unpack_restart(int nlocal, int nth)
    }
 }
 
-/* ---------------------------------------------------------------------- 
-   maxsize of any atom's restart data 
+/* ----------------------------------------------------------------------
+   maxsize of any atom's restart data
 ------------------------------------------------------------------------- */
 
 int FixCMAP::maxsize_restart()
 {
-  return 1 + CMAPMAX*6; 
+  return 1 + CMAPMAX*6;
 }
 
-/* ---------------------------------------------------------------------- 
-   size of atom nlocal's restart data 
+/* ----------------------------------------------------------------------
+   size of atom nlocal's restart data
 ------------------------------------------------------------------------- */
 
 int FixCMAP::size_restart(int nlocal)
@@ -1330,7 +1347,7 @@ int FixCMAP::size_restart(int nlocal)
 ------------------------------------------------------------------------- */
 
 void FixCMAP::grow_arrays(int nmax)
-{ 
+{
   num_crossterm = memory->grow(num_crossterm,nmax,"cmap:num_crossterm");
   crossterm_type = memory->grow(crossterm_type,nmax,CMAPMAX,
                                 "cmap:crossterm_type");
@@ -1382,7 +1399,7 @@ void FixCMAP::set_arrays(int i)
 /* ----------------------------------------------------------------------
    pack values in local atom-based array for exchange with another proc
 ------------------------------------------------------------------------- */
-  
+
 int FixCMAP::pack_exchange(int i, double *buf)
 {
   int n = 0;

src/MOLECULE/fix_cmap.h

@@ -1,4 +1,4 @@
-/* ----------------------------------------------------------------------
+/* -*- c++ -*- ----------------------------------------------------------
    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
    http://lammps.sandia.gov, Sandia National Laboratories
    Steve Plimpton, sjplimp@sandia.gov

src/Make.py

@@ -26,7 +26,8 @@ libclasses = ("atc","awpmd","colvars","cuda","gpu","h5md",
 buildclasses = ("intel","kokkos")
 makeclasses = ("cc","flags","mpi","fft","jpg","png")
 
-setargs = ("gzip","#gzip","ffmpeg","#ffmpeg","smallbig","bigbig","smallsmall","exceptions","#exceptions")
+setargs = ("gzip","#gzip","ffmpeg","#ffmpeg","smallbig","bigbig",
+           "smallsmall","exceptions","#exceptions")
 actionargs = ("lib-all","file","clean","exe")
 
 gpubuildflag = 0
@@ -85,7 +86,9 @@ def switch2str(switches,switch_order):
 def compile_check(compiler,ccflags,warn):
   open("tmpauto.cpp",'w').write("int main(int, char **) {}\n")
   tmp = "%s %s -c tmpauto.cpp" % (compiler,ccflags)
-  txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,shell=True).decode()
+  try: txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,
+                                     shell=True).decode()
+  except subprocess.CalledProcessError as e: txt = e.output
   flag = 1
   if txt or not os.path.isfile("tmpauto.o"):
     flag = 0
@@ -104,7 +107,9 @@ def compile_check(compiler,ccflags,warn):
 def link_check(linker,linkflags,libs,warn):
   open("tmpauto.cpp",'w').write("int main(int, char **) {}\n")
   tmp = "%s %s -o tmpauto tmpauto.cpp %s" % (linker,linkflags,libs)
-  txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,shell=True).decode()
+  try: txt = subprocess.check_output(tmp,stderr=subprocess.STDOUT,
+                                     shell=True).decode()
+  except subprocess.CalledProcessError as e: txt = e.output
   flag = 1
   if txt or not os.path.isfile("tmpauto"):
     flag = 0
@@ -530,7 +535,8 @@ class Actions(object):
     
     if caller == "file" or "file" not in self.alist:
       # make certain that 'MAKE/MINE' folder exists.
-      subprocess.check_output("mkdir -p %s/MAKE/MINE" % dir.src,stderr=subprocess.STDOUT,shell=True)
+      subprocess.check_output("mkdir -p %s/MAKE/MINE" % dir.src,
+                              stderr=subprocess.STDOUT,shell=True)
       make.write("%s/MAKE/MINE/Makefile.auto" % dir.src,1)
       print("Created src/MAKE/MINE/Makefile.auto")
       
@@ -588,7 +594,7 @@ class Actions(object):
     else:
       print(tmp)
       try: subprocess.check_output(tmp,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/lmp_auto" % dir.src):
       error('Unsuccessful "make auto"')
@@ -857,9 +863,12 @@ class Packages(object):
     if self.plist: print("Installing packages ...")
     for one in self.plist:
       if one == "orig": continue
-      subprocess.check_output("cd %s; make %s" % (dir.src,one),stderr=subprocess.STDOUT,shell=True)
+      subprocess.check_output("cd %s; make %s" % (dir.src,one),
+                              stderr=subprocess.STDOUT,shell=True)
     if self.plist and verbose:
-      txt = subprocess.check_output("cd %s; make ps" % dir.src,stderr=subprocess.STDOUT,shell=True).decode()
+      txt = subprocess.check_output("cd %s; make ps" % dir.src,
+                                    stderr=subprocess.STDOUT,
+                                    shell=True).decode()
       print("Package status after installation:")
       print(txt)
       
@@ -869,12 +878,16 @@ class Packages(object):
   def uninstall(self):
     if not self.plist or self.plist[-1] != "orig": return
     print("Restoring packages to original state ...")
-    subprocess.check_output("cd %s; make no-all" % dir.src,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make no-all" % dir.src,
+                            stderr=subprocess.STDOUT,shell=True)
     for one in self.all:
       if self.original[one]:
-        subprocess.check_output("cd %s; make yes-%s" % (dir.src,one),stderr=subprocess.STDOUT,shell=True)
+        subprocess.check_output("cd %s; make yes-%s" % (dir.src,one),
+                                stderr=subprocess.STDOUT,shell=True)
     if verbose:
-      txt = subprocess.check_output("cd %s; make ps" % dir.src,stderr=subprocess.STDOUT,shell=True).decode()
+      txt = subprocess.check_output("cd %s; make ps" % dir.src,
+                                    stderr=subprocess.STDOUT,
+                                    shell=True).decode()
       print("Restored package status:")
       print(txt)
       
@@ -968,7 +981,8 @@ class Settings(object):
   def help(self):
     return """
 -s set1 set2 ...
-  possible settings = gzip #gzip ffmpeg #ffmpeg smallbig bigbig smallsmall exceptions #exceptions
+  possible settings = gzip #gzip ffmpeg #ffmpeg
+                      smallbig bigbig smallsmall exceptions #exceptions
   alter LAMMPS ifdef settings in Makefile.auto
     only happens if new Makefile.auto is created by use of "file" action
   gzip and #gzip turn on/off LAMMPS_GZIP setting
@@ -1057,7 +1071,8 @@ class ATC(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
     
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
     else: txt = "cd %s; make -f Makefile.auto" % libdir
     
@@ -1066,7 +1081,7 @@ class ATC(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libatc.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1108,7 +1123,8 @@ class AWPMD(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
     else: txt = "cd %s; make -f Makefile.auto" % libdir
 
@@ -1117,7 +1133,7 @@ class AWPMD(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
    
     if not os.path.isfile("%s/libawpmd.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1159,7 +1175,8 @@ class COLVARS(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
     else: txt = "cd %s; make -f Makefile.auto" % libdir
 
@@ -1168,7 +1185,7 @@ class COLVARS(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libcolvars.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1211,7 +1228,8 @@ class CUDA(object):
           
   def build(self): 
     libdir = dir.lib + "/cuda"
-    subprocess.check_output("cd %s; make clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make clean" % libdir,
+                            stderr=subprocess.STDOUT,shell=True)
     if self.mode == "double": n = 2
     elif self.mode == "mixed": n = 3
     elif self.mode == "single": n = 1
@@ -1225,7 +1243,7 @@ class CUDA(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/liblammpscuda.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1306,7 +1324,8 @@ class GPU(object):
     make = "make"
     if "shannon" == platform.node(): make = "srun make"
 
-    subprocess.check_output("cd %s; %s -f Makefile.auto clean" % (libdir,make),stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; %s -f Makefile.auto clean" %
+                            (libdir,make),stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; %s -j %d -f Makefile.auto" % (libdir,make,jmake.n)
     else: txt = "cd %s; %s -f Makefile.auto" % (libdir,make)
 
@@ -1315,7 +1334,7 @@ class GPU(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libgpu.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1357,7 +1376,8 @@ class H5MD(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    subprocess.check_output("cd %s; make clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make clean" % libdir,
+                            stderr=subprocess.STDOUT,shell=True)
     txt = "cd %s; make" % libdir
 
     # if verbose, print output as build proceeds, else only print if fails
@@ -1365,7 +1385,7 @@ class H5MD(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libch5md.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1407,7 +1427,8 @@ class MEAM(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
     # do not use -j for MEAM build, parallel build does not work
     txt = "cd %s; make -f Makefile.auto" % libdir
 
@@ -1416,7 +1437,7 @@ class MEAM(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libmeam.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1458,7 +1479,8 @@ class POEMS(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,
+                            stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
     else: txt = "cd %s; make -f Makefile.auto" % libdir
 
@@ -1467,7 +1489,7 @@ class POEMS(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libpoems.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1504,9 +1526,10 @@ class PYTHON(object):
     libdir = dir.lib + "/python"
     if self.lammpsflag:
       subprocess.check_output("cd %s; cp Makefile.lammps.%s Makefile.lammps" %
-                         (libdir,self.lammps))
+                              (libdir,self.lammps))
     if not os.path.isfile("%s/Makefile.lammps.%s" % (libdir,self.lammps)):
-      error("Unsuccessful creation of lib/python/Makefile.lammps.%s file" % self.lammps)
+      error("Unsuccessful creation of lib/python/Makefile.lammps.%s file" %
+            self.lammps)
     else: print("Created lib/python/Makefile.lammps file")
 
 # QMMM lib
@@ -1544,7 +1567,8 @@ class QMMM(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    subprocess.check_output("cd %s; make -f Makefile.auto clean" % libdir,stderr=subprocess.STDOUT,shell=True)
+    subprocess.check_output("cd %s; make -f Makefile.auto clean" %
+                            libdir,stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
     else: txt = "cd %s; make -f Makefile.auto" % libdir
 
@@ -1553,7 +1577,7 @@ class QMMM(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
    
     if not os.path.isfile("%s/libqmmm.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1595,7 +1619,8 @@ class REAX(object):
       make.setvar("EXTRAMAKE","Makefile.lammps.%s" % self.lammps)
     make.write("%s/Makefile.auto" % libdir)
 
-    commands.getoutput("cd %s; make -f Makefile.auto clean" % libdir)
+    cmd = "cd %s; make -f Makefile.auto clean" % libdir
+    subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True)
     if jmake: txt = "cd %s; make -j %d -f Makefile.auto" % (libdir,jmake.n)
     else: txt = "cd %s; make -f Makefile.auto" % libdir
     
@@ -1604,7 +1629,7 @@ class REAX(object):
     if verbose: subprocess.call(txt,shell=True)
     else:
       try: subprocess.check_output(txt,stderr=subprocess.STDOUT,shell=True)
-      except Exception as e: print(e.output)
+      except subprocess.CalledProcessError as e: print(e.output)
 
     if not os.path.isfile("%s/libreax.a" % libdir) or \
           not os.path.isfile("%s/Makefile.lammps" % libdir):
@@ -1641,7 +1666,8 @@ class VORONOI(object):
     if not self.install: return
     libdir = dir.lib + "/voronoi"
     cmd = "cd %s; python install.py %s" % (libdir,self.install)
-    txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,shell=True).decode()
+    txt = subprocess.check_output(cmd,stderr=subprocess.STDOUT,
+                                  shell=True).decode()
     if verbose: print(txt)
     print("Created lib/voronoi library")
 

src/Makefile

@@ -92,7 +92,7 @@ help:
 	@echo 'make package-update (pu) replace src files with updated package files'
 	@echo 'make package-overwrite   replace package files with src files'
 	@echo 'make package-diff (pd)   diff src files against package files'
-	@echo 'make package-purge       purge obsolete copies of package sources'
+	@echo 'make purge               purge obsolete copies of source files'
 	@echo ''
 	@echo 'make machine             build LAMMPS for machine'
 	@echo 'make mode=lib machine    build LAMMPS as static lib for machine'
@@ -314,7 +314,7 @@ package-diff pd:
 	@echo ''
 	@for p in $(PACKUSERUC); do $(SHELL) Package.sh $$p diff; done
 
-package-purge: Purge.list
+purge: Purge.list
 	@echo 'Purging obsolete and auto-generated source files'
 	@for f in `grep -v '#' Purge.list` ;		\
 	    do test -f $$f && rm $$f && echo $$f || : ;		\

src/REPLICA/prd.cpp

@@ -783,8 +783,6 @@ void PRD::log_event()
 
 void PRD::replicate(int ireplica)
 {
-  int nreplica = universe->nworlds;
-  int nprocs_universe = universe->nprocs;
   int i,m;
 
   // -----------------------------------------------------

src/USER-REAXC/pair_reax_c.cpp

@@ -348,7 +348,7 @@ void PairReaxC::init_style( )
 
   int iqeq;
   for (iqeq = 0; iqeq < modify->nfix; iqeq++)
-    if (strcmp(modify->fix[iqeq]->style,"qeq/reax") == 0) break;
+    if (strstr(modify->fix[iqeq]->style,"qeq/reax")) break;
   if (iqeq == modify->nfix && qeqflag == 1)
     error->all(FLERR,"Pair reax/c requires use of fix qeq/reax");
 

src/comm_tiled.cpp

@@ -45,9 +45,6 @@ enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED};    // several files
 
 CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)
 {
-  if (lmp->kokkos)
-    error->all(FLERR,"KOKKOS package does not yet support comm_style tiled");
-
   style = 1;
   layout = LAYOUT_UNIFORM;
   pbc_flag = NULL;
@@ -63,9 +60,6 @@ CommTiled::CommTiled(LAMMPS *lmp) : Comm(lmp)
 
 CommTiled::CommTiled(LAMMPS *lmp, Comm *oldcomm) : Comm(*oldcomm)
 {
-  if (lmp->kokkos)
-    error->all(FLERR,"KOKKOS package does not yet support comm_style tiled");
-
   style = 1;
   layout = oldcomm->layout;
   Comm::copy_arrays(oldcomm);

src/comm_tiled.h

@@ -26,26 +26,26 @@ class CommTiled : public Comm {
 
   void init();
   void setup();                        // setup comm pattern
-  void forward_comm(int dummy = 0);    // forward comm of atom coords
-  void reverse_comm();                 // reverse comm of forces
-  void exchange();                     // move atoms to new procs
-  void borders();                      // setup list of atoms to comm
+  virtual void forward_comm(int dummy = 0);    // forward comm of atom coords
+  virtual void reverse_comm();                 // reverse comm of forces
+  virtual void exchange();                     // move atoms to new procs
+  virtual void borders();                      // setup list of atoms to comm
 
-  void forward_comm_pair(class Pair *);    // forward comm from a Pair
-  void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
+  virtual void forward_comm_pair(class Pair *);    // forward comm from a Pair
+  virtual void reverse_comm_pair(class Pair *);    // reverse comm from a Pair
   virtual void forward_comm_fix(class Fix *, int size=0);
                                                    // forward comm from a Fix
   virtual void reverse_comm_fix(class Fix *, int size=0);
                                                    // reverse comm from a Fix
   virtual void reverse_comm_fix_variable(class Fix *);
                                      // variable size reverse comm from a Fix
-  void forward_comm_compute(class Compute *);  // forward from a Compute
-  void reverse_comm_compute(class Compute *);  // reverse from a Compute
-  void forward_comm_dump(class Dump *);    // forward comm from a Dump
-  void reverse_comm_dump(class Dump *);    // reverse comm from a Dump
+  virtual void forward_comm_compute(class Compute *);  // forward from a Compute
+  virtual void reverse_comm_compute(class Compute *);  // reverse from a Compute
+  virtual void forward_comm_dump(class Dump *);    // forward comm from a Dump
+  virtual void reverse_comm_dump(class Dump *);    // reverse comm from a Dump
 
-  void forward_comm_array(int, double **);          // forward comm of array
-  int exchange_variable(int, double *, double *&);  // exchange on neigh stencil
+  virtual void forward_comm_array(int, double **);          // forward comm of array
+  virtual int exchange_variable(int, double *, double *&);  // exchange on neigh stencil
 
   void coord2proc_setup();
   int coord2proc(double *, int &, int &, int &);

src/fix_addforce.cpp

@@ -203,13 +203,10 @@ void FixAddForce::init()
       update->whichflag == 2 && estyle == NONE)
     error->all(FLERR,"Must use variable energy with fix addforce");
 
-  int max_respa = 0;
-
-  if (strstr(update->integrate_style,"respa"))
-    max_respa = ((Respa *) update->integrate)->nlevels-1;
-
-  if (respa_level >= 0)
-    ilevel_respa = MIN(respa_level,max_respa);
+  if (strstr(update->integrate_style,"respa")) {
+    ilevel_respa = ((Respa *) update->integrate)->nlevels-1;
+    if (respa_level >= 0) ilevel_respa = MIN(respa_level,ilevel_respa);
+  }
 }
 
 /* ---------------------------------------------------------------------- */

src/fix_dt_reset.cpp

@@ -48,6 +48,7 @@ FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) :
   global_freq = 1;
   extscalar = 0;
   extvector = 0;
+  dynamic_group_allow = 1;
 
   nevery = force->inumeric(FLERR,arg[3]);
   if (nevery <= 0) error->all(FLERR,"Illegal fix dt/reset command");

src/fix_group.cpp

@@ -21,6 +21,7 @@
 #include "domain.h"
 #include "region.h"
 #include "modify.h"
+#include "respa.h"
 #include "input.h"
 #include "variable.h"
 #include "memory.h"
@@ -95,6 +96,7 @@ int FixGroup::setmask()
 {
   int mask = 0;
   mask |= POST_INTEGRATE;
+  mask |= POST_INTEGRATE_RESPA;
   return mask;
 }
 
@@ -108,6 +110,9 @@ void FixGroup::init()
   if (group->dynamic[igroup])
     error->all(FLERR,"Group dynamic parent group cannot be dynamic");
 
+  if (strstr(update->integrate_style,"respa"))
+    nlevels_respa = ((Respa *) update->integrate)->nlevels;
+
   // set current indices for region and variable
 
   if (regionflag) {
@@ -169,6 +174,13 @@ void FixGroup::post_integrate()
 
 /* ---------------------------------------------------------------------- */
 
+void FixGroup::post_integrate_respa(int ilevel, int iloop)
+{
+  if (ilevel == nlevels_respa-1) post_integrate();
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixGroup::set_group()
 {
   int nlocal = atom->nlocal;

src/fix_group.h

@@ -32,6 +32,7 @@ class FixGroup : public Fix {
   void init();
   void setup(int);
   void post_integrate();
+  void post_integrate_respa(int,int);
 
  private:
   int gbit,gbitinverse;
@@ -40,6 +41,8 @@ class FixGroup : public Fix {
   char *idregion,*idvar;
   class Region *region;
 
+  int nlevels_respa;
+
   void set_group();
 };
 

src/fix_tmd.cpp

@@ -391,7 +391,7 @@ void FixTMD::readfile(char *file)
 
   char *buffer = new char[CHUNK*MAXLINE];
   char *next,*bufptr;
-  int i,m,nlines,imageflag,ix,iy,iz;
+  int i,m,n,nlines,imageflag,ix,iy,iz;
   tagint itag;
   double x,y,z,xprd,yprd,zprd;
 
@@ -400,7 +400,7 @@ void FixTMD::readfile(char *file)
   char *eof = NULL;
   xprd = yprd = zprd = -1.0;
 
-  while (!eof) {
+  do {
     if (me == 0) {
       m = 0;
       for (nlines = 0; nlines < CHUNK; nlines++) {
@@ -412,7 +412,7 @@ void FixTMD::readfile(char *file)
       m++;
     }
 
-    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    MPI_Bcast(&eof,sizeof(char *)/sizeof(char),MPI_CHAR,0,world);
     MPI_Bcast(&nlines,1,MPI_INT,0,world);
     MPI_Bcast(&m,1,MPI_INT,0,world);
     MPI_Bcast(buffer,m,MPI_CHAR,0,world);
@@ -455,10 +455,17 @@ void FixTMD::readfile(char *file)
       }
 
       if (imageflag)
-        sscanf(bufptr,TAGINT_FORMAT " %lg %lg %lg %d %d %d",
-               &itag,&x,&y,&z,&ix,&iy,&iz);
+        n = sscanf(bufptr,TAGINT_FORMAT " %lg %lg %lg %d %d %d",
+                   &itag,&x,&y,&z,&ix,&iy,&iz);
       else
-        sscanf(bufptr,TAGINT_FORMAT " %lg %lg %lg",&itag,&x,&y,&z);
+        n = sscanf(bufptr,TAGINT_FORMAT " %lg %lg %lg",&itag,&x,&y,&z);
+
+      if (n < 0 && me == 0) {
+        error->warning(FLERR,"Ignoring empty or incorrectly formatted"
+                             " line in target file");
+        bufptr = next + 1;
+        continue;
+      }
 
       m = atom->map(itag);
       if (m >= 0 && m < nlocal && mask[m] & groupbit) {
@@ -473,10 +480,9 @@ void FixTMD::readfile(char *file)
         }
         ncount++;
       }
-
       bufptr = next + 1;
     }
-  }
+  } while (eof != NULL);
 
   // clean up