From 799c7ac3d183824d247482b487ec7fc4c33457ae Mon Sep 17 00:00:00 2001 From: Haohao Fu Date: Sun, 14 Dec 2025 21:24:09 +0800 Subject: [PATCH 1/2] reimplement harmonicDynamics to adapt Jerome's conf --- namd/colvars/src/Makefile.namd | 1 + src/colvar.cpp | 47 +++++++- src/colvar.h | 8 ++ src/colvarbias_restraint.cpp | 136 +++++++++++++++++++---- src/colvarbias_restraint.h | 15 +++ src/colvarcomp.cpp | 6 + src/colvarcomp.h | 14 +-- src/colvarcomp_alchlambda.cpp | 9 ++ src/colvarcomp_harmonicforceconstant.cpp | 70 ++++++++++++ src/colvarcomp_harmonicforceconstant.h | 42 +++++++ 10 files changed, 309 insertions(+), 39 deletions(-) create mode 100644 src/colvarcomp_harmonicforceconstant.cpp create mode 100644 src/colvarcomp_harmonicforceconstant.h diff --git a/namd/colvars/src/Makefile.namd b/namd/colvars/src/Makefile.namd index 9baefa0ea..c14ce25fd 100644 --- a/namd/colvars/src/Makefile.namd +++ b/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 \ diff --git a/src/colvar.cpp b/src/colvar.cpp index 0686e1f5e..b99507d99 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -21,6 +21,7 @@ #include "colvarbias.h" #include "colvars_memstream.h" #include "colvarcomp_torchann.h" +#include "colvarcomp_harmonicforceconstant.h" std::map> colvar::global_cvc_map = std::map>(); @@ -324,6 +325,13 @@ int colvar::init(std::string const &conf) static_cast(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"); add_component_type("alchemical coupling parameter", "alchLambda"); add_component_type("force on alchemical coupling parameter", "alchFLambda"); + add_component_type("force constant as a dynamic variable", "harmonicForceConstant"); add_component_type("arithmetic path collective variables (s)", "aspath"); add_component_type("arithmetic path collective variables (z)", "azpath"); add_component_type("geometrical path collective variables (s)", "gspath"); @@ -1917,8 +1949,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 { diff --git a/src/colvar.h b/src/colvar.h index a8ed2dedc..5b03d13eb 100644 --- a/src/colvar.h +++ b/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(); } diff --git a/src/colvarbias_restraint.cpp b/src/colvarbias_restraint.cpp index 4f33a25d9..bc977f675 100644 --- a/src/colvarbias_restraint.cpp +++ b/src/colvarbias_restraint.cpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include #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; } diff --git a/src/colvarbias_restraint.h b/src/colvarbias_restraint.h index ba29db2ad..e4bfedb85 100644 --- a/src/colvarbias_restraint.h +++ b/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; diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 53854799b..aac1a052b 100644 --- a/src/colvarcomp.cpp +++ b/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 colvar::cvc::cvc_features; diff --git a/src/colvarcomp.h b/src/colvarcomp.h index 277d85e18..e004d4330 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -278,17 +278,8 @@ class colvar::cvc /// \brief Whether or not this CVC will be computed in parallel whenever possible bool b_try_scalable = true; - /// Forcibly set value of CVC - useful for driving an external coordinate, - /// eg. lambda dynamics - inline void set_value(colvarvalue const &new_value, bool now=false) { - x = new_value; - // Cache value to be communicated to back-end between time steps - cvm::proxy->set_alch_lambda(x.real_value); - if (now) { - // If requested (e.g. upon restarting), sync to back-end - cvm::proxy->send_alch_lambda(); - } - } + /// Forcibly set value of CVC - useful for driving an external coordinate + virtual void set_value(colvarvalue const &new_value, bool now=false); protected: @@ -1319,6 +1310,7 @@ class colvar::alch_lambda virtual void calc_force_invgrads(); virtual void calc_Jacobian_derivative(); virtual void apply_force(colvarvalue const &force); + virtual void set_value(colvarvalue const &new_value, bool now=false) override; }; diff --git a/src/colvarcomp_alchlambda.cpp b/src/colvarcomp_alchlambda.cpp index 29168d87c..a358d452e 100644 --- a/src/colvarcomp_alchlambda.cpp +++ b/src/colvarcomp_alchlambda.cpp @@ -84,6 +84,15 @@ void colvar::alch_lambda::apply_force(colvarvalue const & /* force */) // Forces, if any, are applied in colvar::update_extended_Lagrangian() } +void colvar::alch_lambda::set_value(colvarvalue const &new_value, bool now) +{ + x = new_value; + cvm::proxy->set_alch_lambda(x.real_value); + if (now) { + cvm::proxy->send_alch_lambda(); + } +} + colvar::alch_Flambda::alch_Flambda() { diff --git a/src/colvarcomp_harmonicforceconstant.cpp b/src/colvarcomp_harmonicforceconstant.cpp new file mode 100644 index 000000000..e401e3346 --- /dev/null +++ b/src/colvarcomp_harmonicforceconstant.cpp @@ -0,0 +1,70 @@ +// src/colvarcomp_harmonicforceconstant.cpp + +#include "colvarcomp_harmonicforceconstant.h" +#include "colvarmodule.h" + +cvc_harmonicforceconstant::cvc_harmonicforceconstant() + : cvc() // Call base class constructor +{ + set_function_type("harmonicForceConstant"); // Set the type name for this component + // This is a fictitious coordinate, so it does not have gradients with respect to atomic positions. + provide(f_cvc_explicit_gradient, false); + provide(f_cvc_gradient, false); + provide(f_cvc_collect_atom_ids, false); + // This fictitious coordinate does not support total-force calculation + // through inverse gradients / Jacobian machinery. + provide(f_cvc_inv_gradient, false); + provide(f_cvc_Jacobian, false); + + // The colvar object controls the extended Lagrangian dynamics via the f_cv_external flag. + + // This CVC is a scalar defined on the interval [0, 1]. + init_scalar_boundaries(0.0, 1.0); + + x.type(colvarvalue::type_scalar); + x.real_value = 0.0; // Initial value, corresponding to k=0. + + cvm::log("Initializing a harmonicForceConstant component.\n"); +} + +int cvc_harmonicforceconstant::init(std::string const &conf) +{ + // No parameters: behave like alchLambda-style external coordinate + return cvc::init(conf); +} + +void cvc_harmonicforceconstant::calc_force_invgrads() +{ + // Not implemented: this coordinate has no inverse gradients. +} + +void cvc_harmonicforceconstant::calc_value() +{ + // The value of this CV is managed by the parent colvar's extended Lagrangian dynamics. + // This function does nothing. +} + +void cvc_harmonicforceconstant::calc_gradients() +{ + // This is a fictitious coordinate and has no gradients with respect to atomic positions. + // This function does nothing. +} + +void cvc_harmonicforceconstant::calc_Jacobian_derivative() { + // Not implemented + cvm::error("Error: Jacobian derivative is not implemented for harmonicForceConstant.\n", + COLVARS_NOT_IMPLEMENTED); +} + +void cvc_harmonicforceconstant::apply_force(colvarvalue const &force) +{ + // This CVC does not interact with any atoms directly. + // Forces from biases (ABF, Metadynamics) are applied to the extended coordinate + // in colvar::update_forces_energy(). The thermodynamic force from the restraint + // is reported via calc_force_invgrads(). +} + +void cvc_harmonicforceconstant::set_value(colvarvalue const &new_value, bool now) +{ + x = new_value; +} diff --git a/src/colvarcomp_harmonicforceconstant.h b/src/colvarcomp_harmonicforceconstant.h new file mode 100644 index 000000000..1968f393a --- /dev/null +++ b/src/colvarcomp_harmonicforceconstant.h @@ -0,0 +1,42 @@ +// src/colvarcomp_harmonicforceconstant.h + +#ifndef COLVARCOMP_HARMONICFORCECONSTANT_H +#define COLVARCOMP_HARMONICFORCECONSTANT_H + +#include "colvar.h" +#include "colvarcomp.h" + +// Forward declaration to avoid circular dependency with the bias class. +class colvarbias_restraint; + +/// \brief A fictitious scalar coordinate (lambda in [0,1]) intended for extended-Lagrangian dynamics. +/// +/// This CVC does not depend on atomic coordinates and provides no gradients. +/// It is meant to be used as an extended-Lagrangian CV whose value is integrated +/// by colvar::update_extended_Lagrangian(). +/// Forces from biases are applied via colvar::add_bias_force() on the parent colvar. +class cvc_harmonicforceconstant : public colvar::cvc { +public: + cvc_harmonicforceconstant(); + virtual int init(std::string const &conf) override; + virtual ~cvc_harmonicforceconstant() {} + + virtual void calc_value(); + virtual void calc_gradients(); + + /// No total/system force is computed for this fictitious coordinate. + virtual void calc_force_invgrads(); + + /// The force applied to this CVC is handled by the parent colvar's extended Lagrangian dynamics. + virtual void apply_force(colvarvalue const &force); + + /// The Jacobian derivative for this 1D fictitious coordinate is zero. + virtual void calc_Jacobian_derivative(); + + /// Store the value locally; do NOT talk to the MD engine. + virtual void set_value(colvarvalue const &new_value, bool now=false) override; + +protected: +}; + +#endif From 4189bb75e44b5c6613cdadf731d7650a977368fb Mon Sep 17 00:00:00 2001 From: Haohao Fu Date: Sun, 14 Dec 2025 22:06:37 +0800 Subject: [PATCH 2/2] minor fix --- src/colvar.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/colvar.cpp b/src/colvar.cpp index b99507d99..c3bd7c135 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -1727,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; }