First step towards MTS support in lambda-dynamics

namd/src/colvarproxy_namd.C

@@ -1648,15 +1648,15 @@ int colvarproxy_namd::replica_comm_send(char* msg_data, int msg_len,
 }
 
 
-/// Request energy computation every freq steps
+/// Request alchemical energy computation every freq steps
 int colvarproxy_namd::request_alch_energy_freq(int const freq) {
   // This test is only valid for NAMD3
   if (freq % simparams->computeEnergies) {
     cvm::error("computeEnergies must be a divisor of lambda-dynamics period (" + cvm::to_str(freq) + ").\n");
     return COLVARS_INPUT_ERROR;
   }
   if (!simparams->alchOn) {
-    cvm::error("alch must be enabled for lambda-dynamics.\n");
+    cvm::error("alchOn must be enabled for lambda-dynamics.\n");
     return COLVARS_INPUT_ERROR;
   }
   if (!simparams->alchThermIntOn) {

namd/tests/library/Common/da.pdb

@@ -1,15 +1,15 @@
 CRYST1    0.000    0.000    0.000  90.00  90.00  90.00 P 1           1
-ATOM      1  N   ALA B   1       6.652  -0.764  -1.333  1.00  0.00      BH    
-ATOM      2  HT2 ALA B   1       5.760  -0.450  -1.766  1.00  0.00      BH    
-ATOM      3  HT3 ALA B   1       7.342  -1.305  -1.892  1.00  0.00      BH    
-ATOM      4  CA  ALA B   1       7.306   0.408  -0.585  1.00  0.00      BH    
-ATOM      5  HA  ALA B   1       8.253  -0.035  -0.313  1.00  0.00      BH    
-ATOM      6  CB  ALA B   1       7.329   1.530  -1.594  1.00  0.00      BH    
-ATOM      7  HB1 ALA B   1       7.786   2.446  -1.164  1.00  0.00      BH    
-ATOM      8  HB2 ALA B   1       7.862   1.113  -2.475  1.00  0.00      BH    
-ATOM      9  HB3 ALA B   1       6.363   1.945  -1.952  1.00  0.00      BH    
-ATOM     10  C   ALA B   1       6.475   0.890   0.649  1.00  0.00      BH    
-ATOM     11  O   ALA B   1       5.256   1.022   0.506  1.00  0.00      BH    
+ATOM      1  N   ALA B   1       6.652  -0.764  -1.333  1.00 -1.00      BH    
+ATOM      2  HT2 ALA B   1       5.760  -0.450  -1.766  1.00 -1.00      BH    
+ATOM      3  HT3 ALA B   1       7.342  -1.305  -1.892  1.00 -1.00      BH    
+ATOM      4  CA  ALA B   1       7.306   0.408  -0.585  1.00 -1.00      BH    
+ATOM      5  HA  ALA B   1       8.253  -0.035  -0.313  1.00 -1.00      BH    
+ATOM      6  CB  ALA B   1       7.329   1.530  -1.594  1.00 -1.00      BH    
+ATOM      7  HB1 ALA B   1       7.786   2.446  -1.164  1.00 -1.00      BH    
+ATOM      8  HB2 ALA B   1       7.862   1.113  -2.475  1.00 -1.00      BH    
+ATOM      9  HB3 ALA B   1       6.363   1.945  -1.952  1.00 -1.00      BH    
+ATOM     10  C   ALA B   1       6.475   0.890   0.649  1.00 -1.00      BH    
+ATOM     11  O   ALA B   1       5.256   1.022   0.506  1.00 -1.00      BH    
 ATOM     12  N   ALA B   2       7.106   1.140   1.813  1.00  0.00      BH    
 ATOM     13  HN  ALA B   2       8.097   1.077   1.900  1.00  0.00      BH    
 ATOM     14  CA  ALA B   2       6.403   1.216   3.045  1.00  0.00      BH    

src/colvar.cpp

@@ -305,6 +305,8 @@ int colvar::init(std::string const &conf)
   // Detect if we have a single component that is an alchemical lambda
   if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type() == "alchLambda") {
     enable(f_cv_external);
+
+    static_cast<colvar::alch_lambda *>(cvcs[0].get())->init_alchemy(time_step_factor);
   }
 
   // If using scripted biases, any colvar may receive bias forces
@@ -313,6 +315,17 @@ int colvar::init(std::string const &conf)
   }
 
   error_code |= init_extended_Lagrangian(conf);
+
+  // when total atomic forces are obtained from the previous time step,
+  // we cannot (currently) have colvar values and projected total forces for the same timestep
+  // (that would require anticipating the total force request by one timestep)
+  // i.e. the combination of f_cv_total_force_calc and f_cv_multiple_ts requires f_cv_total_force_current_step
+  // Because f_cv_total_force_current_step is static, we can hard-code this, once other features are set
+  // that is f_cv_external and f_cv_extended_Lagrangian
+  if (!is_enabled(f_cv_total_force_current_step)) {
+    exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
+  }
+
   error_code |= init_output_flags(conf);
 
   // Now that the children are defined we can solve dependencies
@@ -1119,6 +1132,8 @@ int colvar::init_dependencies() {
     init_feature(f_cv_total_force, "total_force", f_type_dynamic);
     require_feature_alt(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
 
+    // If this is active, the total force reported to biases (ABF / TI) is from the current step
+    // therefore it does not include Colvars biases -> it is a "system force"
     init_feature(f_cv_total_force_current_step, "total_force_current_step", f_type_dynamic);
 
     // Deps for explicit total force calculation
@@ -1210,15 +1225,6 @@ int colvar::init_dependencies() {
 
     init_feature(f_cv_multiple_ts, "multiple_timestep", f_type_static);
 
-    // when total atomic forces are obtained from the previous time step,
-    // we cannot (currently) have colvar values and projected total forces for the same timestep
-    // TODO this will need refining for driven alchemical parameters
-    // ie. the combination of f_cv_total_force_calc and f_cv_multiple_ts requires f_cv_total_force_current_step
-    // or we need to anticipate the total force request by one timestep
-    if (!cvm::main()->proxy->total_forces_same_step()) {
-      exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
-    }
-
     // check that everything is initialized
     for (i = 0; i < colvardeps::f_cv_ntot; i++) {
       if (is_not_set(i)) {
@@ -1624,29 +1630,31 @@ int colvar::collect_cvc_total_forces()
   if (is_enabled(f_cv_total_force_calc)) {
     ft.reset();
 
-    if (cvm::step_relative() > 0) {
-      // get from the cvcs the total forces from the PREVIOUS step
-      for (size_t i = 0; i < cvcs.size();  i++) {
-        if (!cvcs[i]->is_enabled()) continue;
-            if (cvm::debug())
-            cvm::log("Colvar component no. "+cvm::to_str(i+1)+
-                " within colvar \""+this->name+"\" has total force "+
-                cvm::to_str((cvcs[i])->total_force(),
-                cvm::cv_width, cvm::cv_prec)+".\n");
-        // linear combination is assumed
-        ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
-      }
+    for (size_t i = 0; i < cvcs.size();  i++) {
+      if (!cvcs[i]->is_enabled()) continue;
+          if (cvm::debug())
+          cvm::log("Colvar component no. "+cvm::to_str(i+1)+
+              " within colvar \""+this->name+"\" has total force "+
+              cvm::to_str((cvcs[i])->total_force(),
+              cvm::cv_width, cvm::cv_prec)+".\n");
+      // linear combination is assumed
+      ft += (cvcs[i])->total_force() * (cvcs[i])->sup_coeff / active_cvc_square_norm;
     }
 
     if (!(is_enabled(f_cv_hide_Jacobian) && is_enabled(f_cv_subtract_applied_force))) {
-      // add the Jacobian force to the total force, and don't apply any silent
+      // This is by far the most common case
+      // Add the Jacobian force to the total force, and don't apply any silent
       // correction internally: biases such as colvarbias_abf will handle it
       // If f_cv_hide_Jacobian is enabled, a force of -fj is present in ft due to the
       // Jacobian-compensating force
       ft += fj;
     }
   }
 
+  if (is_enabled(f_cv_total_force_current_step)) {
+   // Report total force value without waiting for calc_colvar_properties()
+    ft_reported = ft;
+  }
   return COLVARS_OK;
 }
 
@@ -1748,9 +1756,12 @@ int colvar::calc_colvar_properties()
     // But we report values at the beginning of the timestep (value at t=0 on the first timestep)
     x_reported = x_ext;
     v_reported = v_ext;
-    // the "total force" with the extended Lagrangian is
-    // calculated in update_forces_energy() below
 
+    // the "total force" for the extended Lagrangian is calculated in update_forces_energy() below
+    // A future improvement could compute a "system force" here, borrowing a part of update_extended_Lagrangian()
+    // this would change the behavior of eABF with respect to other biases
+    // by enabling f_cv_total_force_current_step, and reducing the total force to a system force
+    // giving the behavior of f_cv_subtract_applied_force - this is correct for WTM-eABF etc.
   } else {
 
     if (is_enabled(f_cv_subtract_applied_force) && !cvm::proxy->total_forces_same_step()) {
@@ -1879,14 +1890,18 @@ void colvar::update_extended_Lagrangian()
   }
   f_ext += f_system;
 
-  if (is_enabled(f_cv_subtract_applied_force)) {
-    // Report a "system" force without the biases on this colvar
-    // that is, just the spring force (or alchemical force)
-    ft_reported = f_system;
-  } else {
-    // The total force acting on the extended variable is f_ext
-    // This will be used in the next timestep
-    ft_reported = f_ext;
+  if ( ! is_enabled(f_cv_total_force_current_step)) {
+    if (is_enabled(f_cv_subtract_applied_force)) {
+      // Report a "system" force without the biases on this colvar
+      // that is, just the spring force (or alchemical force)
+      ft_reported = f_system;
+    } else {
+      // The total force acting on the extended variable is f_ext
+      // This will be used in the next timestep
+      ft_reported = f_ext;
+    }
+    // Since biases have already been updated, this ft_reported will only be
+    // communicated to biases at the next timestep
   }
 
   // backup in case we need to revert this integration timestep

src/colvarbias_abf.cpp

@@ -340,11 +340,12 @@ int colvarbias_abf::update()
   if (can_accumulate_data() && is_enabled(f_cvb_history_dependent)) {
 
     if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) {
+      // Note: this will skip step 0 data when available in some cases (extended system),
+      // but not doing so would make the code more complex
       if (samples->index_ok(force_bin)) {
         // Only if requested and within bounds of the grid...
 
-        // get total forces (lagging by 1 timestep) from colvars
-        // and subtract previous ABF force if necessary
+        // get total force and subtract previous ABF force if necessary
         update_system_force();
 
         gradients->acc_force(force_bin, system_force);
@@ -451,14 +452,13 @@ int colvarbias_abf::update_system_force()
   // System force from atomic forces (or extended Lagrangian if applicable)
 
   for (i = 0; i < num_variables(); i++) {
-    if (colvars[i]->is_enabled(f_cv_subtract_applied_force) ||
-       (cvm::proxy->total_forces_same_step() && !colvars[i]->is_enabled(f_cv_external))) {
+    if (colvars[i]->is_enabled(f_cv_subtract_applied_force)
+      || colvars[i]->is_enabled(f_cv_total_force_current_step)) {
       // this colvar is already subtracting the ABF force
-      // or the "total force" is really a system force at current step
-      // (For external parameters, the total force contains biasing forces
-      // unless f_cv_subtract_applied_force is enabled)
+      // or the "total force" is from current step and cannot possibly contain Colvars biases
       system_force[i] = colvars[i]->total_force().real_value;
     } else {
+      // Subtract previous step's bias force from previous step's total force
       system_force[i] = colvars[i]->total_force().real_value
         - colvar_forces[i].real_value;
     }

src/colvarcomp.h

@@ -1212,6 +1212,7 @@ class colvar::alch_lambda
   // No atom groups needed
 public:
   alch_lambda();
+  int init_alchemy(int time_step_factor);
   virtual ~alch_lambda() {}
   virtual void calc_value();
   virtual void calc_force_invgrads();

src/colvarcomp_alchlambda.cpp

@@ -29,15 +29,25 @@ colvar::alch_lambda::alch_lambda()
 
   x.type(colvarvalue::type_scalar);
 
+  // Query initial value from back-end; will be overwritten if restarting from a state file
+  cvm::proxy->get_alch_lambda(&x.real_value);
+}
+
+
+int colvar::alch_lambda::init_alchemy(int factor)
+{
   // We need calculation every time step
   // default in Tinker-HP and NAMD2, must be enforced in NAMD3
   // Also checks back-end settings, ie. that alchemy is enabled
   // (in NAMD3: alchType TI, computeEnergies at the right frequency)
-  cvm::proxy->request_alch_energy_freq(1);
-  // TODO examine how this breaks everything - whereas alchOutFreq seems to work
 
-  // Query initial value from back-end; will be overwritten if restarting from a state file
-  cvm::proxy->get_alch_lambda(&x.real_value);
+  // Forbid MTS until fully implemented
+  if (factor != 1) {
+    return cvm::error("Error: timeStepFactor > 1 is not yet supported for alchemical variables.");
+  }
+  cvm::proxy->request_alch_energy_freq(factor);
+
+  return COLVARS_OK;
 }