git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@8701 f3b2605a-c512-4ea7-a41b-209d697bcdaa

lib/colvars/colvar.cpp

@@ -98,8 +98,8 @@ colvar::colvar (std::string const &conf)
                          "on an axis",       "distanceZ",      distance_z);
   initialize_components ("distance projection "
                          "on a plane",       "distanceXY",     distance_xy);
-  initialize_components ("average distance weighted by inverse sixth power",
+  initialize_components ("average distance weighted by inverse power",
-                         "distance6", distance6);
+                         "distanceInv", distance_inv);
 
   initialize_components ("coordination "
                          "number",           "coordNum",       coordnum);
@@ -132,7 +132,7 @@ colvar::colvar (std::string const &conf)
   initialize_components ("moment of "
                          "inertia",          "inertia",        inertia);
   initialize_components ("moment of inertia around an axis",
-                                             "inertia_z",      inertia_z);
+                                             "inertiaZ",       inertia_z);
   initialize_components ("eigenvector",      "eigenvector",    eigenvector);
 
   if (!cvcs.size())
@@ -318,13 +318,13 @@ colvar::colvar (std::string const &conf)
                             "by enabling a thermostat, or through \"extendedTemp\".\n");
       }
 
-      get_keyval (conf, "extendedFluctuation", tolerance, 0.2*width);
+      get_keyval (conf, "extendedFluctuation", tolerance, width);
       if (tolerance <= 0.0)
         cvm::fatal_error ("Error: \"extendedFluctuation\" must be positive.\n");
       ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance);
       cvm::log ("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2");
 
-      get_keyval (conf, "extendedTimeConstant", period, 40.0 * cvm::dt());
+      get_keyval (conf, "extendedTimeConstant", period, 200.0);
       if (period <= 0.0)
         cvm::fatal_error ("Error: \"extendedTimeConstant\" must be positive.\n");
       ext_mass = (cvm::boltzmann() * temp * period * period)
@@ -339,7 +339,7 @@ colvar::colvar (std::string const &conf)
         }
       }
 
-      get_keyval (conf, "extendedLangevinDamping", ext_gamma, 0.0);
+      get_keyval (conf, "extendedLangevinDamping", ext_gamma, 1.0);
       if (ext_gamma < 0.0)
         cvm::fatal_error ("Error: \"extendedLangevinDamping\" may not be negative.\n");
       if (ext_gamma != 0.0) {

lib/colvars/colvar.h

@@ -471,7 +471,7 @@ public:
   class distance;
   class distance_z;
   class distance_xy;
-  class distance6;
+  class distance_inv;
   class angle;
   class dihedral;
   class coordnum;

lib/colvars/colvaratoms.cpp

@@ -116,8 +116,9 @@ void cvm::atom_group::parse (std::string const &conf,
     std::vector<std::string>::iterator psii = psf_segids.begin();
     while (key_lookup (group_conf, "atomNameResidueRange",
                        range_conf, pos)) {
+      range_count++;
 
-      if (psii > psf_segids.end()) {
+      if (range_count > psf_segids.size()) {
         cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than "
                           "values of \"psfSegID\".\n");
       }
@@ -520,6 +521,7 @@ void cvm::atom_group::read_system_forces()
   }
 }
 
+/* This is deprecated.
 cvm::atom_pos cvm::atom_group::center_of_geometry (cvm::atom_pos const &ref_pos)
 {
   if (b_dummy) {
@@ -535,7 +537,7 @@ cvm::atom_pos cvm::atom_group::center_of_geometry (cvm::atom_pos const &ref_pos)
   }
   cog /= this->size();
   return cog;
-}
+} */
 
 cvm::atom_pos cvm::atom_group::center_of_geometry() const
 {
@@ -551,6 +553,7 @@ cvm::atom_pos cvm::atom_group::center_of_geometry() const
   return cog;
 }
 
+/* This is deprecated.
 cvm::atom_pos cvm::atom_group::center_of_mass (cvm::atom_pos const &ref_pos)
 {
   if (b_dummy) {
@@ -566,7 +569,7 @@ cvm::atom_pos cvm::atom_group::center_of_mass (cvm::atom_pos const &ref_pos)
   }
   com /= this->total_mass;
   return com;
-}
+}*/
 
 cvm::atom_pos cvm::atom_group::center_of_mass() const
 {

lib/colvars/colvaratoms.h

@@ -1,3 +1,5 @@
+// -*- c++ -*-
+
 #ifndef COLVARATOMS_H
 #define COLVARATOMS_H
 
@@ -304,15 +306,8 @@ public:
   /// the colvar has not a scalar value) or the biases require to
   /// micromanage the forces.
   void apply_forces (std::vector<cvm::rvector> const &forces);
+
 };
 
 
-
-
 #endif
-
-
-// Emacs
-// Local Variables:
-// mode: C++
-// End:

lib/colvars/colvarbias_meta.cpp

@@ -1072,8 +1072,12 @@ void colvarbias_meta::read_replica_files()
         }
 
       } else {
-        cvm::fatal_error ("Error: cannot read from file \""+
+        cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+
-                          (replicas[ir])->replica_hills_file+"\".\n");
+                  "\": will try again in "+
+                  cvm::to_str (replica_update_freq)+" steps.\n");
+        (replicas[ir])->update_status++;
+        // cvm::fatal_error ("Error: cannot read from file \""+
+        //                   (replicas[ir])->replica_hills_file+"\".\n");
       }
       is.close();
     }
@@ -1081,12 +1085,12 @@ void colvarbias_meta::read_replica_files()
     size_t const n_flush = (replica_update_freq/new_hill_freq + 1);
     if ((replicas[ir])->update_status > 3*n_flush) {
       // TODO: suspend the calculation?
-      cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+
+      cvm::log ("WARNING: in metadynamics bias \""+this->name+"\""+
-                ": failet do read completely output files from replica \""+
+                " failed to read completely the output of replica \""+
                 (replicas[ir])->replica_id+
                 "\" after more than "+
-                cvm::to_str (n_flush*new_hill_freq)+
+                cvm::to_str ((replicas[ir])->update_status * replica_update_freq)+
-                " steps.  Please check that that simulation is still running.\n");
+                " steps.  Ensure that it is still running.\n");
     }
   }
 }

lib/colvars/colvarcomp.cpp

@@ -71,21 +71,3 @@ void colvar::cvc::calc_Jacobian_derivative()
   cvm::fatal_error ("Error: calculation of inverse gradients is not implemented "
                     "for colvar components of type \""+function_type+"\".\n");
 }
-
-
-colvarvalue colvar::cvc::fdiff_change (cvm::atom_group &group)
-{  
-  colvarvalue change (x.type());
-
-  if (group.old_pos.size()) {
-    for (size_t i = 0; i < group.size(); i++) {
-      cvm::rvector const &pold = group.old_pos[i];
-      cvm::rvector const &p = group[i].pos;
-      change += group[i].grad * (p - pold);
-    }
-  }
-
-  // save for next step
-  group.old_pos = group.positions();
-  return change;
-}

lib/colvars/colvarcomp.h

@@ -111,11 +111,6 @@ public:
   /// report their differences
   bool b_debug_gradients;
 
-  /// \brief When b_debug_gradients is true, this function can be used
-  /// to calculate the estimated change in the value using the change
-  /// in the atomic coordinates times the atomic gradients
-  colvarvalue fdiff_change (cvm::atom_group &group);
-
   /// \brief If this flag is false (default), inverse gradients
   /// (derivatives of atom coordinates with respect to x) are
   /// unavailable; it should be set to true by the constructor of each
@@ -475,16 +470,16 @@ public:
 
 /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power,
 /// as in NMR refinements (colvarvalue::type_scalar type, range (0:*))
-class colvar::distance6
+class colvar::distance_inv
   : public colvar::distance
 {
 protected:
   /// Components of the distance vector orthogonal to the axis
-  cvm::real smoothing;
+  int exponent;
 public:
-  distance6 (std::string const &conf);
+  distance_inv (std::string const &conf);
-  distance6();
+  distance_inv();
-  virtual inline ~distance6() {}
+  virtual inline ~distance_inv() {}
   virtual void calc_value();
   virtual void calc_gradients();
   virtual void apply_force (colvarvalue const &force);
@@ -1301,7 +1296,7 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const
   simple_scalar_dist_functions (distance)
   // NOTE: distance_z has explicit functions, see below 
   simple_scalar_dist_functions (distance_xy)
-  simple_scalar_dist_functions (distance6)
+  simple_scalar_dist_functions (distance_inv)
   simple_scalar_dist_functions (angle)
   simple_scalar_dist_functions (coordnum)
   simple_scalar_dist_functions (selfcoordnum)

lib/colvars/colvarcomp_distances.cpp

@@ -179,7 +179,10 @@ colvar::distance_z::distance_z (std::string const &conf)
     if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
       if (axis.norm2() == 0.0)
         cvm::fatal_error ("Axis vector is zero!");
-      axis = axis.unit();
+      if (axis.norm2() != 1.0) {
+        axis = axis.unit();
+        cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
+      }
     }
     fixed_axis = true;
   }
@@ -451,25 +454,34 @@ void colvar::distance_dir::apply_force (colvarvalue const &force)
 
 
 
-colvar::distance6::distance6 (std::string const &conf)
+colvar::distance_inv::distance_inv (std::string const &conf)
   : distance (conf)
 {
-  function_type = "distance6";
+  function_type = "distance_inv";
+  get_keyval (conf, "exponent", exponent, 6);
+  if (exponent%2) {
+    cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n");
+  }
+  if (exponent <= 0) {
+    cvm::fatal_error ("Error: negative or zero exponent provided.\n");
+  }
+
   b_inverse_gradients = false;
   b_Jacobian_derivative = false;
   x.type (colvarvalue::type_scalar);
 }
 
-colvar::distance6::distance6()
+colvar::distance_inv::distance_inv()
 {
-  function_type = "distance6";
+  function_type = "distance_inv";
+  exponent = 6;
   b_inverse_gradients = false;
   b_Jacobian_derivative = false;
   b_1site_force = false;
   x.type (colvarvalue::type_scalar);
 }
 
-void colvar::distance6::calc_value()
+void colvar::distance_inv::calc_value()
 {
   group1.reset_atoms_data();
   group2.reset_atoms_data();
@@ -483,8 +495,11 @@ void colvar::distance6::calc_value()
       for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
         cvm::rvector const dv = ai2->pos - ai1->pos;
         cvm::real const d2 = dv.norm2();
-        x.real_value += 1.0/(d2*d2*d2);
+        cvm::real dinv = 1.0;
-        cvm::rvector const dsumddv = -6.0/(d2*d2*d2*d2) * dv;
+        for (int ne = 0; ne < exponent/2; ne++)
+          dinv *= 1.0/d2;
+        x.real_value += dinv;
+        cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
         ai1->grad += -1.0 * dsumddv;
         ai2->grad +=        dsumddv;
       }
@@ -494,30 +509,39 @@ void colvar::distance6::calc_value()
       for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
         cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
         cvm::real const d2 = dv.norm2();
-        x.real_value += 1.0/(d2*d2*d2);
+        cvm::real dinv = 1.0;
-        cvm::rvector const dsumddv = -6.0/(d2*d2*d2*d2) * dv;
+        for (int ne = 0; ne < exponent/2; ne++)
+          dinv *= 1.0/d2;
+        x.real_value += dinv;
+        cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
         ai1->grad += -1.0 * dsumddv;
         ai2->grad +=        dsumddv;
       }
     }
   }
 
-  x.real_value = std::pow (x.real_value, -1.0/6.0);
+  x.real_value *= 1.0 / cvm::real (group1.size() * group2.size());
+  x.real_value = std::pow (x.real_value, -1.0/(cvm::real (exponent)));
 }
 
-void colvar::distance6::calc_gradients()
+void colvar::distance_inv::calc_gradients()
 {
+  cvm::real const dxdsum = (-1.0/(cvm::real (exponent))) * std::pow (x.real_value, exponent+1) / cvm::real (group1.size() * group2.size());
+  for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
+    ai1->grad *= dxdsum;
+  }
+  for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
+    ai2->grad *= dxdsum;
+  }
 }
 
-void colvar::distance6::apply_force (colvarvalue const &force)
+void colvar::distance_inv::apply_force (colvarvalue const &force)
 {
-  cvm::real const dxdsum = (-1.0/6.0) * std::pow (x.real_value, -7.0/6.0);
-
   if (!group1.noforce)
-    group1.apply_colvar_force (dxdsum * force.real_value);
+    group1.apply_colvar_force (force.real_value);
 
   if (!group2.noforce)
-    group2.apply_colvar_force (dxdsum * force.real_value);
+    group2.apply_colvar_force (force.real_value);
 }
 
 
@@ -614,21 +638,19 @@ void colvar::inertia::calc_value()
 {
   atoms.reset_atoms_data();
   atoms.read_positions();
-  atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
+  atoms.apply_translation ((-1.0) * atoms.center_of_mass());
 
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     x.real_value += ai->mass * (ai->pos).norm2();
   }
-  x.real_value = std::sqrt (x.real_value / atoms.total_mass);
 }
 
 
 void colvar::inertia::calc_gradients()
 {
-  cvm::real const drdx = 1.0/(atoms.total_mass * x.real_value);
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
-    ai->grad = drdx * ai->mass * ai->pos;
+    ai->grad = 2.0 * ai->mass * ai->pos;
   }
 }
 
@@ -647,7 +669,10 @@ colvar::inertia_z::inertia_z (std::string const &conf)
   if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
     if (axis.norm2() == 0.0)
       cvm::fatal_error ("Axis vector is zero!");
-    axis = axis.unit();
+    if (axis.norm2() != 1.0) {
+      axis = axis.unit();
+      cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
+    }
   }
   x.type (colvarvalue::type_scalar);
 }
@@ -664,22 +689,20 @@ void colvar::inertia_z::calc_value()
 {
   atoms.reset_atoms_data();
   atoms.read_positions();
-  atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
+  atoms.apply_translation ((-1.0) * atoms.center_of_mass());
 
   x.real_value = 0.0;
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
     cvm::real const iprod = ai->pos * axis;
     x.real_value += ai->mass * iprod * iprod;
   }
-  x.real_value = std::sqrt (x.real_value / atoms.total_mass);
 }
 
 
 void colvar::inertia_z::calc_gradients()
 {
-  cvm::real const drdx = 1.0/(atoms.total_mass * x.real_value);
   for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
-    ai->grad = drdx * ai->mass * (ai->pos * axis) * axis;
+    ai->grad = 2.0 * ai->mass * (ai->pos * axis) * axis;
   }
 }
 

lib/colvars/colvarcomp_rotations.cpp

@@ -72,7 +72,6 @@ colvar::orientation::orientation (std::string const &conf)
     rot.request_group2_gradients (atoms.size());
   }
 
-  rot.b_debug_gradients = this->b_debug_gradients;
 }
 
   
@@ -228,38 +227,6 @@ void colvar::tilt::calc_gradients()
       atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
     }
   }
-
-  if (cvm::debug()) {
-
-    std::vector<cvm::atom_pos> pos_test (atoms.positions_shifted (-1.0 * atoms_cog));
-
-    for (size_t comp = 0; comp < 3; comp++) {
-      
-      // correct this cartesian component for the "new" center of geometry
-      for (size_t ia = 0; ia < atoms.size(); ia++) {
-        pos_test[ia][comp] -=
-          cvm::debug_gradients_step_size / cvm::real (atoms.size());
-      }
-
-      for (size_t ia = 0; ia < atoms.size(); ia++) {
-
-        pos_test[ia][comp] += cvm::debug_gradients_step_size;
-        rot.calc_optimal_rotation (ref_pos, pos_test);
-        pos_test[ia][comp] -= cvm::debug_gradients_step_size;
-
-        cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+
-                  cvm::to_str (std::fabs (rot.cos_theta (axis) - x.real_value -
-                                       cvm::debug_gradients_step_size * atoms[ia].grad[comp]) /
-                               std::fabs (rot.cos_theta (axis) - x.real_value),
-                               cvm::cv_width, cvm::cv_prec)+".\n");
-      }
-
-      for (size_t ia = 0; ia < atoms.size(); ia++) {
-        pos_test[ia][comp] +=
-          cvm::debug_gradients_step_size / cvm::real (atoms.size());
-      }
-    }
-  }
 }
 
 
@@ -326,38 +293,6 @@ void colvar::spin_angle::calc_gradients()
       atoms[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
     }
   }
-
-  if (cvm::debug()) {
-
-    std::vector<cvm::atom_pos> pos_test (atoms.positions_shifted (-1.0 * atoms_cog));
-
-    for (size_t comp = 0; comp < 3; comp++) {
-      
-      // correct this cartesian component for the "new" center of geometry
-      for (size_t ia = 0; ia < atoms.size(); ia++) {
-        pos_test[ia][comp] -=
-          cvm::debug_gradients_step_size / cvm::real (atoms.size());
-      }
-
-      for (size_t ia = 0; ia < atoms.size(); ia++) {
-
-        pos_test[ia][comp] += cvm::debug_gradients_step_size;
-        rot.calc_optimal_rotation (ref_pos, pos_test);
-        pos_test[ia][comp] -= cvm::debug_gradients_step_size;
-
-        cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+
-                  cvm::to_str (std::fabs (rot.spin_angle (axis) - x.real_value -
-                                       cvm::debug_gradients_step_size * atoms[ia].grad[comp]) /
-                               std::fabs (rot.spin_angle (axis) - x.real_value),
-                               cvm::cv_width, cvm::cv_prec)+".\n");
-      }
-
-      for (size_t ia = 0; ia < atoms.size(); ia++) {
-        pos_test[ia][comp] +=
-          cvm::debug_gradients_step_size / cvm::real (atoms.size());
-      }
-    }
-  }
 }
 
 

lib/colvars/colvarmodule.cpp

@@ -46,9 +46,8 @@ colvarmodule::colvarmodule (char const  *config_filename,
 
   parse->get_keyval (conf, "analysis", b_analysis, false);
 
-  if (cvm::debug())
+  parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03,
-    parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03,
+                     colvarparse::parse_silent);
-                       colvarparse::parse_silent);
 
   parse->get_keyval (conf, "eigenvalueCrossingThreshold",
                      colvarmodule::rotation::crossing_threshold, 1.0e-04,

lib/colvars/colvarmodule.h

@@ -2,7 +2,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2012-06-20"
+#define COLVARS_VERSION "2012-08-08"
 #endif
 
 #ifndef COLVARS_DEBUG

lib/colvars/colvartypes.h

@@ -810,21 +810,8 @@ public:
   /// associated to this quaternion and another
   inline cvm::real cosine (cvm::quaternion const &q) const
   {
-    if (q0*q0 >= (1.0-1.0E-10)) {
+    cvm::real const iprod = this->inner (q);
-      // null quaternion, return the square of the cosine of the other
+    return 2.0*iprod*iprod - 1.0;
-      // quaternion's rotation angle
-      return (2.0*q.q0*q.q0 - 1.0);
-    } else if (q.q0*q.q0 >= (1.0-1.0E-10)) {
-      return (2.0*q0*q0 - 1.0);
-    } else {
-      // get a vector orthogonal to both axes of rotation
-      cvm::rvector const op = (cvm::rvector::outer (this->get_vector(), q.get_vector()));
-      cvm::real const  opl2 = op.norm2();
-      // rotate it with both quaternions and get the normalized inner product
-      return ( (opl2 > 0.0) ?
-               ((this->rotate (op)) * (q.rotate (op))) / opl2 :
-               1.0 );
-    }
   }
 
   /// \brief Square distance from another quaternion on the