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

lib/colvars/colvar.cpp

@@ -203,7 +203,7 @@ colvar::colvar(std::string const &conf)
       cvm::error("Could not parse scripted colvar type.");
       return;
     }
-    x_reported.type(x.type());
+
     cvm::log(std::string("Expecting colvar value of type ")
       + colvarvalue::type_desc(x.type()));
 
@@ -216,6 +216,8 @@ colvar::colvar(std::string const &conf)
       x.vector1d_value.resize(size);
     }
 
+    x_reported.type(x);
+
     // Sort array of cvcs based on their names
     // Note: default CVC names are in input order for same type of CVC
     std::sort(cvcs.begin(), cvcs.end(), compare);
@@ -335,6 +337,7 @@ colvar::colvar(std::string const &conf)
 
   // at this point, the colvar's type is defined
   f.type(value());
+  f_accumulated.type(value());
   fb.type(value());
 
   get_keyval(conf, "width", width, 1.0);
@@ -419,6 +422,8 @@ colvar::colvar(std::string const &conf)
                INPUT_ERROR);
   }
 
+  get_keyval(conf, "timeStepFactor", time_step_factor, 1);
+
   {
     bool b_extended_Lagrangian;
     get_keyval(conf, "extendedLagrangian", b_extended_Lagrangian, false);
@@ -432,6 +437,7 @@ colvar::colvar(std::string const &conf)
       // Make feature available only on user request
       provide(f_cv_extended_Lagrangian);
       enable(f_cv_extended_Lagrangian);
+      provide(f_cv_Langevin);
 
       xr.type(value());
       vr.type(value());
@@ -1187,6 +1193,7 @@ cvm::real colvar::update_forces_energy()
     }
   }
 
+  f_accumulated += f;
 
   if (is_enabled(f_cv_fdiff_velocity)) {
     // set it for the next step
@@ -1229,14 +1236,14 @@ void colvar::communicate_forces()
       if (!cvcs[i]->is_enabled()) continue;
       // cvc force is colvar force times colvar/cvc Jacobian
       // (vector-matrix product)
-      (cvcs[i])->apply_force(colvarvalue(f.as_vector() * func_grads[grad_index++],
+      (cvcs[i])->apply_force(colvarvalue(f_accumulated.as_vector() * func_grads[grad_index++],
                              cvcs[i]->value().type()));
     }
   } else if (x.type() == colvarvalue::type_scalar) {
 
     for (i = 0; i < cvcs.size(); i++) {
       if (!cvcs[i]->is_enabled()) continue;
-      (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff *
+      (cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff *
                               cvm::real((cvcs[i])->sup_np) *
                               (std::pow((cvcs[i])->value().real_value,
                                       (cvcs[i])->sup_np-1)) );
@@ -1246,10 +1253,14 @@ void colvar::communicate_forces()
 
     for (i = 0; i < cvcs.size(); i++) {
       if (!cvcs[i]->is_enabled()) continue;
-      (cvcs[i])->apply_force(f * (cvcs[i])->sup_coeff);
+      (cvcs[i])->apply_force(f_accumulated * (cvcs[i])->sup_coeff);
     }
   }
 
+  // Accumulated forces have been applied, impulse-style
+  // Reset to start accumulating again
+  f_accumulated.reset();
+
   if (cvm::debug())
     cvm::log("Done communicating forces from colvar \""+this->name+"\".\n");
 }

lib/colvars/colvar.h

@@ -296,10 +296,24 @@ protected:
   /// Sum of square coefficients for active cvcs
   cvm::real active_cvc_square_norm;
 
+  /// Time step multiplier (for coarse-time-step colvars)
+  /// Colvar will only be calculated at those times; biases may ignore the information and
+  /// always update their own forces (which is typically inexpensive) especially if
+  /// they rely on other colvars. In this case, the colvar will accumulate forces applied between
+  /// colvar updates. Alternately they may use it to calculate "impulse" biasing
+  /// forces at longer intervals. Impulse forces must be multiplied by the timestep factor.
+  int   time_step_factor;
+
+  /// Biasing force collected between updates, to be applied at next update for coarse-time-step colvars
+  colvarvalue f_accumulated;
+
 public:
   /// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
   inline size_t num_active_cvcs() const { return n_active_cvcs; }
 
+  /// \brief returns time_step_factor
+  inline int get_time_step_factor() const {return time_step_factor;}
+
   /// \brief Use the internal metrics (as from \link cvc
   /// \endlink objects) to calculate square distances and gradients
   ///
@@ -352,7 +366,6 @@ public:
 
 
 protected:
-
   /// Previous value (to calculate velocities during analysis)
   colvarvalue            x_old;
 

lib/colvars/colvardeps.cpp

@@ -263,6 +263,7 @@ void cvm::deps::init_cv_requires() {
     f_description(f_cv_extended_Lagrangian, "extended Lagrangian");
 
     f_description(f_cv_Langevin, "Langevin dynamics");
+    f_req_self(f_cv_Langevin, f_cv_extended_Lagrangian);
 
     f_description(f_cv_linear, "linear");
 

lib/colvars/colvarmodule.cpp

@@ -509,6 +509,11 @@ int colvarmodule::calc_colvars()
   int error_code = COLVARS_OK;
   std::vector<colvar *>::iterator cvi;
 
+  // Determine which colvars are active at this time step
+  for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+    (*cvi)->feature_states[cvm::deps::f_cv_active]->enabled = (step_absolute() % (*cvi)->get_time_step_factor() == 0);
+  }
+
   // if SMP support is available, split up the work
   if (proxy->smp_enabled() == COLVARS_OK) {
 
@@ -525,6 +530,7 @@ int colvarmodule::calc_colvars()
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
 
+      if (!(*cvi)->is_enabled()) continue;
       combine_errors(error_code, (*cvi)->update_cvc_flags());
 
       size_t num_items = (*cvi)->num_active_cvcs();
@@ -544,6 +550,7 @@ int colvarmodule::calc_colvars()
 
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+      if (!(*cvi)->is_enabled()) continue;
       combine_errors(error_code, (*cvi)->collect_cvc_data());
     }
     cvm::decrease_depth();
@@ -553,6 +560,7 @@ int colvarmodule::calc_colvars()
     // calculate colvars one at a time
     cvm::increase_depth();
     for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+      if (!(*cvi)->is_enabled()) continue;
       combine_errors(error_code, (*cvi)->calc());
       if (cvm::get_error()) {
         return COLVARS_ERROR;
@@ -645,6 +653,8 @@ int colvarmodule::update_colvar_forces()
              "of colvars (if they have any).\n");
   cvm::increase_depth();
   for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
+    // Here we call even inactive colvars, so they accumulate biasing forces
+    // as well as update their extended-system dynamics
     total_colvar_energy += (*cvi)->update_forces_energy();
     if (cvm::get_error()) {
       return COLVARS_ERROR;
@@ -660,6 +670,7 @@ int colvarmodule::update_colvar_forces()
   cvm::increase_depth();
   for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
     if ((*cvi)->is_enabled(cvm::deps::f_cv_gradient)) {
+      if (!(*cvi)->is_enabled()) continue;
       (*cvi)->communicate_forces();
       if (cvm::get_error()) {
         return COLVARS_ERROR;

lib/colvars/colvarmodule.h

@@ -4,7 +4,7 @@
 #define COLVARMODULE_H
 
 #ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2016-05-10"
+#define COLVARS_VERSION "2016-05-23"
 #endif
 
 #ifndef COLVARS_DEBUG