Extended dynamics for harmonic restraints

namd/colvars/src/Makefile.namd

@@ -24,6 +24,7 @@ COLVARSLIB = \
 	$(DSTDIR)/colvarcomp_combination.o \
 	$(DSTDIR)/colvarcomp_neuralnetwork.o \
 	$(DSTDIR)/colvarcomp_torchann.o \
+	$(DSTDIR)/colvarcomp_harmonicforceconstant.o \
 	$(DSTDIR)/colvar_neuralnetworkcompute.o \
 	$(DSTDIR)/colvardeps.o \
 	$(DSTDIR)/colvargrid.o \

src/colvar.cpp

@@ -21,6 +21,7 @@
 #include "colvarbias.h"
 #include "colvars_memstream.h"
 #include "colvarcomp_torchann.h"
+#include "colvarcomp_harmonicforceconstant.h"
 
 std::map<std::string, std::function<colvar::cvc *()>> colvar::global_cvc_map =
     std::map<std::string, std::function<colvar::cvc *()>>();
@@ -324,6 +325,13 @@ int colvar::init(std::string const &conf)
 
     static_cast<colvar::alch_lambda *>(cvcs[0].get())->init_alchemy(time_step_factor);
   }
+  // harmonicForceConstant: fictitious external coordinate integrated internally
+  // Needs to behave as an external parameter so that f_cv_apply_force can be enabled
+  // (it has no atomic gradients).
+  if (is_enabled(f_cv_single_cvc) && (cvcs[0]->function_type() == "harmonicForceConstant")) {
+    enable(f_cv_external); 
+    cvm::log("Enabled external parameter mode for harmonicForceConstant to allow extended Lagrangian dynamics.\n");
+  }
 
   // If using scripted biases, any colvar may receive bias forces
   if (cvm::scripted_forces()) {
@@ -356,6 +364,24 @@ int colvar::init(std::string const &conf)
 }
 
 
+void colvar::add_system_force(colvarvalue const &force)
+{
+  if (cvm::debug()) {
+    cvm::log("Adding system force "+cvm::to_str(force)+" to colvar \""+name+"\".\n");
+  }
+  system_force_accumulator += force;
+
+  if (is_enabled(f_cv_external)) {
+    if (cvcs.size() > 0 && cvcs[0]->function_type() == "harmonicForceConstant") {
+        if (ft_reported.type() != value().type()) {
+            ft_reported.type(value());
+        }
+        ft_reported = system_force_accumulator;
+    }
+  }
+}
+
+
 #ifdef LEPTON
 int colvar::init_custom_function(std::string const &conf)
 {
@@ -721,8 +747,13 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
       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);
+
+      if (is_enabled(f_cv_single_cvc) && (cvcs[0]->function_type() == "harmonicForceConstant")) {
+        // bypass the requirement of jacobian
+      } else {
+        // 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 && proxy->simulation_running()) { // Then a finite temperature is required
@@ -926,6 +957,7 @@ void colvar::define_component_types()
   add_component_type<eigenvector>("eigenvector", "eigenvector");
   add_component_type<alch_lambda>("alchemical coupling parameter", "alchLambda");
   add_component_type<alch_Flambda>("force on alchemical coupling parameter", "alchFLambda");
+  add_component_type<cvc_harmonicforceconstant>("force constant as a dynamic variable", "harmonicForceConstant");
   add_component_type<aspath>("arithmetic path collective variables (s)", "aspath");
   add_component_type<azpath>("arithmetic path collective variables (z)", "azpath");
   add_component_type<gspath>("geometrical path collective variables (s)", "gspath");
@@ -1695,12 +1727,13 @@ int colvar::collect_cvc_total_forces()
       // 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;
+    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;
 }
 
@@ -1917,8 +1950,15 @@ void colvar::update_extended_Lagrangian()
   colvarvalue f_system(fr.type()); // force exterted by the system on the extended DOF
 
   if (is_enabled(f_cv_external)) {
-    // Add "alchemical" force from external variable
-    f_system = cvcs[0]->total_force();
+    // External parameter case:
+    // - alchLambda obtains system force from backend (dE/dlambda)
+    // - harmonicForceConstant is a fictitious internal coordinate: no backend system force
+    if (is_enabled(f_cv_single_cvc) && (cvcs[0]->function_type() == "harmonicForceConstant")) {
+      f_system = system_force_accumulator;
+    } else {
+      // Add "alchemical" force from external variable
+      f_system = cvcs[0]->total_force();
+    }
     // f is now irrelevant because we are not applying atomic forces in the simulation
     // just driving the external variable lambda
   } else {

src/colvar.h

@@ -351,6 +351,12 @@ class colvar : public colvarparse, public colvardeps {
   /// Apply a force to the actual value (only meaningful with extended Lagrangian)
   void add_bias_force_actual_value(colvarvalue const &force);
 
+  /// \brief Accumulator for system forces determined internally (e.g. by harmonic restraints on lambda)
+  colvarvalue system_force_accumulator;
+
+  /// \brief Add to the system force accumulator
+  void add_system_force(colvarvalue const &force);
+
   /// \brief Collect all forces on this colvar, integrate internal
   /// equations of motion of internal degrees of freedom; see also
   /// colvar::communicate_forces()
@@ -793,6 +799,8 @@ inline void colvar::reset_bias_force() {
   fb.reset();
   fb_actual.type(value());
   fb_actual.reset();
+  system_force_accumulator.type(value());
+  system_force_accumulator.reset();
 }
 
 

src/colvarbias_restraint.cpp

@@ -10,6 +10,8 @@
 #include <fstream>
 #include <iostream>
 #include <iomanip>
+#include <cmath>
+#include <limits>
 
 #include "colvarmodule.h"
 #include "colvarproxy.h"
@@ -270,6 +272,12 @@ int colvarbias_restraint_moving::update() {
     if (b_chg_walls) update_walls(lambda);
   }
 
+  // For dynamicForceConstantLambda mode: always update k based on the lambda CV value
+  // This is needed because the above conditions may not be satisfied when targetNumSteps is not set
+  if (b_dynamic_force_k) {
+    update_k(0.0);  // The argument is ignored in update_k() when b_dynamic_force_k is true
+  }
+
   return COLVARS_OK;
 }
 
@@ -562,41 +570,99 @@ int colvarbias_restraint_k_moving::init(std::string const &conf)
     starting_force_k = force_k;
   }
 
+  // --- dynamicForceConstantLambda: drive k using an external lambda colvar ---
+  // Note: this enables b_chg_force_k so that colvarbias_restraint_moving::update()
+  // will call update_k() every step.
+  force_k_max = force_k; // store k_max
+  if (get_keyval(conf, "dynamicForceConstantLambda",
+                 dynamic_force_k_lambda_cv_name,
+                 std::string(""))) {
+    b_dynamic_force_k = true;
+    b_chg_force_k = true;  // IMPORTANT: ensure moving::update() calls update_k()
+    starting_force_k = force_k; // not strictly needed, but keeps state consistent
+    get_keyval(conf, "dynamicForceConstantExponent",
+               dynamic_force_k_exponent, 1.0);
+    if (dynamic_force_k_exponent < 1.0) {
+      return cvm::error("Error: dynamicForceConstantExponent must be >= 1.0.\n",
+                        COLVARS_INPUT_ERROR);
+    }
+    dynamic_force_k_lambda_cv = cvm::colvar_by_name(dynamic_force_k_lambda_cv_name);
+    if (!dynamic_force_k_lambda_cv) {
+      return cvm::error("Error: dynamicForceConstantLambda \"" +
+                        dynamic_force_k_lambda_cv_name +
+                        "\" does not match any defined colvar.\n",
+                        COLVARS_INPUT_ERROR);
+    }
+    if (!dynamic_force_k_lambda_cv->is_enabled(f_cv_extended_Lagrangian)) {
+      return cvm::error("Error: dynamicForceConstantLambda colvar \"" +
+                        dynamic_force_k_lambda_cv_name +
+                        "\" must have extendedLagrangian enabled.\n",
+                        COLVARS_INPUT_ERROR);
+    }
+    if (!dynamic_force_k_lambda_cv->is_enabled(f_cv_scalar)) {
+      return cvm::error("Error: dynamicForceConstantLambda colvar \"" +
+                        dynamic_force_k_lambda_cv_name +
+                        "\" must be scalar.\n",
+                        COLVARS_INPUT_ERROR);
+    }
+    // Disallow legacy time-dependent k schedule options
+    if (b_decoupling || (target_force_k > 0.0)) {
+      return cvm::error("Error: dynamicForceConstantLambda is not compatible with "
+                        "targetForceConstant/decoupling in this version.\n",
+                        COLVARS_INPUT_ERROR);
+    }
+    cvm::log("Restraint \"" + name + "\" will use dynamic force constant: "
+             "k_eff = k_max * lambda^n, with lambda colvar \"" +
+             dynamic_force_k_lambda_cv_name + "\" and n=" +
+             cvm::to_str(dynamic_force_k_exponent) + ".\n");
+  }
+
   if (!b_chg_force_k) {
     return COLVARS_OK;
   }
 
-  // parse moving restraint options
-  colvarbias_restraint_moving::init(conf);
+  // parse moving restraint options (only for time-dependent k schedules)
+  if (!b_dynamic_force_k) {
+    colvarbias_restraint_moving::init(conf);
 
-  if (get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule) &&
-      target_nstages > 0) {
-    cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", COLVARS_INPUT_ERROR);
-    return cvm::get_error();
-  }
+    if (get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule) &&
+        target_nstages > 0) {
+      cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", COLVARS_INPUT_ERROR);
+      return cvm::get_error();
+    }
 
-  if (lambda_schedule.size()) {
-    // There is one more lambda-point than stages
-    target_nstages = lambda_schedule.size() - 1;
-  }
+    if (lambda_schedule.size()) {
+      // There is one more lambda-point than stages
+      target_nstages = lambda_schedule.size() - 1;
+    }
 
-  if ((get_keyval(conf, "targetForceExponent", lambda_exp, lambda_exp, parse_deprecated)
-    || get_keyval(conf, "lambdaExponent", lambda_exp, lambda_exp))
-    && !b_chg_force_k) {
-    cvm::error("Error: cannot set lambdaExponent unless a changing force constant is active.\n", COLVARS_INPUT_ERROR);
-  }
-  if (lambda_exp < 1.0) {
-    cvm::log("Warning: for all practical purposes, lambdaExponent should be 1.0 or greater.\n");
+    if ((get_keyval(conf, "targetForceExponent", lambda_exp, lambda_exp, parse_deprecated)
+      || get_keyval(conf, "lambdaExponent", lambda_exp, lambda_exp))
+      && !b_chg_force_k) {
+      cvm::error("Error: cannot set lambdaExponent unless a changing force constant is active.\n", COLVARS_INPUT_ERROR);
+    }
+    if (lambda_exp < 1.0) {
+      cvm::log("Warning: for all practical purposes, lambdaExponent should be 1.0 or greater.\n");
+    }
   }
 
   return COLVARS_OK;
 }
 
 
-void colvarbias_restraint_k_moving::update_k(cvm::real lambda) {
+void colvarbias_restraint_k_moving::update_k(cvm::real lambda_in) {
+  if (b_dynamic_force_k) {
+    cvm::real lambda_dyn = dynamic_force_k_lambda_cv->value().real_value;
+    if (lambda_dyn < 0.0) lambda_dyn = 0.0;
+    if (lambda_dyn > 1.0) lambda_dyn = 1.0;
+    force_k = force_k_max * cvm::pow(lambda_dyn, dynamic_force_k_exponent);
+    force_k_incr = 0.0; // dynamic k doesn't accumulate work in same way
+    return;
+  }
+  // original moving-k behavior
   cvm::real const force_k_old = force_k;
 
-  force_k = starting_force_k + (target_force_k - starting_force_k) * cvm::pow(lambda, lambda_exp);
+  force_k = starting_force_k + (target_force_k - starting_force_k) * cvm::pow(lambda_in, lambda_exp);
   force_k_incr = force_k - force_k_old;
   if (!target_nstages && (cvm::step_absolute() > first_step + target_nsteps)) {
     force_k_incr = 0.0;
@@ -605,6 +671,26 @@ void colvarbias_restraint_k_moving::update_k(cvm::real lambda) {
            this->name+"\": "+cvm::to_str(force_k)+".\n");
 }
 
+void colvarbias_restraint_k_moving::apply_dynamic_lambda_force()
+{
+  if (!b_dynamic_force_k) return;
+  cvm::real lambda = dynamic_force_k_lambda_cv->value().real_value;
+  if (lambda < 0.0) lambda = 0.0;
+  if (lambda > 1.0) lambda = 1.0;
+  cvm::real dU_dk = 0.0;
+  for (size_t i = 0; i < num_variables(); i++) {
+    dU_dk += d_restraint_potential_dk(i);
+  }
+  cvm::real const dk_dlambda =
+    force_k_max * dynamic_force_k_exponent *
+    ((dynamic_force_k_exponent == 1.0) ? 1.0 : cvm::pow(lambda, dynamic_force_k_exponent - 1.0));
+  cvm::real const F_lambda = - dU_dk * dk_dlambda;
+  colvarvalue f_l;
+  f_l.type(colvarvalue::type_scalar);
+  f_l.real_value = F_lambda;
+  dynamic_force_k_lambda_cv->add_system_force(f_l);
+}
+
 
 cvm::real colvarbias_restraint_k_moving::dU_dlambda_k() const {
   if (force_k == 0.0) return 0.0;
@@ -735,16 +821,18 @@ int colvarbias_restraint_harmonic::update()
   // update the TI estimator (if defined)
   error_code |= colvarbias_ti::update();
 
-  // update parameters (centers or force constant)
+  // update parameters (this is where moving restraints are implemented)
   error_code |= colvarbias_restraint_moving::update();
-
+  
   // update restraint energy and forces
   error_code |= colvarbias_restraint::update();
-
-  // update accumulated work using the current forces
+  
+  // Update accumulated work using the current forces (if moving restraints are active)
   error_code |= colvarbias_restraint_centers_moving::update_acc_work();
   error_code |= colvarbias_restraint_k_moving::update_acc_work();
 
+  // Apply thermodynamic force to dynamic lambda CV (if active)
+  colvarbias_restraint_k_moving::apply_dynamic_lambda_force();
   return error_code;
 }
 

src/colvarbias_restraint.h

@@ -87,6 +87,7 @@ class colvarbias_restraint_k
   colvarbias_restraint_k(char const *key);
   virtual int init(std::string const &conf);
   virtual int change_configuration(std::string const &conf);
+  cvm::real get_force_k() const { return force_k; }
 
 protected:
 
@@ -128,6 +129,10 @@ class colvarbias_restraint_moving
   /// \brief Changing wall locations?
   bool b_chg_walls = false;
 
+  /// \brief Dynamic force constant driven by external lambda CV?
+  /// If true, update_k() should be called every step
+  bool b_dynamic_force_k = false;
+
   /// @brief Update the force constant by interpolating between initial and target
   virtual void update_k(cvm::real /* lambda */) {}
   /// @brief Update the centers by interpolating between initial and target
@@ -256,6 +261,16 @@ class colvarbias_restraint_k_moving
   /// \brief Exponent for varying the force constant
   cvm::real lambda_exp;
 
+  // --- Dynamic k driven by an external lambda colvar (extended-Lagrangian CV) ---
+  // Note: b_dynamic_force_k is defined in base class colvarbias_restraint_moving
+  std::string dynamic_force_k_lambda_cv_name;
+  colvar *dynamic_force_k_lambda_cv = nullptr;
+  cvm::real dynamic_force_k_exponent = 1.0;
+  cvm::real force_k_max = 0.0; // store k_max = user "forceConstant"
+
+  // Apply thermodynamic force to lambda CV: F_lambda = - dU/dlambda
+  void apply_dynamic_lambda_force();
+
   /// \brief Increment of the force constant at each step
   cvm::real force_k_incr;
 

src/colvarcomp.cpp

@@ -1148,6 +1148,12 @@ void colvar::cvc::wrap(colvarvalue &x_unwrapped) const
 }
 
 
+void colvar::cvc::set_value(colvarvalue const &new_value, bool /*now*/)
+{
+  x = new_value;
+}
+
+
 // Static members
 
 std::vector<colvardeps::feature *> colvar::cvc::cvc_features;