From 2a9be4275850a8ed653e6b5ae055c6a5db81534b Mon Sep 17 00:00:00 2001 From: Giacomo Fiorin Date: Tue, 3 Aug 2021 18:03:09 -0400 Subject: [PATCH] Update Colvars to version 2021-08-03 --- lib/colvars/Makefile.common | 1 + lib/colvars/Makefile.deps | 19 +- lib/colvars/colvar.cpp | 147 +++++-- lib/colvars/colvar.h | 13 +- lib/colvars/colvar_arithmeticpath.h | 12 + lib/colvars/colvaratoms.cpp | 83 ++-- lib/colvars/colvaratoms.h | 34 +- lib/colvars/colvarbias_abf.cpp | 13 +- lib/colvars/colvarbias_abf.h | 2 + lib/colvars/colvarbias_histogram.cpp | 4 +- lib/colvars/colvarbias_histogram.h | 1 - lib/colvars/colvarbias_meta.cpp | 199 +++++---- lib/colvars/colvarbias_meta.h | 6 +- lib/colvars/colvarcomp.cpp | 12 +- lib/colvars/colvarcomp.h | 106 ++++- lib/colvars/colvarcomp_alchlambda.cpp | 56 +++ lib/colvars/colvarcomp_apath.cpp | 42 +- lib/colvars/colvarcomp_distances.cpp | 35 +- lib/colvars/colvarcomp_gpath.cpp | 8 +- lib/colvars/colvarcomp_rotations.cpp | 332 +++++++++++++++ lib/colvars/colvarcomp_volmaps.cpp | 93 ++++- lib/colvars/colvardeps.cpp | 34 +- lib/colvars/colvardeps.h | 5 + lib/colvars/colvargrid.cpp | 101 ++++- lib/colvars/colvargrid.h | 29 +- lib/colvars/colvarmodule.cpp | 111 +++-- lib/colvars/colvarmodule.h | 33 +- lib/colvars/colvarmodule_utils.h | 80 ++++ lib/colvars/colvarparse.cpp | 24 ++ lib/colvars/colvarparse.h | 10 +- lib/colvars/colvarproxy.cpp | 173 +++++++- lib/colvars/colvarproxy.h | 164 +++++++- lib/colvars/colvarproxy_tcl.cpp | 2 + lib/colvars/colvarproxy_tcl.h | 8 +- lib/colvars/colvarproxy_volmaps.cpp | 73 +++- lib/colvars/colvarproxy_volmaps.h | 58 ++- lib/colvars/colvars_version.h | 2 +- lib/colvars/colvarscript.cpp | 447 +++++++++++++++++++-- lib/colvars/colvarscript.h | 135 +++++-- lib/colvars/colvarscript_commands.cpp | 48 +++ lib/colvars/colvarscript_commands.h | 169 +++++++- lib/colvars/colvarscript_commands_bias.h | 40 +- lib/colvars/colvarscript_commands_colvar.h | 89 ++-- lib/colvars/colvartypes.cpp | 3 + lib/colvars/colvarvalue.cpp | 14 +- lib/colvars/colvarvalue.h | 6 +- src/COLVARS/colvarproxy_lammps_version.h | 2 +- src/COLVARS/fix_colvars.cpp | 6 +- 48 files changed, 2557 insertions(+), 527 deletions(-) create mode 100644 lib/colvars/colvarcomp_alchlambda.cpp create mode 100644 lib/colvars/colvarmodule_utils.h diff --git a/lib/colvars/Makefile.common b/lib/colvars/Makefile.common index ce501dc5f3..5743507507 100644 --- a/lib/colvars/Makefile.common +++ b/lib/colvars/Makefile.common @@ -30,6 +30,7 @@ COLVARS_SRCS = \ colvarbias_histogram.cpp \ colvarbias_meta.cpp \ colvarbias_restraint.cpp \ + colvarcomp_alchlambda.cpp \ colvarcomp_angles.cpp \ colvarcomp_apath.cpp \ colvarcomp_coordnums.cpp \ diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps index 78e8768edd..cdf948a471 100644 --- a/lib/colvars/Makefile.deps +++ b/lib/colvars/Makefile.deps @@ -86,6 +86,20 @@ $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \ lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ lepton/include/lepton/Exception.h \ lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h +$(COLVARS_OBJ_DIR)colvarcomp_alchlambda.o: colvarcomp_alchlambda.cpp \ + colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \ + colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \ + lepton/include/lepton/CompiledExpression.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/windowsIncludes.h \ + lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/ExpressionProgram.h \ + lepton/include/lepton/ExpressionTreeNode.h \ + lepton/include/lepton/Operation.h lepton/include/lepton/CustomFunction.h \ + lepton/include/lepton/Exception.h \ + lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ + colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \ colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \ colvarparse.h colvarparams.h colvardeps.h lepton/include/Lepton.h \ @@ -285,7 +299,7 @@ $(COLVARS_OBJ_DIR)colvarproxy.o: colvarproxy.cpp colvarmodule.h \ lepton/include/lepton/Exception.h \ lepton/include/lepton/ParsedExpression.h lepton/include/lepton/Parser.h \ colvarscript_commands.h colvarscript_commands_colvar.h \ - colvarscript_commands_bias.h colvaratoms.h + colvarscript_commands_bias.h colvaratoms.h colvarmodule_utils.h $(COLVARS_OBJ_DIR)colvarproxy_replicas.o: colvarproxy_replicas.cpp \ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h @@ -294,7 +308,8 @@ $(COLVARS_OBJ_DIR)colvarproxy_tcl.o: colvarproxy_tcl.cpp colvarmodule.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvaratoms.h colvarparse.h \ colvarparams.h colvardeps.h $(COLVARS_OBJ_DIR)colvarproxy_volmaps.o: colvarproxy_volmaps.cpp \ - colvarmodule.h colvars_version.h colvarproxy_volmaps.h + colvarmodule.h colvars_version.h colvarproxy_volmaps.h \ + colvarmodule_utils.h $(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarproxy.h \ colvarmodule.h colvars_version.h colvartypes.h colvarvalue.h \ colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \ diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index ed6fec1ac5..54e11aaa75 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -96,7 +96,7 @@ int colvar::init(std::string const &conf) "", colvarparse::parse_silent)) { enable(f_cv_scripted); - cvm::log("This colvar uses scripted function \"" + scripted_function + "\"."); + cvm::log("This colvar uses scripted function \"" + scripted_function + "\".\n"); std::string type_str; get_keyval(conf, "scriptedFunctionType", type_str, "scalar"); @@ -134,7 +134,7 @@ int colvar::init(std::string const &conf) std::sort(cvcs.begin(), cvcs.end(), compare); if(cvcs.size() > 1) { - cvm::log("Sorted list of components for this scripted colvar:"); + cvm::log("Sorted list of components for this scripted colvar:\n"); for (i = 0; i < cvcs.size(); i++) { cvm::log(cvm::to_str(i+1) + " " + cvcs[i]->name); } @@ -298,6 +298,11 @@ int colvar::init(std::string const &conf) error_code |= init_extended_Lagrangian(conf); error_code |= init_output_flags(conf); + // Detect if we have one component that is an alchemical lambda + if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type == "alch_lambda") { + enable(f_cv_external); + } + // Now that the children are defined we can solve dependencies enable(f_cv_active); @@ -636,7 +641,6 @@ int colvar::init_extended_Lagrangian(std::string const &conf) x_ext.type(colvarvalue::type_notset); v_ext.type(value()); fr.type(value()); - const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); if (temp <= 0.0) { if (found) @@ -654,7 +658,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf) return INPUT_ERROR; } ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); - cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2"); + cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2\n"); get_keyval(conf, "extendedTimeConstant", extended_period, 200.0); if (extended_period <= 0.0) { @@ -662,7 +666,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf) } ext_mass = (cvm::boltzmann() * temp * extended_period * extended_period) / (4.0 * PI * PI * tolerance * tolerance); - cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)"); + cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)\n"); { bool b_output_energy; @@ -770,7 +774,9 @@ template int colvar::init_components_type(std::string c if ( (cvcp->function_type != std::string("distance_z")) && (cvcp->function_type != std::string("dihedral")) && (cvcp->function_type != std::string("polar_phi")) && - (cvcp->function_type != std::string("spin_angle")) ) { + (cvcp->function_type != std::string("spin_angle")) && + (cvcp->function_type != std::string("euler_phi")) && + (cvcp->function_type != std::string("euler_psi"))) { cvm::error("Error: invalid use of period and/or " "wrapAround in a \""+ std::string(def_config_key)+ @@ -860,6 +866,7 @@ int colvar::init_components(std::string const &conf) "inertia", "inertia"); error_code |= init_components_type(conf, "moment of inertia around an axis", "inertiaZ"); error_code |= init_components_type(conf, "eigenvector", "eigenvector"); + error_code |= init_components_type(conf, "alchemical coupling parameter", "alchLambda"); error_code |= init_components_type(conf, "geometrical path collective variables (s)", "gspath"); error_code |= init_components_type(conf, "geometrical path collective variables (z)", "gzpath"); error_code |= init_components_type(conf, "linear combination of other collective variables", "linearCombination"); @@ -867,6 +874,9 @@ int colvar::init_components(std::string const &conf) error_code |= init_components_type(conf, "geometrical path collective variables (z) for other CVs", "gzpathCV"); error_code |= init_components_type(conf, "arithmetic path collective variables (s) for other CVs", "aspathCV"); error_code |= init_components_type(conf, "arithmetic path collective variables (s) for other CVs", "azpathCV"); + error_code |= init_components_type(conf, "euler phi angle of the optimal orientation", "eulerPhi"); + error_code |= init_components_type(conf, "euler psi angle of the optimal orientation", "eulerPsi"); + error_code |= init_components_type(conf, "euler theta angle of the optimal orientation", "eulerTheta"); error_code |= init_components_type(conf, "total value of atomic map", "mapTotal"); @@ -1074,6 +1084,7 @@ int colvar::init_dependencies() { init_feature(f_cv_hide_Jacobian, "hide_Jacobian_force", f_type_user); require_feature_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated + exclude_feature_self(f_cv_hide_Jacobian, f_cv_extended_Lagrangian); init_feature(f_cv_extended_Lagrangian, "extended_Lagrangian", f_type_user); require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar); @@ -1082,6 +1093,9 @@ int colvar::init_dependencies() { 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); + require_feature_self(f_cv_external, f_cv_single_cvc); + init_feature(f_cv_single_cvc, "single_component", f_type_static); init_feature(f_cv_linear, "linear", f_type_static); @@ -1195,6 +1209,21 @@ std::vector > colvar::get_atom_lists() } +std::vector const &colvar::get_volmap_ids() +{ + volmap_ids_.resize(cvcs.size()); + for (size_t i = 0; i < cvcs.size(); i++) { + if (cvcs[i]->param_exists("mapID") == COLVARS_OK) { + volmap_ids_[i] = + *(reinterpret_cast(cvcs[i]->get_param_ptr("mapID"))); + } else { + volmap_ids_[i] = -1; + } + } + return volmap_ids_; +} + + colvar::~colvar() { // There is no need to call free_children_deps() here @@ -1552,9 +1581,11 @@ int colvar::collect_cvc_total_forces() } } - if (!is_enabled(f_cv_hide_Jacobian)) { + 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 // 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; } } @@ -1687,14 +1718,16 @@ cvm::real colvar::update_forces_energy() // add the biases' force, which at this point should already have // been summed over each bias using this colvar + // fb is already multiplied by the relevant time step factor for each bias f += fb; if (is_enabled(f_cv_Jacobian)) { // the instantaneous Jacobian force was not included in the reported total force; // instead, it is subtracted from the applied force (silent Jacobian correction) // This requires the Jacobian term for the *current* timestep + // Need to scale it for impulse MTS if (is_enabled(f_cv_hide_Jacobian)) - f -= fj; + f -= fj * cvm::real(time_step_factor); } // At this point f is the force f from external biases that will be applied to the @@ -1727,26 +1760,44 @@ cvm::real colvar::update_forces_energy() colvarvalue f_ext(fr.type()); // force acting on the extended variable f_ext.reset(); - // the total force is applied to the fictitious mass, while the - // atoms only feel the harmonic force + wall force - // fr: bias force on extended variable (without harmonic spring), for output in trajectory - // f_ext: total force on extended variable (including harmonic spring) - // f: - initially, external biasing force - // - after this code block, colvar force to be applied to atomic coordinates - // ie. spring force (fb_actual will be added just below) + if (is_enabled(f_cv_external)) { + // There are no forces on the "actual colvar" bc there is no gradient wrt atomic coordinates + // So we apply this to the extended DOF + f += fb_actual; + } + fr = f; // External force has been scaled for a 1-timestep impulse, scale it back because we will // integrate it with the colvar's own timestep factor f_ext = f / cvm::real(time_step_factor); - f_ext += (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x); - f = (-0.5 * ext_force_k) * this->dist2_rgrad(x_ext, x); - // Coupling force is a slow force, to be applied to atomic coords impulse-style - f *= cvm::real(time_step_factor); + + 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(); + // f is now irrelevant because we are not applying atomic forces in the simulation + // just driving the external variable lambda + } else { + // the total force is applied to the fictitious mass, while the + // atoms only feel the harmonic force + wall force + // fr: bias force on extended variable (without harmonic spring), for output in trajectory + // f_ext: total force on extended variable (including harmonic spring) + // f: - initially, external biasing force + // - after this code block, colvar force to be applied to atomic coordinates + // ie. spring force (fb_actual will be added just below) + f_system = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x); + f = -1.0 * f_system; + // Coupling force is a slow force, to be applied to atomic coords impulse-style + // over a single MD timestep + f *= cvm::real(time_step_factor); + } + 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 - ft_reported = (-0.5 * ext_force_k) * this->dist2_lgrad(x_ext, x); + // 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 @@ -1778,10 +1829,19 @@ cvm::real colvar::update_forces_energy() (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) { x_ext -= 2.0 * delta; v_ext *= -1.0; + if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) || + (is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) { + cvm::error("Error: extended coordinate value " + cvm::to_str(x_ext) + " is still outside boundaries after reflection.\n"); + } } x_ext.apply_constraints(); this->wrap(x_ext); + if (is_enabled(f_cv_external)) { + // Colvar value is constrained to the extended value + x = x_ext; + cvcs[0]->set_value(x_ext); + } } else { // If this is a postprocessing run (eg. in VMD), the extended DOF // is equal to the actual coordinate @@ -1789,9 +1849,11 @@ cvm::real colvar::update_forces_energy() } } - // Now adding the force on the actual colvar (for those biases that - // bypass the extended Lagrangian mass) - f += fb_actual; + if (!is_enabled(f_cv_external)) { + // Now adding the force on the actual colvar (for those biases that + // bypass the extended Lagrangian mass) + f += fb_actual; + } if (cvm::debug()) cvm::log("Done updating colvar \""+this->name+"\".\n"); @@ -1954,7 +2016,7 @@ int colvar::update_cvc_flags() int colvar::update_cvc_config(std::vector const &confs) { - cvm::log("Updating configuration for colvar \""+name+"\""); + cvm::log("Updating configuration for colvar \""+name+"\"\n"); if (confs.size() != cvcs.size()) { return cvm::error("Error: Wrong number of CVC config strings. " @@ -2076,6 +2138,15 @@ bool colvar::periodic_boundaries() const cvm::real colvar::dist2(colvarvalue const &x1, colvarvalue const &x2) const { + if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) { + if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) { + cvm::real diff = x1.real_value - x2.real_value; + const cvm::real period_lower_boundary = wrap_center - period / 2.0; + const cvm::real period_upper_boundary = wrap_center + period / 2.0; + diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff)); + return diff * diff; + } + } if (is_enabled(f_cv_homogeneous)) { return (cvcs[0])->dist2(x1, x2); } else { @@ -2086,6 +2157,15 @@ cvm::real colvar::dist2(colvarvalue const &x1, colvarvalue colvar::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { + if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) { + if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) { + cvm::real diff = x1.real_value - x2.real_value; + const cvm::real period_lower_boundary = wrap_center - period / 2.0; + const cvm::real period_upper_boundary = wrap_center + period / 2.0; + diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff)); + return 2.0 * diff; + } + } if (is_enabled(f_cv_homogeneous)) { return (cvcs[0])->dist2_lgrad(x1, x2); } else { @@ -2096,6 +2176,15 @@ colvarvalue colvar::dist2_lgrad(colvarvalue const &x1, colvarvalue colvar::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { + if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) { + if (is_enabled(f_cv_periodic) && is_enabled(f_cv_scalar)) { + cvm::real diff = x1.real_value - x2.real_value; + const cvm::real period_lower_boundary = wrap_center - period / 2.0; + const cvm::real period_upper_boundary = wrap_center + period / 2.0; + diff = (diff < period_lower_boundary ? diff + period : ( diff > period_upper_boundary ? diff - period : diff)); + return (-2.0) * diff; + } + } if (is_enabled(f_cv_homogeneous)) { return (cvcs[0])->dist2_rgrad(x1, x2); } else { @@ -2302,7 +2391,7 @@ std::ostream & colvar::write_traj_label(std::ostream & os) os << " " << cvm::wrap_string(this->name, this_cv_width); - if (is_enabled(f_cv_extended_Lagrangian)) { + if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) { // extended DOF os << " r_" << cvm::wrap_string(this->name, this_cv_width-2); @@ -2314,7 +2403,7 @@ std::ostream & colvar::write_traj_label(std::ostream & os) os << " v_" << cvm::wrap_string(this->name, this_cv_width-2); - if (is_enabled(f_cv_extended_Lagrangian)) { + if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) { // extended DOF os << " vr_" << cvm::wrap_string(this->name, this_cv_width-3); @@ -2347,7 +2436,7 @@ std::ostream & colvar::write_traj(std::ostream &os) os << " "; if (is_enabled(f_cv_output_value)) { - if (is_enabled(f_cv_extended_Lagrangian)) { + if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << x; @@ -2360,7 +2449,7 @@ std::ostream & colvar::write_traj(std::ostream &os) if (is_enabled(f_cv_output_velocity)) { - if (is_enabled(f_cv_extended_Lagrangian)) { + if (is_enabled(f_cv_extended_Lagrangian) && !is_enabled(f_cv_external)) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) << v_fdiff; diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 5cf0077a39..2856170727 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -192,7 +192,7 @@ protected: /// Amplitude of Gaussian white noise for Langevin extended dynamics cvm::real ext_sigma; - /// \brief Harmonic restraint force + /// \brief Applied force on extended DOF, for output (unscaled if using MTS) colvarvalue fr; /// \brief Jacobian force, when Jacobian_force is enabled @@ -591,6 +591,7 @@ public: class alpha_dihedrals; class alpha_angles; class dihedPC; + class alch_lambda; class componentDisabled; class CartesianBasedPath; class gspath; @@ -601,6 +602,9 @@ public: class gzpathCV; class aspathCV; class azpathCV; + class euler_phi; + class euler_psi; + class euler_theta; // non-scalar components class distance_vec; @@ -658,6 +662,9 @@ protected: static std::map> global_cvc_map; #endif + /// Volmap numeric IDs, one for each CVC (-1 if not available) + std::vector volmap_ids_; + public: /// \brief Sorted array of (zero-based) IDs for all atoms involved @@ -670,6 +677,10 @@ public: /// \brief Get vector of vectors of atom IDs for all atom groups virtual std::vector > get_atom_lists(); + + /// Volmap numeric IDs, one for each CVC (-1 if not available) + std::vector const &get_volmap_ids(); + }; diff --git a/lib/colvars/colvar_arithmeticpath.h b/lib/colvars/colvar_arithmeticpath.h index 552edf23e3..cf613af389 100644 --- a/lib/colvars/colvar_arithmeticpath.h +++ b/lib/colvars/colvar_arithmeticpath.h @@ -7,6 +7,7 @@ #include #include #include +#include namespace ArithmeticPathCV { @@ -24,6 +25,7 @@ public: virtual void computeValue(); virtual void computeDerivatives(); virtual void compute(); + virtual void reComputeLambda(const vector& rmsd_between_refs); protected: scalar_type lambda; vector weights; @@ -124,6 +126,16 @@ void ArithmeticPathBase::computeDerivative } } +template +void ArithmeticPathBase::reComputeLambda(const vector& rmsd_between_refs) { + scalar_type mean_square_displacements = 0.0; + for (size_t i_frame = 1; i_frame < total_frames; ++i_frame) { + cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n")); + mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1]; + } + mean_square_displacements /= scalar_type(total_frames - 1); + lambda = 1.0 / mean_square_displacements; +} } #endif // ARITHMETICPATHCV_H diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 1a4280a202..9bbd91cf05 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -37,6 +37,7 @@ cvm::atom::atom(int atom_number) } id = p->get_atom_id(index); update_mass(); + update_charge(); reset_data(); } @@ -53,6 +54,7 @@ cvm::atom::atom(cvm::residue_id const &residue, } id = p->get_atom_id(index); update_mass(); + update_charge(); reset_data(); } @@ -62,6 +64,7 @@ cvm::atom::atom(atom const &a) { id = (cvm::proxy)->get_atom_id(index); update_mass(); + update_charge(); reset_data(); } @@ -214,8 +217,6 @@ int cvm::atom_group::init() index = -1; b_dummy = false; - b_center = false; - b_rotate = false; b_user_defined_fit = false; fitting_group = NULL; @@ -240,8 +241,9 @@ int cvm::atom_group::init_dependencies() { } init_feature(f_ag_active, "active", f_type_dynamic); - init_feature(f_ag_center, "translational_fit", f_type_static); - init_feature(f_ag_rotate, "rotational_fit", f_type_static); + init_feature(f_ag_center, "center_to_reference", f_type_user); + init_feature(f_ag_center_origin, "center_to_origin", f_type_user); + init_feature(f_ag_rotate, "rotate_to_origin", f_type_user); init_feature(f_ag_fitting_group, "fitting_group", f_type_static); init_feature(f_ag_explicit_gradient, "explicit_atom_gradient", f_type_dynamic); init_feature(f_ag_fit_gradients, "fit_gradients", f_type_user); @@ -272,6 +274,10 @@ int cvm::atom_group::init_dependencies() { // Features that are implemented (or not) by all atom groups feature_states[f_ag_active].available = true; + feature_states[f_ag_center].available = true; + feature_states[f_ag_center_origin].available = true; + feature_states[f_ag_rotate].available = true; + // f_ag_scalable_com is provided by the CVC iff it is COM-based feature_states[f_ag_scalable_com].available = false; // TODO make f_ag_scalable depend on f_ag_scalable_com (or something else) @@ -317,6 +323,9 @@ void cvm::atom_group::update_total_mass() total_mass += ai->mass; } } + if (total_mass < 1e-15) { + cvm::error("ERROR: " + description + " has zero total mass.\n"); + } } @@ -381,12 +390,22 @@ int cvm::atom_group::parse(std::string const &group_conf) // We need to know about fitting to decide whether the group is scalable // and we need to know about scalability before adding atoms - bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false); - bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false); + bool b_defined_center = get_keyval_feature(this, group_conf, "centerToOrigin", f_ag_center_origin, false); + // Legacy alias + b_defined_center |= get_keyval_feature(this, group_conf, "centerReference", f_ag_center, is_enabled(f_ag_center_origin), parse_deprecated); + b_defined_center |= get_keyval_feature(this, group_conf, "centerToReference", f_ag_center, is_enabled(f_ag_center)); + + if (is_enabled(f_ag_center_origin) && ! is_enabled(f_ag_center)) { + return cvm::error("centerToReference may not be disabled if centerToOrigin is enabled.\n"); + } + // Legacy alias + bool b_defined_rotate = get_keyval_feature(this, group_conf, "rotateReference", f_ag_rotate, false, parse_deprecated); + b_defined_rotate |= get_keyval_feature(this, group_conf, "rotateToReference", f_ag_rotate, is_enabled(f_ag_rotate)); + // is the user setting explicit options? b_user_defined_fit = b_defined_center || b_defined_rotate; - if (is_available(f_ag_scalable_com) && !b_rotate && !b_center) { + if (is_available(f_ag_scalable_com) && !is_enabled(f_ag_rotate) && !is_enabled(f_ag_center)) { enable(f_ag_scalable_com); enable(f_ag_scalable); } @@ -507,12 +526,8 @@ int cvm::atom_group::parse(std::string const &group_conf) // whether these atoms will ever receive forces or not bool enable_forces = true; - // disableForces is deprecated - if (get_keyval(group_conf, "enableForces", enable_forces, true)) { - noforce = !enable_forces; - } else { - get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent); - } + get_keyval(group_conf, "enableForces", enable_forces, true, colvarparse::parse_silent); + noforce = !enable_forces; } // Now that atoms are defined we can parse the detailed fitting options @@ -760,10 +775,10 @@ std::string const cvm::atom_group::print_atom_ids() const int cvm::atom_group::parse_fitting_options(std::string const &group_conf) { - if (b_center || b_rotate) { + if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) { if (b_dummy) - cvm::error("Error: centerReference or rotateReference " + cvm::error("Error: centerToReference or rotateToReference " "cannot be defined for a dummy atom.\n"); bool b_ref_pos_group = false; @@ -827,7 +842,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf) if (ref_pos.size()) { - if (b_rotate) { + if (is_enabled(f_ag_rotate)) { if (ref_pos.size() != group_for_fit->size()) cvm::error("Error: the number of reference positions provided("+ cvm::to_str(ref_pos.size())+ @@ -846,7 +861,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf) return COLVARS_ERROR; } - if (b_rotate && !noforce) { + if (is_enabled(f_ag_rotate) && !noforce) { cvm::log("Warning: atom group \""+key+ "\" will be aligned to a fixed orientation given by the reference positions provided. " "If the internal structure of the group changes too much (i.e. its RMSD is comparable " @@ -865,7 +880,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf) bool b_fit_gradients; get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true); - if (b_fit_gradients && (b_center || b_rotate)) { + if (b_fit_gradients && (is_enabled(f_ag_center) || is_enabled(f_ag_rotate))) { enable(f_ag_fit_gradients); } } @@ -879,7 +894,7 @@ void cvm::atom_group::do_feature_side_effects(int id) // If enabled features are changed upstream, the features below should be refreshed switch (id) { case f_ag_fit_gradients: - if (b_center || b_rotate) { + if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) { atom_group *group_for_fit = fitting_group ? fitting_group : this; group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0)); rot.request_group1_gradients(group_for_fit->size()); @@ -971,7 +986,7 @@ int cvm::atom_group::calc_required_properties() calc_center_of_geometry(); if (!is_enabled(f_ag_scalable)) { - if (b_center || b_rotate) { + if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) { if (fitting_group) { fitting_group->calc_center_of_geometry(); } @@ -1001,7 +1016,7 @@ void cvm::atom_group::calc_apply_roto_translation() fitting_group->cog_orig = fitting_group->center_of_geometry(); } - if (b_center) { + if (is_enabled(f_ag_center)) { // center on the origin first cvm::atom_pos const rpg_cog = fitting_group ? fitting_group->center_of_geometry() : this->center_of_geometry(); @@ -1011,9 +1026,9 @@ void cvm::atom_group::calc_apply_roto_translation() } } - if (b_rotate) { - // rotate the group (around the center of geometry if b_center is - // true, around the origin otherwise) + if (is_enabled(f_ag_rotate)) { + // rotate the group (around the center of geometry if f_ag_center is + // enabled, around the origin otherwise) rot.calc_optimal_rotation(fitting_group ? fitting_group->positions() : this->positions(), @@ -1030,7 +1045,7 @@ void cvm::atom_group::calc_apply_roto_translation() } } - if (b_center) { + if (is_enabled(f_ag_center) && !is_enabled(f_ag_center_origin)) { // align with the center of geometry of ref_pos apply_translation(ref_pos_cog); if (fitting_group) { @@ -1063,7 +1078,7 @@ void cvm::atom_group::read_velocities() { if (b_dummy) return; - if (b_rotate) { + if (is_enabled(f_ag_rotate)) { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_velocity(); @@ -1084,7 +1099,7 @@ void cvm::atom_group::read_total_forces() { if (b_dummy) return; - if (b_rotate) { + if (is_enabled(f_ag_rotate)) { for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_total_force(); @@ -1172,14 +1187,14 @@ void cvm::atom_group::calc_fit_gradients() cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this; - if (b_center) { + if (is_enabled(f_ag_center)) { // add the center of geometry contribution to the gradients cvm::rvector atom_grad; for (size_t i = 0; i < this->size(); i++) { atom_grad += atoms[i].grad; } - if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad); + if (is_enabled(f_ag_rotate)) atom_grad = (rot.inverse()).rotate(atom_grad); atom_grad *= (-1.0)/(cvm::real(group_for_fit->size())); for (size_t j = 0; j < group_for_fit->size(); j++) { @@ -1187,7 +1202,7 @@ void cvm::atom_group::calc_fit_gradients() } } - if (b_rotate) { + if (is_enabled(f_ag_rotate)) { // add the rotation matrix contribution to the gradients cvm::rotation const rot_inv = rot.inverse(); @@ -1196,7 +1211,7 @@ void cvm::atom_group::calc_fit_gradients() // compute centered, unrotated position cvm::atom_pos const pos_orig = - rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos))); + rot_inv.rotate((is_enabled(f_ag_center) ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos))); // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i cvm::quaternion const dxdq = @@ -1340,7 +1355,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force) return; } - if (b_rotate) { + if (is_enabled(f_ag_rotate)) { // rotate forces back to the original frame cvm::rotation const rot_inv = rot.inverse(); @@ -1355,7 +1370,7 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force) } } - if ((b_center || b_rotate) && is_enabled(f_ag_fit_gradients)) { + if ((is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) && is_enabled(f_ag_fit_gradients)) { atom_group *group_for_fit = fitting_group ? fitting_group : this; @@ -1386,7 +1401,7 @@ void cvm::atom_group::apply_force(cvm::rvector const &force) return; } - if (b_rotate) { + if (is_enabled(f_ag_rotate)) { cvm::rotation const rot_inv = rot.inverse(); for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index af6a529f8a..490b008d06 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -101,18 +101,14 @@ public: inline void update_mass() { colvarproxy *p = cvm::proxy; - if (p->updated_masses()) { - mass = p->get_atom_mass(index); - } + mass = p->get_atom_mass(index); } /// Get the latest value of the charge inline void update_charge() { colvarproxy *p = cvm::proxy; - if (p->updated_charges()) { - charge = p->get_atom_charge(index); - } + charge = p->get_atom_charge(index); } /// Get the current position @@ -328,35 +324,23 @@ public: /// If yes, returns 1-based number of a common atom; else, returns 0 static int overlap(const atom_group &g1, const atom_group &g2); - /// \brief When updating atomic coordinates, translate them to align with the - /// center of mass of the reference coordinates - bool b_center; - - /// \brief When updating atom coordinates (and after - /// centering them if b_center is set), rotate the group to - /// align with the reference coordinates. - /// - /// Note: gradients will be calculated in the rotated frame: when - /// forces will be applied, they will rotated back to the original - /// frame - bool b_rotate; - /// The rotation calculated automatically if b_rotate is defined + /// The rotation calculated automatically if f_ag_rotate is defined cvm::rotation rot; - /// \brief Indicates that the user has explicitly set centerReference or + /// \brief Indicates that the user has explicitly set centerToReference or /// rotateReference, and the corresponding reference: /// cvc's (eg rmsd, eigenvector) will not override the user's choice bool b_user_defined_fit; - /// \brief use reference coordinates for b_center or b_rotate + /// \brief use reference coordinates for f_ag_center or f_ag_rotate std::vector ref_pos; /// \brief Center of geometry of the reference coordinates; regardless - /// of whether b_center is true, ref_pos is centered to zero at + /// of whether f_ag_center is true, ref_pos is centered to zero at /// initialization, and ref_pos_cog serves to center the positions cvm::atom_pos ref_pos_cog; - /// \brief If b_center or b_rotate is true, use this group to + /// \brief If f_ag_center or f_ag_rotate is true, use this group to /// define the transformation (default: this group itself) atom_group *fitting_group; @@ -395,12 +379,12 @@ public: void apply_translation(cvm::rvector const &t); /// \brief Get the current velocities; this must be called always - /// *after* read_positions(); if b_rotate is defined, the same + /// *after* read_positions(); if f_ag_rotate is defined, the same /// rotation applied to the coordinates will be used void read_velocities(); /// \brief Get the current total_forces; this must be called always - /// *after* read_positions(); if b_rotate is defined, the same + /// *after* read_positions(); if f_ag_rotate is defined, the same /// rotation applied to the coordinates will be used void read_total_forces(); diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index a84ec358dd..4893ceb46d 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -136,10 +136,11 @@ int colvarbias_abf::init(std::string const &conf) colvars[i]->enable(f_cv_hide_Jacobian); } - // If any colvar is extended-system, we need to collect the extended - // system gradient - if (colvars[i]->is_enabled(f_cv_extended_Lagrangian)) + // If any colvar is extended-system (restrained style, not external with constraint), we are running eABF + if (colvars[i]->is_enabled(f_cv_extended_Lagrangian) + && !colvars[i]->is_enabled(f_cv_external)) { b_extended = true; + } // Cannot mix and match coarse time steps with ABF because it gives // wrong total force averages - total force needs to be averaged over @@ -475,7 +476,7 @@ int colvarbias_abf::update() last_gradients->copy_grid(*gradients); last_samples->copy_grid(*samples); shared_last_step = cvm::step_absolute(); - cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+"."); + cvm::log("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".\n"); } // update UI estimator every step @@ -812,8 +813,10 @@ int colvarbias_abf::write_output_files() cvm::log("ABF bias trying to write gradients and samples to disk"); } - if (shared_on && cvm::main()->proxy->replica_index() > 0) { + if (shared_on && cvm::main()->proxy->replica_index() > 0 + && ! (b_CZAR_estimator || b_UI_estimator) ) { // No need to report the same data as replica 0, let it do the I/O job + // except if using an eABF FE estimator return COLVARS_OK; } diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index 15c81e9466..3a3120a058 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -121,6 +121,8 @@ private: } else { system_force[i] = colvars[i]->total_force().real_value - colvar_forces[i].real_value; + // If hideJacobian is active then total_force has an extra term of -fj + // which is the Jacobian-compensating force at the colvar level } if (cvm::debug()) cvm::log("ABF System force calc: cv " + cvm::to_str(i) + diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index 3efe0b0acc..92cc112845 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -129,14 +129,14 @@ int colvarbias_histogram::update() // output_prefix is unset during the constructor if (cvm::step_relative() == 0) { out_name = cvm::output_prefix() + "." + this->name + ".dat"; - cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\""); + cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\"\n"); } } if (out_name_dx.size() == 0) { if (cvm::step_relative() == 0) { out_name_dx = cvm::output_prefix() + "." + this->name + ".dx"; - cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\""); + cvm::log("Histogram " + this->name + " will be written to file \"" + out_name_dx + "\"\n"); } } diff --git a/lib/colvars/colvarbias_histogram.h b/lib/colvars/colvarbias_histogram.h index 2e6c6884fb..6044fa189e 100644 --- a/lib/colvars/colvarbias_histogram.h +++ b/lib/colvars/colvarbias_histogram.h @@ -35,7 +35,6 @@ protected: colvar_grid_scalar *grid; std::vector bin; std::string out_name, out_name_dx; - size_t output_freq; /// If one or more of the variables are \link colvarvalue::type_vector \endlink, treat them as arrays of this length size_t colvar_array_size; diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index d448b1f2e4..ccbd6c406f 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -50,6 +50,7 @@ colvarbias_meta::colvarbias_meta(char const *key) dump_fes = true; keep_hills = false; + restart_keep_hills = false; dump_fes_save = false; dump_replica_fes = false; @@ -902,7 +903,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, count++) { size_t i; for (i = 0; i < num_variables(); i++) { - new_colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i); + new_colvar_values[i] = he->bin_to_value_scalar(he_ix[i], i); } // loop over the hills and increment the energy grid locally @@ -934,40 +935,12 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, } } else { - - // TODO delete this (never used) - - // simpler version, with just the energy - - for ( ; (he->index_ok(he_ix)); ) { - - for (size_t i = 0; i < num_variables(); i++) { - new_colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i); - } - - hills_energy_here = 0.0; - calc_hills(h_first, h_last, hills_energy_here, &new_colvar_values); - he->acc_value(he_ix, hills_energy_here); - - he->incr(he_ix); - - count++; - if ((count % print_frequency) == 0) { - if (print_progress) { - cvm::real const progress = cvm::real(count) / cvm::real(he->number_of_points()); - std::ostringstream os; - os.setf(std::ios::fixed, std::ios::floatfield); - os << std::setw(6) << std::setprecision(2) - << 100.0 * progress - << "% done."; - cvm::log(os.str()); - } - } - } + cvm::error("No grid object provided in metadynamics::project_hills()\n", + BUG_ERROR); } if (print_progress) { - cvm::log("100.00% done."); + cvm::log("100.00% done.\n"); } if (! keep_hills) { @@ -1281,9 +1254,25 @@ int colvarbias_meta::set_state_params(std::string const &state_conf) return error_code; } + colvarparse::get_keyval(state_conf, "keepHills", restart_keep_hills, false, + colvarparse::parse_restart); + + if ((!restart_keep_hills) && (cvm::main()->restart_version_number() < 20210604)) { + if (keep_hills) { + cvm::log("Warning: could not ensure that keepHills was enabled when " + "this state file was written; because it is enabled now, " + "it is assumed that it was also then, but please verify.\n"); + restart_keep_hills = true; + } + } else { + if (restart_keep_hills) { + cvm::log("This state file/stream contains explicit hills.\n"); + } + } + std::string check_replica = ""; if (colvarparse::get_keyval(state_conf, "replicaID", check_replica, - std::string(""), colvarparse::parse_silent) && + std::string(""), colvarparse::parse_restart) && (check_replica != this->replica_id)) { return cvm::error("Error: in the state file , the " "\"metadynamics\" block has a different replicaID ("+ @@ -1297,8 +1286,6 @@ int colvarbias_meta::set_state_params(std::string const &state_conf) std::istream & colvarbias_meta::read_state_data(std::istream& is) { - bool grids_from_restart_file = use_grids; - if (use_grids) { if (expand_grids) { @@ -1343,7 +1330,6 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) !(hills_energy->read_restart(is))) { is.clear(); is.seekg(hills_energy_pos, std::ios::beg); - grids_from_restart_file = false; if (!rebin_grids) { if (hills_energy_backup == NULL) cvm::fatal_error("Error: couldn't read the free energy grid for metadynamics bias \""+ @@ -1381,7 +1367,6 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) !(hills_energy_gradients->read_restart(is))) { is.clear(); is.seekg(hills_energy_gradients_pos, std::ios::beg); - grids_from_restart_file = false; if (!rebin_grids) { if (hills_energy_backup == NULL) cvm::fatal_error("Error: couldn't read the free energy gradients grid for metadynamics bias \""+ @@ -1419,28 +1404,29 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) } } + // Save references to the end of the list of existing hills, so that it can + // be cleared if hills are read successfully state bool const existing_hills = !hills.empty(); size_t const old_hills_size = hills.size(); hill_iter old_hills_end = hills.end(); hill_iter old_hills_off_grid_end = hills_off_grid.end(); - // read the hills explicitly written (if there are any) + // Read any hills following the grid data (if any) while (read_hill(is)) { - if (cvm::debug()) + if (cvm::debug()) { cvm::log("Read a previously saved hill under the " "metadynamics bias \""+ this->name+"\", created at step "+ cvm::to_str((hills.back()).it)+".\n"); + } } is.clear(); new_hills_begin = hills.end(); - if (grids_from_restart_file) { - if (hills.size() > old_hills_size) - cvm::log("Read "+cvm::to_str(hills.size())+ - " hills in addition to the grids.\n"); - } else { - if (!hills.empty()) - cvm::log("Read "+cvm::to_str(hills.size())+" hills.\n"); + cvm::log("Read "+cvm::to_str(hills.size() - old_hills_size)+" hills.\n"); + + if (existing_hills) { + hills.erase(hills.begin(), old_hills_end); + hills_off_grid.erase(hills_off_grid.begin(), old_hills_off_grid_end); } if (rebin_grids) { @@ -1454,7 +1440,16 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) colvar_grid_gradient *new_hills_energy_gradients = new colvar_grid_gradient(colvars); - if (!grids_from_restart_file || (keep_hills && !hills.empty())) { + if (cvm::debug()) { + std::ostringstream tmp_os; + tmp_os << "hills_energy parameters:\n"; + hills_energy->write_params(tmp_os); + tmp_os << "new_hills_energy parameters:\n"; + new_hills_energy->write_params(tmp_os); + cvm::log(tmp_os.str()); + } + + if (restart_keep_hills && !hills.empty()) { // if there are hills, recompute the new grids from them cvm::log("Rebinning the energy and forces grids from "+ cvm::to_str(hills.size())+" hills (this may take a while)...\n"); @@ -1494,11 +1489,6 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) if (cvm::debug()) cvm::log("colvarbias_meta::read_restart() done\n"); - if (existing_hills) { - hills.erase(hills.begin(), old_hills_end); - hills_off_grid.erase(hills_off_grid.begin(), old_hills_off_grid_end); - } - has_data = true; if (comm != single_replica) { @@ -1509,6 +1499,15 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) } +inline std::istream & reset_istream(std::istream &is, size_t start_pos) +{ + is.clear(); + is.seekg(start_pos, std::ios::beg); + is.setstate(std::ios::failbit); + return is; +} + + std::istream & colvarbias_meta::read_hill(std::istream &is) { if (!is) return is; // do nothing if failbit is set @@ -1518,45 +1517,72 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) std::string data; if ( !(is >> read_block("hill", &data)) ) { - is.clear(); - is.seekg(start_pos, std::ios::beg); - is.setstate(std::ios::failbit); - return is; + return reset_istream(is, start_pos); } + std::istringstream data_is(data); + cvm::step_number h_it = 0L; - get_keyval(data, "step", h_it, h_it, parse_restart); - if (h_it <= state_file_step) { - if (cvm::debug()) - cvm::log("Skipping a hill older than the state file for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); - return is; - } - cvm::real h_weight; - get_keyval(data, "weight", h_weight, hill_weight, parse_restart); - std::vector h_centers(num_variables()); for (i = 0; i < num_variables(); i++) { h_centers[i].type(variables(i)->value()); } - get_keyval(data, "centers", h_centers, h_centers, parse_restart); - std::vector h_sigmas(num_variables()); - get_keyval(data, "widths", h_sigmas, h_sigmas, parse_restart); - for (i = 0; i < num_variables(); i++) { - // For backward compatibility, read the widths instead of the sigmas - h_sigmas[i] /= 2.0; - } + std::string h_replica; - std::string h_replica = ""; - if (comm != single_replica) { - get_keyval(data, "replicaID", h_replica, replica_id, parse_restart); - if (h_replica != replica_id) - cvm::fatal_error("Error: trying to read a hill created by replica \""+h_replica+ - "\" for replica \""+replica_id+ - "\"; did you swap output files?\n"); + std::string keyword; + while (data_is >> keyword) { + + if (keyword == "step") { + if ( !(data_is >> h_it)) { + return reset_istream(is, start_pos); + } + if ((h_it <= state_file_step) && !restart_keep_hills) { + if (cvm::debug()) + cvm::log("Skipping a hill older than the state file for metadynamics bias \""+ + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); + return is; + } + } + + if (keyword == "weight") { + if ( !(data_is >> h_weight)) { + return reset_istream(is, start_pos); + } + } + + if (keyword == "centers") { + for (i = 0; i < num_variables(); i++) { + if ( !(data_is >> h_centers[i])) { + return reset_istream(is, start_pos); + } + } + } + + if (keyword == "widths") { + for (i = 0; i < num_variables(); i++) { + if ( !(data_is >> h_sigmas[i])) { + return reset_istream(is, start_pos); + } + // For backward compatibility, read the widths instead of the sigmas + h_sigmas[i] /= 2.0; + } + } + + if (comm != single_replica) { + if (keyword == "replicaID") { + if ( !(data_is >> h_replica)) { + return reset_istream(is, start_pos); + } + if (h_replica != replica_id) { + cvm::error("Error: trying to read a hill created by replica \""+ + h_replica+"\" for replica \""+replica_id+ + "\"; did you swap output files?\n", INPUT_ERROR); + } + } + } } hill_iter const hills_end = hills.end(); @@ -1699,8 +1725,12 @@ std::string const colvarbias_meta::hills_traj_file_name() const std::string const colvarbias_meta::get_state_params() const { std::ostringstream os; - if (this->comm != single_replica) + if (keep_hills) { + os << "keepHills on" << "\n"; + } + if (this->comm != single_replica) { os << "replicaID " << this->replica_id << "\n"; + } return (colvarbias::get_state_params() + os.str()); } @@ -1944,6 +1974,7 @@ colvarbias_meta::hill::hill(cvm::step_number it_in, sigmas(cv_values.size()), replica(replica_in) { + hill_value = 0.0; for (size_t i = 0; i < cv_values.size(); i++) { centers[i].type(cv_values[i]); centers[i] = cv_values[i]; @@ -1967,7 +1998,9 @@ colvarbias_meta::hill::hill(colvarbias_meta::hill const &h) centers(h.centers), sigmas(h.sigmas), replica(h.replica) -{} +{ + hill_value = 0.0; +} colvarbias_meta::hill::~hill() diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 2f34abcc77..bf1710170e 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -155,10 +155,12 @@ protected: /// \brief How often the hills should be projected onto the grids size_t grids_freq; - /// \brief Whether to keep the hills in the restart file (e.g. to do - /// meaningful accurate rebinning afterwards) + /// Keep hills in the restart file (e.g. to accurately rebin later) bool keep_hills; + /// value of keepHills saved in the most recent restart file + bool restart_keep_hills; + /// \brief Dump the free energy surface (.pmf file) every restartFrequency bool dump_fes; diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index a1542d05bf..8938e64a0e 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -51,7 +51,7 @@ int colvar::cvc::init(std::string const &conf) std::string const old_name(name); if (name.size() > 0) { - cvm::log("Updating configuration for component \""+name+"\""); + cvm::log("Updating configuration for component \""+name+"\"\n"); } if (get_keyval(conf, "name", name, name)) { @@ -112,7 +112,7 @@ int colvar::cvc::init_total_force_params(std::string const &conf) } if (get_keyval_feature(this, conf, "oneSiteTotalForce", f_cvc_one_site_total_force, is_enabled(f_cvc_one_site_total_force))) { - cvm::log("Computing total force on group 1 only"); + cvm::log("Computing total force on group 1 only\n"); } if (! is_enabled(f_cvc_one_site_total_force)) { @@ -426,7 +426,7 @@ void colvar::cvc::collect_gradients(std::vector const &atom_ids, std::vecto // If necessary, apply inverse rotation to get atomic // gradient in the laboratory frame - if (ag.b_rotate) { + if (ag.is_enabled(f_ag_rotate)) { cvm::rotation const rot_inv = ag.rot.inverse(); for (size_t k = 0; k < ag.size(); k++) { @@ -505,7 +505,7 @@ void colvar::cvc::debug_gradients() cvm::atom_pos fit_gradient_sum, gradient_sum; // print the values of the fit gradients - if (group->b_rotate || group->b_center) { + if (group->is_enabled(f_ag_center) || group->is_enabled(f_ag_rotate)) { if (group->is_enabled(f_ag_fit_gradients)) { size_t j; @@ -514,7 +514,7 @@ void colvar::cvc::debug_gradients() for (j = 0; j < group_for_fit->fit_gradients.size(); j++) { cvm::log((group->fitting_group ? std::string("refPosGroup") : group->key) + "[" + cvm::to_str(j) + "] = " + - (group->b_rotate ? + (group->is_enabled(f_ag_rotate) ? cvm::to_str(rot_0.rotate(group_for_fit->fit_gradients[j])) : cvm::to_str(group_for_fit->fit_gradients[j]))); } @@ -525,7 +525,7 @@ void colvar::cvc::debug_gradients() for (size_t ia = 0; ia < group->size(); ia++) { // tests are best conducted in the unrotated (simulation) frame - cvm::rvector const atom_grad = (group->b_rotate ? + cvm::rvector const atom_grad = (group->is_enabled(f_ag_rotate) ? rot_inv.rotate((*group)[ia].grad) : (*group)[ia].grad); gradient_sum += atom_grad; diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 57df68ab6e..2f68ac6abc 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -271,6 +271,12 @@ public: /// \brief Whether or not this CVC will be computed in parallel whenever possible bool b_try_scalable; + /// Forcibly set value of CVC - useful for driving an external coordinate, + /// eg. lambda dynamics + inline void set_value(colvarvalue const &new_value) { + x = new_value; + } + protected: /// \brief Cached value @@ -1344,6 +1350,71 @@ public: }; +class colvar::euler_phi + : public colvar::orientation +{ +public: + euler_phi(std::string const &conf); + euler_phi(); + virtual int init(std::string const &conf); + virtual ~euler_phi() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to handle the 2*PI periodicity + virtual void wrap(colvarvalue &x_unwrapped) const; +}; + + +class colvar::euler_psi + : public colvar::orientation +{ +public: + euler_psi(std::string const &conf); + euler_psi(); + virtual int init(std::string const &conf); + virtual ~euler_psi() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + /// Redefined to handle the 2*PI periodicity + virtual void wrap(colvarvalue &x_unwrapped) const; +}; + + +class colvar::euler_theta + : public colvar::orientation +{ +public: + euler_theta(std::string const &conf); + euler_theta(); + virtual int init(std::string const &conf); + virtual ~euler_theta() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); + // theta angle is a scalar variable and not periodic + // we need to override the virtual functions from orientation + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; +}; + /// \brief Colvar component: root mean square deviation (RMSD) of a /// group with respect to a set of reference coordinates; uses \link @@ -1406,6 +1477,28 @@ public: }; +// \brief Colvar component: alch_lambda +// To communicate value with back-end in lambda-dynamics +class colvar::alch_lambda + : public colvar::cvc +{ +protected: + // No atom groups needed +public: + alch_lambda(std::string const &conf); + alch_lambda(); + virtual ~alch_lambda() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); + virtual cvm::real dist2(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const; +}; + class colvar::componentDisabled : public colvar::cvc @@ -1681,11 +1774,20 @@ public: protected: - /// Identifier of the map object (as used by the simulation engine) - std::string map_name; + /// String identifier of the map object (as used by the simulation engine) + std::string volmap_name; + + /// Numeric identifier of the map object (as used by the simulation engine) + int volmap_id; /// Index of the map objet in the proxy arrays int volmap_index; + + /// Group of atoms selected internally (optional) + cvm::atom_group *atoms; + + /// Weights assigned to each atom (default: uniform weights) + std::vector atom_weights; }; diff --git a/lib/colvars/colvarcomp_alchlambda.cpp b/lib/colvars/colvarcomp_alchlambda.cpp new file mode 100644 index 0000000000..5e10a1dab5 --- /dev/null +++ b/lib/colvars/colvarcomp_alchlambda.cpp @@ -0,0 +1,56 @@ +// -*- c++ -*- + +// This file is part of the Collective Variables module (Colvars). +// The original version of Colvars and its updates are located at: +// https://github.com/Colvars/colvars +// Please update all Colvars source files before making any changes. +// If you wish to distribute your changes, please submit them to the +// Colvars repository at GitHub. + + +#include + +#include "colvarmodule.h" +#include "colvarvalue.h" +#include "colvarparse.h" +#include "colvar.h" +#include "colvarcomp.h" + + +colvar::alch_lambda::alch_lambda(std::string const &conf) + : cvc(conf) +{ + function_type = "alch_lambda"; + + disable(f_cvc_explicit_gradient); + disable(f_cvc_gradient); + + x.type(colvarvalue::type_scalar); + // Query initial value from back-end + cvm::proxy->get_alch_lambda(&x.real_value); +} + + +void colvar::alch_lambda::calc_value() +{ + // Special workflow: + // at the beginning of the timestep we get a force instead of calculating the value + + cvm::proxy->get_dE_dLambda(&ft.real_value); + ft.real_value *= -1.0; // Energy derivative to force +} + + +void colvar::alch_lambda::calc_gradients() +{ +} + + +void colvar::alch_lambda::apply_force(colvarvalue const &force) +{ + // Special workflow: + // at the end of the time step we send a new value + cvm::proxy->set_alch_lambda(&x.real_value); +} + +simple_scalar_dist_functions(alch_lambda) diff --git a/lib/colvars/colvarcomp_apath.cpp b/lib/colvars/colvarcomp_apath.cpp index 591d8d0012..a4d1df9e50 100644 --- a/lib/colvars/colvarcomp_apath.cpp +++ b/lib/colvars/colvarcomp_apath.cpp @@ -19,17 +19,8 @@ colvar::aspathCV::aspathCV(std::string const &conf): CVBasedPath(conf) { get_keyval(conf, "weights", p_weights, std::vector(cv.size(), 1.0)); x.type(colvarvalue::type_scalar); use_explicit_gradients = true; - std::vector rmsd_between_refs(total_reference_frames - 1, 0.0); - computeDistanceBetweenReferenceFrames(rmsd_between_refs); - cvm::real mean_square_displacements = 0.0; - for (size_t i_frame = 1; i_frame < total_reference_frames; ++i_frame) { - cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n")); - mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1]; - } - mean_square_displacements /= cvm::real(total_reference_frames - 1); - cvm::real suggested_lambda = 1.0 / mean_square_displacements; cvm::real p_lambda; - get_keyval(conf, "lambda", p_lambda, suggested_lambda); + get_keyval(conf, "lambda", p_lambda, -1.0); ArithmeticPathCV::ArithmeticPathBase::initialize(cv.size(), total_reference_frames, p_lambda, ref_cv[0], p_weights); cvm::log(std::string("Lambda is ") + cvm::to_str(lambda) + std::string("\n")); for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { @@ -58,6 +49,16 @@ void colvar::aspathCV::updateDistanceToReferenceFrames() { } void colvar::aspathCV::calc_value() { + if (lambda < 0) { + // this implies that the user may not set a valid lambda value + // so recompute it by the suggested value in Parrinello's paper + cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n"); + cvm::log("This component (aspathCV) will recompute a value for lambda following the suggestion in the origin paper.\n"); + std::vector rmsd_between_refs(total_reference_frames - 1, 0.0); + computeDistanceBetweenReferenceFrames(rmsd_between_refs); + reComputeLambda(rmsd_between_refs); + cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(lambda)); + } computeValue(); x = s; } @@ -107,17 +108,8 @@ colvar::azpathCV::azpathCV(std::string const &conf): CVBasedPath(conf) { get_keyval(conf, "weights", p_weights, std::vector(cv.size(), 1.0)); x.type(colvarvalue::type_scalar); use_explicit_gradients = true; - std::vector rmsd_between_refs(total_reference_frames - 1, 0.0); - computeDistanceBetweenReferenceFrames(rmsd_between_refs); - cvm::real mean_square_displacements = 0.0; - for (size_t i_frame = 1; i_frame < total_reference_frames; ++i_frame) { - cvm::log(std::string("Distance between frame ") + cvm::to_str(i_frame) + " and " + cvm::to_str(i_frame + 1) + " is " + cvm::to_str(rmsd_between_refs[i_frame - 1]) + std::string("\n")); - mean_square_displacements += rmsd_between_refs[i_frame - 1] * rmsd_between_refs[i_frame - 1]; - } - mean_square_displacements /= cvm::real(total_reference_frames - 1); - cvm::real suggested_lambda = 1.0 / mean_square_displacements; cvm::real p_lambda; - get_keyval(conf, "lambda", p_lambda, suggested_lambda); + get_keyval(conf, "lambda", p_lambda, -1.0); ArithmeticPathCV::ArithmeticPathBase::initialize(cv.size(), total_reference_frames, p_lambda, ref_cv[0], p_weights); cvm::log(std::string("Lambda is ") + cvm::to_str(lambda) + std::string("\n")); for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { @@ -146,6 +138,16 @@ void colvar::azpathCV::updateDistanceToReferenceFrames() { } void colvar::azpathCV::calc_value() { + if (lambda < 0) { + // this implies that the user may not set a valid lambda value + // so recompute it by the suggested value in Parrinello's paper + cvm::log("A non-positive value of lambda is detected, which implies that it may not set in the configuration.\n"); + cvm::log("This component (azpathCV) will recompute a value for lambda following the suggestion in the origin paper.\n"); + std::vector rmsd_between_refs(total_reference_frames - 1, 0.0); + computeDistanceBetweenReferenceFrames(rmsd_between_refs); + reComputeLambda(rmsd_between_refs); + cvm::log("Ok, the value of lambda is updated to " + cvm::to_str(lambda)); + } computeValue(); x = z; } diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 1b1b295ff1..c044da86b3 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -194,10 +194,10 @@ colvar::distance_z::distance_z(std::string const &conf) ref2 = parse_group(conf, "ref2", true); if ( ref2 ) { - cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\""); + cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"\n"); fixed_axis = false; if (key_lookup(conf, "axis")) - cvm::log("Warning: explicit axis definition will be ignored!"); + cvm::log("Warning: explicit axis definition will be ignored!\n"); } else { if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) { if (axis.norm2() == 0.0) { @@ -808,9 +808,9 @@ colvar::gyration::gyration(std::string const &conf) atoms = parse_group(conf, "atoms"); if (atoms->b_user_defined_fit) { - cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); + cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); } else { - atoms->b_center = true; + atoms->enable(f_ag_center); atoms->ref_pos.assign(1, cvm::atom_pos(0.0, 0.0, 0.0)); atoms->fit_gradients.assign(atoms->size(), cvm::rvector(0.0, 0.0, 0.0)); } @@ -1025,13 +1025,13 @@ colvar::rmsd::rmsd(std::string const &conf) } if (atoms->b_user_defined_fit) { - cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); + cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); } else { // Default: fit everything - cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating it as a variable: " + cvm::log("Enabling \"centerToReference\" and \"rotateToReference\", to minimize RMSD before calculating it as a variable: " "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); - atoms->b_center = true; - atoms->b_rotate = true; + atoms->enable(f_ag_center); + atoms->enable(f_ag_rotate); // default case: reference positions for calculating the rmsd are also those used // for fitting atoms->ref_pos = ref_pos; @@ -1156,7 +1156,7 @@ void colvar::rmsd::calc_Jacobian_derivative() cvm::real rotation_term = 0.0; // The rotation term only applies is coordinates are rotated - if (atoms->b_rotate) { + if (atoms->is_enabled(f_ag_rotate)) { // gradient of the rotation matrix cvm::matrix2d grad_rot_mat(3, 3); @@ -1202,7 +1202,7 @@ void colvar::rmsd::calc_Jacobian_derivative() } // The translation term only applies is coordinates are centered - cvm::real translation_term = atoms->b_center ? 3.0 : 0.0; + cvm::real translation_term = atoms->is_enabled(f_ag_center) ? 3.0 : 0.0; jd.real_value = x.real_value > 0.0 ? (3.0 * atoms->size() - 1.0 - translation_term - rotation_term) / x.real_value : @@ -1284,10 +1284,10 @@ colvar::eigenvector::eigenvector(std::string const &conf) cvm::log("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n"); } else { // default: fit everything - cvm::log("Enabling \"centerReference\" and \"rotateReference\", to minimize RMSD before calculating the vector projection: " + cvm::log("Enabling \"centerToReference\" and \"rotateToReference\", to minimize RMSD before calculating the vector projection: " "if this is not the desired behavior, disable them explicitly within the \"atoms\" block.\n"); - atoms->b_center = true; - atoms->b_rotate = true; + atoms->enable(f_ag_center); + atoms->enable(f_ag_rotate); atoms->ref_pos = ref_pos; atoms->center_ref_pos(); atoms->disable(f_ag_fit_gradients); // cancel out if group is fitted on itself @@ -1355,14 +1355,14 @@ colvar::eigenvector::eigenvector(std::string const &conf) if (b_difference_vector) { - if (atoms->b_center) { + if (atoms->is_enabled(f_ag_center)) { // both sets should be centered on the origin for fitting for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] -= eig_center; ref_pos[i] -= ref_pos_center; } } - if (atoms->b_rotate) { + if (atoms->is_enabled(f_ag_rotate)) { atoms->rot.calc_optimal_rotation(eigenvec, ref_pos); for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] = atoms->rot.rotate(eigenvec[i]); @@ -1372,7 +1372,7 @@ colvar::eigenvector::eigenvector(std::string const &conf) for (size_t i = 0; i < atoms->size(); i++) { eigenvec[i] -= ref_pos[i]; } - if (atoms->b_center) { + if (atoms->is_enabled(f_ag_center)) { // bring back the ref positions to where they were for (size_t i = 0; i < atoms->size(); i++) { ref_pos[i] += ref_pos_center; @@ -1521,7 +1521,8 @@ colvar::cartesian::cartesian(std::string const &conf) x.type(colvarvalue::type_vector); disable(f_cvc_explicit_gradient); - x.vector1d_value.resize(atoms->size() * axes.size()); + // Don't try to access atoms if creation of the atom group failed + if (atoms != NULL) x.vector1d_value.resize(atoms->size() * axes.size()); } diff --git a/lib/colvars/colvarcomp_gpath.cpp b/lib/colvars/colvarcomp_gpath.cpp index 8849f9dbe4..884dfdaee9 100644 --- a/lib/colvars/colvarcomp_gpath.cpp +++ b/lib/colvars/colvarcomp_gpath.cpp @@ -65,8 +65,8 @@ colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(con cvm::atom_group* tmp_atoms = parse_group(conf, "atoms"); if (!has_user_defined_fitting) { // Swipe from the rmsd class - tmp_atoms->b_center = true; - tmp_atoms->b_rotate = true; + tmp_atoms->enable(f_ag_center); + tmp_atoms->enable(f_ag_rotate); tmp_atoms->ref_pos = reference_frames[i_frame]; tmp_atoms->center_ref_pos(); tmp_atoms->enable(f_ag_fit_gradients); @@ -87,8 +87,8 @@ colvar::CartesianBasedPath::CartesianBasedPath(std::string const &conf): cvc(con std::vector reference_fitting_position(tmp_fitting_atoms->size()); cvm::load_coords(reference_position_filename.c_str(), &reference_fitting_position, tmp_fitting_atoms, reference_column, reference_column_value); // setup the atom group for calculating - tmp_atoms->b_center = true; - tmp_atoms->b_rotate = true; + tmp_atoms->enable(f_ag_center); + tmp_atoms->enable(f_ag_rotate); tmp_atoms->b_user_defined_fit = true; tmp_atoms->disable(f_ag_scalable); tmp_atoms->disable(f_ag_scalable_com); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index c7a414b006..fbf572b3f3 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -442,3 +442,335 @@ void colvar::spin_angle::wrap(colvarvalue &x_unwrapped) const return; } + + +colvar::euler_phi::euler_phi(std::string const &conf) + : orientation() +{ + function_type = "euler_phi"; + period = 360.0; + enable(f_cvc_periodic); + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); + init(conf); +} + + +colvar::euler_phi::euler_phi() + : orientation() +{ + function_type = "euler_phi"; + period = 360.0; + enable(f_cvc_periodic); + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); +} + + +int colvar::euler_phi::init(std::string const &conf) +{ + int error_code = COLVARS_OK; + error_code |= orientation::init(conf); + return error_code; +} + + +void colvar::euler_phi::calc_value() +{ + atoms_cog = atoms->center_of_geometry(); + + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); + + const cvm::real& q0 = rot.q.q0; + const cvm::real& q1 = rot.q.q1; + const cvm::real& q2 = rot.q.q2; + const cvm::real& q3 = rot.q.q3; + const cvm::real tmp_y = 2 * (q0 * q1 + q2 * q3); + const cvm::real tmp_x = 1 - 2 * (q1 * q1 + q2 * q2); + x.real_value = cvm::atan2(tmp_y, tmp_x) * (180.0/PI); +} + + +void colvar::euler_phi::calc_gradients() +{ + const cvm::real& q0 = rot.q.q0; + const cvm::real& q1 = rot.q.q1; + const cvm::real& q2 = rot.q.q2; + const cvm::real& q3 = rot.q.q3; + const cvm::real denominator = (2 * q0 * q1 + 2 * q2 * q3) * (2 * q0 * q1 + 2 * q2 * q3) + (-2 * q1 * q1 - 2 * q2 * q2 + 1) * (-2 * q1 * q1 - 2 * q2 * q2 + 1); + const cvm::real dxdq0 = (180.0/PI) * 2 * q1 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) / denominator; + const cvm::real dxdq1 = (180.0/PI) * (2 * q0 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) - 4 * q1 * (-2 * q0 * q1 - 2 * q2 * q3)) / denominator; + const cvm::real dxdq2 = (180.0/PI) * (-4 * q2 * (-2 * q0 * q1 - 2 * q2 * q3) + 2 * q3 * (-2 * q1 * q1 - 2 * q2 * q2 + 1)) / denominator; + const cvm::real dxdq3 = (180.0/PI) * 2 * q2 * (-2 * q1 * q1 - 2 * q2 * q2 + 1) / denominator; + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) + + (dxdq1 * (rot.dQ0_2[ia])[1]) + + (dxdq2 * (rot.dQ0_2[ia])[2]) + + (dxdq3 * (rot.dQ0_2[ia])[3]); + } +} + + +void colvar::euler_phi::apply_force(colvarvalue const &force) +{ + cvm::real const &fw = force.real_value; + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); + } +} + + +cvm::real colvar::euler_phi::dist2(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); + return diff * diff; +} + + +colvarvalue colvar::euler_phi::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); + return 2.0 * diff; +} + + +colvarvalue colvar::euler_phi::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); + return (-2.0) * diff; +} + + +void colvar::euler_phi::wrap(colvarvalue &x_unwrapped) const +{ + if ((x_unwrapped.real_value - wrap_center) >= 180.0) { + x_unwrapped.real_value -= 360.0; + return; + } + + if ((x_unwrapped.real_value - wrap_center) < -180.0) { + x_unwrapped.real_value += 360.0; + return; + } + + return; +} + + +colvar::euler_psi::euler_psi(std::string const &conf) + : orientation() +{ + function_type = "euler_psi"; + period = 360.0; + enable(f_cvc_periodic); + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); + init(conf); +} + + +colvar::euler_psi::euler_psi() + : orientation() +{ + function_type = "euler_psi"; + period = 360.0; + enable(f_cvc_periodic); + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); +} + + +int colvar::euler_psi::init(std::string const &conf) +{ + int error_code = COLVARS_OK; + error_code |= orientation::init(conf); + return error_code; +} + + +void colvar::euler_psi::calc_value() +{ + atoms_cog = atoms->center_of_geometry(); + + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); + + const cvm::real& q0 = rot.q.q0; + const cvm::real& q1 = rot.q.q1; + const cvm::real& q2 = rot.q.q2; + const cvm::real& q3 = rot.q.q3; + const cvm::real tmp_y = 2 * (q0 * q3 + q1 * q2); + const cvm::real tmp_x = 1 - 2 * (q2 * q2 + q3 * q3); + x.real_value = cvm::atan2(tmp_y, tmp_x) * (180.0/PI); +} + + +void colvar::euler_psi::calc_gradients() +{ + const cvm::real& q0 = rot.q.q0; + const cvm::real& q1 = rot.q.q1; + const cvm::real& q2 = rot.q.q2; + const cvm::real& q3 = rot.q.q3; + const cvm::real denominator = (2 * q0 * q3 + 2 * q1 * q2) * (2 * q0 * q3 + 2 * q1 * q2) + (-2 * q2 * q2 - 2 * q3 * q3 + 1) * (-2 * q2 * q2 - 2 * q3 * q3 + 1); + const cvm::real dxdq0 = (180.0/PI) * 2 * q3 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) / denominator; + const cvm::real dxdq1 = (180.0/PI) * 2 * q2 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) / denominator; + const cvm::real dxdq2 = (180.0/PI) * (2 * q1 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) - 4 * q2 * (-2 * q0 * q3 - 2 * q1 * q2)) / denominator; + const cvm::real dxdq3 = (180.0/PI) * (2 * q0 * (-2 * q2 * q2 - 2 * q3 * q3 + 1) - 4 * q3 * (-2 * q0 * q3 - 2 * q1 * q2)) / denominator; + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) + + (dxdq1 * (rot.dQ0_2[ia])[1]) + + (dxdq2 * (rot.dQ0_2[ia])[2]) + + (dxdq3 * (rot.dQ0_2[ia])[3]); + } +} + + +void colvar::euler_psi::apply_force(colvarvalue const &force) +{ + cvm::real const &fw = force.real_value; + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); + } +} + + +cvm::real colvar::euler_psi::dist2(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); + return diff * diff; +} + + +colvarvalue colvar::euler_psi::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); + return 2.0 * diff; +} + + +colvarvalue colvar::euler_psi::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + cvm::real diff = x1.real_value - x2.real_value; + diff = (diff < -180.0 ? diff + 360.0 : (diff > 180.0 ? diff - 360.0 : diff)); + return (-2.0) * diff; +} + + +void colvar::euler_psi::wrap(colvarvalue &x_unwrapped) const +{ + if ((x_unwrapped.real_value - wrap_center) >= 180.0) { + x_unwrapped.real_value -= 360.0; + return; + } + + if ((x_unwrapped.real_value - wrap_center) < -180.0) { + x_unwrapped.real_value += 360.0; + return; + } + + return; +} + + +colvar::euler_theta::euler_theta(std::string const &conf) + : orientation() +{ + function_type = "euler_theta"; + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); + init(conf); +} + + +colvar::euler_theta::euler_theta() + : orientation() +{ + function_type = "euler_theta"; + enable(f_cvc_explicit_gradient); + x.type(colvarvalue::type_scalar); +} + + +int colvar::euler_theta::init(std::string const &conf) +{ + int error_code = COLVARS_OK; + error_code |= orientation::init(conf); + return error_code; +} + + +void colvar::euler_theta::calc_value() +{ + atoms_cog = atoms->center_of_geometry(); + + rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog)); + + const cvm::real& q0 = rot.q.q0; + const cvm::real& q1 = rot.q.q1; + const cvm::real& q2 = rot.q.q2; + const cvm::real& q3 = rot.q.q3; + x.real_value = cvm::asin(2 * (q0 * q2 - q3 * q1)) * (180.0/PI); +} + + +void colvar::euler_theta::calc_gradients() +{ + const cvm::real& q0 = rot.q.q0; + const cvm::real& q1 = rot.q.q1; + const cvm::real& q2 = rot.q.q2; + const cvm::real& q3 = rot.q.q3; + const cvm::real denominator = cvm::sqrt(1 - (2 * q0 * q2 - 2 * q1 * q3) * (2 * q0 * q2 - 2 * q1 * q3)); + const cvm::real dxdq0 = (180.0/PI) * 2 * q2 / denominator; + const cvm::real dxdq1 = (180.0/PI) * -2 * q3 / denominator; + const cvm::real dxdq2 = (180.0/PI) * 2 * q0 / denominator; + const cvm::real dxdq3 = (180.0/PI) * -2 * q1 / denominator; + for (size_t ia = 0; ia < atoms->size(); ia++) { + (*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]) + + (dxdq1 * (rot.dQ0_2[ia])[1]) + + (dxdq2 * (rot.dQ0_2[ia])[2]) + + (dxdq3 * (rot.dQ0_2[ia])[3]); + } +} + + +void colvar::euler_theta::apply_force(colvarvalue const &force) +{ + cvm::real const &fw = force.real_value; + if (!atoms->noforce) { + atoms->apply_colvar_force(fw); + } +} + + +cvm::real colvar::euler_theta::dist2(colvarvalue const &x1, + colvarvalue const &x2) const +{ + // theta angle is not periodic + return cvc::dist2(x1, x2); +} + + +colvarvalue colvar::euler_theta::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + // theta angle is not periodic + return cvc::dist2_lgrad(x1, x2); +} + + +colvarvalue colvar::euler_theta::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + // theta angle is not periodic + return cvc::dist2_rgrad(x1, x2); +} diff --git a/lib/colvars/colvarcomp_volmaps.cpp b/lib/colvars/colvarcomp_volmaps.cpp index 8ea2b76de4..5f80e46717 100644 --- a/lib/colvars/colvarcomp_volmaps.cpp +++ b/lib/colvars/colvarcomp_volmaps.cpp @@ -19,6 +19,9 @@ colvar::map_total::map_total() : cvc(), volmap_index(-1) { function_type = "map_total"; + volmap_id = -1; + volmap_index = -1; + atoms = NULL; x.type(colvarvalue::type_scalar); } @@ -27,6 +30,9 @@ colvar::map_total::map_total(std::string const &conf) : cvc(), volmap_index(-1) { function_type = "map_total"; + volmap_id = -1; + volmap_index = -1; + atoms = NULL; x.type(colvarvalue::type_scalar); map_total::init(conf); } @@ -35,26 +41,101 @@ colvar::map_total::map_total(std::string const &conf) int colvar::map_total::init(std::string const &conf) { int error_code = cvc::init(conf); - get_keyval(conf, "mapName", map_name, map_name); - volmap_index = (cvm::proxy)->init_volmap(map_name); - error_code |= volmap_index > 0 ? COLVARS_OK : INPUT_ERROR; + colvarproxy *proxy = cvm::main()->proxy; + get_keyval(conf, "mapName", volmap_name, volmap_name); + get_keyval(conf, "mapID", volmap_id, volmap_id); + register_param("mapID", reinterpret_cast(&volmap_id)); + + if ((volmap_name.size() > 0) && (volmap_id >= 0)) { + error_code |= + cvm::error("Error: mapName and mapID are mutually exclusive.\n"); + } + + // Parse optional group + atoms = parse_group(conf, "atoms", true); + if (atoms != NULL) { + + // Using internal selection + if (volmap_name.size()) { + error_code |= proxy->check_volmap_by_name(volmap_name); + } + if (volmap_id >= 0) { + error_code |= proxy->check_volmap_by_id(volmap_id); + } + + } else { + + // Using selection from the MD engine + if (volmap_name.size()) { + volmap_index = proxy->init_volmap_by_name(volmap_name); + } + if (volmap_id >= 0) { + volmap_index = proxy->init_volmap_by_id(volmap_id); + } + error_code |= volmap_index > 0 ? COLVARS_OK : INPUT_ERROR; + } + + if (get_keyval(conf, "atomWeights", atom_weights, atom_weights)) { + if (atoms == NULL) { + error_code |= cvm::error("Error: weights can only be assigned when atoms " + "are selected explicitly in Colvars.\n", + INPUT_ERROR); + } else { + if (atoms->size() != atom_weights.size()) { + error_code |= cvm::error("Error: if defined, the number of weights ("+ + cvm::to_str(atom_weights.size())+ + ") must equal the number of atoms ("+ + cvm::to_str(atoms->size())+ + ").\n", INPUT_ERROR); + } + } + } + + if (volmap_name.size() > 0) { + volmap_id = proxy->get_volmap_id_from_name(volmap_name.c_str()); + } + return error_code; } void colvar::map_total::calc_value() { - x.real_value = (cvm::proxy)->get_volmap_value(volmap_index); + colvarproxy *proxy = cvm::main()->proxy; + int flags = is_enabled(f_cvc_gradient) ? colvarproxy::volmap_flag_gradients : + colvarproxy::volmap_flag_gradients; + + if (atoms != NULL) { + // Compute the map inside Colvars + x.real_value = 0.0; + + cvm::real *w = NULL; + if (atom_weights.size() > 0) { + flags |= colvarproxy::volmap_flag_use_atom_field; + w = &(atom_weights[0]); + } + proxy->compute_volmap(flags, volmap_id, atoms->begin(), atoms->end(), + &(x.real_value), w); + } else { + // Get the externally computed value + x.real_value = proxy->get_volmap_value(volmap_index); + } } void colvar::map_total::calc_gradients() { - // Atomic coordinates are not available here + // Computed in calc_value() or by the MD engine } void colvar::map_total::apply_force(colvarvalue const &force) { - (cvm::proxy)->apply_volmap_force(volmap_index, force.real_value); + colvarproxy *proxy = cvm::main()->proxy; + if (atoms) { + if (!atoms->noforce) + atoms->apply_colvar_force(force.real_value); + } else { + proxy->apply_volmap_force(volmap_index, force.real_value); + } } diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp index 94c5cfbda7..c240ce5a96 100644 --- a/lib/colvars/colvardeps.cpp +++ b/lib/colvars/colvardeps.cpp @@ -49,7 +49,7 @@ void colvardeps::free_children_deps() { // Cannot be in the base class destructor because it needs the derived class features() size_t i,j,fid; - if (cvm::debug()) cvm::log("DEPS: freeing children deps for " + description); + if (cvm::debug()) cvm::log("DEPS: freeing children deps for " + description + "\n"); cvm::increase_depth(); for (fid = 0; fid < feature_states.size(); fid++) { @@ -58,7 +58,7 @@ void colvardeps::free_children_deps() { int g = features()[fid]->requires_children[i]; for (j=0; jfeatures()[g]->description); + + children[j]->features()[g]->description + "\n"); children[j]->decr_ref_count(g); } } @@ -80,7 +80,7 @@ void colvardeps::restore_children_deps() { int g = features()[fid]->requires_children[i]; for (j=0; jfeatures()[g]->description); + + children[j]->features()[g]->description + "\n"); children[j]->enable(g, false, false); } } @@ -135,7 +135,7 @@ int colvardeps::enable(int feature_id, if (cvm::debug()) { cvm::log("DEPS: " + description + (dry_run ? " testing " : " enabling ") + - "\"" + f->description +"\""); + "\"" + f->description +"\"\n"); } if (fs->enabled) { @@ -144,7 +144,7 @@ int colvardeps::enable(int feature_id, // as requirement is enabled fs->ref_count++; if (cvm::debug()) - cvm::log("DEPS: bumping ref_count to " + cvm::to_str(fs->ref_count)); + cvm::log("DEPS: bumping ref_count to " + cvm::to_str(fs->ref_count) + "\n"); } // Do not try to further resolve deps return COLVARS_OK; @@ -243,7 +243,7 @@ int colvardeps::enable(int feature_id, enable(g, false, false); // Just for printing error output } cvm::decrease_depth(); - cvm::log("-----------------------------------------"); + cvm::log("-----------------------------------------\n"); if (toplevel) { cvm::error("Error: Failed dependency in " + description + "."); } @@ -285,7 +285,7 @@ int colvardeps::enable(int feature_id, do_feature_side_effects(feature_id); if (cvm::debug()) cvm::log("DEPS: feature \"" + f->description + "\" in " - + description + " enabled, ref_count = 1."); + + description + " enabled, ref_count = 1." + "\n"); } return COLVARS_OK; } @@ -297,7 +297,7 @@ int colvardeps::disable(int feature_id) { feature_state *fs = &feature_states[feature_id]; if (cvm::debug()) cvm::log("DEPS: disabling feature \"" - + f->description + "\" in " + description); + + f->description + "\" in " + description + "\n"); if (fs->enabled == false) { return COLVARS_OK; @@ -313,14 +313,14 @@ int colvardeps::disable(int feature_id) { // internal deps (self) for (i=0; irequires_self.size(); i++) { if (cvm::debug()) cvm::log("DEPS: dereferencing self " - + features()[f->requires_self[i]]->description); + + features()[f->requires_self[i]]->description + "\n"); decr_ref_count(f->requires_self[i]); } // alternates for (i=0; ialternate_refs.size(); i++) { if (cvm::debug()) cvm::log("DEPS: dereferencing alt " - + features()[fs->alternate_refs[i]]->description); + + features()[fs->alternate_refs[i]]->description + "\n"); decr_ref_count(fs->alternate_refs[i]); } // Forget these, now that they are dereferenced @@ -337,7 +337,7 @@ int colvardeps::disable(int feature_id) { int g = f->requires_children[i]; for (j=0; jfeatures()[g]->description); + + children[j]->features()[g]->description + "\n"); children[j]->decr_ref_count(g); } } @@ -430,11 +430,13 @@ void colvardeps::require_feature_alt(int f, int g, int h, int i, int j) { void colvardeps::print_state() { size_t i; - cvm::log("Features of \"" + description + "\" ON/OFF (refcount)"); + cvm::log("Features of \"" + description + "\" (refcount)\n"); for (i = 0; i < feature_states.size(); i++) { - std::string onoff = is_enabled(i) ? "ON" : "OFF"; - cvm::log("- " + features()[i]->description + " " + onoff + " (" - + cvm::to_str(feature_states[i].ref_count) + ")"); + std::string onoff = is_enabled(i) ? "ON " : " "; + // Only display refcount if non-zero for less clutter + std::string refcount = feature_states[i].ref_count != 0 ? + " (" + cvm::to_str(feature_states[i].ref_count) + ") " : ""; + cvm::log("- " + onoff + features()[i]->description + refcount + "\n"); } cvm::increase_depth(); for (i=0; irequires_children.size(); i++) { int g = features()[fid]->requires_children[i]; if (cvm::debug()) cvm::log("DEPS: re-enabling children's " - + child->features()[g]->description); + + child->features()[g]->description + "\n"); child->enable(g, false, false); } } diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h index f4f543f336..9fd441ffec 100644 --- a/lib/colvars/colvardeps.h +++ b/lib/colvars/colvardeps.h @@ -283,6 +283,10 @@ public: /// center with fictitious mass; bias forces will be applied to /// the center f_cv_extended_Lagrangian, + /// \brief An extended variable that sets an external variable in the + /// back-end (eg. an alchemical coupling parameter for lambda-dynamics) + /// Can have a single component + f_cv_external, /// \brief The extended system coordinate undergoes Langevin dynamics f_cv_Langevin, /// \brief Output the potential and kinetic energies @@ -375,6 +379,7 @@ public: enum features_atomgroup { f_ag_active, f_ag_center, + f_ag_center_origin, f_ag_rotate, f_ag_fitting_group, /// Perform a standard minimum msd fit for given atoms diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index df122b1be4..2074c97aca 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -122,17 +122,83 @@ cvm::real colvar_grid_scalar::entropy() const colvar_grid_gradient::colvar_grid_gradient() - : colvar_grid(), samples(NULL) + : colvar_grid(), + samples(NULL), + weights(NULL) {} colvar_grid_gradient::colvar_grid_gradient(std::vector const &nx_i) - : colvar_grid(nx_i, 0.0, nx_i.size()), samples(NULL) + : colvar_grid(nx_i, 0.0, nx_i.size()), + samples(NULL), + weights(NULL) {} colvar_grid_gradient::colvar_grid_gradient(std::vector &colvars) - : colvar_grid(colvars, 0.0, colvars.size()), samples(NULL) + : colvar_grid(colvars, 0.0, colvars.size()), + samples(NULL), + weights(NULL) {} + +colvar_grid_gradient::colvar_grid_gradient(std::string &filename) + : colvar_grid(), + samples(NULL), + weights(NULL) +{ + std::ifstream is; + is.open(filename.c_str()); + if (!is.is_open()) { + cvm::error("Error opening multicol gradient file " + filename + " for reading.\n"); + return; + } + + // Data in the header: nColvars, then for each + // xiMin, dXi, nPoints, periodic flag + + std::string hash; + size_t i; + + if ( !(is >> hash) || (hash != "#") ) { + cvm::error("Error reading grid at position "+ + cvm::to_str(static_cast(is.tellg()))+ + " in stream(read \"" + hash + "\")\n"); + return; + } + + is >> nd; + mult = nd; + std::vector lower_in(nd), widths_in(nd); + std::vector nx_in(nd); + std::vector periodic_in(nd); + + for (i = 0; i < nd; i++ ) { + if ( !(is >> hash) || (hash != "#") ) { + cvm::error("Error reading grid at position "+ + cvm::to_str(static_cast(is.tellg()))+ + " in stream(read \"" + hash + "\")\n"); + return; + } + + is >> lower_in[i] >> widths_in[i] >> nx_in[i] >> periodic_in[i]; + } + + this->setup(nx_in, 0., mult); + + widths = widths_in; + + for (i = 0; i < nd; i++ ) { + lower_boundaries.push_back(colvarvalue(lower_in[i])); + periodic.push_back(static_cast(periodic_in[i])); + } + + // Reset the istream for read_multicol, which expects the whole file + is.clear(); + is.seekg(0); + read_multicol(is); + is.close(); +} + + void colvar_grid_gradient::write_1D_integral(std::ostream &os) { cvm::real bin, min, integral; @@ -202,7 +268,7 @@ integrate_potential::integrate_potential(std::vector &colvars, colvar_ // Compute inverse of Laplacian diagonal for Jacobi preconditioning // For now all code related to preconditioning is commented out // until a method better than Jacobi is implemented -// cvm::log("Preparing inverse diagonal for preconditioning..."); +// cvm::log("Preparing inverse diagonal for preconditioning...\n"); // inv_lap_diag.resize(nt); // std::vector id(nt), lap_col(nt); // for (int i = 0; i < nt; i++) { @@ -213,7 +279,30 @@ integrate_potential::integrate_potential(std::vector &colvars, colvar_ // id[i] = 0.; // inv_lap_diag[i] = 1. / lap_col[i]; // } -// cvm::log("Done."); +// cvm::log("Done.\n"); + } +} + + +integrate_potential::integrate_potential(colvar_grid_gradient * gradients) + : gradients(gradients) +{ + nd = gradients->num_variables(); + nx = gradients->number_of_points_vec(); + widths = gradients->widths; + periodic = gradients->periodic; + + // Expand grid by 1 bin in non-periodic dimensions + for (size_t i = 0; i < nd; i++ ) { + if (!periodic[i]) nx[i]++; + // Shift the grid by half the bin width (values at edges instead of center of bins) + lower_boundaries.push_back(gradients->lower_boundaries[i].real_value - 0.5 * widths[i]); + } + + setup(nx); + + if (nd > 1) { + divergence.resize(nt); } } @@ -246,7 +335,7 @@ int integrate_potential::integrate(const int itmax, const cvm::real &tol, cvm::r } else if (nd <= 3) { nr_linbcg_sym(divergence, data, tol, itmax, iter, err); - cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err)); + cvm::log("Integrated in " + cvm::to_str(iter) + " steps, error: " + cvm::to_str(err) + "\n"); } else { cvm::error("Cannot integrate PMF in dimension > 3\n"); diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 17d049a6eb..3642ae387c 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -102,6 +102,12 @@ public: return nd; } + /// Return the numbers of points in all dimensions + inline std::vector const &number_of_points_vec() const + { + return nx; + } + /// Return the number of points in the i-th direction, if provided, or /// the total number inline size_t number_of_points(int const icv = -1) const @@ -199,6 +205,7 @@ public: { nd = nt = 0; mult = 1; + has_parent_data = false; this->setup(); } @@ -222,9 +229,9 @@ public: hard_lower_boundaries(g.hard_lower_boundaries), hard_upper_boundaries(g.hard_upper_boundaries), widths(g.widths), + has_parent_data(false), has_data(false) - { - } + {} /// \brief Constructor from explicit grid sizes \param nx_i Number /// of grid points along each dimension \param t Initial value for @@ -233,7 +240,7 @@ public: colvar_grid(std::vector const &nx_i, T const &t = T(), size_t mult_i = 1) - : has_data(false) + : has_parent_data(false), has_data(false) { this->setup(nx_i, t, mult_i); } @@ -245,7 +252,7 @@ public: T const &t = T(), size_t mult_i = 1, bool add_extra_bin = false) - : has_data(false) + : has_parent_data(false), has_data(false) { this->init_from_colvars(colvars, t, mult_i, add_extra_bin); } @@ -1066,8 +1073,8 @@ public: std::vector nx_read; std::vector bin; - if ( cv.size() != nd ) { - cvm::error("Cannot read grid file: missing reference to colvars."); + if ( cv.size() > 0 && cv.size() != nd ) { + cvm::error("Cannot read grid file: number of variables in file differs from number referenced by grid.\n"); return is; } @@ -1525,6 +1532,9 @@ public: /// Constructor from a vector of colvars colvar_grid_gradient(std::vector &colvars); + /// Constructor from a multicol file + colvar_grid_gradient(std::string &filename); + /// \brief Get a vector with the binned value(s) indexed by ix, normalized if applicable inline void vector_value(std::vector const &ix, std::vector &v) const { @@ -1658,10 +1668,13 @@ class integrate_potential : public colvar_grid_scalar {} /// Constructor from a vector of colvars + gradient grid - integrate_potential (std::vector &colvars, colvar_grid_gradient * gradients); + integrate_potential(std::vector &colvars, colvar_grid_gradient * gradients); + + /// Constructor from a gradient grid (for processing grid files without a Colvars config) + integrate_potential(colvar_grid_gradient * gradients); /// \brief Calculate potential from divergence (in 2D); return number of steps - int integrate (const int itmax, const cvm::real & tol, cvm::real & err); + int integrate(const int itmax, const cvm::real & tol, cvm::real & err); /// \brief Update matrix containing divergence and boundary conditions /// based on new gradient point value, in neighboring bins diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 95d42560fc..405c68244b 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -34,6 +34,9 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) xyz_reader_use_count = 0; + restart_version_str.clear(); + restart_version_int = 0; + if (proxy == NULL) { proxy = proxy_in; // Pointer to the proxy object parse = new colvarparse(); // Parsing object for global options @@ -48,7 +51,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) cvm::log(cvm::line_marker); cvm::log("Initializing the collective variables module, version "+ - cvm::to_str(COLVARS_VERSION)+".\n"); + version()+".\n"); cvm::log("Please cite Fiorin et al, Mol Phys 2013:\n " "https://dx.doi.org/10.1080/00268976.2013.813594\n" "in any publication based on this calculation.\n"); @@ -58,7 +61,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) } #if (__cplusplus >= 201103L) - cvm::log("This version was built with the C++11 standard or higher."); + cvm::log("This version was built with the C++11 standard or higher.\n"); #else cvm::log("This version was built without the C++11 standard: some features are disabled.\n" "Please see the following link for details:\n" @@ -186,6 +189,7 @@ std::istream & colvarmodule::getline(std::istream &is, std::string &line) size_t const sz = l.size(); if (sz > 0) { if (l[sz-1] == '\r' ) { + // Replace Windows newlines with Unix newlines line = l.substr(0, sz-1); } else { line = l; @@ -200,6 +204,7 @@ std::istream & colvarmodule::getline(std::istream &is, std::string &line) int colvarmodule::parse_config(std::string &conf) { + // Auto-generated additional configuration extra_conf.clear(); // Check that the input has matching braces @@ -208,6 +213,9 @@ int colvarmodule::parse_config(std::string &conf) INPUT_ERROR); } + // Check that the input has only ASCII characters, and warn otherwise + colvarparse::check_ascii(conf); + // Parse global options if (catch_input_errors(parse_global_params(conf))) { return get_error(); @@ -472,7 +480,7 @@ int colvarmodule::parse_biases(std::string const &conf) if (use_scripted_forces) { cvm::log(cvm::line_marker); cvm::increase_depth(); - cvm::log("User forces script will be run at each bias update."); + cvm::log("User forces script will be run at each bias update.\n"); cvm::decrease_depth(); } @@ -754,6 +762,9 @@ int colvarmodule::calc() error_code |= end_of_step(); + // TODO move this to a base-class proxy method that calls this function + error_code |= proxy->end_of_step(); + return error_code; } @@ -1311,21 +1322,23 @@ std::istream & colvarmodule::read_restart(std::istream &is) colvarparse::parse_restart); it = it_restart; - std::string restart_version; - int restart_version_int = 0; + restart_version_str.clear(); + restart_version_int = 0; parse->get_keyval(restart_conf, "version", - restart_version, std::string(""), + restart_version_str, std::string(""), colvarparse::parse_restart); - if (restart_version.size()) { - if (restart_version != std::string(COLVARS_VERSION)) { - cvm::log("This state file was generated with version "+ - restart_version+"\n"); - } + if (restart_version_str.size()) { + // Initialize integer version number of this restart file restart_version_int = - proxy->get_version_from_string(restart_version.c_str()); + proxy->get_version_from_string(restart_version_str.c_str()); } - if (restart_version_int < 20160810) { + if (restart_version() != version()) { + cvm::log("This state file was generated with version "+ + restart_version()+"\n"); + } + + if (restart_version_number() < 20160810) { // check for total force change if (proxy->total_forces_enabled()) { warn_total_forces = true; @@ -1769,6 +1782,8 @@ int cvm::read_index_file(char const *filename) cvm::error("Error: in opening index file \""+ std::string(filename)+"\".\n", FILE_ERROR); + } else { + index_file_names.push_back(std::string(filename)); } while (is.good()) { @@ -1861,6 +1876,7 @@ int colvarmodule::reset_index_groups() } index_group_names.clear(); index_groups.clear(); + index_file_names.clear(); return COLVARS_OK; } @@ -1924,48 +1940,75 @@ int cvm::load_coords_xyz(char const *filename, std::string line; cvm::real x = 0.0, y = 0.0, z = 0.0; + std::string const error_msg("Error: cannot parse XYZ file \""+ + std::string(filename)+"\".\n"); + if ( ! (xyz_is >> natoms) ) { - cvm::error("Error: cannot parse XYZ file " - + std::string(filename) + ".\n", INPUT_ERROR); + return cvm::error(error_msg, INPUT_ERROR); } ++xyz_reader_use_count; if (xyz_reader_use_count < 2) { - cvm::log("Warning: beginning from 2019-11-26 the XYZ file reader assumes Angstrom units."); + cvm::log("Warning: beginning from 2019-11-26 the XYZ file reader assumes Angstrom units.\n"); } - // skip comment line - cvm::getline(xyz_is, line); - cvm::getline(xyz_is, line); - xyz_is.width(255); - std::vector::iterator pos_i = pos->begin(); + if (xyz_is.good()) { + // skip comment line + cvm::getline(xyz_is, line); + cvm::getline(xyz_is, line); + xyz_is.width(255); + } else { + return cvm::error(error_msg, INPUT_ERROR); + } + std::vector::iterator pos_i = pos->begin(); + size_t xyz_natoms = 0; if (pos->size() != natoms) { // Use specified indices int next = 0; // indices are zero-based std::vector::const_iterator index = atoms->sorted_ids().begin(); - for ( ; pos_i != pos->end() ; pos_i++, index++) { + for ( ; pos_i != pos->end() ; pos_i++, index++) { while ( next < *index ) { cvm::getline(xyz_is, line); next++; } - xyz_is >> symbol; - xyz_is >> x >> y >> z; - // XYZ files are assumed to be in Angstrom (as eg. VMD will) - (*pos_i)[0] = proxy->angstrom_to_internal(x); - (*pos_i)[1] = proxy->angstrom_to_internal(y); - (*pos_i)[2] = proxy->angstrom_to_internal(z); + if (xyz_is.good()) { + xyz_is >> symbol; + xyz_is >> x >> y >> z; + // XYZ files are assumed to be in Angstrom (as eg. VMD will) + (*pos_i)[0] = proxy->angstrom_to_internal(x); + (*pos_i)[1] = proxy->angstrom_to_internal(y); + (*pos_i)[2] = proxy->angstrom_to_internal(z); + xyz_natoms++; + } else { + return cvm::error(error_msg, INPUT_ERROR); + } } + } else { // Use all positions + for ( ; pos_i != pos->end() ; pos_i++) { - xyz_is >> symbol; - xyz_is >> x >> y >> z; - (*pos_i)[0] = proxy->angstrom_to_internal(x); - (*pos_i)[1] = proxy->angstrom_to_internal(y); - (*pos_i)[2] = proxy->angstrom_to_internal(z); + if (xyz_is.good()) { + xyz_is >> symbol; + xyz_is >> x >> y >> z; + (*pos_i)[0] = proxy->angstrom_to_internal(x); + (*pos_i)[1] = proxy->angstrom_to_internal(y); + (*pos_i)[2] = proxy->angstrom_to_internal(z); + xyz_natoms++; + } else { + return cvm::error(error_msg, INPUT_ERROR); + } } } - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + + if (xyz_natoms != pos->size()) { + return cvm::error("Error: The number of positions read from file \""+ + std::string(filename)+"\" does not match the number of "+ + "positions required: "+cvm::to_str(xyz_natoms)+" vs. "+ + cvm::to_str(pos->size())+".\n", INPUT_ERROR); + } + + return COLVARS_OK; } diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 1cad4195d6..3d4296a4c9 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -81,6 +81,12 @@ private: public: + /// Get the version string (YYYY-MM-DD format) + std::string version() const + { + return std::string(COLVARS_VERSION); + } + /// Get the version number (higher = more recent) int version_number() const { @@ -150,6 +156,12 @@ public: return ::cos(static_cast(x)); } + /// Reimplemented to work around MS compiler issues + static inline real asin(real const &x) + { + return ::asin(static_cast(x)); + } + /// Reimplemented to work around MS compiler issues static inline real acos(real const &x) { @@ -685,6 +697,9 @@ public: static rvector position_distance(atom_pos const &pos1, atom_pos const &pos2); + /// \brief Names of .ndx files that have been loaded + std::vector index_file_names; + /// \brief Names of groups from one or more Gromacs .ndx files std::vector index_group_names; @@ -758,7 +773,11 @@ protected: /// Write labels at the next iteration bool cv_traj_write_labels; -private: + /// Version of the most recent state file read + std::string restart_version_str; + + /// Integer version of the most recent state file read + int restart_version_int; /// Counter for the current depth in the object hierarchy (useg e.g. in output) size_t depth_s; @@ -771,6 +790,18 @@ private: public: + /// Version of the most recent state file read + inline std::string restart_version() const + { + return restart_version_str; + } + + /// Integer version of the most recent state file read + inline int restart_version_number() const + { + return restart_version_int; + } + /// Get the current object depth in the hierarchy static size_t & depth(); diff --git a/lib/colvars/colvarmodule_utils.h b/lib/colvars/colvarmodule_utils.h new file mode 100644 index 0000000000..a7004edd92 --- /dev/null +++ b/lib/colvars/colvarmodule_utils.h @@ -0,0 +1,80 @@ +// -*- c++ -*- + +// This file is part of the Collective Variables module (Colvars). +// The original version of Colvars and its updates are located at: +// https://github.com/Colvars/colvars +// Please update all Colvars source files before making any changes. +// If you wish to distribute your changes, please submit them to the +// Colvars repository at GitHub. + + +#ifndef COLVARMODULE_UTILS_H +#define COLVARMODULE_UTILS_H + + +#include "colvarmodule.h" + + +template +cvm::real get_force_norm2(T const &x) +{ + return x.norm2(); +} + + +template <> +inline cvm::real get_force_norm2(cvm::real const &x) +{ + return x*x; +} + + +template +cvm::real compute_norm2_stats(std::vector const &v, + int *minmax_index = NULL) +{ + cvm::real result = 0.0; + if (flag == -1) { + // Initialize for minimum search, using approx. largest float32 value + result = 1.0e38; + } + + typename std::vector::const_iterator xi = v.begin(); + size_t i = 0; + + if (get_index) *minmax_index = -1; // Let's not assume minmax_index is initialized to -1 + + for ( ; xi != v.end(); xi++, i++) { + cvm::real const norm2 = get_force_norm2(*xi); + if (flag == 0) { + result += norm2; + } + if (flag == 1) { + if (norm2 > result) { + result = norm2; + if (get_index) *minmax_index = i; + } + } + if (flag == -1) { + if (norm2 < result) { + result = norm2; + if (get_index) *minmax_index = i; + } + } + } + + size_t const n = v.size(); + + if (flag == 0) { + if (n > 0) { + result /= cvm::real(n); + } + } + + result = cvm::sqrt(result); + + return result; +} + + +#endif diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index cbe19c1914..3449a681f3 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -125,6 +125,10 @@ void colvarparse::mark_key_set_user(std::string const &key_str, cvm::log("# "+key_str+" = "+cvm::to_str(value)+"\n", cvm::log_user_params()); } + if (parse_mode & parse_deprecation_warning) { + cvm::log("Warning: keyword "+key_str+ + " is deprecated. Check the documentation for the current equivalent.\n"); + } } @@ -919,6 +923,26 @@ int colvarparse::check_braces(std::string const &conf, return (brace_count != 0) ? INPUT_ERROR : COLVARS_OK; } + +int colvarparse::check_ascii(std::string const &conf) +{ + // Check for non-ASCII characters + std::string line; + std::istringstream is(conf); + while (cvm::getline(is, line)) { + unsigned char const * const uchars = + reinterpret_cast(line.c_str()); + for (size_t i = 0; i < line.size(); i++) { + if (uchars[i] & 0x80U) { + cvm::log("Warning: non-ASCII character detected in this line: \""+ + line+"\".\n"); + } + } + } + return COLVARS_OK; +} + + void colvarparse::split_string(const std::string& data, const std::string& delim, std::vector& dest) { size_t index = 0, new_index = 0; std::string tmpstr; diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index b7d42fdffa..8e35896f89 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -56,6 +56,8 @@ public: parse_echo = (1<<1), /// Print the default value of a keyword, if it is NOT given parse_echo_default = (1<<2), + /// Print a deprecation warning if the keyword is given + parse_deprecation_warning = (1<<3), /// Do not print the keyword parse_silent = 0, /// Raise error if the keyword is not provided @@ -66,7 +68,9 @@ public: /// The call is being executed from a read_restart() function parse_restart = (1<<18), /// Alias for old default behavior (should be phased out) - parse_normal = (1<<2) | (1<<1) | (1<<17) + parse_normal = (1<<2) | (1<<1) | (1<<17), + /// Settings for a deprecated keyword + parse_deprecated = (1<<1) | (1<<3) | (1<<17) }; /// \brief Check that all the keywords within "conf" are in the list @@ -317,6 +321,10 @@ public: /// from this position static int check_braces(std::string const &conf, size_t const start_pos); + /// \brief Check that a config string contains non-ASCII characters + /// \param conf The configuration string + static int check_ascii(std::string const &conf); + /// \brief Split a string with a specified delimiter into a vector /// \param data The string to be splitted /// \param delim A delimiter diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp index 24f833f857..d0f83c70a7 100644 --- a/lib/colvars/colvarproxy.cpp +++ b/lib/colvars/colvarproxy.cpp @@ -24,12 +24,15 @@ #include "colvarproxy.h" #include "colvarscript.h" #include "colvaratoms.h" +#include "colvarmodule_utils.h" colvarproxy_system::colvarproxy_system() { angstrom_value = 0.0; + kcal_mol_value = 0.0; + boundaries_type = boundaries_unsupported; total_force_requested = false; reset_pbc_lattice(); } @@ -38,6 +41,46 @@ colvarproxy_system::colvarproxy_system() colvarproxy_system::~colvarproxy_system() {} +int colvarproxy_system::set_unit_system(std::string const & /* units */, + bool /* check_only */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +cvm::real colvarproxy_system::backend_angstrom_value() +{ + return 1.0; +} + + +cvm::real colvarproxy_system::boltzmann() +{ + return 0.001987191; +} + + +cvm::real colvarproxy_system::temperature() +{ + // TODO define, document and implement a user method to set the value of this + return 300.0; +} + + +cvm::real colvarproxy_system::dt() +{ + // TODO define, document and implement a user method to set the value of this + return 1.0; +} + + +cvm::real colvarproxy_system::rand_gaussian() +{ + // TODO define, document and implement a user method to set the value of this + return 0.0; +} + + void colvarproxy_system::add_energy(cvm::real /* energy */) {} @@ -139,9 +182,31 @@ int colvarproxy_system::get_molid(int &) } +int colvarproxy_system::get_alch_lambda(cvm::real* lambda) +{ + return cvm::error("Error in get_alch_lambda: alchemical lambda dynamics is not supported by this build.", + COLVARS_NOT_IMPLEMENTED); +} + + +int colvarproxy_system::set_alch_lambda(cvm::real* lambda) +{ + return cvm::error("Error in set_alch_lambda: alchemical lambda dynamics is not supported by this build.", + COLVARS_NOT_IMPLEMENTED); +} + + +int colvarproxy_system::get_dE_dLambda(cvm::real* force) +{ + return cvm::error("Error in get_dE_dLambda: alchemical lambda dynamics is not supported by this build.", + COLVARS_NOT_IMPLEMENTED); +} + colvarproxy_atoms::colvarproxy_atoms() { + atoms_rms_applied_force_ = atoms_max_applied_force_ = 0.0; + atoms_max_applied_force_id_ = -1; updated_masses_ = updated_charges_ = false; } @@ -178,6 +243,18 @@ int colvarproxy_atoms::add_atom_slot(int atom_id) } +int colvarproxy_atoms::init_atom(int /* atom_number */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_atoms::check_atom_id(int /* atom_number */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + int colvarproxy_atoms::init_atom(cvm::residue_id const & /* residue */, std::string const & /* atom_name */, std::string const & /* segment_id */) @@ -232,8 +309,39 @@ int colvarproxy_atoms::load_coords(char const * /* filename */, } +void colvarproxy_atoms::compute_rms_atoms_applied_force() +{ + atoms_rms_applied_force_ = + compute_norm2_stats(atoms_new_colvar_forces); +} -colvarproxy_atom_groups::colvarproxy_atom_groups() {} + +void colvarproxy_atoms::compute_max_atoms_applied_force() +{ + int minmax_index = -1; + size_t const n_atoms_ids = atoms_ids.size(); + if ((n_atoms_ids > 0) && (n_atoms_ids == atoms_new_colvar_forces.size())) { + atoms_max_applied_force_ = + compute_norm2_stats(atoms_new_colvar_forces, + &minmax_index); + if (minmax_index >= 0) { + atoms_max_applied_force_id_ = atoms_ids[minmax_index]; + } else { + atoms_max_applied_force_id_ = -1; + } + } else { + atoms_max_applied_force_ = + compute_norm2_stats(atoms_new_colvar_forces); + atoms_max_applied_force_id_ = -1; + } +} + + + +colvarproxy_atom_groups::colvarproxy_atom_groups() +{ + atom_groups_rms_applied_force_ = atom_groups_max_applied_force_ = 0.0; +} colvarproxy_atom_groups::~colvarproxy_atom_groups() @@ -296,6 +404,20 @@ void colvarproxy_atom_groups::clear_atom_group(int index) } +void colvarproxy_atom_groups::compute_rms_atom_groups_applied_force() +{ + atom_groups_rms_applied_force_ = + compute_norm2_stats(atom_groups_new_colvar_forces); +} + + +void colvarproxy_atom_groups::compute_max_atom_groups_applied_force() +{ + atom_groups_max_applied_force_ = + compute_norm2_stats(atom_groups_new_colvar_forces); +} + + colvarproxy_smp::colvarproxy_smp() { @@ -464,28 +586,14 @@ int colvarproxy_smp::smp_unlock() colvarproxy_script::colvarproxy_script() { script = NULL; + force_script_defined = false; + have_scripts = false; } colvarproxy_script::~colvarproxy_script() {} -char const *colvarproxy_script::script_obj_to_str(unsigned char *obj) -{ - cvm::error("Error: trying to print a script object without a scripting " - "language interface.\n", BUG_ERROR); - return reinterpret_cast(obj); -} - - -std::vector colvarproxy_script::script_obj_to_str_vector(unsigned char * /* obj */) -{ - cvm::error("Error: trying to print a script object without a scripting " - "language interface.\n", BUG_ERROR); - return std::vector(); -} - - int colvarproxy_script::run_force_callback() { return COLVARS_NOT_IMPLEMENTED; @@ -512,6 +620,7 @@ int colvarproxy_script::run_colvar_gradient_callback(std::string const & /* name colvarproxy_io::colvarproxy_io() { input_buffer_ = NULL; + restart_frequency_engine = 0; } @@ -660,6 +769,23 @@ int colvarproxy::update_output() } +int colvarproxy::end_of_step() +{ + // Disable flags that Colvars doesn't need any more + updated_masses_ = updated_charges_ = false; + + // Compute force statistics + compute_rms_atoms_applied_force(); + compute_max_atoms_applied_force(); + compute_rms_atom_groups_applied_force(); + compute_max_atom_groups_applied_force(); + compute_rms_volmaps_applied_force(); + compute_max_volmaps_applied_force(); + + return COLVARS_OK; +} + + int colvarproxy::post_run() { int error_code = COLVARS_OK; @@ -672,6 +798,19 @@ int colvarproxy::post_run() } +void colvarproxy::log(std::string const &message) +{ + fprintf(stdout, "colvars: %s", message.c_str()); +} + + +void colvarproxy::error(std::string const &message) +{ + // TODO handle errors? + colvarproxy::log(message); +} + + void colvarproxy::add_error_msg(std::string const &message) { std::istringstream is(message); diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 7a60292092..c9841ebdf4 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -59,7 +59,7 @@ public: std::string units; /// \brief Request to set the units used internally by Colvars - virtual int set_unit_system(std::string const &units, bool check_only) = 0; + virtual int set_unit_system(std::string const &units, bool check_only); /// \brief Value of 1 Angstrom in the internal (front-end) Colvars unit for atomic coordinates /// * defaults to 0. in the base class; derived proxy classes must set it @@ -68,7 +68,7 @@ public: cvm::real angstrom_value; /// \brief Value of 1 Angstrom in the backend's unit for atomic coordinates - virtual cvm::real backend_angstrom_value() = 0; + virtual cvm::real backend_angstrom_value(); /// \brief Value of 1 kcal/mol in the internal Colvars unit for energy cvm::real kcal_mol_value; @@ -79,6 +79,12 @@ public: return l * angstrom_value; } + /// \brief Convert a length from internal to Angstrom + inline cvm::real internal_to_angstrom(cvm::real l) const + { + return l / angstrom_value; + } + // /// \brief Convert a length from back-end unit to internal // inline cvm::real back_end_to_internal_unit(cvm::real l) { // if (angstrom_value == 0.) { @@ -88,19 +94,19 @@ public: // } /// \brief Boltzmann constant in internal Colvars units - virtual cvm::real boltzmann() = 0; + virtual cvm::real boltzmann(); /// \brief Target temperature of the simulation (K units) - virtual cvm::real temperature() = 0; + virtual cvm::real temperature(); /// \brief Time step of the simulation (fs) - virtual cvm::real dt() = 0; + virtual cvm::real dt(); /// \brief Pseudo-random number with Gaussian distribution - virtual cvm::real rand_gaussian(void) = 0; + virtual cvm::real rand_gaussian(void); /// Pass restraint energy value for current timestep to MD engine - virtual void add_energy(cvm::real energy) = 0; + virtual void add_energy(cvm::real energy); /// \brief Get the PBC-aware distance vector between two positions virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, @@ -126,6 +132,15 @@ public: /// \param molid Set this argument equal to the current VMD molid virtual int get_molid(int &molid); + /// Get value of alchemical lambda parameter from back-end (if available) + virtual int get_alch_lambda(cvm::real* lambda); + + /// Set value of alchemical lambda parameter in back-end (if available) + virtual int set_alch_lambda(cvm::real* lambda); + + /// Get energy derivative with respect to lambda (if available) + virtual int get_dE_dLambda(cvm::real* force); + protected: /// Whether the total forces have been requested @@ -167,11 +182,11 @@ public: /// Prepare this atom for collective variables calculation, selecting it by /// numeric index (1-based) - virtual int init_atom(int atom_number) = 0; + virtual int init_atom(int atom_number); /// Check that this atom number is valid, but do not initialize the /// corresponding atom yet - virtual int check_atom_id(int atom_number) = 0; + virtual int check_atom_id(int atom_number); /// Select this atom for collective variables calculation, using name and /// residue number. Not all programs support this: leave this function as @@ -262,11 +277,16 @@ public: return cvm::rvector(0.0); } - inline std::vector *modify_atom_ids() + inline std::vector const *get_atom_ids() const { return &atoms_ids; } + inline std::vector const *get_atom_masses() const + { + return &atoms_masses; + } + inline std::vector *modify_atom_masses() { // assume that we are requesting masses to change them @@ -274,6 +294,11 @@ public: return &atoms_masses; } + inline std::vector const *get_atom_charges() + { + return &atoms_charges; + } + inline std::vector *modify_atom_charges() { // assume that we are requesting charges to change them @@ -281,21 +306,60 @@ public: return &atoms_charges; } + inline std::vector const *get_atom_positions() const + { + return &atoms_positions; + } + inline std::vector *modify_atom_positions() { return &atoms_positions; } + inline std::vector const *get_atom_total_forces() const + { + return &atoms_total_forces; + } + inline std::vector *modify_atom_total_forces() { return &atoms_total_forces; } - inline std::vector *modify_atom_new_colvar_forces() + inline std::vector const *get_atom_applied_forces() const { return &atoms_new_colvar_forces; } + inline std::vector *modify_atom_applied_forces() + { + return &atoms_new_colvar_forces; + } + + /// Compute the root-mean-square of the applied forces + void compute_rms_atoms_applied_force(); + + /// Compute the maximum norm among all applied forces + void compute_max_atoms_applied_force(); + + /// Get the root-mean-square of the applied forces + inline cvm::real rms_atoms_applied_force() const + { + return atoms_rms_applied_force_; + } + + /// Get the maximum norm among all applied forces + inline cvm::real max_atoms_applied_force() const + { + return atoms_max_applied_force_; + } + + /// Get the atom ID with the largest applied force + inline int max_atoms_applied_force_id() const + { + return atoms_max_applied_force_id_; + } + /// Record whether masses have been updated inline bool updated_masses() const { @@ -326,6 +390,15 @@ protected: /// \brief Forces applied from colvars, to be communicated to the MD integrator std::vector atoms_new_colvar_forces; + /// Root-mean-square of the applied forces + cvm::real atoms_rms_applied_force_; + + /// Maximum norm among all applied forces + cvm::real atoms_max_applied_force_; + + /// ID of the atom with the maximum norm among all applied forces + int atoms_max_applied_force_id_; + /// Whether the masses and charges have been updated from the host code bool updated_masses_, updated_charges_; @@ -404,6 +477,56 @@ public: return cvm::rvector(0.0); } + inline std::vector const *get_atom_group_ids() const + { + return &atom_groups_ids; + } + + inline std::vector *modify_atom_group_masses() + { + // TODO updated_masses + return &atom_groups_masses; + } + + inline std::vector *modify_atom_group_charges() + { + // TODO updated masses + return &atom_groups_charges; + } + + inline std::vector *modify_atom_group_positions() + { + return &atom_groups_coms; + } + + inline std::vector *modify_atom_group_total_forces() + { + return &atom_groups_total_forces; + } + + inline std::vector *modify_atom_group_applied_forces() + { + return &atom_groups_new_colvar_forces; + } + + /// Compute the root-mean-square of the applied forces + void compute_rms_atom_groups_applied_force(); + + /// Compute the maximum norm among all applied forces + void compute_max_atom_groups_applied_force(); + + /// Get the root-mean-square of the applied forces + inline cvm::real rms_atom_groups_applied_force() const + { + return atom_groups_rms_applied_force_; + } + + /// Get the maximum norm among all applied forces + inline cvm::real max_atom_groups_applied_force() const + { + return atom_groups_max_applied_force_; + } + protected: /// \brief Array of 0-based integers used to uniquely associate atom groups @@ -422,6 +545,12 @@ protected: /// \brief Forces applied from colvars, to be communicated to the MD integrator std::vector atom_groups_new_colvar_forces; + /// Root-mean-square of the applied group forces + cvm::real atom_groups_rms_applied_force_; + + /// Maximum norm among all applied group forces + cvm::real atom_groups_max_applied_force_; + /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet int add_atom_group_slot(int atom_group_id); }; @@ -519,12 +648,6 @@ public: /// Destructor virtual ~colvarproxy_script(); - /// Convert a script object (Tcl or Python call argument) to a C string - virtual char const *script_obj_to_str(unsigned char *obj); - - /// Convert a script object (Tcl or Python call argument) to a vector of strings - virtual std::vector script_obj_to_str_vector(unsigned char *obj); - /// Pointer to the scripting interface object /// (does not need to be allocated in a new interface) colvarscript *script; @@ -706,11 +829,14 @@ public: /// \brief Update data based from the results of a module update (e.g. send forces) virtual int update_output(); + /// Carry out operations needed before next step is run + int end_of_step(); + /// Print a message to the main log - virtual void log(std::string const &message) = 0; + virtual void log(std::string const &message); /// Print a message to the main log and/or let the host code know about it - virtual void error(std::string const &message) = 0; + virtual void error(std::string const &message); /// Record error message (used by VMD to collect them after a script call) void add_error_msg(std::string const &message); diff --git a/lib/colvars/colvarproxy_tcl.cpp b/lib/colvars/colvarproxy_tcl.cpp index a799bead7d..0a5473cae9 100644 --- a/lib/colvars/colvarproxy_tcl.cpp +++ b/lib/colvars/colvarproxy_tcl.cpp @@ -11,6 +11,8 @@ #if defined(NAMD_TCL) || defined(VMDTCL) #define COLVARS_TCL +#endif +#ifdef COLVARS_TCL #include #endif diff --git a/lib/colvars/colvarproxy_tcl.h b/lib/colvars/colvarproxy_tcl.h index 371b3c0224..badb62f900 100644 --- a/lib/colvars/colvarproxy_tcl.h +++ b/lib/colvars/colvarproxy_tcl.h @@ -27,7 +27,7 @@ public: /// Is Tcl available? (trigger initialization if needed) int tcl_available(); - /// Tcl implementation of script_obj_to_str() + /// Get a string representation of the Tcl object pointed to by obj char const *tcl_get_str(void *obj); /// Tcl implementation of run_force_callback() @@ -51,6 +51,12 @@ public: return tcl_interp_; } + /// Set the pointer to the Tcl interpreter + inline void set_tcl_interp(void *interp) + { + tcl_interp_ = interp; + } + protected: /// Pointer to Tcl interpreter object diff --git a/lib/colvars/colvarproxy_volmaps.cpp b/lib/colvars/colvarproxy_volmaps.cpp index fc665eec99..03e5d303d1 100644 --- a/lib/colvars/colvarproxy_volmaps.cpp +++ b/lib/colvars/colvarproxy_volmaps.cpp @@ -9,9 +9,13 @@ #include "colvarmodule.h" #include "colvarproxy_volmaps.h" +#include "colvarmodule_utils.h" -colvarproxy_volmaps::colvarproxy_volmaps() {} +colvarproxy_volmaps::colvarproxy_volmaps() +{ + volmaps_rms_applied_force_ = volmaps_max_applied_force_ = 0.0; +} colvarproxy_volmaps::~colvarproxy_volmaps() {} @@ -46,25 +50,41 @@ int colvarproxy_volmaps::add_volmap_slot(int volmap_id) } -int colvarproxy_volmaps::init_volmap(int volmap_id) +int colvarproxy_volmaps::check_volmap_by_id(int /* volmap_id */) { - return cvm::error("Error: access to volumetric maps is unavailable " - "in this build.\n", + return cvm::error("Error: selecting volumetric maps is not available.\n", COLVARS_NOT_IMPLEMENTED); } -int colvarproxy_volmaps::init_volmap(const char *volmap_name) +int colvarproxy_volmaps::check_volmap_by_name(const char * /* volmap_name */) { - return cvm::error("Error: access to volumetric maps is unavailable " - "in this build.\n", - COLVARS_NOT_IMPLEMENTED); + return cvm::error("Error: selecting volumetric maps by name is not " + "available.\n", COLVARS_NOT_IMPLEMENTED); } -int colvarproxy_volmaps::init_volmap(const std::string &volmap_name) +int colvarproxy_volmaps::init_volmap_by_name(char const *volmap_name) { - return init_volmap(volmap_name.c_str()); + return -1; +} + + +int colvarproxy_volmaps::init_volmap_by_id(int volmap_id) +{ + return -1; +} + + +int colvarproxy_volmaps::init_volmap_by_name(std::string const &volmap_name) +{ + return init_volmap_by_name(volmap_name.c_str()); +} + + +int colvarproxy_volmaps::check_volmap_by_name(std::string const &volmap_name) +{ + return check_volmap_by_name(volmap_name.c_str()); } @@ -79,3 +99,36 @@ void colvarproxy_volmaps::clear_volmap(int index) volmaps_ncopies[index] -= 1; } } + + +int colvarproxy_volmaps::get_volmap_id_from_name(char const *volmap_name) +{ + // Raise error + colvarproxy_volmaps::check_volmap_by_name(volmap_name); + return -1; +} + + +int colvarproxy_volmaps::compute_volmap(int /* flags */, + int /* volmap_id */, + cvm::atom_iter /* atom_begin */, + cvm::atom_iter /* atom_end */, + cvm::real * /* value */, + cvm::real * /* atom_field */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +void colvarproxy_volmaps::compute_rms_volmaps_applied_force() +{ + volmaps_rms_applied_force_ = + compute_norm2_stats(volmaps_new_colvar_forces); +} + + +void colvarproxy_volmaps::compute_max_volmaps_applied_force() +{ + volmaps_max_applied_force_ = + compute_norm2_stats(volmaps_new_colvar_forces); +} diff --git a/lib/colvars/colvarproxy_volmaps.h b/lib/colvars/colvarproxy_volmaps.h index 67cbb2cd2d..6e88ee83f9 100644 --- a/lib/colvars/colvarproxy_volmaps.h +++ b/lib/colvars/colvarproxy_volmaps.h @@ -25,17 +25,37 @@ public: int add_volmap_slot(int volmap_id); /// Request and prepare this volumetric map for use by Colvars - virtual int init_volmap(int volmap_id); + /// \param volmap_id Numeric ID used by the MD engine + /// \returns Index of the map in the colvarproxy arrays + virtual int init_volmap_by_id(int volmap_id); /// Request and prepare this volumetric map for use by Colvars - virtual int init_volmap(char const *volmap_name); + /// \param volmap_name Name used by the MD engine + /// \returns Index of the map in the colvarproxy arrays + virtual int init_volmap_by_name(char const *volmap_name); + + /// Check that the given volmap ID is valid (return COLVARS_OK if it is) + /// \param volmap_id Numeric ID used by the MD engine + /// \returns Error code + virtual int check_volmap_by_id(int volmap_id); + + /// Check that the given volmap name is valid (return COLVARS_OK if it is) + /// \param volmap_name Name used by the MD engine + /// \returns Error code + virtual int check_volmap_by_name(char const *volmap_name); /// Request and prepare this volumetric map for use by Colvars - int init_volmap(std::string const &volmap_name); + int init_volmap_by_name(std::string const &volmap_name); + + /// Check that the given volmap name is valid (return COLVARS_OK if it is) + int check_volmap_by_name(std::string const &volmap_name); /// \brief Used by the CVC destructors virtual void clear_volmap(int index); + /// Get the numeric ID of the given volumetric map (for the MD program) + virtual int get_volmap_id_from_name(char const *volmap_name); + /// Get the numeric ID of the given volumetric map (for the MD program) inline int get_volmap_id(int index) const { @@ -54,6 +74,32 @@ public: volmaps_new_colvar_forces[index] += new_force; } + /// Re-weigh an atomic field (e.g. a colvar) by the value of a volumetric map + /// \param flags Combination of flags + /// \param volmap_id Numeric index of the map (no need to request it) + /// \param atom_begin Iterator pointing to first atom + /// \param atom_end Iterator pointing past the last atom + /// \param value Pointer to location of total to increment + /// \param atom_field Array of atomic field values (if NULL, ones are used) + virtual int compute_volmap(int flags, + int volmap_id, + cvm::atom_iter atom_begin, + cvm::atom_iter atom_end, + cvm::real *value, + cvm::real *atom_field); + + /// Flags controlling what computation is done on the map + enum { + volmap_flag_null = 0, + volmap_flag_gradients = 1, + volmap_flag_use_atom_field = (1<<8) + }; + + /// Compute the root-mean-square of the applied forces + void compute_rms_volmaps_applied_force(); + + /// Compute the maximum norm among all applied forces + void compute_max_volmaps_applied_force(); protected: @@ -70,6 +116,12 @@ protected: /// \brief Forces applied from colvars, to be communicated to the MD /// integrator std::vector volmaps_new_colvar_forces; + + /// Root-mean-square of the the applied forces + cvm::real volmaps_rms_applied_force_; + + /// Maximum norm among all applied forces + cvm::real volmaps_max_applied_force_; }; diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index 33e05d72d1..7847bd4aba 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -1,3 +1,3 @@ #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2020-09-17" +#define COLVARS_VERSION "2021-08-03" #endif diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index ebd52b10ad..490ff6e81c 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -13,6 +13,8 @@ #if defined(NAMD_TCL) || defined(VMDTCL) #define COLVARS_TCL +#endif +#ifdef COLVARS_TCL #include #endif @@ -23,6 +25,19 @@ +#ifdef COLVARS_TCL +/// Run the script API via Tcl command-line interface +/// \param clientData Not used +/// \param my_interp Pointer to Tcl_Interp object (read from Colvars if NULL) +/// \param objc Number of Tcl command parameters +/// \param objv Array of command parameters +/// \return Result of the script command +extern "C" int tcl_run_colvarscript_command(ClientData clientData, + Tcl_Interp *interp_in, + int objc, Tcl_Obj *const objv[]); +#endif + + colvarscript::colvarscript(colvarproxy *p) : proxy_(p), colvars(p->colvars), @@ -57,9 +72,11 @@ int colvarscript::init_commands() } cmd_help.resize(colvarscript::cv_n_commands); + cmd_rethelp.resize(colvarscript::cv_n_commands); cmd_n_args_min.resize(colvarscript::cv_n_commands); cmd_n_args_max.resize(colvarscript::cv_n_commands); cmd_arghelp.resize(colvarscript::cv_n_commands); + cmd_full_help.resize(colvarscript::cv_n_commands); cmd_fns.resize(colvarscript::cv_n_commands); if (cmd_names) { @@ -95,24 +112,60 @@ int colvarscript::init_command(colvarscript::command const &comm, { cmd_str_map[std::string(name)] = comm; cmd_names[comm] = name; - cmd_help[comm] = help; + + // Initialize short help string and return-value help string (if present) + { + std::string const help_str(help); + std::istringstream is(help_str); + std::string line; + std::getline(is, line); + cmd_help[comm] = line; + cmd_rethelp[comm] = ""; + while (std::getline(is, line)) { + cmd_rethelp[comm] += line + "\n"; + } + } + + // Initialize arguments' help strings cmd_n_args_min[comm] = n_args_min; cmd_n_args_max[comm] = n_args_max; - std::string const arghelp_str(arghelp); - std::istringstream is(arghelp_str); - std::string line; - for (int iarg = 0; iarg < n_args_max; iarg++) { - if (! std::getline(is, line)) { - return cvm::error("Error: could not initialize help string for scripting " - "command \""+std::string(name)+"\".\n", BUG_ERROR); + { + std::string const arghelp_str(arghelp); + std::istringstream is(arghelp_str); + std::string line; + for (int iarg = 0; iarg < n_args_max; iarg++) { + if (! std::getline(is, line)) { + return cvm::error("Error: could not initialize help string for scripting " + "command \""+std::string(name)+"\".\n", BUG_ERROR); + } + cmd_arghelp[comm].push_back(line); } - cmd_arghelp[comm].push_back(line); } + + cmd_full_help[comm] = cmd_help[comm]+"\n"; + if (cmd_n_args_min[comm] > 0) { + cmd_full_help[comm] += "\nParameters\n"; + cmd_full_help[comm] += "----------\n\n"; + size_t i; + for (i = 0; i < cmd_n_args_min[comm]; i++) { + cmd_full_help[comm] += cmd_arghelp[comm][i]+"\n"; + } + for (i = cmd_n_args_min[comm]; i < cmd_n_args_max[comm]; i++) { + cmd_full_help[comm] += cmd_arghelp[comm][i]+" (optional)\n"; + } + } + if (cmd_rethelp[comm].size() > 0) { + cmd_full_help[comm] += "\nReturns\n"; + cmd_full_help[comm] += "-------\n\n"; + cmd_full_help[comm] += cmd_rethelp[comm]+"\n"; + } + cmd_fns[comm] = fn; if (cvm::debug()) { cvm::log("Defined command \""+std::string(name)+"\", with help string:\n"); - cvm::log(get_command_help(name)); + cvm::log(get_command_full_help(name)); } + return COLVARS_OK; } @@ -133,27 +186,76 @@ std::string colvarscript::get_cmd_prefix(colvarscript::Object_type t) } -std::string colvarscript::get_command_help(char const *cmd) + +char const *colvarscript::get_command_help(char const *cmd) { if (cmd_str_map.count(cmd) > 0) { colvarscript::command const c = cmd_str_map[std::string(cmd)]; - std::string new_result(cmd_help[c]+"\n"); - if (cmd_n_args_max[c] == 0) return new_result; - new_result += "\nParameters\n"; - new_result += "----------\n\n"; - size_t i; - for (i = 0; i < cmd_n_args_min[c]; i++) { - new_result += cmd_arghelp[c][i]+"\n"; - } - for (i = cmd_n_args_min[c]; i < cmd_n_args_max[c]; i++) { - new_result += cmd_arghelp[c][i]+" (optional)\n"; - } - return new_result; + return cmd_help[c].c_str(); } - cvm::error("Error: command "+std::string(cmd)+ " is not implemented.\n", INPUT_ERROR); - return std::string(""); + return NULL; +} + + +char const *colvarscript::get_command_rethelp(char const *cmd) +{ + if (cmd_str_map.count(cmd) > 0) { + colvarscript::command const c = cmd_str_map[std::string(cmd)]; + return cmd_rethelp[c].c_str(); + } + cvm::error("Error: command "+std::string(cmd)+ + " is not implemented.\n", INPUT_ERROR); + return NULL; +} + + +char const *colvarscript::get_command_arghelp(char const *cmd, int i) +{ + if (cmd_str_map.count(cmd) > 0) { + colvarscript::command const c = cmd_str_map[std::string(cmd)]; + return cmd_arghelp[c][i].c_str(); + } + cvm::error("Error: command "+std::string(cmd)+ + " is not implemented.\n", INPUT_ERROR); + return NULL; +} + + +int colvarscript::get_command_n_args_min(char const *cmd) +{ + if (cmd_str_map.count(cmd) > 0) { + colvarscript::command const c = cmd_str_map[std::string(cmd)]; + return cmd_n_args_min[c]; + } + cvm::error("Error: command "+std::string(cmd)+ + " is not implemented.\n", INPUT_ERROR); + return -1; +} + + +int colvarscript::get_command_n_args_max(char const *cmd) +{ + if (cmd_str_map.count(cmd) > 0) { + colvarscript::command const c = cmd_str_map[std::string(cmd)]; + return cmd_n_args_max[c]; + } + cvm::error("Error: command "+std::string(cmd)+ + " is not implemented.\n", INPUT_ERROR); + return -1; +} + + +char const *colvarscript::get_command_full_help(char const *cmd) +{ + if (cmd_str_map.count(cmd) > 0) { + colvarscript::command const c = cmd_str_map[std::string(cmd)]; + return cmd_full_help[c].c_str(); + } + cvm::error("Error: command "+std::string(cmd)+ + " is not implemented.\n", INPUT_ERROR); + return NULL; } @@ -234,7 +336,7 @@ std::string colvarscript::get_command_cmdline_help(colvarscript::Object_type t, if (cmd_str_map.count(cmdkey) > 0) { command const c = cmd_str_map[cmdkey]; return get_command_cmdline_syntax(t, c)+"\n\n"+ - get_command_help(cmd_names[c]); + get_command_full_help(cmd_names[c]); } cvm::error("Error: could not find scripting command \""+cmd+"\".", INPUT_ERROR); @@ -244,7 +346,7 @@ std::string colvarscript::get_command_cmdline_help(colvarscript::Object_type t, int colvarscript::run(int objc, unsigned char *const objv[]) { - result.clear(); + clear_str_result(); if (cvm::debug()) { cvm::log("Called script run with " + cvm::to_str(objc) + " args:"); @@ -346,6 +448,60 @@ int colvarscript::run(int objc, unsigned char *const objv[]) } +char *colvarscript::obj_to_str(unsigned char *obj) +{ + char *strobj = reinterpret_cast(obj); + if (cvm::debug()) { + cvm::log("Using simple-cast script::obj_to_str(): result = \"" + + (strobj ? std::string(strobj) : std::string("(null)")) + "\""); + } + return strobj; +} + + +std::vector colvarscript::obj_to_str_vector(unsigned char *obj) +{ + if (cvm::debug()) { + cvm::log("Using simple-cast colvarscript::obj_to_str_vector().\n"); + } + + std::vector new_result; + std::string const str(reinterpret_cast(obj)); + + // TODO get rid of this once colvarscript can handle both fix_modify and Tcl? + // LAMMPS has a nicer function in the utils class + + for (size_t i = 0; i < str.length(); i++) { + char const c = str[i]; + if (c == '\"') { + i++; + if (i >= str.length()) { + cvm::error("Error: could not split the following string:\n"+ + str+"\n", INPUT_ERROR); + break; + } + new_result.push_back(std::string("")); + while (str[i] != '\"') { + new_result.back().append(1, str[i]); + if (i >= str.length()) { + cvm::error("Error: could not split the following string:\n"+ + str+"\n", INPUT_ERROR); + break; + } else { + i++; + } + } + } + } + + if (cvm::debug()) { + cvm::log("result = "+cvm::to_str(new_result)+".\n"); + } + + return new_result; +} + + int colvarscript::proc_features(colvardeps *obj, int objc, unsigned char *const objv[]) { @@ -428,9 +584,9 @@ int colvarscript::set_result_str(std::string const &s) { if (cvm::get_error() != COLVARS_OK) { // Avoid overwriting the error message - result += s; + modify_str_result() += s; } else { - result = s; + modify_str_result() = s; } return COLVARS_OK; } @@ -438,17 +594,17 @@ int colvarscript::set_result_str(std::string const &s) void colvarscript::add_error_msg(std::string const &s) { - result += s; + modify_str_result() += s; // Ensure terminating newlines if (s[s.size()-1] != '\n') { - result += "\n"; + modify_str_result() += "\n"; } } int colvarscript::clear_str_result() { - result.clear(); + modify_str_result().clear(); return COLVARS_OK; } @@ -487,32 +643,69 @@ const char * get_colvarscript_result() int tcl_colvars_vmd_init(Tcl_Interp *interp, int molid); #endif -extern "C" -int tcl_run_colvarscript_command(ClientData /* clientData */, - Tcl_Interp *my_interp, - int objc, Tcl_Obj *const objv[]) +#if !defined(VMDTCL) && !defined(NAMD_TCL) +extern "C" { + int Colvars_Init(Tcl_Interp *interp) { + colvarproxy *proxy = new colvarproxy(); + colvarmodule *colvars = new colvarmodule(proxy); + proxy->set_tcl_interp(reinterpret_cast(interp)); + proxy->colvars = colvars; + proxy->script = new colvarscript(proxy); + Tcl_CreateObjCommand(interp, "cv", tcl_run_colvarscript_command, + (ClientData *) NULL, (Tcl_CmdDeleteProc *) NULL); + Tcl_EvalEx(interp, "package provide colvars", -1, 0); + return TCL_OK; + } +} +#endif + + +extern "C" int tcl_run_colvarscript_command(ClientData /* clientData */, + Tcl_Interp *my_interp, + int objc, Tcl_Obj *const objv[]) { colvarmodule *colvars = cvm::main(); if (!colvars) { #if defined(VMDTCL) + + if (objc == 2) { + if (!strcmp(Tcl_GetString(objv[1]), "molid")) { + // return invalid molid + Tcl_SetResult(my_interp, (char *) "-1", TCL_STATIC); + } + if (!strcmp(Tcl_GetString(objv[1]), "delete") || + !strcmp(Tcl_GetString(objv[1]), "reset")) { + // nothing to delete or reset + Tcl_SetResult(my_interp, NULL, TCL_STATIC); + } + if (!strcmp(Tcl_GetString(objv[1]), "help")) { + // print message + Tcl_SetResult(my_interp, + (char *) "First, setup the Colvars module with: " + "cv molid |top", TCL_STATIC); + } + return TCL_OK; + } + if (objc >= 3) { // require a molid to create the module if (!strcmp(Tcl_GetString(objv[1]), "molid")) { - int molid = -1; + int molid = -(1<<16); // This value is used to indicate "top" if (strcmp(Tcl_GetString(objv[2]), "top")) { // If this is not "top", get the integer value Tcl_GetIntFromObj(my_interp, objv[2], &molid); } return tcl_colvars_vmd_init(my_interp, molid); } else { - // TODO allow calling cv help after this - Tcl_SetResult(my_interp, (char *) "Syntax error.", TCL_STATIC); + Tcl_SetResult(my_interp, (char *) "Syntax error. First, setup the Colvars module with cv molid |top", TCL_STATIC); return TCL_ERROR; } } + Tcl_SetResult(my_interp, (char *) "First, setup the Colvars module with: " - "cv molid ", TCL_STATIC); + "cv molid |top", TCL_STATIC); + #else Tcl_SetResult(my_interp, const_cast("Error: Colvars module not yet initialized"), @@ -534,10 +727,19 @@ int tcl_run_colvarscript_command(ClientData /* clientData */, cvm::clear_error(); - int retval = script->run(objc, - reinterpret_cast(objv)); + unsigned char * arg_pointers_[100]; + if (objc > 100) { + std::string const errstr = "Too many positional arguments ("+ + cvm::to_str(objc)+") passed to the \"cv\" command.\n"; + Tcl_SetResult(interp, const_cast(errstr.c_str()), TCL_VOLATILE); + return TCL_ERROR; + } + for (int i = 0; i < objc; i++) { + arg_pointers_[i] = reinterpret_cast(const_cast(proxy->tcl_get_str(objv[i]))); + } + int retval = script->run(objc, arg_pointers_); - std::string result = proxy->get_error_msgs() + script->result; + std::string result = proxy->get_error_msgs() + script->str_result(); Tcl_SetResult(interp, const_cast(result.c_str()), TCL_VOLATILE); @@ -558,3 +760,162 @@ int tcl_run_colvarscript_command(ClientData /* clientData */, } #endif // #if defined(COLVARS_TCL) + + + + +int colvarscript::set_result_text_from_str(std::string const &x_str, + unsigned char *obj) { + if (obj) { + strcpy(reinterpret_cast(obj), x_str.c_str()); + } else { + set_result_str(x_str); + } + return COLVARS_OK; +} + +// Template to convert everything to string and use the above + +template +int colvarscript::set_result_text(T const &x, unsigned char *obj) { + std::string const x_str = x.to_simple_string(); + return set_result_text_from_str(x_str, obj); +} + + +template +int colvarscript::pack_vector_elements_text(std::vector const &x, + std::string &x_str) { + x_str.clear(); + for (size_t i = 0; i < x.size(); ++i) { + if (i > 0) x_str.append(1, ' '); + x_str += cvm::to_str(x[i]); + } + return COLVARS_OK; +} + + +// Specializations for plain old data types that don't have a stringifier member + +template <> +int colvarscript::set_result_text(int const &x, unsigned char *obj) { + std::string const x_str = cvm::to_str(x); + return set_result_text_from_str(x_str, obj); +} + +template <> +int colvarscript::set_result_text(std::vector const &x, + unsigned char *obj) { + std::string x_str(""); + pack_vector_elements_text(x, x_str); + return set_result_text_from_str(x_str, obj); +} + + +template <> +int colvarscript::set_result_text(long int const &x, unsigned char *obj) { + std::string const x_str = cvm::to_str(x); + return set_result_text_from_str(x_str, obj); +} + +template <> +int colvarscript::set_result_text(std::vector const &x, + unsigned char *obj) { + std::string x_str(""); + pack_vector_elements_text(x, x_str); + return set_result_text_from_str(x_str, obj); +} + + +template <> +int colvarscript::set_result_text(cvm::real const &x, unsigned char *obj) { + std::string const x_str = cvm::to_str(x); + return set_result_text_from_str(x_str, obj); +} + +template <> +int colvarscript::set_result_text(std::vector const &x, + unsigned char *obj) { + std::string x_str(""); + pack_vector_elements_text(x, x_str); + return set_result_text_from_str(x_str, obj); +} + + +// TODO these can be removed after the Tcl backend is ready (otherwise, the +// default template syntax may break scripts or the Dashboard) + +template <> +int colvarscript::set_result_text(std::vector const &x, + unsigned char *obj) { + std::string x_str(""); + for (size_t i = 0; i < x.size(); i++) { + if (i > 0) x_str.append(1, ' '); + x_str += "{ "+x[i].to_simple_string()+" }"; + } + return set_result_text_from_str(x_str, obj); +} + +template <> +int colvarscript::set_result_text(std::vector const &x, + unsigned char *obj) { + std::string x_str(""); + for (size_t i = 0; i < x.size(); i++) { + if (i > 0) x_str.append(1, ' '); + x_str += "{ "+x[i].to_simple_string()+" }"; + } + return set_result_text_from_str(x_str, obj); +} + + +// Member functions to set script results for each typexc + +int colvarscript::set_result_int(int const &x, unsigned char *obj) { + return set_result_text(x, obj); +} + +int colvarscript::set_result_int_vec(std::vector const &x, + unsigned char *obj) { + return set_result_text< std::vector >(x, obj); +} + + +int colvarscript::set_result_long_int(long int const &x, unsigned char *obj) { + return set_result_text(x, obj); +} + +int colvarscript::set_result_long_int_vec(std::vector const &x, + unsigned char *obj) { + return set_result_text< std::vector >(x, obj); +} + + +int colvarscript::set_result_real(cvm::real const &x, unsigned char *obj) { + return set_result_text(x, obj); +} + +int colvarscript::set_result_real_vec(std::vector const &x, + unsigned char *obj) { + return set_result_text< std::vector >(x, obj); +} + + +int colvarscript::set_result_rvector(cvm::rvector const &x, unsigned char *obj) { + return set_result_text(x, obj); +} + +int colvarscript::set_result_rvector_vec(std::vector const &x, + unsigned char *obj) { + return set_result_text< std::vector >(x, obj); +} + + +int colvarscript::set_result_colvarvalue(colvarvalue const &x, + unsigned char *obj) { + return set_result_text(x, obj); +} + +int colvarscript::set_result_colvarvalue_vec(std::vector const &x, + unsigned char *obj) { + return set_result_text< std::vector >(x, obj); +} diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h index d6f77668e6..7eac15c09d 100644 --- a/lib/colvars/colvarscript.h +++ b/lib/colvars/colvarscript.h @@ -46,9 +46,8 @@ public: /// COLVARSCRIPT_ERROR int proxy_error; - /// If an error is returned by one of the methods, it should set this to the - /// error message - std::string result; + /// String representation of the result of a script call + std::string str_result_; /// Run a script command with space-separated positional arguments (objects) int run(int objc, unsigned char *const objv[]); @@ -56,13 +55,13 @@ public: /// Get the string result of the current scripting call inline std::string const &str_result() const { - return result; + return str_result_; } /// Modify the string result of the current scripting call inline std::string &modify_str_result() { - return result; + return str_result_; } /// Set the return value to the given string @@ -137,21 +136,38 @@ public: template int cmd_arg_shift(); - /// Use scripting language to get the string representation of an object - inline char const *obj_to_str(unsigned char *const obj) - { - return (obj == NULL ? NULL : proxy_->script_obj_to_str(obj)); - } - /// Get names of all commands inline char const **get_command_names() const { return cmd_names; } + /// Get one-line help summary for a command + /// \param cmd Name of the command's function (e.g. "cv_units") + char const *get_command_help(char const *cmd); + + /// Get description of the return value of a command + /// \param cmd Name of the command's function (e.g. "cv_units") + char const *get_command_rethelp(char const *cmd); + + /// Get description of the argument of a command (excluding prefix) + /// \param cmd Name of the command's function (e.g. "cv_units") + /// \param i Index of the argument; 0 is the first argument after the + /// prefix, e.g. "value" has an index of 0 in the array of arguments: + /// { "cv", "colvar", "xi", "value" } + char const *get_command_arghelp(char const *cmd, int i); + + /// Get number of required arguments (excluding prefix) + /// \param cmd Name of the command's function (e.g. "cv_units") + int get_command_n_args_min(char const *cmd); + + /// Get number of total arguments (excluding prefix) + /// \param cmd Name of the command's function (e.g. "cv_units") + int get_command_n_args_max(char const *cmd); + /// Get help string for a command (does not specify how it is launched) /// \param cmd Name of the command's function (e.g. "cv_units") - std::string get_command_help(char const *cmd); + char const *get_command_full_help(char const *cmd); /// Get summary of command line syntax for all commands of a given context /// \param t One of use_module, use_colvar or use_bias @@ -182,6 +198,53 @@ public: return this->proxy_; } + // Input functions - get the string reps of script argument objects + + /// Get the string representation of an object (by default, a simple cast) + char *obj_to_str(unsigned char *obj); + + /// Get a list of strings from an object (does not work with a simple cast) + std::vector obj_to_str_vector(unsigned char *obj); + + + // Output functions - convert internal objects to representations suitable + // for use in the scripting language. At the moment only conversion to C + // strings is supported, and obj is assumed to be a char * pointer. + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_int(int const &x, unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_int_vec(std::vector const &x, unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_long_int(long int const &x, unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_long_int_vec(std::vector const &x, + unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_real(cvm::real const &x, unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_real_vec(std::vector const &x, + unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_rvector(cvm::rvector const &x, unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_rvector_vec(std::vector const &x, + unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_colvarvalue(colvarvalue const &x, unsigned char *obj = NULL); + + /// Copy x into obj if not NULL, or into the script object's result otherwise + int set_result_colvarvalue_vec(std::vector const &x, + unsigned char *obj = NULL); + private: /// Set up all script API functions @@ -193,14 +256,6 @@ private: int n_args_min, int n_args_max, char const *arghelp, int (*fn)(void *, int, unsigned char * const *)); - /// Execute a script command - inline int exec_command(command c, - void *pobj, - int objc, unsigned char * const *objv) - { - return (*(cmd_fns[c]))(pobj, objc, objv); - } - public: // TODO this function will be removed soon /// Run subcommands on base colvardeps object (colvar, bias, ...) @@ -218,6 +273,9 @@ private: // TODO /// Help strings for each command std::vector cmd_help; + /// Description of the return values of each command (may be empty) + std::vector cmd_rethelp; + /// Minimum number of arguments for each command std::vector cmd_n_args_min; @@ -227,6 +285,9 @@ private: // TODO /// Help strings for each command argument std::vector< std::vector > cmd_arghelp; + /// Full help strings for each command + std::vector cmd_full_help; + /// Implementations of each command std::vector cmd_fns; @@ -241,6 +302,18 @@ private: // TODO return NULL; } + /// Set obj equal to x, using its string representation + template + int set_result_text(T const &x, unsigned char *obj); + + /// Code reused by instances of set_result_text() + template + int pack_vector_elements_text(std::vector const &x, std::string &x_str); + + /// Code reused by all instances of set_result_text() + int set_result_text_from_str(std::string const &x_str, unsigned char *obj); + + }; @@ -305,13 +378,15 @@ int colvarscript::check_cmd_nargs(char const *cmd, { int const shift = cmd_arg_shift(); if (objc < shift+n_args_min) { - add_error_msg("Missing arguments for script function \""+std::string(cmd)+ - "\":\n"+get_command_help(cmd)); + add_error_msg("Insufficient number of arguments ("+cvm::to_str(objc)+ + ") for script function \""+std::string(cmd)+ + "\":\n"+get_command_full_help(cmd)); return COLVARSCRIPT_ERROR; } if (objc > shift+n_args_max) { - add_error_msg("Too many arguments for script function \""+std::string(cmd)+ - "\":\n"+get_command_help(cmd)); + add_error_msg("Too many arguments ("+cvm::to_str(objc)+ + ") for script function \""+std::string(cmd)+ + "\":\n"+get_command_full_help(cmd)); return COLVARSCRIPT_ERROR; } return COLVARSCRIPT_OK; @@ -364,18 +439,6 @@ int colvarscript::cmd_arg_shift() extern "C" { -#if defined(COLVARS_TCL) - /// Run the script API via Tcl command-line interface - /// \param clientData Not used - /// \param my_interp Pointer to Tcl_Interp object (read from Colvars if NULL) - /// \param objc Number of Tcl command parameters - /// \param objv Array of command parameters - /// \return Result of the script command - int tcl_run_colvarscript_command(ClientData clientData, - Tcl_Interp *interp_in, - int objc, Tcl_Obj *const objv[]); -#endif // #if defined(COLVARS_TCL) - /// Generic wrapper for string-based scripting int run_colvarscript_command(int objc, unsigned char *const objv[]); diff --git a/lib/colvars/colvarscript_commands.cpp b/lib/colvars/colvarscript_commands.cpp index c74663d2fd..0029979912 100644 --- a/lib/colvars/colvarscript_commands.cpp +++ b/lib/colvars/colvarscript_commands.cpp @@ -33,6 +33,54 @@ char const **cvscript_command_names() } +extern "C" +char const *cvscript_command_help(char const *c) +{ + colvarscript *script = colvarscript_obj(); + return script->get_command_help(c); +} + + +extern "C" +char const *cvscript_command_rethelp(char const *c) +{ + colvarscript *script = colvarscript_obj(); + return script->get_command_rethelp(c); +} + + +extern "C" +char const *cvscript_command_arghelp(char const *c, int i) +{ + colvarscript *script = colvarscript_obj(); + return script->get_command_arghelp(c, i); +} + + +extern "C" +char const *cvscript_command_full_help(char const *c) +{ + colvarscript *script = colvarscript_obj(); + return script->get_command_full_help(c); +} + + +extern "C" +int cvscript_command_n_args_min(char const *c) +{ + colvarscript *script = colvarscript_obj(); + return script->get_command_n_args_min(c); +} + + +extern "C" +int cvscript_command_n_args_max(char const *c) +{ + colvarscript *script = colvarscript_obj(); + return script->get_command_n_args_max(c); +} + + // Instantiate the body of all script commands #define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \ diff --git a/lib/colvars/colvarscript_commands.h b/lib/colvars/colvarscript_commands.h index d90a3ac422..021aae1cc6 100644 --- a/lib/colvars/colvarscript_commands.h +++ b/lib/colvars/colvarscript_commands.h @@ -25,7 +25,8 @@ // COMM = the id of the command (must be a member of colvarscript::command) -// HELP = a one-line description (C string literal) for the command +// HELP = short description (C string literal) for the command; the second line +// is optional, and documents the return value (if any) // N_ARGS_MIN = the lowest number of arguments allowed @@ -68,6 +69,33 @@ extern "C" { /// Get the names of all commands (array of strings) char const ** cvscript_command_names(); + /// Get the help summary of the given command + /// \param cmd Name of the command's function (e.g. "cv_units") + char const *cvscript_command_help(char const *cmd); + + /// Get description of the return value of a command + /// \param cmd Name of the command's function (e.g. "cv_units") + char const *cvscript_command_rethelp(char const *cmd); + + /// Get description of the arguments of a command (excluding prefix) + /// \param cmd Name of the command's function (e.g. "cv_units") + /// \param i Index of the argument; 0 is the first argument after the + /// prefix, e.g. "value" has an index of 0 in the array of arguments: + /// { "cv", "colvar", "xi", "value" } + char const *cvscript_command_arghelp(char const *cmd, int i); + + /// Get the full help string of a command + /// \param cmd Name of the command's function (e.g. "cv_units") + char const *cvscript_command_full_help(char const *cmd); + + /// Get number of required arguments (excluding prefix) + /// \param cmd Name of the command's function (e.g. "cv_units") + int cvscript_command_n_args_min(char const *cmd); + + /// Get number of total arguments (excluding prefix) + /// \param cmd Name of the command's function (e.g. "cv_units") + int cvscript_command_n_args_max(char const *cmd); + } #endif @@ -135,7 +163,8 @@ CVSCRIPT(cv_delete, ) CVSCRIPT(cv_frame, - "Get or set current frame number (VMD only)", + "Get or set current frame number (VMD only)\n" + "frame : integer - Frame number", 0, 1, "frame : integer - Frame number", char const *arg = @@ -143,7 +172,7 @@ CVSCRIPT(cv_frame, if (arg == NULL) { long int f = -1; if (script->proxy()->get_frame(f) == COLVARS_OK) { - script->set_result_str(cvm::to_str(f)); + script->set_result_long_int(f); return COLVARS_OK; } else { script->add_error_msg("Frame number is not available"); @@ -161,8 +190,90 @@ CVSCRIPT(cv_frame, return COLVARS_OK; ) +CVSCRIPT(cv_getatomappliedforces, + "Get the list of forces applied by Colvars to atoms\n" + "forces : array of arrays of floats - Atomic forces", + 0, 0, + "", + script->set_result_rvector_vec(*(script->proxy()->get_atom_applied_forces())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatomappliedforcesmax, + "Get the maximum norm of forces applied by Colvars to atoms\n" + "force : float - Maximum atomic force", + 0, 0, + "", + script->set_result_real(script->proxy()->max_atoms_applied_force()); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatomappliedforcesmaxid, + "Get the atom ID with the largest applied force\n" + "id : int - ID of the atom with the maximum atomic force", + 0, 0, + "", + script->set_result_int(script->proxy()->max_atoms_applied_force_id()); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatomappliedforcesrms, + "Get the root-mean-square norm of forces applied by Colvars to atoms\n" + "force : float - RMS atomic force", + 0, 0, + "", + script->set_result_real(script->proxy()->rms_atoms_applied_force()); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatomids, + "Get the list of indices of atoms used in Colvars\n" + "indices : array of ints - Atomic indices", + 0, 0, + "", + script->set_result_int_vec(*(script->proxy()->get_atom_ids())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatomcharges, + "Get the list of charges of atoms used in Colvars\n" + "charges : array of floats - Atomic charges", + 0, 0, + "", + script->set_result_real_vec(*(script->proxy()->get_atom_charges())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatommasses, + "Get the list of masses of atoms used in Colvars\n" + "masses : array of floats - Atomic masses", + 0, 0, + "", + script->set_result_real_vec(*(script->proxy()->get_atom_masses())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatompositions, + "Get the list of cached positions of atoms used in Colvars\n" + "positions : array of arrays of floats - Atomic positions", + 0, 0, + "", + script->set_result_rvector_vec(*(script->proxy()->get_atom_positions())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getatomtotalforces, + "Get the list of cached total forces of atoms used in Colvars\n" + "forces : array of arrays of floats - Atomic total foces", + 0, 0, + "", + script->set_result_rvector_vec(*(script->proxy()->get_atom_total_forces())); + return COLVARS_OK; + ) + CVSCRIPT(cv_getconfig, - "Get the module's configuration string read so far", + "Get the module's configuration string read so far\n" + "conf : string - Current configuration string", 0, 0, "", script->set_result_str(cvm::main()->get_config()); @@ -170,15 +281,17 @@ CVSCRIPT(cv_getconfig, ) CVSCRIPT(cv_getenergy, - "Get the current Colvars energy", + "Get the current Colvars energy\n" + "E : float - Amount of energy (internal units)", 0, 0, "", - script->set_result_str(cvm::to_str(cvm::main()->total_bias_energy)); + script->set_result_real(cvm::main()->total_bias_energy); return COLVARS_OK; ) CVSCRIPT(cv_help, - "Get the help string of the Colvars scripting interface", + "Get the help string of the Colvars scripting interface\n" + "help : string - Help string", 0, 1, "command : string - Get the help string of this specific command", unsigned char *const cmdobj = @@ -205,7 +318,8 @@ CVSCRIPT(cv_help, ) CVSCRIPT(cv_list, - "Return a list of all variables or biases", + "Return a list of all variables or biases\n" + "list : sequence of strings - List of elements", 0, 1, "param : string - \"colvars\" or \"biases\"; default is \"colvars\"", std::string res; @@ -235,7 +349,8 @@ CVSCRIPT(cv_list, ) CVSCRIPT(cv_listcommands, - "Get the list of script functions, prefixed with \"cv_\", \"colvar_\" or \"bias_\"", + "Get the list of script functions, prefixed with \"cv_\", \"colvar_\" or \"bias_\"\n" + "list : sequence of strings - List of commands", 0, 0, "", int const n_commands = cvscript_n_commands(); @@ -249,6 +364,20 @@ CVSCRIPT(cv_listcommands, return COLVARS_OK; ) +CVSCRIPT(cv_listindexfiles, + "Get a list of the index files loaded in this session", + 0, 0, + "", + int const n_files = script->module()->index_file_names.size(); + std::string result; + for (int i = 0; i < n_files; i++) { + if (i > 0) result.append(1, ' '); + result.append(script->module()->index_file_names[i]); + } + script->set_result_str(result); + return COLVARS_OK; + ) + CVSCRIPT(cv_load, "Load data from a state file into all matching colvars and biases", 1, 1, @@ -280,15 +409,16 @@ CVSCRIPT(cv_loadfromstring, ) CVSCRIPT(cv_molid, - "Get or set the molecule ID on which Colvars is defined (VMD only)", + "Get or set the molecule ID on which Colvars is defined (VMD only)\n" + "molid : integer - Current molecule ID", 0, 1, - "molid : integer - Molecule ID; -1 means undefined", + "molid : integer - New molecule ID; -1 means undefined", char const *arg = script->obj_to_str(script->get_module_cmd_arg(0, objc, objv)); if (arg == NULL) { int molid = -1; script->proxy()->get_molid(molid); - script->set_result_str(cvm::to_str(molid)); + script->set_result_int(molid); return COLVARS_OK; } else { script->add_error_msg("Error: To change the molecule ID in VMD, use cv delete first."); @@ -297,7 +427,8 @@ CVSCRIPT(cv_molid, ) CVSCRIPT(cv_printframe, - "Return the values that would be written to colvars.traj", + "Return the values that would be written to colvars.traj\n" + "values : string - The values\n", 0, 0, "", std::ostringstream os; @@ -307,7 +438,8 @@ CVSCRIPT(cv_printframe, ) CVSCRIPT(cv_printframelabels, - "Return the labels that would be written to colvars.traj", + "Return the labels that would be written to colvars.traj\n" + "Labels : string - The labels", 0, 0, "", std::ostringstream os; @@ -348,14 +480,16 @@ CVSCRIPT(cv_save, ) CVSCRIPT(cv_savetostring, - "Write the Colvars state to a string and return it", + "Write the Colvars state to a string and return it\n" + "state : string - The saved state", 0, 0, "", return script->module()->write_restart_string(script->modify_str_result()); ) CVSCRIPT(cv_units, - "Get or set the current Colvars unit system", + "Get or set the current Colvars unit system\n" + "units : string - The current unit system", 0, 1, "units : string - The new unit system", char const *argstr = @@ -390,7 +524,8 @@ CVSCRIPT(cv_update, ) CVSCRIPT(cv_version, - "Get the Colvars Module version number", + "Get the Colvars Module version number\n" + "version : string - Colvars version", 0, 0, "", script->set_result_str(COLVARS_VERSION); diff --git a/lib/colvars/colvarscript_commands_bias.h b/lib/colvars/colvarscript_commands_bias.h index 990902e239..f83b3422ff 100644 --- a/lib/colvars/colvarscript_commands_bias.h +++ b/lib/colvars/colvarscript_commands_bias.h @@ -9,15 +9,17 @@ CVSCRIPT(bias_bin, - "Get the current grid bin index (1D ABF only for now)", + "Get the current grid bin index (1D ABF only for now)\n" + "bin : integer - Bin index", 0, 0, "", - script->set_result_str(cvm::to_str(this_bias->current_bin())); + script->set_result_int(this_bias->current_bin()); return COLVARS_OK; ) CVSCRIPT(bias_bincount, - "Get the number of samples at the given grid bin (1D ABF only for now)", + "Get the number of samples at the given grid bin (1D ABF only for now)\n" + "samples : integer - Number of samples", 0, 1, "index : integer - Grid index; defaults to current bin", int index = this_bias->current_bin(); @@ -30,12 +32,13 @@ CVSCRIPT(bias_bincount, return COLVARSCRIPT_ERROR; } } - script->set_result_str(cvm::to_str(this_bias->bin_count(index))); + script->set_result_int(this_bias->bin_count(index)); return COLVARS_OK; ) CVSCRIPT(bias_binnum, - "Get the total number of grid points of this bias (1D ABF only for now)", + "Get the total number of grid points of this bias (1D ABF only for now)\n" + "Bins : integer - Number of grid points", 0, 0, "", int r = this_bias->bin_num(); @@ -44,7 +47,7 @@ CVSCRIPT(bias_binnum, this_bias->name); return COLVARSCRIPT_ERROR; } - script->set_result_str(cvm::to_str(r)); + script->set_result_int(r); return COLVARS_OK; ) @@ -57,22 +60,25 @@ CVSCRIPT(bias_delete, ) CVSCRIPT(bias_energy, - "Get the current energy of this bias", + "Get the current energy of this bias\n" + "E : float - Energy value", 0, 0, "", - script->set_result_str(cvm::to_str(this_bias->get_energy())); + script->set_result_real(this_bias->get_energy()); return COLVARS_OK; ) CVSCRIPT(bias_get, - "Get the value of the given feature for this bias", + "Get the value of the given feature for this bias\n" + "state : 1/0 - State of the given feature", 1, 1, "feature : string - Name of the feature", return script->proc_features(this_bias, objc, objv); ) CVSCRIPT(bias_getconfig, - "Return the configuration string of this bias", + "Return the configuration string of this bias\n" + "conf : string - Current configuration string", 0, 0, "", script->set_result_str(this_bias->get_config()); @@ -80,7 +86,8 @@ CVSCRIPT(bias_getconfig, ) CVSCRIPT(bias_help, - "Get a help summary or the help string of one bias subcommand", + "Get a help summary or the help string of one bias subcommand\n" + "help : string - Help string", 0, 1, "command : string - Get the help string of this specific command", unsigned char *const cmdobj = @@ -129,7 +136,8 @@ CVSCRIPT(bias_save, ) CVSCRIPT(bias_savetostring, - "Save data from this bias into a string and return it", + "Save data from this bias into a string and return it\n" + "state : string - The bias state", 0, 0, "", return this_bias->write_state_string(script->modify_str_result()); @@ -156,7 +164,8 @@ CVSCRIPT(bias_share, ) CVSCRIPT(bias_state, - "Print a string representation of the feature state of this bias", + "Print a string representation of the feature state of this bias\n" + "state : string - String representation of the bias features", 0, 0, "", this_bias->print_state(); @@ -164,10 +173,11 @@ CVSCRIPT(bias_state, ) CVSCRIPT(bias_update, - "Recompute this bias and return its up-to-date energy", + "Recompute this bias and return its up-to-date energy\n" + "E : float - Energy value", 0, 0, "", this_bias->update(); - script->set_result_str(cvm::to_str(this_bias->get_energy())); + script->set_result_colvarvalue(this_bias->get_energy()); return COLVARS_OK; ) diff --git a/lib/colvars/colvarscript_commands_colvar.h b/lib/colvars/colvarscript_commands_colvar.h index b880b5b8be..7c4c2d67f7 100644 --- a/lib/colvars/colvarscript_commands_colvar.h +++ b/lib/colvars/colvarscript_commands_colvar.h @@ -9,7 +9,8 @@ CVSCRIPT(colvar_addforce, - "Apply the given force onto this colvar and return the same", + "Apply the given force onto this colvar and return the same\n" + "force : float or array - Applied force; matches colvar dimensionality", 1, 1, "force : float or array - Applied force; must match colvar dimensionality", std::string const f_str(script->obj_to_str(script->get_colvar_cmd_arg(0, objc, objv))); @@ -23,7 +24,7 @@ CVSCRIPT(colvar_addforce, return COLVARSCRIPT_ERROR; } this_colvar->add_bias_force(force); - script->set_result_str(force.to_simple_string()); + script->set_result_colvarvalue(force); return COLVARS_OK; ) @@ -56,22 +57,25 @@ CVSCRIPT(colvar_delete, ) CVSCRIPT(colvar_get, - "Get the value of the given feature for this colvar", + "Get the value of the given feature for this colvar\n" + "state : 1/0 - State of the given feature", 1, 1, "feature : string - Name of the feature", return script->proc_features(this_colvar, objc, objv); ) CVSCRIPT(colvar_getappliedforce, - "Return the total of the forces applied to this colvar", + "Return the total of the forces applied to this colvar\n" + "force : float - Applied force; matches the colvar dimensionality", 0, 0, "", - script->set_result_str((this_colvar->applied_force()).to_simple_string()); + script->set_result_colvarvalue(this_colvar->applied_force()); return COLVARS_OK; ) CVSCRIPT(colvar_getatomgroups, - "Return the atom indices used by this colvar as a list of lists", + "Return the atom indices used by this colvar as a list of lists\n" + "groups : array of arrays of ints - Atom indices", 0, 0, "", std::string result; @@ -91,21 +95,17 @@ CVSCRIPT(colvar_getatomgroups, ) CVSCRIPT(colvar_getatomids, - "Return the list of atom indices used by this colvar", + "Return the list of atom indices used by this colvar\n" + "indices : array of ints - Atom indices", 0, 0, "", - std::string result; - std::vector::iterator li = this_colvar->atom_ids.begin(); - for ( ; li != this_colvar->atom_ids.end(); ++li) { - result += cvm::to_str(*li); - result += " "; - } - script->set_result_str(result); + script->set_result_int_vec(this_colvar->atom_ids); return COLVARS_OK; ) CVSCRIPT(colvar_getconfig, - "Return the configuration string of this colvar", + "Return the configuration string of this colvar\n" + "conf : string - Current configuration string", 0, 0, "", script->set_result_str(this_colvar->get_config()); @@ -113,35 +113,34 @@ CVSCRIPT(colvar_getconfig, ) CVSCRIPT(colvar_getgradients, - "Return the atomic gradients of this colvar", + "Return the atomic gradients of this colvar\n" + "gradients : array of arrays of floats - Atomic gradients", 0, 0, "", - std::string result; - std::vector::iterator li = - this_colvar->atomic_gradients.begin(); - for ( ; li != this_colvar->atomic_gradients.end(); ++li) { - result += "{"; - int j; - for (j = 0; j < 3; ++j) { - result += cvm::to_str((*li)[j]); - result += " "; - } - result += "} "; - } - script->set_result_str(result); + script->set_result_rvector_vec(this_colvar->atomic_gradients); return COLVARS_OK; ) CVSCRIPT(colvar_gettotalforce, - "Return the sum of internal and external forces to this colvar", + "Return the sum of internal and external forces to this colvar\n" + "force : float - Total force; matches the colvar dimensionality", 0, 0, "", - script->set_result_str((this_colvar->total_force()).to_simple_string()); + script->set_result_colvarvalue(this_colvar->total_force()); + return COLVARS_OK; + ) + +CVSCRIPT(colvar_getvolmapids, + "Return the list of volumetric map indices used by this colvar", + 0, 0, + "", + script->set_result_int_vec(this_colvar->get_volmap_ids()); return COLVARS_OK; ) CVSCRIPT(colvar_help, - "Get a help summary or the help string of one colvar subcommand", + "Get a help summary or the help string of one colvar subcommand\n" + "help : string - Help string", 0, 1, "command : string - Get the help string of this specific command", unsigned char *const cmdobj = @@ -167,7 +166,7 @@ CVSCRIPT(colvar_modifycvcs, "Modify configuration of individual components by passing string arguments", 1, 1, "confs : sequence of strings - New configurations; empty strings are skipped", - std::vector const confs(script->proxy()->script_obj_to_str_vector(script->get_colvar_cmd_arg(0, objc, objv))); + std::vector const confs(script->obj_to_str_vector(script->get_colvar_cmd_arg(0, objc, objv))); cvm::increase_depth(); int res = this_colvar->update_cvc_config(confs); cvm::decrease_depth(); @@ -180,10 +179,11 @@ CVSCRIPT(colvar_modifycvcs, ) CVSCRIPT(colvar_run_ave, - "Get the current running average of the value of this colvar", + "Get the current running average of the value of this colvar\n" + "value : float or array - Averaged value; matches the colvar dimensionality", 0, 0, "", - script->set_result_str(this_colvar->run_ave().to_simple_string()); + script->set_result_colvarvalue(this_colvar->run_ave()); return COLVARS_OK; ) @@ -196,7 +196,8 @@ CVSCRIPT(colvar_set, ) CVSCRIPT(colvar_state, - "Print a string representation of the feature state of this colvar", + "Print a string representation of the feature state of this colvar\n" + "state : string - The feature state", 0, 0, "", this_colvar->print_state(); @@ -204,7 +205,8 @@ CVSCRIPT(colvar_state, ) CVSCRIPT(colvar_type, - "Get the type description of this colvar", + "Get the type description of this colvar\n" + "type : string - Type description", 0, 0, "", script->set_result_str(this_colvar->value().type_desc(this_colvar->value().value_type)); @@ -212,25 +214,28 @@ CVSCRIPT(colvar_type, ) CVSCRIPT(colvar_update, - "Recompute this colvar and return its up-to-date value", + "Recompute this colvar and return its up-to-date value\n" + "value : float or array - Current value; matches the colvar dimensionality", 0, 0, "", this_colvar->calc(); this_colvar->update_forces_energy(); - script->set_result_str((this_colvar->value()).to_simple_string()); + script->set_result_colvarvalue(this_colvar->value()); return COLVARS_OK; ) CVSCRIPT(colvar_value, - "Get the current value of this colvar", + "Get the current value of this colvar\n" + "value : float or array - Current value; matches the colvar dimensionality", 0, 0, "", - script->set_result_str(this_colvar->value().to_simple_string()); + script->set_result_colvarvalue(this_colvar->value()); return COLVARS_OK; ) CVSCRIPT(colvar_width, - "Get the width of this colvar", + "Get the width of this colvar\n" + "width : float - Value of the width", 0, 0, "", script->set_result_str(cvm::to_str(this_colvar->width, 0, diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index 3cbaed63a8..4117a7a68f 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -256,6 +256,7 @@ namespace { colvarmodule::rotation::rotation() { b_debug_gradients = false; + lambda = 0.0; #ifdef COLVARS_LAMMPS jacobi = new_Jacobi_solver(4); #else @@ -268,6 +269,7 @@ colvarmodule::rotation::rotation(cvm::quaternion const &qi) : q(qi) { b_debug_gradients = false; + lambda = 0.0; #ifdef COLVARS_LAMMPS jacobi = new_Jacobi_solver(4); #else @@ -283,6 +285,7 @@ colvarmodule::rotation::rotation(cvm::real angle, cvm::rvector const &axis) cvm::real const sina = cvm::sin(angle/2.0); q = cvm::quaternion(cvm::cos(angle/2.0), sina * axis_n.x, sina * axis_n.y, sina * axis_n.z); + lambda = 0.0; #ifdef COLVARS_LAMMPS jacobi = new_Jacobi_solver(4); #else diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index 7ab617bc44..24a2790f6e 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -224,7 +224,7 @@ void colvarvalue::is_derivative() colvarvalue::colvarvalue(colvarvalue const &x) - : value_type(x.type()) + : value_type(x.type()), real_value(0.0) { switch (x.type()) { case type_scalar: @@ -252,6 +252,7 @@ colvarvalue::colvarvalue(colvarvalue const &x) colvarvalue::colvarvalue(cvm::vector1d const &v, Type vti) + : real_value(0.0) { if ((vti != type_vector) && (v.size() != num_dimensions(vti))) { cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+ @@ -579,16 +580,7 @@ colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const cvm::rvector const &v1 = this->rvector_value; cvm::rvector const &v2 = x2.rvector_value; cvm::real const cos_t = v1 * v2; - cvm::real const sin_t = cvm::sqrt(1.0 - cos_t*cos_t); - return colvarvalue( 2.0 * sin_t * - cvm::rvector((-1.0) * sin_t * v2.x + - cos_t/sin_t * (v1.x - cos_t*v2.x), - (-1.0) * sin_t * v2.y + - cos_t/sin_t * (v1.y - cos_t*v2.y), - (-1.0) * sin_t * v2.z + - cos_t/sin_t * (v1.z - cos_t*v2.z) - ), - colvarvalue::type_unit3vectorderiv ); + return colvarvalue(2.0 * (cos_t * v1 - v2), colvarvalue::type_unit3vectorderiv); } case colvarvalue::type_quaternion: case colvarvalue::type_quaternionderiv: diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index ca367dd43f..3f26c35df4 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -124,7 +124,7 @@ public: /// Constructor from a type specification inline colvarvalue(Type const &vti) - : value_type(vti) + : value_type(vti), real_value(0.0) { reset(); } @@ -138,12 +138,12 @@ public: /// by default a type \link type_3vector \endlink , if you want a /// \link type_unit3vector \endlink you must set it explicitly) inline colvarvalue(cvm::rvector const &v, Type vti = type_3vector) - : value_type(vti), rvector_value(v) + : value_type(vti), real_value(0.0), rvector_value(v) {} /// \brief Copy constructor from quaternion base type inline colvarvalue(cvm::quaternion const &q, Type vti = type_quaternion) - : value_type(vti), quaternion_value(q) + : value_type(vti), real_value(0.0), quaternion_value(q) {} /// Copy constructor from vector1d base type diff --git a/src/COLVARS/colvarproxy_lammps_version.h b/src/COLVARS/colvarproxy_lammps_version.h index d9f2955233..0399595533 100644 --- a/src/COLVARS/colvarproxy_lammps_version.h +++ b/src/COLVARS/colvarproxy_lammps_version.h @@ -1,3 +1,3 @@ #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2020-04-07" +#define COLVARPROXY_VERSION "2021-03-02" #endif diff --git a/src/COLVARS/fix_colvars.cpp b/src/COLVARS/fix_colvars.cpp index f4f6974c70..45e12951bb 100644 --- a/src/COLVARS/fix_colvars.cpp +++ b/src/COLVARS/fix_colvars.cpp @@ -480,7 +480,7 @@ void FixColvars::one_time_init() memory->create(force_buf,3*num_coords,"colvars:force_buf"); if (me == 0) { - std::vector &tl = *(proxy->modify_atom_ids()); + std::vector const &tl = *(proxy->get_atom_ids()); inthash_t *hashtable=new inthash_t; inthash_init(hashtable, num_coords); idmap = (void *)hashtable; @@ -563,7 +563,7 @@ void FixColvars::setup(int vflag) if (me == 0) { - std::vector &id = *(proxy->modify_atom_ids()); + std::vector const &id = *(proxy->get_atom_ids()); std::vector &tp = *(proxy->modify_atom_types()); std::vector &cd = *(proxy->modify_atom_positions()); std::vector &of = *(proxy->modify_atom_total_forces()); @@ -836,7 +836,7 @@ void FixColvars::post_force(int /*vflag*/) if (me == 0) { - std::vector &fo = *(proxy->modify_atom_new_colvar_forces()); + std::vector &fo = *(proxy->modify_atom_applied_forces()); double *fbuf = force_buf; for (int j=0; j < num_coords; ++j) {