Lambda-ABF implementation for the NAMD interface

colvartools/abmd.tcl

@@ -53,7 +53,7 @@ proc calc_z_max_gradient { args } { return 0 }
 proc calc_z_min_gradient { args } { return 0 }
 
 proc setup_ABMD { colvar force_k z_stop {direction up} } {
-  # cv config "scriptedColvarForces on"
+  cv config "scriptedColvarForces on"
 
   namespace eval ::ABMD {}
   set ::ABMD::cvname $colvar

namd/src/colvarproxy_namd.C

@@ -345,6 +345,8 @@ void colvarproxy_namd::calculate()
     colvars->update_engine_parameters();
     colvars->setup_input();
     colvars->setup_output();
+    // Controller is only available after full startup phase, so now
+    controller = &Node::Object()->state->getController();
 
     first_timestep = false;
 
@@ -373,7 +375,7 @@ void colvarproxy_namd::calculate()
   }
 
   previous_NAMD_step = step;
-  update_accelMD_info();
+  if (accelMDOn) update_accelMD_info();
 
   {
     Vector const a = lattice->a();
@@ -605,14 +607,8 @@ void colvarproxy_namd::calculate()
 }
 
 void colvarproxy_namd::update_accelMD_info() {
-  if (accelMDOn == false) {
-    return;
-  }
-  const Controller& c = Node::Object()->state->getController();
   // This aMD factor is from previous step!
-  amd_weight_factor = std::exp(c.accelMDdV /
-                               (target_temperature() * boltzmann()));
-//   std::cout << "Step: " << cvm::to_str(colvars->it) << " accelMD dV in colvars: " << c.accelMDdV << std::endl;
+  amd_weight_factor = std::exp(controller->accelMDdV / (target_temperature() * boltzmann()));
 }
 
 
@@ -1651,3 +1647,48 @@ int colvarproxy_namd::replica_comm_send(char* msg_data, int msg_len,
   replica_send(msg_data, msg_len, dest_rep, CkMyPe());
   return msg_len;
 }
+
+
+/// Request 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");
+    return COLVARS_INPUT_ERROR;
+  }
+  if (!simparams->alchThermIntOn) {
+    cvm::error("alchType must be set to TI for lambda-dynamics.\n");
+    return COLVARS_INPUT_ERROR;
+  }
+  return COLVARS_OK;
+}
+
+
+/// Get value of alchemical lambda parameter from back-end
+int colvarproxy_namd::get_alch_lambda(cvm::real* lambda) {
+  *lambda = simparams->alchLambda;
+  return COLVARS_OK;
+}
+
+
+/// Set value of alchemical lambda parameter in back-end
+int colvarproxy_namd::send_alch_lambda(void) {
+  simparams->alchLambda = cached_alch_lambda;
+  return COLVARS_OK;
+}
+
+
+/// Get energy derivative with respect to lambda
+int colvarproxy_namd::get_dE_dlambda(cvm::real* dE_dlambda) {
+  // Force data at step zero is garbage in NAMD3, zero in NAMD2
+  if (cvm::step_relative() > 0) {
+    *dE_dlambda = controller->getTIderivative();
+  } else {
+    *dE_dlambda = 0.0;
+  }
+  return COLVARS_OK;
+}

namd/src/colvarproxy_namd.h

@@ -46,6 +46,9 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster {
   /// Pointer to the NAMD simulation input object
   SimParameters *simparams;
 
+  /// Pointer to Controller object
+  Controller const *controller;
+
   /// NAMD-style PRNG object
   Random random;
 
@@ -266,6 +269,18 @@ class colvarproxy_namd : public colvarproxy, public GlobalMaster {
 
   int backup_file(char const *filename) override;
 
+  /// Get value of alchemical lambda parameter from back-end
+  int get_alch_lambda(cvm::real* lambda);
+
+  /// Set value of alchemical lambda parameter in back-end
+  int send_alch_lambda(void);
+
+  /// Request energy computation every freq steps
+  int request_alch_energy_freq(int const freq);
+
+  /// Get energy derivative with respect to lambda
+  int get_dE_dlambda(cvm::real* dE_dlambda);
+
 };
 
 

src/colvar.cpp

@@ -182,12 +182,6 @@ int colvar::init(std::string const &conf)
 
   set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
 
-  // If using scripted biases, any colvar may receive bias forces
-  // and will need its gradient
-  if (cvm::scripted_forces()) {
-    enable(f_cv_gradient);
-  }
-
   // check for linear combinations
   {
     bool lin = !(is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function));
@@ -320,6 +314,11 @@ int colvar::init(std::string const &conf)
     enable(f_cv_external);
   }
 
+  // If using scripted biases, any colvar may receive bias forces
+  if (cvm::scripted_forces()) {
+    enable(f_cv_apply_force);
+  }
+
   error_code |= init_extended_Lagrangian(conf);
   error_code |= init_output_flags(conf);
 
@@ -693,14 +692,15 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
     x_ext.type(colvarvalue::type_notset);
     v_ext.type(value());
     fr.type(value());
-    const bool temp_provided = get_keyval(conf, "extendedTemp", temp,
-                                          proxy->target_temperature());
+    const bool temp_provided = get_keyval(conf, "extendedTemp", temp, proxy->target_temperature());
     if (is_enabled(f_cv_external)) {
-      // In the case of an "external" coordinate, there is no coupling potential:
+      // In the case of a driven external parameter in the back-end, there is no coupling potential:
       // only the fictitious mass is meaningful
       get_keyval(conf, "extendedMass", ext_mass);
       // Ensure that the computed restraint energy term is zero
       ext_force_k = 0.0;
+      // Then we need forces from the back-end
+      enable(f_cv_total_force_calc);
     } else {
       // Standard case of coupling to a geometric colvar
       if (temp <= 0.0) { // Then a finite temperature is required
@@ -1135,6 +1135,9 @@ int colvar::init_dependencies() {
     init_feature(f_cv_gradient, "gradient", f_type_dynamic);
     require_feature_children(f_cv_gradient, f_cvc_gradient);
 
+    init_feature(f_cv_apply_force, "apply_force", f_type_dynamic);
+    require_feature_alt(f_cv_apply_force, f_cv_gradient, f_cv_external);
+
     init_feature(f_cv_collect_gradient, "collect_gradient", f_type_dynamic);
     require_feature_self(f_cv_collect_gradient, f_cv_gradient);
     require_feature_self(f_cv_collect_gradient, f_cv_scalar);
@@ -1153,6 +1156,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);
 
+    init_feature(f_cv_total_force_current_step, "total_force_current_step", f_type_dynamic);
+
     // Deps for explicit total force calculation
     init_feature(f_cv_total_force_calc, "total_force_calculation", f_type_dynamic);
     require_feature_self(f_cv_total_force_calc, f_cv_scalar);
@@ -1171,13 +1176,15 @@ int colvar::init_dependencies() {
 
     init_feature(f_cv_extended_Lagrangian, "extended_Lagrangian", f_type_user);
     require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar);
-    require_feature_self(f_cv_extended_Lagrangian, f_cv_gradient);
+    require_feature_self(f_cv_extended_Lagrangian, f_cv_apply_force);
 
     init_feature(f_cv_Langevin, "Langevin_dynamics", f_type_user);
     require_feature_self(f_cv_Langevin, f_cv_extended_Lagrangian);
 
-    init_feature(f_cv_external, "external", f_type_user);
+    init_feature(f_cv_external, "external_parameter", f_type_static);
     require_feature_self(f_cv_external, f_cv_single_cvc);
+    // External parameters always report the total force for current step
+    require_feature_self(f_cv_external, f_cv_total_force_current_step);
 
     init_feature(f_cv_single_cvc, "single_component", f_type_static);
 
@@ -1238,10 +1245,16 @@ int colvar::init_dependencies() {
     init_feature(f_cv_linear, "linear", f_type_static);
     init_feature(f_cv_homogeneous, "homogeneous", f_type_static);
 
-    // because total forces are obtained from the previous time step,
-    // we cannot (currently) have colvar values and total forces for the same timestep
     init_feature(f_cv_multiple_ts, "multiple_timestep", f_type_static);
-    exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
+
+    // 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++) {
@@ -1262,6 +1275,10 @@ int colvar::init_dependencies() {
   feature_states[f_cv_fdiff_velocity].available =
     cvm::main()->proxy->simulation_running();
 
+  // Some back-ends report current total forces for all colvars
+  if (cvm::main()->proxy->total_forces_same_step())
+    enable(f_cv_total_force_current_step);
+
   return COLVARS_OK;
 }
 
@@ -1388,23 +1405,22 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
     cvm::log("Calculating colvar \""+this->name+"\", components "+
              cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
 
-  colvarproxy *proxy = cvm::main()->proxy;
   int error_code = COLVARS_OK;
 
   error_code |= check_cvc_range(first_cvc, num_cvcs);
   if (error_code != COLVARS_OK) {
     return error_code;
   }
 
-  if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
+  if ((cvm::step_relative() > 0) && (!is_enabled(f_cv_total_force_current_step))){
     // Use Jacobian derivative from previous timestep
     error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
   }
   // atom coordinates are updated by the next line
   error_code |= calc_cvc_values(first_cvc, num_cvcs);
   error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
   error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
-  if (proxy->total_forces_same_step()){
+  if (is_enabled(f_cv_total_force_current_step)){
     // Use Jacobian derivative from this timestep
     error_code |= calc_cvc_total_force(first_cvc, num_cvcs);
   }
@@ -1421,18 +1437,17 @@ int colvar::collect_cvc_data()
   if (cvm::debug())
     cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
 
-  colvarproxy *proxy = cvm::main()->proxy;
   int error_code = COLVARS_OK;
 
-  if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
+  if ((cvm::step_relative() > 0) && (!is_enabled(f_cv_total_force_current_step))){
     // Total force depends on Jacobian derivative from previous timestep
     // collect_cvc_total_forces() uses the previous value of jd
     error_code |= collect_cvc_total_forces();
   }
   error_code |= collect_cvc_values();
   error_code |= collect_cvc_gradients();
   error_code |= collect_cvc_Jacobians();
-  if (proxy->total_forces_same_step()){
+  if (is_enabled(f_cv_total_force_current_step)){
     // Use Jacobian derivative from this timestep
     error_code |= collect_cvc_total_forces();
   }
@@ -1871,13 +1886,13 @@ void colvar::update_extended_Lagrangian()
     f += fb_actual;
   }
 
-  // fr: bias force on extended variable (without harmonic spring), for output in trajectory
-  fr = f;
-
   // External force has been scaled for an inner-timestep impulse (for the back-end integrator)
   // here we scale it back because this integrator uses only the outer (long) timestep
   f_ext = f / cvm::real(time_step_factor);
 
+  // fr: bias force on extended variable (without harmonic spring), for output in trajectory
+  fr = f_ext;
+
   colvarvalue f_system(fr.type()); // force exterted by the system on the extended DOF
 
   if (is_enabled(f_cv_external)) {
@@ -2382,6 +2397,10 @@ int colvar::set_state_params(std::string const &conf)
              cvm::to_str(x)+"\n");
     x_restart = x;
     after_restart = true;
+    // Externally driven cv (e.g. alchemical lambda) is imposed by restart value
+    if (is_enabled(f_cv_external) && is_enabled(f_cv_extended_Lagrangian)) {
+      cvcs[0]->set_value(x);
+    }
   }
 
   if (is_enabled(f_cv_extended_Lagrangian)) {
@@ -2524,8 +2543,14 @@ std::string const colvar::get_state_params() const
   os << "  name " << name << "\n"
      << "  x "
      << std::setprecision(cvm::cv_prec)
-     << std::setw(cvm::cv_width)
-     << x << "\n";
+     << std::setw(cvm::cv_width);
+  if (is_enabled(f_cv_external)) {
+    // For an external colvar, x is one timestep in the future after integration
+    // write x at beginning of timestep
+    os << x_reported << "\n";
+  } else {
+    os << x << "\n";
+  }
 
   if (is_enabled(f_cv_output_velocity)) {
     os << "  v "

src/colvar.h

@@ -760,7 +760,7 @@ inline colvarvalue const & colvar::total_force() const
 
 inline void colvar::add_bias_force(colvarvalue const &force)
 {
-  check_enabled(f_cv_gradient,
+  check_enabled(f_cv_apply_force,
                 std::string("applying a force to the variable \""+name+"\""));
   if (cvm::debug()) {
     cvm::log("Adding biasing force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");

src/colvarbias.cpp

@@ -93,11 +93,18 @@ int colvarbias::init(std::string const &conf)
     cvm::log("Reinitializing bias \""+name+"\".\n");
   }
 
+  feature_states[f_cvb_step_zero_data].available = true;
+
   colvar_values.resize(num_variables());
   for (i = 0; i < num_variables(); i++) {
     colvar_values[i].type(colvars[i]->value().type());
     colvar_forces[i].type(colvar_values[i].type());
     previous_colvar_forces[i].type(colvar_values[i].type());
+    if (!colvars[i]->is_enabled(f_cv_total_force_current_step)) {
+      // If any colvar does not have current-step total force, then
+      // we can't do step 0 data
+      feature_states[f_cvb_step_zero_data].available = false;
+    }
   }
 
   output_prefix = cvm::output_prefix();
@@ -157,7 +164,7 @@ int colvarbias::init_dependencies() {
     init_feature(f_cvb_step_zero_data, "step_zero_data", f_type_user);
 
     init_feature(f_cvb_apply_force, "apply_force", f_type_user);
-    require_feature_children(f_cvb_apply_force, f_cv_gradient);
+    require_feature_children(f_cvb_apply_force, f_cv_apply_force);
 
     init_feature(f_cvb_bypass_ext_lagrangian, "bypass_extended_Lagrangian_coordinates", f_type_user);
 
@@ -221,8 +228,6 @@ int colvarbias::init_dependencies() {
   // The feature f_cvb_bypass_ext_lagrangian is only implemented by some derived classes
   // (initially, harmonicWalls)
   feature_states[f_cvb_bypass_ext_lagrangian].available = false;
-  // disabled by default; can be changed by derived classes that implement it
-  feature_states[f_cvb_bypass_ext_lagrangian].enabled = false;
 
   return COLVARS_OK;
 }