diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index 14deceeb87..dd028ae65e 100644 Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ diff --git a/lib/colvars/Makefile.common b/lib/colvars/Makefile.common index c555c534d8..0482cff74a 100644 --- a/lib/colvars/Makefile.common +++ b/lib/colvars/Makefile.common @@ -31,6 +31,7 @@ COLVARS_SRCS = \ colvarbias_meta.cpp \ colvarbias_restraint.cpp \ colvarcomp_angles.cpp \ + colvarcomp_apath.cpp \ colvarcomp_coordnums.cpp \ colvarcomp.cpp \ colvarcomp_distances.cpp \ @@ -41,8 +42,10 @@ COLVARS_SRCS = \ colvardeps.cpp \ colvargrid.cpp \ colvarmodule.cpp \ + colvarparams.cpp \ colvarparse.cpp \ colvarproxy.cpp \ + colvarproxy_replicas.cpp \ colvarscript.cpp \ colvartypes.cpp \ colvarvalue.cpp diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps index a0d8515bc1..33a9dd10ed 100644 --- a/lib/colvars/Makefile.deps +++ b/lib/colvars/Makefile.deps @@ -1,235 +1,90 @@ $(COLVARS_OBJ_DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ - colvarparse.h colvaratoms.h colvardeps.h + colvarparse.h colvarparams.h colvaratoms.h colvardeps.h $(COLVARS_OBJ_DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvar.h \ - colvarparse.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 \ - colvarbias_abf.h colvarbias.h colvargrid.h colvar_UIestimator.h + colvarparse.h colvarparams.h colvardeps.h colvarbias_abf.h colvarbias.h \ + colvargrid.h colvar_UIestimator.h $(COLVARS_OBJ_DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h \ colvars_version.h colvarbias.h colvar.h colvarvalue.h colvartypes.h \ - colvarparse.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 \ - colvarbias_alb.h + colvarparse.h colvarparams.h colvardeps.h colvarbias_alb.h $(COLVARS_OBJ_DIR)colvarbias.o: colvarbias.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvarbias.h \ - colvar.h colvarparse.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 \ - colvargrid.h + colvar.h colvarparse.h colvarparams.h colvardeps.h colvargrid.h $(COLVARS_OBJ_DIR)colvarbias_histogram.o: colvarbias_histogram.cpp \ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ - colvarvalue.h colvar.h colvarparse.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 \ + colvarvalue.h colvar.h colvarparse.h colvarparams.h colvardeps.h \ colvarbias_histogram.h colvarbias.h colvargrid.h $(COLVARS_OBJ_DIR)colvarbias_meta.o: colvarbias_meta.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvar.h \ - colvarparse.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 \ - colvarbias_meta.h colvarbias.h colvargrid.h + colvarparse.h colvarparams.h colvardeps.h colvarbias_meta.h colvarbias.h \ + colvargrid.h $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ colvarvalue.h colvarbias_restraint.h colvarbias.h colvar.h colvarparse.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 + colvarparams.h colvardeps.h $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \ colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \ - colvarparse.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 + colvarparse.h colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvar_arithmeticpath.h +$(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp $(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \ colvarmodule.h colvars_version.h colvarparse.h colvarvalue.h \ - colvartypes.h colvaratoms.h colvarproxy.h colvardeps.h colvar.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 + colvartypes.h colvarparams.h colvaratoms.h colvarproxy.h colvardeps.h \ + colvar.h colvarcomp.h colvar_arithmeticpath.h $(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvar.h colvarparse.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 + colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \ + colvar_arithmeticpath.h $(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvar.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 + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvar_arithmeticpath.h +$(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp $(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvar.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 + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvar_arithmeticpath.h $(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvar.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 + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvar_arithmeticpath.h $(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \ - colvarvalue.h colvartypes.h colvarparse.h colvar.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 colvarscript.h colvarbias.h + colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \ + colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \ + colvar_arithmeticpath.h colvarscript.h colvarbias.h $(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvardeps.h \ - colvarparse.h + colvarparse.h colvarparams.h $(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.h \ - colvars_version.h colvarvalue.h colvartypes.h colvarparse.h colvar.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 colvargrid.h + colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ + colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvar_arithmeticpath.h colvargrid.h $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \ colvars_version.h colvarparse.h colvarvalue.h colvartypes.h \ - colvarproxy.h colvar.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 \ - colvarbias.h colvarbias_abf.h colvargrid.h colvar_UIestimator.h \ - colvarbias_alb.h colvarbias_histogram.h colvarbias_meta.h \ - colvarbias_restraint.h colvarscript.h colvaratoms.h colvarcomp.h + colvarparams.h colvarproxy.h colvar.h colvardeps.h colvarbias.h \ + colvarbias_abf.h colvargrid.h colvar_UIestimator.h colvarbias_alb.h \ + colvarbias_histogram.h colvarbias_meta.h colvarbias_restraint.h \ + colvarscript.h colvaratoms.h colvarcomp.h colvar_arithmeticpath.h +$(COLVARS_OBJ_DIR)colvarparams.o: colvarparams.cpp colvarmodule.h \ + colvars_version.h colvarvalue.h colvartypes.h colvarparams.h $(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \ - colvars_version.h colvarvalue.h colvartypes.h colvarparse.h + colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ + colvarparams.h $(COLVARS_OBJ_DIR)colvarproxy.o: colvarproxy.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ - colvarscript.h colvarbias.h colvar.h colvarparse.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 \ - colvaratoms.h + colvarscript.h colvarbias.h colvar.h colvarparse.h colvarparams.h \ + colvardeps.h colvaratoms.h +$(COLVARS_OBJ_DIR)colvarproxy_replicas.o: colvarproxy_replicas.cpp \ + colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ + colvarvalue.h $(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarscript.h \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarbias.h colvar.h colvarparse.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 \ + colvarbias.h colvar.h colvarparse.h colvarparams.h colvardeps.h \ colvarproxy.h $(COLVARS_OBJ_DIR)colvartypes.o: colvartypes.cpp colvarmodule.h \ - colvars_version.h colvartypes.h colvarparse.h colvarvalue.h + colvars_version.h colvartypes.h colvarparse.h colvarvalue.h \ + colvarparams.h $(COLVARS_OBJ_DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index e3676084ac..b853aaa1b8 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -29,6 +29,10 @@ colvar::colvar() kinetic_energy = 0.0; potential_energy = 0.0; +#ifdef LEPTON + dev_null = 0.0; +#endif + description = "uninitialized colvar"; init_dependencies(); } @@ -194,23 +198,32 @@ int colvar::init(std::string const &conf) { bool homogeneous = is_enabled(f_cv_linear); for (i = 0; i < cvcs.size(); i++) { - if ((cvm::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) { + if (cvm::fabs(cvm::fabs(cvcs[i]->sup_coeff) - 1.0) > 1.0e-10) { homogeneous = false; } } set_enabled(f_cv_homogeneous, homogeneous); } + // A single-component variable almost concides with its CVC object + if ((cvcs.size() == 1) && is_enabled(f_cv_homogeneous)) { + if ( !is_enabled(f_cv_scripted) && !is_enabled(f_cv_custom_function) && + (cvm::fabs(cvcs[0]->sup_coeff - 1.0) < 1.0e-10) && + (cvcs[0]->sup_np == 1) ) { + enable(f_cv_single_cvc); + } + } + // Colvar is deemed periodic if: // - it is homogeneous // - all cvcs are periodic // - all cvcs have the same period - if (is_enabled(f_cv_homogeneous) && cvcs[0]->b_periodic) { // TODO make this a CVC feature + if (is_enabled(f_cv_homogeneous) && cvcs[0]->is_enabled(f_cvc_periodic)) { bool b_periodic = true; period = cvcs[0]->period; wrap_center = cvcs[0]->wrap_center; for (i = 1; i < cvcs.size(); i++) { - if (!cvcs[i]->b_periodic || cvcs[i]->period != period) { + if (!cvcs[i]->is_enabled(f_cvc_periodic) || cvcs[i]->period != period) { b_periodic = false; period = 0.0; cvm::log("Warning: although one component is periodic, this colvar will " @@ -296,13 +309,13 @@ int colvar::init(std::string const &conf) #ifdef LEPTON int colvar::init_custom_function(std::string const &conf) { - std::string expr; + std::string expr, expr_in; // expr_in is a buffer to remember expr after unsuccessful parsing std::vector pexprs; Lepton::ParsedExpression pexpr; size_t pos = 0; // current position in config string double *ref; - if (!key_lookup(conf, "customFunction", &expr, &pos)) { + if (!key_lookup(conf, "customFunction", &expr_in, &pos)) { return COLVARS_OK; } @@ -310,6 +323,7 @@ int colvar::init_custom_function(std::string const &conf) cvm::log("This colvar uses a custom function.\n"); do { + expr = expr_in; if (cvm::debug()) cvm::log("Parsing expression \"" + expr + "\".\n"); try { @@ -337,8 +351,7 @@ int colvar::init_custom_function(std::string const &conf) // To keep the same workflow, we use a pointer to a double here // that will receive CVC values - even though none was allocated by Lepton ref = &dev_null; - if (cvm::debug()) - cvm::log("Variable " + vn + " is absent from expression \"" + expr + "\".\n"); + cvm::log("Warning: Variable " + vn + " is absent from expression \"" + expr + "\".\n"); } value_eval_var_refs.push_back(ref); } @@ -348,7 +361,7 @@ int colvar::init_custom_function(std::string const &conf) cvm::error("Error compiling expression \"" + expr + "\".\n", INPUT_ERROR); return INPUT_ERROR; } - } while (key_lookup(conf, "customFunction", &expr, &pos)); + } while (key_lookup(conf, "customFunction", &expr_in, &pos)); // Now define derivative with respect to each scalar sub-component @@ -373,8 +386,7 @@ int colvar::init_custom_function(std::string const &conf) catch (...) { // Variable is absent from derivative // To keep the same workflow, we use a pointer to a double here // that will receive CVC values - even though none was allocated by Lepton - if (cvm::debug()) - cvm::log("Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n"); + cvm::log("Warning: Variable " + vvn + " is absent from derivative of \"" + expr + "\" wrt " + vn + ".\n"); ref = &dev_null; } grad_eval_var_refs.push_back(ref); @@ -439,6 +451,21 @@ int colvar::init_custom_function(std::string const &conf) int colvar::init_custom_function(std::string const &conf) { + + std::string expr; + size_t pos = 0; + if (key_lookup(conf, "customFunction", &expr, &pos)) { + std::string msg("Error: customFunction requires the Lepton library."); +#if (__cplusplus < 201103L) + // NOTE: this is not ideal; testing for the Lepton library's version would + // be more accurate, but also less portable + msg += + std::string(" Note also that recent versions of Lepton require C++11: " + "please see https://colvars.github.io/README-c++11.html."); +#endif + return cvm::error(msg, COLVARS_NOT_IMPLEMENTED); + } + return COLVARS_OK; } @@ -449,7 +476,19 @@ int colvar::init_grid_parameters(std::string const &conf) { colvarmodule *cv = cvm::main(); - get_keyval(conf, "width", width, 1.0); + cvm::real default_width = width; + + if (!key_already_set("width")) { + // The first time, check if the CVC has a width to provide + default_width = 1.0; + if (is_enabled(f_cv_single_cvc) && cvcs[0]->is_enabled(f_cvc_width)) { + cvm::real const cvc_width = cvcs[0]->get_param("width"); + default_width = cvc_width; + } + } + + get_keyval(conf, "width", width, default_width); + if (width <= 0.0) { cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); return INPUT_ERROR; @@ -460,12 +499,32 @@ int colvar::init_grid_parameters(std::string const &conf) if (is_enabled(f_cv_scalar)) { + if (is_enabled(f_cv_single_cvc)) { + // Get the default boundaries from the component + if (cvcs[0]->is_enabled(f_cvc_lower_boundary)) { + enable(f_cv_lower_boundary); + enable(f_cv_hard_lower_boundary); + lower_boundary = + *(reinterpret_cast(cvcs[0]->get_param_ptr("lowerBoundary"))); + } + if (cvcs[0]->is_enabled(f_cvc_upper_boundary)) { + enable(f_cv_upper_boundary); + enable(f_cv_hard_upper_boundary); + upper_boundary = + *(reinterpret_cast(cvcs[0]->get_param_ptr("upperBoundary"))); + } + } + if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) { enable(f_cv_lower_boundary); + // Because this is the user's choice, we cannot assume it is a true + // physical boundary + disable(f_cv_hard_lower_boundary); } if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) { enable(f_cv_upper_boundary); + disable(f_cv_hard_upper_boundary); } std::string lw_conf, uw_conf; @@ -520,12 +579,11 @@ harmonicWalls {\n\ } } - if (is_enabled(f_cv_lower_boundary)) { - get_keyval(conf, "hardLowerBoundary", hard_lower_boundary, false); - } - if (is_enabled(f_cv_upper_boundary)) { - get_keyval(conf, "hardUpperBoundary", hard_upper_boundary, false); - } + get_keyval_feature(this, conf, "hardLowerBoundary", f_cv_hard_lower_boundary, + is_enabled(f_cv_hard_lower_boundary)); + + get_keyval_feature(this, conf, "hardUpperBoundary", f_cv_hard_upper_boundary, + is_enabled(f_cv_hard_upper_boundary)); // consistency checks for boundaries and walls if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) { @@ -546,7 +604,8 @@ harmonicWalls {\n\ INPUT_ERROR); return INPUT_ERROR; } - if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) { + if (expand_boundaries && is_enabled(f_cv_hard_lower_boundary) && + is_enabled(f_cv_hard_upper_boundary)) { cvm::error("Error: inconsistent configuration " "(trying to expand boundaries with both " "hardLowerBoundary and hardUpperBoundary enabled).\n", @@ -568,7 +627,8 @@ int colvar::init_extended_Lagrangian(std::string const &conf) cvm::log("Enabling the extended Lagrangian term for colvar \""+ this->name+"\".\n"); - x_ext.type(value()); + // Mark x_ext as uninitialized so we can initialize it to the colvar value when updating + x_ext.type(colvarvalue::type_notset); v_ext.type(value()); fr.type(value()); @@ -646,7 +706,7 @@ int colvar::init_output_flags(std::string const &conf) bool temp; if (get_keyval(conf, "outputSystemForce", temp, false, colvarparse::parse_silent)) { cvm::error("Option outputSystemForce is deprecated: only outputTotalForce is supported instead.\n" - "The two are NOT identical: see http://colvars.github.io/totalforce.html.\n", INPUT_ERROR); + "The two are NOT identical: see https://colvars.github.io/totalforce.html.\n", INPUT_ERROR); return INPUT_ERROR; } } @@ -664,7 +724,7 @@ int colvar::init_output_flags(std::string const &conf) // read the configuration and set up corresponding instances, for // each type of component implemented template int colvar::init_components_type(std::string const &conf, - char const *def_desc, + char const * /* def_desc */, char const *def_config_key) { size_t def_count = 0; @@ -744,6 +804,7 @@ template int colvar::init_components_type(std::string c int colvar::init_components(std::string const &conf) { int error_code = COLVARS_OK; + size_t i = 0, j = 0; error_code |= init_components_type(conf, "distance", "distance"); error_code |= init_components_type(conf, "distance vector", "distanceVec"); @@ -793,9 +854,11 @@ int colvar::init_components(std::string const &conf) error_code |= init_components_type(conf, "eigenvector", "eigenvector"); 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", "subColvar"); + error_code |= init_components_type(conf, "linear combination of other collective variables", "linearCombination"); error_code |= init_components_type(conf, "geometrical path collective variables (s) for other CVs", "gspathCV"); 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"); if (!cvcs.size() || (error_code != COLVARS_OK)) { cvm::error("Error: no valid components were provided " @@ -804,15 +867,26 @@ int colvar::init_components(std::string const &conf) return INPUT_ERROR; } + // Check for uniqueness of CVC names (esp. if user-provided) + for (i = 0; i < cvcs.size(); i++) { + for (j = i+1; j < cvcs.size(); j++) { + if (cvcs[i]->name == cvcs[j]->name) { + cvm::error("Components " + cvm::to_str(i) + " and " + cvm::to_str(j) +\ + " cannot have the same name \"" + cvcs[i]->name+ "\".\n", INPUT_ERROR); + return INPUT_ERROR; + } + } + } + n_active_cvcs = cvcs.size(); - cvm::log("All components initialized.\n"); - // Store list of children cvcs for dependency checking purposes - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { add_child(cvcs[i]); } + cvm::log("All components initialized.\n"); + return COLVARS_OK; } @@ -921,7 +995,7 @@ int colvar::parse_analysis(std::string const &conf) return cvm::error("Error: collective variable \""+acf_colvar_name+ "\" is not defined at this time.\n", INPUT_ERROR); } - cv2->enable(f_cv_fdiff_velocity); + cv2->enable(f_cv_fdiff_velocity); // Manual dependency to object of same type } else if (acf_type_str == to_lower_cppstr(std::string("coordinate_p2"))) { acf_type = acf_p2coor; } else { @@ -998,6 +1072,8 @@ 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_single_cvc, "single component", f_type_static); + init_feature(f_cv_linear, "linear", f_type_static); init_feature(f_cv_scalar, "scalar", f_type_static); @@ -1023,6 +1099,12 @@ int colvar::init_dependencies() { init_feature(f_cv_upper_boundary, "upper boundary", f_type_user); require_feature_self(f_cv_upper_boundary, f_cv_scalar); + init_feature(f_cv_hard_lower_boundary, "hard lower boundary", f_type_user); + require_feature_self(f_cv_hard_lower_boundary, f_cv_lower_boundary); + + init_feature(f_cv_hard_upper_boundary, "hard upper boundary", f_type_user); + require_feature_self(f_cv_hard_upper_boundary, f_cv_upper_boundary); + init_feature(f_cv_grid, "grid", f_type_dynamic); require_feature_self(f_cv_grid, f_cv_lower_boundary); require_feature_self(f_cv_grid, f_cv_upper_boundary); @@ -1227,7 +1309,7 @@ int colvar::collect_cvc_data() } -int colvar::check_cvc_range(int first_cvc, size_t num_cvcs) +int colvar::check_cvc_range(int first_cvc, size_t /* num_cvcs */) { if ((first_cvc < 0) || (first_cvc >= ((int) cvcs.size()))) { cvm::error("Error: trying to address a component outside the " @@ -1512,7 +1594,7 @@ int colvar::calc_colvar_properties() // initialize the restraint center in the first step to the value // just calculated from the cvcs - if (cvm::step_relative() == 0 && !after_restart) { + if ((cvm::step_relative() == 0 && !after_restart) || x_ext.type() == colvarvalue::type_notset) { x_ext = x; v_ext.reset(); // (already 0; added for clarity) } @@ -1859,6 +1941,59 @@ int colvar::update_cvc_config(std::vector const &confs) } +int colvar::cvc_param_exists(std::string const ¶m_name) +{ + if (is_enabled(f_cv_single_cvc)) { + return cvcs[0]->param_exists(param_name); + } + return cvm::error("Error: calling colvar::cvc_param_exists() for a variable " + "with more than one component.\n", COLVARS_NOT_IMPLEMENTED); +} + + +cvm::real colvar::get_cvc_param(std::string const ¶m_name) +{ + if (is_enabled(f_cv_single_cvc)) { + return cvcs[0]->get_param(param_name); + } + cvm::error("Error: calling colvar::get_cvc_param() for a variable " + "with more than one component.\n", COLVARS_NOT_IMPLEMENTED); + return 0.0; +} + + +void const *colvar::get_cvc_param_ptr(std::string const ¶m_name) +{ + if (is_enabled(f_cv_single_cvc)) { + return cvcs[0]->get_param_ptr(param_name); + } + cvm::error("Error: calling colvar::get_cvc_param() for a variable " + "with more than one component.\n", COLVARS_NOT_IMPLEMENTED); + return NULL; +} + + +colvarvalue const *colvar::get_cvc_param_grad(std::string const ¶m_name) +{ + if (is_enabled(f_cv_single_cvc)) { + return cvcs[0]->get_param_grad(param_name); + } + cvm::error("Error: calling colvar::get_cvc_param_grad() for a variable " + "with more than one component.\n", COLVARS_NOT_IMPLEMENTED); + return NULL; +} + + +int colvar::set_cvc_param(std::string const ¶m_name, void const *new_value) +{ + if (is_enabled(f_cv_single_cvc)) { + return cvcs[0]->set_param(param_name, new_value); + } + return cvm::error("Error: calling colvar::set_cvc_param() for a variable " + "with more than one component.\n", COLVARS_NOT_IMPLEMENTED); +} + + // ******************** METRIC FUNCTIONS ******************** // Use the metrics defined by \link colvar::cvc \endlink objects diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index e2ab0f3c1d..f7d5a55627 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -217,11 +217,11 @@ public: /// should equal the system force plus \link f \endlink colvarvalue ft; - /// Period, if this variable is periodic cvm::real period; - cvm::real wrap_center; + /// Center of wrapping, if this variable is periodic + cvm::real wrap_center; /// \brief Expand the boundaries of multiples of width, to keep the /// value always within range @@ -233,8 +233,6 @@ public: colvarvalue lower_wall; /// \brief Force constant for the lower boundary potential (|x-xb|^2) cvm::real lower_wall_k; - /// \brief Whether this colvar has a hard lower boundary - bool hard_lower_boundary; /// \brief Location of the upper boundary colvarvalue upper_boundary; @@ -242,8 +240,6 @@ public: colvarvalue upper_wall; /// \brief Force constant for the upper boundary potential (|x-xb|^2) cvm::real upper_wall_k; - /// \brief Whether this colvar has a hard upper boundary - bool hard_upper_boundary; /// \brief Is the interval defined by the two boundaries periodic? bool periodic_boundaries() const; @@ -356,7 +352,7 @@ public: /// colvar::update()) to the external degrees of freedom void communicate_forces(); - /// \brief Enables and disables individual CVCs based on flags + /// \brief Enables and disables individual CVCs based on the given array int set_cvc_flags(std::vector const &flags); /// \brief Updates the flags in the CVC objects, and their number @@ -365,7 +361,23 @@ public: /// \brief Modify the configuration of CVCs (currently, only base class data) int update_cvc_config(std::vector const &confs); + /// Whether this named parameter exists (in the first and only component) + int cvc_param_exists(std::string const ¶m_name); + + /// Get the value of the named parameter (from the first and only component) + cvm::real get_cvc_param(std::string const ¶m_name); + + /// Get a pointer to the named parameter (from the first and only component) + void const *get_cvc_param_ptr(std::string const ¶m_name); + + /// Pointer to the derivative of the variable with respect to param_name + colvarvalue const *get_cvc_param_grad(std::string const ¶m_name); + + /// Set the named parameter in the first and only component to the given value + int set_cvc_param(std::string const ¶m_name, void const *new_value); + protected: + /// \brief Number of CVC objects with an active flag size_t n_active_cvcs; @@ -380,11 +392,24 @@ protected: public: - /// \brief Return the number of CVC objects defined - inline size_t num_cvcs() const { return cvcs.size(); } + /// \brief Number of dimensions of the value of this colvar + inline size_t num_dimensions() const + { + return value().size(); + } - /// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags) - inline size_t num_active_cvcs() const { return n_active_cvcs; } + /// \brief Number of CVC objects defined + inline size_t num_cvcs() const + { + return cvcs.size(); + } + + /// \brief number of CVC objects with an active flag (as set by + /// update_cvc_flags) + inline size_t num_active_cvcs() const + { + return n_active_cvcs; + } /// \brief Use the internal metrics (as from \link colvar::cvc /// \endlink objects) to calculate square distances and gradients @@ -413,22 +438,20 @@ public: /// Handles correctly symmetries and periodic boundary conditions void wrap(colvarvalue &x_unwrapped) const; - /// Read the analysis tasks int parse_analysis(std::string const &conf); /// Perform analysis tasks int analyze(); - /// Read the value from a collective variable trajectory file std::istream & read_traj(std::istream &is); + /// Output formatted values to the trajectory file std::ostream & write_traj(std::ostream &os); /// Write a label to the trajectory file (comment line) std::ostream & write_traj_label(std::ostream &os); - /// Read the collective variable from a restart file std::istream & read_restart(std::istream &is); /// Write the collective variable to a restart file @@ -437,8 +460,8 @@ public: /// Write output files (if defined, e.g. in analysis mode) int write_output_files(); - protected: + /// Previous value (to calculate velocities during analysis) colvarvalue x_old; @@ -530,11 +553,12 @@ protected: /// Calculate the running average and its standard deviation int calc_runave(); - /// If extended Lagrangian active: colvar energies (kinetic and harmonic potential) + /// If extended Lagrangian active: colvar kinetic energy cvm::real kinetic_energy; + /// If extended Lagrangian active: colvar harmonic potential cvm::real potential_energy; -public: +public: // collective variable component base class class cvc; @@ -577,6 +601,8 @@ public: class CVBasedPath; class gspathCV; class gzpathCV; + class aspathCV; + class azpathCV; // non-scalar components class distance_vec; @@ -596,7 +622,6 @@ protected: /// in all cvcs (called when enabling f_cv_collect_gradients) void build_atom_list(void); -private: /// Name of scripted function to be used std::string scripted_function; @@ -620,6 +645,7 @@ private: #endif public: + /// \brief Sorted array of (zero-based) IDs for all atoms involved std::vector atom_ids; @@ -628,34 +654,35 @@ public: /// For scalar variables only! std::vector atomic_gradients; - inline size_t n_components() const { - return cvcs.size(); - } - /// \brief Get vector of vectors of atom IDs for all atom groups virtual std::vector > get_atom_lists(); }; + inline cvm::real const & colvar::force_constant() const { return ext_force_k; } + inline colvarvalue const & colvar::value() const { return x_reported; } + inline colvarvalue const & colvar::actual_value() const { return x; } + inline colvarvalue const & colvar::run_ave() const { return runave; } + inline colvarvalue const & colvar::velocity() const { return v_reported; diff --git a/lib/colvars/colvar_UIestimator.h b/lib/colvars/colvar_UIestimator.h index 365f46148a..cb4c7ed57d 100644 --- a/lib/colvars/colvar_UIestimator.h +++ b/lib/colvars/colvar_UIestimator.h @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvar_arithmeticpath.h b/lib/colvars/colvar_arithmeticpath.h new file mode 100644 index 0000000000..97d42a1bbc --- /dev/null +++ b/lib/colvars/colvar_arithmeticpath.h @@ -0,0 +1,129 @@ +#ifndef ARITHMETICPATHCV_H +#define ARITHMETICPATHCV_H + +#include "colvarmodule.h" + +#include +#include +#include +#include + +namespace ArithmeticPathCV { + +using std::vector; + +enum path_sz {S, Z}; + +template +class ArithmeticPathBase { +public: + ArithmeticPathBase() {} + virtual ~ArithmeticPathBase() {} + virtual void initialize(size_t p_num_elements, size_t p_total_frames, double p_lambda, const vector& p_element, const vector& p_weights); + virtual void updateDistanceToReferenceFrames() = 0; + virtual void computeValue(); + virtual void computeDerivatives(); + virtual void compute(); +protected: + scalar_type lambda; + vector weights; + size_t num_elements; + size_t total_frames; + vector< vector > frame_element_distances; + scalar_type s; + scalar_type z; + vector dsdx; + vector dzdx; +private: + // intermediate variables + vector s_numerator_frame; + vector s_denominator_frame; + scalar_type numerator_s; + scalar_type denominator_s; + scalar_type normalization_factor; +}; + +template +void ArithmeticPathBase::initialize(size_t p_num_elements, size_t p_total_frames, double p_lambda, const vector& p_element, const vector& p_weights) { + lambda = p_lambda; + weights = p_weights; + num_elements = p_num_elements; + total_frames = p_total_frames; + frame_element_distances.resize(total_frames, p_element); + for (size_t i_frame = 0; i_frame < frame_element_distances.size(); ++i_frame) { + for (size_t j_elem = 0; j_elem < num_elements; ++j_elem) { + frame_element_distances[i_frame][j_elem].reset(); + } + } + s = scalar_type(0); + z = scalar_type(0); + dsdx = p_element; + dzdx = p_element; + s_numerator_frame.resize(total_frames, scalar_type(0)); + s_denominator_frame.resize(total_frames, scalar_type(0)); + numerator_s = scalar_type(0); + denominator_s = scalar_type(0); + normalization_factor = 1.0 / static_cast(total_frames - 1); +} + +template +void ArithmeticPathBase::computeValue() { + updateDistanceToReferenceFrames(); + numerator_s = scalar_type(0); + denominator_s = scalar_type(0); + for (size_t i_frame = 0; i_frame < frame_element_distances.size(); ++i_frame) { + scalar_type exponent_tmp = scalar_type(0); + for (size_t j_elem = 0; j_elem < num_elements; ++j_elem) { + exponent_tmp += weights[j_elem] * frame_element_distances[i_frame][j_elem] * weights[j_elem] * frame_element_distances[i_frame][j_elem]; + } + exponent_tmp = exponent_tmp * -1.0 * lambda; + // prevent underflow if the argument of cvm::exp is less than -708.4 + if (exponent_tmp > -708.4) { + exponent_tmp = cvm::exp(exponent_tmp); + } else { + exponent_tmp = 0; + } + numerator_s += static_cast(i_frame) * exponent_tmp; + denominator_s += exponent_tmp; + s_numerator_frame[i_frame] = static_cast(i_frame) * exponent_tmp; + s_denominator_frame[i_frame] = exponent_tmp; + } + s = numerator_s / denominator_s * normalization_factor; + z = -1.0 / lambda * cvm::logn(denominator_s); +} + +template +void ArithmeticPathBase::compute() { + computeValue(); + computeDerivatives(); +} + +template +void ArithmeticPathBase::computeDerivatives() { + for (size_t j_elem = 0; j_elem < num_elements; ++j_elem) { + element_type dsdxj_numerator_part1(dsdx[j_elem]); + element_type dsdxj_numerator_part2(dsdx[j_elem]); + element_type dzdxj_numerator(dsdx[j_elem]); + dsdxj_numerator_part1.reset(); + dsdxj_numerator_part2.reset(); + dzdxj_numerator.reset(); + for (size_t i_frame = 0; i_frame < frame_element_distances.size(); ++i_frame) { + element_type derivative_tmp = -2.0 * lambda * weights[j_elem] * weights[j_elem] * frame_element_distances[i_frame][j_elem]; + dsdxj_numerator_part1 += s_numerator_frame[i_frame] * derivative_tmp; + dsdxj_numerator_part2 += s_denominator_frame[i_frame] * derivative_tmp; + dzdxj_numerator += s_denominator_frame[i_frame] * derivative_tmp; + } + dsdxj_numerator_part1 *= denominator_s; + dsdxj_numerator_part2 *= numerator_s; + if ((dsdxj_numerator_part1 - dsdxj_numerator_part2).norm() < std::numeric_limits::min()) { + dsdx[j_elem] = 0; + } else { + dsdx[j_elem] = (dsdxj_numerator_part1 - dsdxj_numerator_part2) / (denominator_s * denominator_s) * normalization_factor; + } + dzdx[j_elem] = -1.0 / lambda * dzdxj_numerator / denominator_s; + } +} + +} + +#endif // ARITHMETICPATHCV_H diff --git a/lib/colvars/colvar_geometricpath.h b/lib/colvars/colvar_geometricpath.h index a7ef7f337c..7f8acfd233 100644 --- a/lib/colvars/colvar_geometricpath.h +++ b/lib/colvars/colvar_geometricpath.h @@ -2,12 +2,14 @@ #define GEOMETRICPATHCV_H // 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 +// 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 "colvarmodule.h" + #include #include #include @@ -67,8 +69,8 @@ public: virtual ~GeometricPathBase() {} virtual void initialize(size_t vector_size, const element_type& element = element_type(), size_t total_frames = 1, bool p_use_second_closest_frame = true, bool p_use_third_closest_frame = false, bool p_use_z_square = false); virtual void initialize(size_t vector_size, const std::vector& elements, size_t total_frames = 1, bool p_use_second_closest_frame = true, bool p_use_third_closest_frame = false, bool p_use_z_square = false); - virtual void prepareVectors(); - virtual void updateReferenceDistances(); + virtual void prepareVectors() = 0; + virtual void updateDistanceToReferenceFrames() = 0; virtual void compute(); virtual void determineClosestFrames(); virtual void computeValue(); @@ -119,7 +121,7 @@ void GeometricPathBase::initialize(size_t } template -void GeometricPathBase::initialize(size_t vector_size, const std::vector& elements, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) { +void GeometricPathBase::initialize(size_t /* vector_size */, const std::vector& elements, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) { v1v1 = scalar_type(); v2v2 = scalar_type(); v3v3 = scalar_type(); @@ -151,18 +153,6 @@ void GeometricPathBase::initialize(size_t m = static_cast(1.0); } -template -void GeometricPathBase::prepareVectors() { - std::cout << "Warning: you should not call the prepareVectors() in base class!\n"; - std::cout << std::flush; -} - -template -void GeometricPathBase::updateReferenceDistances() { - std::cout << "Warning: you should not call the updateReferenceDistances() in base class!\n"; - std::cout << std::flush; -} - template void GeometricPathBase::compute() { computeValue(); @@ -182,7 +172,7 @@ void GeometricPathBase::determineClosestFr // sigma(z) is on the right side of the closest frame sign = -1; } - if (std::abs(static_cast(frame_index[0]) - static_cast(frame_index[1])) > 1) { + if (cvm::fabs(static_cast(frame_index[0]) - static_cast(frame_index[1])) > 1) { std::cout << "Warning: Geometrical pathCV relies on the assumption that the second closest frame is the neighbouring frame\n"; std::cout << " Please check your configuration or increase restraint on z(σ)\n"; for (size_t i_frame = 0; i_frame < frame_index.size(); ++i_frame) { @@ -197,7 +187,7 @@ void GeometricPathBase::determineClosestFr template void GeometricPathBase::computeValue() { - updateReferenceDistances(); + updateDistanceToReferenceFrames(); determineClosestFrames(); prepareVectors(); v1v1 = scalar_type(); @@ -218,14 +208,14 @@ void GeometricPathBase::computeValue() { v4v4 += v4[i_elem] * v4[i_elem]; } } - f = (std::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)) - v1v3) / v3v3; + f = (cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2)) - v1v3) / v3v3; if (path_type == Z) { dx = 0.5 * (f - 1); zz = v1v1 + 2 * dx * v1v4 + dx * dx * v4v4; if (use_z_square) { z = zz; } else { - z = std::sqrt(std::fabs(zz)); + z = cvm::sqrt(cvm::fabs(zz)); } } if (path_type == S) { @@ -235,7 +225,7 @@ void GeometricPathBase::computeValue() { template void GeometricPathBase::computeDerivatives() { - const scalar_type factor1 = 1.0 / (2.0 * v3v3 * std::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2))); + const scalar_type factor1 = 1.0 / (2.0 * v3v3 * cvm::sqrt(v1v3 * v1v3 - v3v3 * (v1v1 - v2v2))); const scalar_type factor2 = 1.0 / v3v3; for (size_t i_elem = 0; i_elem < v1.size(); ++i_elem) { // Compute the derivative of f with vector v1 diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index af17d8df66..180fc69a85 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -628,39 +628,42 @@ int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf) int cvm::atom_group::add_index_group(std::string const &index_group_name) { - colvarmodule *cv = cvm::main(); + std::vector const &index_group_names = + cvm::main()->index_group_names; + std::vector *> const &index_groups = + cvm::main()->index_groups; - std::list::iterator names_i = cv->index_group_names.begin(); - std::list >::iterator index_groups_i = cv->index_groups.begin(); - for ( ; names_i != cv->index_group_names.end() ; ++names_i, ++index_groups_i) { - if (*names_i == index_group_name) + size_t i_group = 0; + for ( ; i_group < index_groups.size(); i_group++) { + if (index_group_names[i_group] == index_group_name) break; } - if (names_i == cv->index_group_names.end()) { - cvm::error("Error: could not find index group "+ - index_group_name+" among those provided by the index file.\n", - INPUT_ERROR); - return COLVARS_ERROR; + if (i_group >= index_group_names.size()) { + return cvm::error("Error: could not find index group "+ + index_group_name+" among those already provided.\n", + INPUT_ERROR); } - atoms_ids.reserve(atoms_ids.size()+index_groups_i->size()); + int error_code = COLVARS_OK; + + std::vector const &index_group = *(index_groups[i_group]); + + atoms_ids.reserve(atoms_ids.size()+index_group.size()); if (is_enabled(f_ag_scalable)) { - for (size_t i = 0; i < index_groups_i->size(); i++) { - add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i])); + for (size_t i = 0; i < index_group.size(); i++) { + error_code |= + add_atom_id((cvm::proxy)->check_atom_id(index_group[i])); } } else { - atoms.reserve(atoms.size()+index_groups_i->size()); - for (size_t i = 0; i < index_groups_i->size(); i++) { - add_atom(cvm::atom((*index_groups_i)[i])); + atoms.reserve(atoms.size()+index_group.size()); + for (size_t i = 0; i < index_group.size(); i++) { + error_code |= add_atom(cvm::atom(index_group[i])); } } - if (cvm::get_error()) - return COLVARS_ERROR; - - return COLVARS_OK; + return error_code; } diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index c8b6ac4541..af6a529f8a 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 724326d3b4..ff7bbad19b 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -233,25 +233,18 @@ int colvarbias::add_colvar(std::string const &cv_name) } colvars.push_back(cv); + cv->biases.push_back(this); // add back-reference to this bias to colvar + + // Add dependency link. All biases need at least the value of each colvar + // although possibly not at all timesteps + add_child(cv); colvar_forces.push_back(colvarvalue()); colvar_forces.back().type(cv->value()); // make sure each force is initialized to zero colvar_forces.back().is_derivative(); // colvar constraints are not applied to the force colvar_forces.back().reset(); - previous_colvar_forces.push_back(colvar_forces.back()); - cv->biases.push_back(this); // add back-reference to this bias to colvar - - if (is_enabled(f_cvb_apply_force)) { - cv->enable(f_cv_gradient); - } - - // Add dependency link. - // All biases need at least the value of each colvar - // although possibly not at all timesteps - add_child(cv); - } else { cvm::error("Error: cannot find a colvar named \""+ cv_name+"\".\n", INPUT_ERROR); @@ -268,13 +261,29 @@ int colvarbias::update() cvm::log("Updating the "+bias_type+" bias \""+this->name+"\".\n"); } + int error_code = COLVARS_OK; + has_data = true; + error_code |= calc_energy(NULL); + error_code |= calc_forces(NULL); + + return error_code; +} + + +int colvarbias::calc_energy(std::vector const *) +{ bias_energy = 0.0; + return COLVARS_OK; +} + + +int colvarbias::calc_forces(std::vector const *) +{ for (size_t ir = 0; ir < num_variables(); ir++) { colvar_forces[ir].reset(); } - return COLVARS_OK; } @@ -309,7 +318,7 @@ int colvarbias::end_of_step() } -int colvarbias::change_configuration(std::string const &conf) +int colvarbias::change_configuration(std::string const & /* conf */) { cvm::error("Error: change_configuration() not implemented.\n", COLVARS_NOT_IMPLEMENTED); @@ -317,7 +326,7 @@ int colvarbias::change_configuration(std::string const &conf) } -cvm::real colvarbias::energy_difference(std::string const &conf) +cvm::real colvarbias::energy_difference(std::string const & /* conf */) { cvm::error("Error: energy_difference() not implemented.\n", COLVARS_NOT_IMPLEMENTED); @@ -336,7 +345,7 @@ int colvarbias::current_bin() cvm::error("Error: current_bin() not implemented.\n"); return COLVARS_NOT_IMPLEMENTED; } -int colvarbias::bin_count(int bin_index) +int colvarbias::bin_count(int /* bin_index */) { cvm::error("Error: bin_count() not implemented.\n"); return COLVARS_NOT_IMPLEMENTED; @@ -410,6 +419,8 @@ std::istream & colvarbias::read_state(std::istream &is) !(is >> brace) || !(brace == "{") || !(is >> colvarparse::read_block("configuration", conf)) || (set_state_params(conf) != COLVARS_OK) ) { + if (key != bias_type) + cvm::log("Found key \"" + key + "\" instead of \"" + bias_type + "\"\n"); cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+ this->name+"\" at position "+ cvm::to_str(static_cast(is.tellg()))+ @@ -646,7 +657,7 @@ std::string const colvarbias_ti::get_state_params() const } -int colvarbias_ti::set_state_params(std::string const &state_conf) +int colvarbias_ti::set_state_params(std::string const & /* state_conf */) { return COLVARS_OK; } diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index eac88a7f18..ed72c6063d 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -51,17 +51,18 @@ public: } /// Retrieve colvar values and calculate their biasing forces - /// Return bias energy + /// Some implementations may use calc_energy() and calc_forces() virtual int update(); - /// \brief Compute the energy of the bias with alternative values of the - /// collective variables (suitable for bias exchange) - virtual int calc_energy(std::vector const &values = - std::vector(0)) - { - cvm::error("Error: calc_energy() not implemented.\n", COLVARS_NOT_IMPLEMENTED); - return COLVARS_NOT_IMPLEMENTED; - } + /// Compute the energy of the bias + /// Uses the vector of colvar values provided if not NULL, and the values + /// currently cached in the bias instance otherwise + virtual int calc_energy(std::vector const *values); + + /// Compute the forces due to the bias + /// Uses the vector of colvar values provided if not NULL, and the values + /// currently cached in the bias instance otherwise + virtual int calc_forces(std::vector const *values); /// Send forces to the collective variables virtual void communicate_forces(); diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index bac021be99..bcd532e680 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -35,6 +35,8 @@ int colvarbias_abf::init(std::string const &conf) { colvarbias::init(conf); + colvarproxy *proxy = cvm::main()->proxy; + enable(f_cvb_scalar_variables); enable(f_cvb_calc_pmf); @@ -76,11 +78,13 @@ int colvarbias_abf::init(std::string const &conf) // shared ABF get_keyval(conf, "shared", shared_on, false); if (shared_on) { - if (!cvm::replica_enabled() || cvm::replica_num() <= 1) { - cvm::error("Error: shared ABF requires more than one replica."); - return COLVARS_ERROR; + if ((proxy->replica_enabled() != COLVARS_OK) || + (proxy->num_replicas() <= 1)) { + return cvm::error("Error: shared ABF requires more than one replica.", + INPUT_ERROR); } - cvm::log("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n"); + cvm::log("shared ABF will be applied among "+ + cvm::to_str(proxy->num_replicas()) + " replicas.\n"); if (cvm::proxy->smp_enabled() == COLVARS_OK) { cvm::error("Error: shared ABF is currently not available with SMP parallelism; " "please set \"SMP off\" at the top of the Colvars configuration file.\n", @@ -483,13 +487,18 @@ int colvarbias_abf::update() eabf_UI.update(cvm::step_absolute(), x, y); } + /// Add the bias energy for 1D ABF + bias_energy = calc_energy(NULL); + return COLVARS_OK; } int colvarbias_abf::replica_share() { - if ( !cvm::replica_enabled() ) { + colvarproxy *proxy = cvm::main()->proxy; + + if (proxy->replica_enabled() != COLVARS_OK) { cvm::error("Error: shared ABF: No replicas.\n"); return COLVARS_ERROR; } @@ -508,12 +517,12 @@ int colvarbias_abf::replica_share() { size_t msg_total = data_n*sizeof(size_t) + samp_start; char* msg_data = new char[msg_total]; - if (cvm::replica_index() == 0) { + if (proxy->replica_index() == 0) { int p; // Replica 0 collects the delta gradient and count from the others. - for (p = 1; p < cvm::replica_num(); p++) { + for (p = 1; p < proxy->num_replicas(); p++) { // Receive the deltas. - cvm::replica_comm_recv(msg_data, msg_total, p); + proxy->replica_comm_recv(msg_data, msg_total, p); // Map the deltas from the others into the grids. last_gradients->raw_data_in((cvm::real*)(&msg_data[0])); @@ -528,8 +537,8 @@ int colvarbias_abf::replica_share() { // Now we must send the combined gradient to the other replicas. gradients->raw_data_out((cvm::real*)(&msg_data[0])); samples->raw_data_out((size_t*)(&msg_data[samp_start])); - for (p = 1; p < cvm::replica_num(); p++) { - cvm::replica_comm_send(msg_data, msg_total, p); + for (p = 1; p < proxy->num_replicas(); p++) { + proxy->replica_comm_send(msg_data, msg_total, p); } } else { @@ -541,10 +550,10 @@ int colvarbias_abf::replica_share() { // Cast the raw char data to the gradient and samples. last_gradients->raw_data_out((cvm::real*)(&msg_data[0])); last_samples->raw_data_out((size_t*)(&msg_data[samp_start])); - cvm::replica_comm_send(msg_data, msg_total, 0); + proxy->replica_comm_send(msg_data, msg_total, 0); // We now receive the combined gradient from Replica 0. - cvm::replica_comm_recv(msg_data, msg_total, 0); + proxy->replica_comm_recv(msg_data, msg_total, 0); // We sync to the combined gradient computed by Replica 0. gradients->raw_data_in((cvm::real*)(&msg_data[0])); samples->raw_data_in((size_t*)(&msg_data[samp_start])); @@ -552,7 +561,7 @@ int colvarbias_abf::replica_share() { // Without a barrier it's possible that one replica starts // share 2 when other replicas haven't finished share 1. - cvm::replica_comm_barrier(); + proxy->replica_comm_barrier(); // Done syncing the replicas. delete[] msg_data; @@ -832,3 +841,50 @@ int colvarbias_abf::write_output_files() write_gradients_samples(output_prefix); return COLVARS_OK; } + +int colvarbias_abf::calc_energy(std::vector const *values) +{ + if (values) { + return cvm::error("colvarbias_abf::calc_energy() with an argument " + "is currently not implemented.\n", + COLVARS_NOT_IMPLEMENTED); + } + + if (num_variables() != 1) return 0.0; + + // Get the home bin. + int home0 = gradients->current_bin_scalar(0); + if (home0 < 0) return 0.0; + int gradient_len = (int)(gradients->number_of_points(0)); + int home = (home0 < gradient_len) ? home0 : (gradient_len-1); + + // Integrate the gradient up to the home bin. + cvm::real sum = 0.0; + for (int i = 0; i < home; i++) { + std::vector ix(1,i); + + // Include the full_samples factor if necessary. + unsigned int count = samples->value(ix); + cvm::real fact = 1.0; + if ( count < full_samples ) { + fact = (count < min_samples) ? 0.0 : + (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples)); + } + if (count > 0) sum += fact*gradients->value(ix)/count*gradients->widths[0]; + } + + // Integrate the gradient up to the current position in the home interval, a fractional portion of a bin. + std::vector ix(1,home); + cvm::real frac = gradients->current_bin_scalar_fraction(0); + unsigned int count = samples->value(ix); + cvm::real fact = 1.0; + if ( count < full_samples ) { + fact = (count < min_samples) ? 0.0 : + (cvm::real(count - min_samples)) / (cvm::real(full_samples - min_samples)); + } + if (count > 0) + sum += fact*gradients->value(ix)/count*gradients->widths[0]*frac; + + // The applied potential is the negative integral of force samples. + return -sum; +} diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index 4bcc149da5..5a1a1c8128 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -159,7 +159,9 @@ private: virtual std::istream& read_state_data(std::istream&); virtual std::ostream& write_state_data(std::ostream&); virtual int write_output_files(); + + /// Calculate the bias energy for 1D ABF + virtual int calc_energy(std::vector const *values); }; #endif - diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 187ecc363a..1d20245fb2 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -255,14 +255,12 @@ int colvarbias_alb::set_state_params(std::string const &conf) if (!get_keyval(conf, "maxCouplingRange", max_coupling_range)) cvm::fatal_error("Error: maxCouplingRange is missing from the restart.\n"); - if (!get_keyval(conf, "couplingRate", coupling_rate)) cvm::fatal_error("Error: current setCoupling is missing from the restart.\n"); if (!get_keyval(conf, "couplingAccum", coupling_accum)) cvm::fatal_error("Error: couplingAccum is missing from the restart.\n"); - if (!get_keyval(conf, "mean", means)) cvm::fatal_error("Error: current mean is missing from the restart.\n"); @@ -388,7 +386,7 @@ std::ostream & colvarbias_alb::write_traj(std::ostream &os) for (size_t i = 0; i < means.size(); i++) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) - << -2. * (means[i] / (static_cast (colvar_centers[i])) - 1) * ssd[i] / (fmax(update_calls,2) - 1); + << -2.0 * (means[i] / (static_cast(colvar_centers[i])) - 1) * ssd[i] / (fmax(update_calls, 2.0) - 1); } @@ -405,8 +403,8 @@ cvm::real colvarbias_alb::restraint_potential(cvm::real k, colvarvalue colvarbias_alb::restraint_force(cvm::real k, - colvar const *x, - colvarvalue const &xcenter) const + colvar const * /* x */, + colvarvalue const & /* xcenter */) const { return k; } diff --git a/lib/colvars/colvarbias_alb.h b/lib/colvars/colvarbias_alb.h index d003c244e9..4d16a4e7e2 100644 --- a/lib/colvars/colvarbias_alb.h +++ b/lib/colvars/colvarbias_alb.h @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index 329b1d9dc0..b59a63f6e3 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -78,7 +78,7 @@ int colvarbias_histogram::init(std::string const &conf) } for (i = 0; i < num_variables(); i++) { - colvars[i]->enable(f_cv_grid); + colvars[i]->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature } grid = new colvar_grid_scalar(); diff --git a/lib/colvars/colvarbias_histogram.h b/lib/colvars/colvarbias_histogram.h index 23565caa5c..2e6c6884fb 100644 --- a/lib/colvars/colvarbias_histogram.h +++ b/lib/colvars/colvarbias_histogram.h @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 8540e4a945..267d556be6 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -38,14 +38,29 @@ colvarbias_meta::colvarbias_meta(char const *key) new_hills_begin = hills.end(); hills_traj_os = NULL; + use_grids = true; + rebin_grids = false; + hills_energy = NULL; + hills_energy_gradients = NULL; + + dump_fes = true; + keep_hills = false; + dump_fes_save = false; + dump_replica_fes = false; + ebmeta_equil_steps = 0L; + + replica_update_freq = 0; + replica_id.clear(); } int colvarbias_meta::init(std::string const &conf) { - colvarbias::init(conf); - colvarbias_ti::init(conf); + int error_code = COLVARS_OK; + + error_code |= colvarbias::init(conf); + error_code |= colvarbias_ti::init(conf); enable(f_cvb_calc_pmf); @@ -77,21 +92,16 @@ int colvarbias_meta::init(std::string const &conf) comm = single_replica; } - // in all cases, the first replica is this bias itself - if (replicas.size() == 0) { - replicas.push_back(this); - } - - get_keyval(conf, "useGrids", use_grids, true); + get_keyval(conf, "useGrids", use_grids, use_grids); if (use_grids) { get_keyval(conf, "gridsUpdateFrequency", grids_freq, new_hill_freq); - get_keyval(conf, "rebinGrids", rebin_grids, false); + get_keyval(conf, "rebinGrids", rebin_grids, rebin_grids); expand_grids = false; size_t i; for (i = 0; i < num_variables(); i++) { - variables(i)->enable(f_cv_grid); + variables(i)->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature if (variables(i)->expand_boundaries) { expand_grids = true; cvm::log("Metadynamics bias \""+this->name+"\""+ @@ -101,67 +111,24 @@ int colvarbias_meta::init(std::string const &conf) } } - get_keyval(conf, "keepHills", keep_hills, false); - if (! get_keyval(conf, "writeFreeEnergyFile", dump_fes, true)) - get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent); - if (get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false, colvarparse::parse_silent)) { - cvm::log("Option \"saveFreeEnergyFile\" is deprecated, " - "please use \"keepFreeEnergyFiles\" instead."); - } + get_keyval(conf, "writeFreeEnergyFile", dump_fes, dump_fes); + + get_keyval(conf, "keepHills", keep_hills, keep_hills); get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save); hills_energy = new colvar_grid_scalar(colvars); hills_energy_gradients = new colvar_grid_gradient(colvars); + } else { - rebin_grids = false; - keep_hills = false; + dump_fes = false; - dump_fes_save = false; - dump_replica_fes = false; - - hills_energy = NULL; - hills_energy_gradients = NULL; - } - - if (comm != single_replica) { - - if (expand_grids) - cvm::fatal_error("Error: expandBoundaries is not supported when " - "using more than one replicas; please allocate " - "wide enough boundaries for each colvar" - "ahead of time.\n"); - - if (get_keyval(conf, "dumpPartialFreeEnergyFile", dump_replica_fes, false)) { - if (dump_replica_fes && (! dump_fes)) { - cvm::log("Enabling \"dumpFreeEnergyFile\".\n"); - } - } - - get_keyval(conf, "replicaID", replica_id, std::string("")); - if (!replica_id.size()) - cvm::error("Error: replicaID must be defined " - "when using more than one replica.\n", INPUT_ERROR); - - get_keyval(conf, "replicasRegistry", - replicas_registry_file, - (this->name+".replicas.registry.txt")); - - get_keyval(conf, "replicaUpdateFrequency", - replica_update_freq, new_hill_freq); - - if (keep_hills) - cvm::log("Warning: in metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": keepHills with more than one replica can lead to a very " - "large amount of input/output and slow down your calculations. " - "Please consider disabling it.\n"); - } get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false); - init_well_tempered_params(conf); - init_ebmeta_params(conf); + error_code |= init_replicas_params(conf); + error_code |= init_well_tempered_params(conf); + error_code |= init_ebmeta_params(conf); if (cvm::debug()) cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+ @@ -171,6 +138,72 @@ int colvarbias_meta::init(std::string const &conf) } +int colvarbias_meta::init_replicas_params(std::string const &conf) +{ + colvarproxy *proxy = cvm::main()->proxy; + + // in all cases, the first replica is this bias itself + if (replicas.size() == 0) { + replicas.push_back(this); + } + + if (comm != single_replica) { + + if (!get_keyval(conf, "writePartialFreeEnergyFile", + dump_replica_fes, dump_replica_fes)) { + get_keyval(conf, "dumpPartialFreeEnergyFile", dump_replica_fes, + dump_replica_fes, colvarparse::parse_silent); + } + + if (dump_replica_fes && (! dump_fes)) { + dump_fes = true; + cvm::log("Enabling \"writeFreeEnergyFile\".\n"); + } + + get_keyval(conf, "replicaID", replica_id, replica_id); + if (!replica_id.size()) { + if (proxy->replica_enabled() == COLVARS_OK) { + // Obtain replicaID from the communicator + replica_id = cvm::to_str(proxy->replica_index()); + cvm::log("Setting replicaID from communication layer: replicaID = "+ + replica_id+".\n"); + } else { + return cvm::error("Error: using more than one replica, but replicaID " + "could not be obtained.\n", INPUT_ERROR); + } + } + + get_keyval(conf, "replicasRegistry", replicas_registry_file, + replicas_registry_file); + if (!replicas_registry_file.size()) { + return cvm::error("Error: the name of the \"replicasRegistry\" file " + "must be provided.\n", INPUT_ERROR); + } + + get_keyval(conf, "replicaUpdateFrequency", + replica_update_freq, replica_update_freq); + if (replica_update_freq == 0) { + return cvm::error("Error: replicaUpdateFrequency must be positive.\n", + INPUT_ERROR); + } + + if (expand_grids) { + return cvm::error("Error: expandBoundaries is not supported when " + "using more than one replicas; please allocate " + "wide enough boundaries for each colvar" + "ahead of time.\n", INPUT_ERROR); + } + + if (keep_hills) { + return cvm::error("Error: multipleReplicas and keepHills are not " + "supported together.\n", INPUT_ERROR); + } + } + + return COLVARS_OK; +} + + int colvarbias_meta::init_well_tempered_params(std::string const &conf) { // for well-tempered metadynamics @@ -379,8 +412,8 @@ int colvarbias_meta::update() error_code |= replica_share(); } - error_code |= calc_energy(); - error_code |= calc_forces(); + error_code |= calc_energy(NULL); + error_code |= calc_forces(NULL); return error_code; } @@ -418,7 +451,7 @@ int colvarbias_meta::update_grid_params() int &new_size = new_sizes[i]; bool changed_lb = false, changed_ub = false; - if (!variables(i)->hard_lower_boundary) + if (!variables(i)->is_enabled(f_cv_hard_lower_boundary)) if (curr_bin[i] < min_buffer) { int const extra_points = (min_buffer - curr_bin[i]); new_lb -= extra_points * variables(i)->width; @@ -435,7 +468,7 @@ int colvarbias_meta::update_grid_params() cvm::to_str(new_lower_boundaries[i])+".\n"); } - if (!variables(i)->hard_upper_boundary) + if (!variables(i)->is_enabled(f_cv_hard_upper_boundary)) if (curr_bin[i] > new_size - min_buffer - 1) { int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer); new_ub += extra_points * variables(i)->width; @@ -521,7 +554,7 @@ int colvarbias_meta::update_bias() std::vector curr_bin = hills_energy->get_colvars_index(); hills_energy_sum_here = hills_energy->value(curr_bin); } else { - calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here); + calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here, NULL); } hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); } @@ -578,7 +611,7 @@ int colvarbias_meta::update_grid_data() } -int colvarbias_meta::calc_energy(std::vector const &values) +int colvarbias_meta::calc_energy(std::vector const *values) { size_t ir = 0; @@ -586,8 +619,8 @@ int colvarbias_meta::calc_energy(std::vector const &values) replicas[ir]->bias_energy = 0.0; } - std::vector const curr_bin = values.size() ? - hills_energy->get_colvars_index(values) : + std::vector const curr_bin = values ? + hills_energy->get_colvars_index(*values) : hills_energy->get_colvars_index(); if (hills_energy->index_ok(curr_bin)) { @@ -619,7 +652,8 @@ int colvarbias_meta::calc_energy(std::vector const &values) for (ir = 0; ir < replicas.size(); ir++) { calc_hills(replicas[ir]->new_hills_begin, replicas[ir]->hills.end(), - bias_energy); + bias_energy, + values); if (cvm::debug()) { cvm::log("Hills energy = "+cvm::to_str(bias_energy)+".\n"); } @@ -629,7 +663,7 @@ int colvarbias_meta::calc_energy(std::vector const &values) } -int colvarbias_meta::calc_forces(std::vector const &values) +int colvarbias_meta::calc_forces(std::vector const *values) { size_t ir = 0, ic = 0; for (ir = 0; ir < replicas.size(); ir++) { @@ -638,8 +672,8 @@ int colvarbias_meta::calc_forces(std::vector const &values) } } - std::vector const curr_bin = values.size() ? - hills_energy->get_colvars_index(values) : + std::vector const curr_bin = values ? + hills_energy->get_colvars_index(*values) : hills_energy->get_colvars_index(); if (hills_energy->index_ok(curr_bin)) { @@ -693,7 +727,7 @@ int colvarbias_meta::calc_forces(std::vector const &values) void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, cvm::real &energy, - std::vector const &colvar_values) + std::vector const *values) { size_t i = 0; std::vector curr_values(num_variables()); @@ -701,9 +735,9 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, curr_values[i].type(variables(i)->value()); } - if (colvar_values.size()) { + if (values) { for (i = 0; i < num_variables(); i++) { - curr_values[i] = colvar_values[i]; + curr_values[i] = (*values)[i]; } } else { for (i = 0; i < num_variables(); i++) { @@ -738,10 +772,10 @@ void colvarbias_meta::calc_hills_force(size_t const &i, colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, std::vector &forces, - std::vector const &values) + std::vector const *values) { // Retrieve the value of the colvar - colvarvalue const x(values.size() ? values[i] : variables(i)->value()); + colvarvalue const x(values ? (*values)[i] : variables(i)->value()); // do the type check only once (all colvarvalues in the hills series // were already saved with their types matching those in the @@ -846,12 +880,12 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, // loop over the hills and increment the energy grid locally hills_energy_here = 0.0; - calc_hills(h_first, h_last, hills_energy_here, colvar_values); + calc_hills(h_first, h_last, hills_energy_here, &colvar_values); he->acc_value(he_ix, hills_energy_here); for (i = 0; i < num_variables(); i++) { hills_forces_here[i].reset(); - calc_hills_force(i, h_first, h_last, hills_forces_here, colvar_values); + calc_hills_force(i, h_first, h_last, hills_forces_here, &colvar_values); colvar_forces_scalar[i] = hills_forces_here[i].real_value; } hg->acc_force(hg_ix, &(colvar_forces_scalar.front())); @@ -883,7 +917,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, } hills_energy_here = 0.0; - calc_hills(h_first, h_last, hills_energy_here, colvar_values); + calc_hills(h_first, h_last, hills_energy_here, &colvar_values); he->acc_value(he_ix, hills_energy_here); he->incr(he_ix); @@ -915,7 +949,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first, colvarbias_meta::hill_iter h_last, - colvar_grid_scalar *he) + colvar_grid_scalar * /* he */) { hills_off_grid.clear(); @@ -985,11 +1019,6 @@ void colvarbias_meta::update_replicas_registry() if (new_replica == this->replica_id) { // this is the record for this same replica, skip it - if (cvm::debug()) - cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": skipping this replica's own record: \""+ - new_replica+"\", \""+new_replica_file+"\"\n"); new_replica_file.clear(); new_replica.clear(); continue; @@ -1039,11 +1068,12 @@ void colvarbias_meta::update_replicas_registry() (replicas.back())->enable(f_cvb_calc_ti_samples); (replicas.back())->colvarbias_ti::init_grids(); } + (replicas.back())->update_status = 1; } } } else { - cvm::fatal_error("Error: cannot read the replicas registry file \""+ - replicas_registry+"\".\n"); + cvm::error("Error: cannot read the replicas registry file \""+ + replicas_registry+"\".\n", FILE_ERROR); } // now (re)read the list file of each replica @@ -1068,7 +1098,6 @@ void colvarbias_meta::update_replicas_registry() cvm::to_str(replica_update_freq)+" steps.\n"); (replicas[ir])->update_status++; } else { - (replicas[ir])->update_status = 0; if (new_state_file != (replicas[ir])->replica_state_file) { cvm::log("Metadynamics bias \""+this->name+"\""+ ": replica \""+(replicas[ir])->replica_id+ @@ -1092,27 +1121,41 @@ void colvarbias_meta::read_replica_files() // Note: we start from the 2nd replica. for (size_t ir = 1; ir < replicas.size(); ir++) { - if (! (replicas[ir])->replica_state_file_in_sync) { - // if a new state file is being read, the hills file is also new - (replicas[ir])->replica_hills_file_pos = 0; - } - // (re)read the state file if necessary if ( (! (replicas[ir])->has_data) || (! (replicas[ir])->replica_state_file_in_sync) ) { - - cvm::log("Metadynamics bias \""+this->name+"\""+ - ": reading the state of replica \""+ - (replicas[ir])->replica_id+"\" from file \""+ - (replicas[ir])->replica_state_file+"\".\n"); - - std::ifstream is((replicas[ir])->replica_state_file.c_str()); - if ((replicas[ir])->read_state(is)) { - // state file has been read successfully - (replicas[ir])->replica_state_file_in_sync = true; - (replicas[ir])->update_status = 0; + if ((replicas[ir])->replica_state_file.size()) { + cvm::log("Metadynamics bias \""+this->name+"\""+ + ": reading the state of replica \""+ + (replicas[ir])->replica_id+"\" from file \""+ + (replicas[ir])->replica_state_file+"\".\n"); + std::ifstream is((replicas[ir])->replica_state_file.c_str()); + if ((replicas[ir])->read_state(is)) { + // state file has been read successfully + (replicas[ir])->replica_state_file_in_sync = true; + (replicas[ir])->update_status = 0; + } else { + cvm::log("Failed to read the file \""+ + (replicas[ir])->replica_state_file+ + "\": will try again in "+ + cvm::to_str(replica_update_freq)+" steps.\n"); + (replicas[ir])->replica_state_file_in_sync = false; + (replicas[ir])->update_status++; + } + is.close(); + } else { + cvm::log("Metadynamics bias \""+this->name+"\""+ + ": the state file of replica \""+ + (replicas[ir])->replica_id+"\" is currently undefined: " + "will try again after "+ + cvm::to_str(replica_update_freq)+" steps.\n"); + (replicas[ir])->update_status++; } - is.close(); + } + + if (! (replicas[ir])->replica_state_file_in_sync) { + // if a new state file is being read, the hills file is also new + (replicas[ir])->replica_hills_file_pos = 0; } // now read the hills added after writing the state file @@ -1124,14 +1167,16 @@ void colvarbias_meta::read_replica_files() (replicas[ir])->replica_id+"\" in the file \""+ (replicas[ir])->replica_hills_file+"\".\n"); - // read hills from the other replicas' files; for each file, resume - // the position recorded previously + // read hills from the other replicas' files std::ifstream is((replicas[ir])->replica_hills_file.c_str()); if (is.is_open()) { - // try to resume the previous position - is.seekg((replicas[ir])->replica_hills_file_pos, std::ios::beg); + // try to resume the previous position (if not the beginning) + if ((replicas[ir])->replica_hills_file_pos > 0) { + is.seekg((replicas[ir])->replica_hills_file_pos, std::ios::beg); + } + if (!is.is_open()){ // if fail (the file may have been overwritten), reset this // position @@ -1149,7 +1194,6 @@ void colvarbias_meta::read_replica_files() } else { while ((replicas[ir])->read_hill(is)) { - // if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ ": received a hill from replica \""+ (replicas[ir])->replica_id+ @@ -1159,11 +1203,13 @@ void colvarbias_meta::read_replica_files() is.clear(); // store the position for the next read (replicas[ir])->replica_hills_file_pos = is.tellg(); - if (cvm::debug()) + if (cvm::debug()) { cvm::log("Metadynamics bias \""+this->name+"\""+ - ": stopped reading file \""+(replicas[ir])->replica_hills_file+ + ": stopped reading file \""+ + (replicas[ir])->replica_hills_file+ "\" at position "+ cvm::to_str((replicas[ir])->replica_hills_file_pos)+".\n"); + } // test whether this is the end of the file is.seekg(0, std::ios::end); @@ -1175,12 +1221,11 @@ void colvarbias_meta::read_replica_files() } } else { - cvm::log("Failed to read the file \""+(replicas[ir])->replica_hills_file+ + cvm::log("Failed to read the file \""+ + (replicas[ir])->replica_hills_file+ "\": will try again in "+ cvm::to_str(replica_update_freq)+" steps.\n"); (replicas[ir])->update_status++; - // cvm::fatal_error ("Error: cannot read from file \""+ - // (replicas[ir])->replica_hills_file+"\".\n"); } is.close(); } @@ -1188,8 +1233,8 @@ void colvarbias_meta::read_replica_files() size_t const n_flush = (replica_update_freq/new_hill_freq + 1); if ((replicas[ir])->update_status > 3*n_flush) { // TODO: suspend the calculation? - cvm::log("WARNING: in metadynamics bias \""+this->name+"\""+ - " failed to read completely the output of replica \""+ + cvm::log("WARNING: metadynamics bias \""+this->name+"\""+ + " could not read information from replica \""+ (replicas[ir])->replica_id+ "\" after more than "+ cvm::to_str((replicas[ir])->update_status * replica_update_freq)+ @@ -1712,7 +1757,7 @@ void colvarbias_meta::write_pmf() cvm:: real target_val=target_dist->value(i); if (target_val>0) { pmf_val=pmf->value(i); - pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * std::log(target_val); + pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val); } pmf->set_value(i,pmf_val); } @@ -1753,7 +1798,7 @@ void colvarbias_meta::write_pmf() cvm:: real target_val=target_dist->value(i); if (target_val>0) { pmf_val=pmf->value(i); - pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * std::log(target_val); + pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val); } pmf->set_value(i,pmf_val); } diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 8e98274b2d..abf1417c32 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -41,6 +41,7 @@ public: virtual ~colvarbias_meta(); virtual int init(std::string const &conf); + virtual int init_replicas_params(std::string const &conf); virtual int init_well_tempered_params(std::string const &conf); virtual int init_ebmeta_params(std::string const &conf); @@ -52,10 +53,8 @@ public: virtual int update_grid_data(); virtual int replica_share(); - virtual int calc_energy(std::vector const &values = - std::vector(0)); - virtual int calc_forces(std::vector const &values = - std::vector(0)); + virtual int calc_energy(std::vector const *values); + virtual int calc_forces(std::vector const *values); virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &state_conf); @@ -125,7 +124,7 @@ protected: virtual void calc_hills(hill_iter h_first, hill_iter h_last, cvm::real &energy, - std::vector const &values = std::vector(0)); + std::vector const *values); /// \brief Calculate the forces acting on the i-th colvar, /// incrementing colvar_forces[i]; must be called after calc_hills @@ -134,7 +133,7 @@ protected: hill_iter h_first, hill_iter h_last, std::vector &forces, - std::vector const &values = std::vector(0)); + std::vector const *values); /// Height of new hills diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index ab02820cf0..1fd9c84e98 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -174,7 +174,7 @@ int colvarbias_restraint_k::change_configuration(std::string const &conf) -colvarbias_restraint_moving::colvarbias_restraint_moving(char const *key) +colvarbias_restraint_moving::colvarbias_restraint_moving(char const * /* key */) { target_nstages = 0; target_nsteps = 0L; diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index 6493f7f16b..644835f129 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -26,10 +26,10 @@ public: virtual int update(); /// Load new configuration - force constant and/or centers only - virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; } + virtual int change_configuration(std::string const & /* conf */) { return COLVARS_NOT_IMPLEMENTED; } /// Calculate change in energy from using alternate configuration - virtual cvm::real energy_difference(std::string const &conf) { return 0.0; } + virtual cvm::real energy_difference(std::string const & /* conf */) { return 0.0; } virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &conf); @@ -107,7 +107,7 @@ public: // Note: despite the diamond inheritance, most of this function gets only executed once virtual int init(std::string const &conf); virtual int update() { return COLVARS_OK; } - virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; } + virtual int change_configuration(std::string const & /* conf */) { return COLVARS_NOT_IMPLEMENTED; } virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &conf); @@ -149,7 +149,7 @@ public: colvarbias_restraint_centers_moving(char const *key); virtual int init(std::string const &conf); virtual int update(); - virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; } + virtual int change_configuration(std::string const & /* conf */) { return COLVARS_NOT_IMPLEMENTED; } virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &conf); @@ -188,7 +188,7 @@ public: colvarbias_restraint_k_moving(char const *key); virtual int init(std::string const &conf); virtual int update(); - virtual int change_configuration(std::string const &conf) { return COLVARS_NOT_IMPLEMENTED; } + virtual int change_configuration(std::string const & /* conf */) { return COLVARS_NOT_IMPLEMENTED; } virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &conf); diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index 431884a877..05d1195d58 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -17,30 +17,28 @@ colvar::cvc::cvc() - : sup_coeff(1.0), - sup_np(1), - b_periodic(false), - b_try_scalable(true) { description = "uninitialized colvar component"; - init_dependencies(); + b_try_scalable = true; sup_coeff = 1.0; + sup_np = 1; period = 0.0; wrap_center = 0.0; + width = 0.0; + init_dependencies(); } colvar::cvc::cvc(std::string const &conf) - : sup_coeff(1.0), - sup_np(1), - b_periodic(false), - b_try_scalable(true) { description = "uninitialized colvar component"; - init_dependencies(); + b_try_scalable = true; sup_coeff = 1.0; + sup_np = 1; period = 0.0; wrap_center = 0.0; + width = 0.0; + init_dependencies(); init(conf); } @@ -72,11 +70,17 @@ int colvar::cvc::init(std::string const &conf) get_keyval(conf, "componentCoeff", sup_coeff, sup_coeff); get_keyval(conf, "componentExp", sup_np, sup_np); + // TODO these could be condensed into get_keyval() + register_param("componentCoeff", reinterpret_cast(&sup_coeff)); + register_param("componentExp", reinterpret_cast(&sup_np)); get_keyval(conf, "period", period, period); get_keyval(conf, "wrapAround", wrap_center, wrap_center); + // TODO when init() is called after all constructors, check periodic flag + register_param("period", reinterpret_cast(&period)); + register_param("wrapAround", reinterpret_cast(&wrap_center)); - get_keyval_feature(dynamic_cast(this), conf, "debugGradients", + get_keyval_feature(this, conf, "debugGradients", f_cvc_debug_gradient, false, parse_silent); bool b_no_PBC = !is_enabled(f_cvc_pbc_minimum_image); // Enabled by default @@ -194,6 +198,14 @@ int colvar::cvc::init_dependencies() { init_feature(f_cvc_scalar, "scalar", f_type_static); + init_feature(f_cvc_periodic, "periodic", f_type_static); + + init_feature(f_cvc_width, "defined width", f_type_static); + + init_feature(f_cvc_lower_boundary, "defined lower boundary", f_type_static); + + init_feature(f_cvc_upper_boundary, "defined upper boundary", f_type_static); + init_feature(f_cvc_gradient, "gradient", f_type_dynamic); init_feature(f_cvc_explicit_gradient, "explicit gradient", f_type_static); @@ -209,15 +221,15 @@ int colvar::cvc::init_dependencies() { init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic); require_feature_self(f_cvc_Jacobian, f_cvc_inv_gradient); - init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static); - - init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user); - // Compute total force on first site only to avoid unwanted // coupling to other colvars (see e.g. Ciccotti et al., 2005) init_feature(f_cvc_one_site_total_force, "compute total force from one group", f_type_user); require_feature_self(f_cvc_one_site_total_force, f_cvc_com_based); + init_feature(f_cvc_com_based, "depends on group centers of mass", f_type_static); + + init_feature(f_cvc_pbc_minimum_image, "use minimum-image distances with PBCs", f_type_user); + init_feature(f_cvc_scalable, "scalable calculation", f_type_static); require_feature_self(f_cvc_scalable, f_cvc_scalable_com); @@ -291,6 +303,77 @@ colvar::cvc::~cvc() } +void colvar::cvc::init_as_distance() +{ + x.type(colvarvalue::type_scalar); + enable(f_cvc_lower_boundary); + lower_boundary.type(colvarvalue::type_scalar); + lower_boundary.real_value = 0.0; + register_param("lowerBoundary", reinterpret_cast(&lower_boundary)); +} + + +void colvar::cvc::init_as_angle() +{ + x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, 180); +} + + +void colvar::cvc::init_scalar_boundaries(cvm::real lb, cvm::real ub) +{ + enable(f_cvc_lower_boundary); + lower_boundary.type(colvarvalue::type_scalar); + lower_boundary.real_value = lb; + enable(f_cvc_upper_boundary); + upper_boundary.type(colvarvalue::type_scalar); + upper_boundary.real_value = ub; + register_param("lowerBoundary", reinterpret_cast(&lower_boundary)); + register_param("upperBoundary", reinterpret_cast(&upper_boundary)); +} + + +void colvar::cvc::register_atom_group(cvm::atom_group *ag) +{ + atom_groups.push_back(ag); + add_child(ag); +} + + +colvarvalue const *colvar::cvc::get_param_grad(std::string const ¶m_name) +{ + colvarvalue const *ptr = + reinterpret_cast(get_param_grad_ptr(param_name)); + return ptr != NULL ? ptr : NULL; +} + + +int colvar::cvc::set_param(std::string const ¶m_name, + void const *new_value) +{ + if (param_map.count(param_name) > 0) { + + // TODO When we can use C++11, make this a proper function map + if (param_name.compare("componentCoeff") == 0) { + sup_coeff = *(reinterpret_cast(new_value)); + } + if (param_name.compare("componentExp") == 0) { + sup_np = *(reinterpret_cast(new_value)); + } + if (is_enabled(f_cvc_periodic)) { + if (param_name.compare("period") == 0) { + period = *(reinterpret_cast(new_value)); + } + if (param_name.compare("wrapAround") == 0) { + wrap_center = *(reinterpret_cast(new_value)); + } + } + } + + return colvarparams::set_param(param_name, new_value); +} + + void colvar::cvc::read_data() { size_t ig; @@ -537,7 +620,7 @@ colvarvalue colvar::cvc::dist2_rgrad(colvarvalue const &x1, } -void colvar::cvc::wrap(colvarvalue &x_unwrapped) const +void colvar::cvc::wrap(colvarvalue & /* x_unwrapped */) const { return; } diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 1d95cbe600..b01c7e48e1 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -23,11 +23,13 @@ #include "colvarmodule.h" #include "colvar.h" #include "colvaratoms.h" +#include "colvar_arithmeticpath.h" #if (__cplusplus >= 201103L) +// C++11-only functions #include "colvar_geometricpath.h" #include -#endif // C++11 checking +#endif #include @@ -97,10 +99,7 @@ public: /// \brief Exponent in the polynomial combination (default: 1) int sup_np; - /// \brief Is this a periodic component? - bool b_periodic; - - /// \brief Period of this cvc value, (default: 0.0, non periodic) + /// \brief Period of the values of this CVC (default: 0.0, non periodic) cvm::real period; /// \brief If the component is periodic, wrap around this value (default: 0.0) @@ -261,10 +260,13 @@ public: std::vector atom_groups; /// \brief Store a pointer to new atom group, and list as child for dependencies - inline void register_atom_group(cvm::atom_group *ag) { - atom_groups.push_back(ag); - add_child(ag); - } + void register_atom_group(cvm::atom_group *ag); + + /// Pointer to the gradient of parameter param_name + virtual colvarvalue const *get_param_grad(std::string const ¶m_name); + + /// Set the named parameter to the given value + virtual int set_param(std::string const ¶m_name, void const *new_value); /// \brief Whether or not this CVC will be computed in parallel whenever possible bool b_try_scalable; @@ -286,6 +288,24 @@ protected: /// \brief Calculated Jacobian derivative (divergence of the inverse /// gradients): serves to calculate the phase space correction colvarvalue jd; + + /// \brief Set data types for a scalar distance (convenience function) + void init_as_distance(); + + /// \brief Set data types for a bounded angle (convenience function) + void init_as_angle(); + + /// \brief Set two scalar boundaries (convenience function) + void init_scalar_boundaries(cvm::real lb, cvm::real ub); + + /// \brief Location of the lower boundary (not defined by user choice) + colvarvalue lower_boundary; + + /// \brief Location of the upper boundary (not defined by user choice) + colvarvalue upper_boundary; + + /// \brief CVC-specific default colvar width + cvm::real width; }; @@ -538,7 +558,6 @@ protected: bool b_no_PBC; public: distance_inv(std::string const &conf); - distance_inv(); virtual ~distance_inv() {} virtual void calc_value(); virtual void calc_gradients(); @@ -614,9 +633,7 @@ protected: /// Atoms involved cvm::atom_group *atoms; public: - /// Constructor gyration(std::string const &conf); - gyration(); virtual ~gyration() {} virtual void calc_value(); virtual void calc_gradients(); @@ -752,7 +769,6 @@ public: angle(std::string const &conf); /// \brief Initialize the three groups after three atoms angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3); - angle(); virtual ~angle() {} virtual void calc_value(); virtual void calc_gradients(); @@ -902,7 +918,6 @@ protected: public: coordnum(std::string const &conf); - coordnum(); ~coordnum(); virtual void calc_value(); @@ -972,7 +987,6 @@ protected: public: selfcoordnum(std::string const &conf); - selfcoordnum(); ~selfcoordnum(); virtual void calc_value(); virtual void calc_gradients(); @@ -1011,7 +1025,6 @@ protected: public: /// Constructor groupcoordnum(std::string const &conf); - groupcoordnum(); virtual ~groupcoordnum() {} virtual void calc_value(); virtual void calc_gradients(); @@ -1232,7 +1245,6 @@ class colvar::orientation_angle public: orientation_angle(std::string const &conf); - orientation_angle(); virtual int init(std::string const &conf); virtual ~orientation_angle() {} virtual void calc_value(); @@ -1285,7 +1297,6 @@ protected: public: tilt(std::string const &conf); - tilt(); virtual int init(std::string const &conf); virtual ~tilt() {} virtual void calc_value(); @@ -1394,13 +1405,13 @@ class colvar::componentDisabled : public colvar::cvc { public: - componentDisabled(std::string const &conf) { + componentDisabled(std::string const & /* conf */) { cvm::error("Error: this component is not enabled in the current build; please see https://colvars.github.io/README-c++11.html"); } virtual ~componentDisabled() {} virtual void calc_value() {} virtual void calc_gradients() {} - virtual void apply_force(colvarvalue const &force) {} + virtual void apply_force(colvarvalue const & /* force */) {} }; @@ -1410,7 +1421,7 @@ class colvar::CartesianBasedPath : public colvar::cvc { protected: - virtual void computeReferenceDistance(std::vector& result); + virtual void computeDistanceToReferenceFrames(std::vector& result); /// Selected atoms cvm::atom_group *atoms; /// Fitting options @@ -1440,7 +1451,7 @@ private: cvm::rotation rot_v3; protected: virtual void prepareVectors(); - virtual void updateReferenceDistances(); + virtual void updateDistanceToReferenceFrames(); public: gspath(std::string const &conf); virtual ~gspath() {} @@ -1462,7 +1473,7 @@ private: cvm::rotation rot_v4; protected: virtual void prepareVectors(); - virtual void updateReferenceDistances(); + virtual void updateDistanceToReferenceFrames(); public: gzpath(std::string const &conf); virtual ~gzpath() {} @@ -1508,7 +1519,9 @@ protected: /// Total number of reference frames size_t total_reference_frames; protected: - virtual void computeReferenceDistance(std::vector& result); + virtual void computeDistanceToReferenceFrames(std::vector& result); + /// Helper function to determine the distance between reference frames + virtual void computeDistanceBetweenReferenceFrames(std::vector& result) const; cvm::real getPolynomialFactorOfCVGradient(size_t i_cv) const; public: CVBasedPath(std::string const &conf); @@ -1526,7 +1539,7 @@ class colvar::gspathCV : public colvar::CVBasedPath, public GeometricPathCV::GeometricPathBase { protected: - virtual void updateReferenceDistances(); + virtual void updateDistanceToReferenceFrames(); virtual void prepareVectors(); public: gspathCV(std::string const &conf); @@ -1542,7 +1555,7 @@ class colvar::gzpathCV : public colvar::CVBasedPath, public GeometricPathCV::GeometricPathBase { protected: - virtual void updateReferenceDistances(); + virtual void updateDistanceToReferenceFrames(); virtual void prepareVectors(); public: gzpathCV(std::string const &conf); @@ -1552,6 +1565,35 @@ public: virtual void apply_force(colvarvalue const &force); }; + + +class colvar::aspathCV + : public colvar::CVBasedPath, public ArithmeticPathCV::ArithmeticPathBase +{ +protected: + virtual void updateDistanceToReferenceFrames(); +public: + aspathCV(std::string const &conf); + virtual ~aspathCV(); + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); +}; + + +class colvar::azpathCV + : public colvar::CVBasedPath, public ArithmeticPathCV::ArithmeticPathBase +{ +protected: + virtual void updateDistanceToReferenceFrames(); +public: + azpathCV(std::string const &conf); + virtual ~azpathCV(); + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); +}; + #else // if the compiler doesn't support C++11 class colvar::linearCombination @@ -1603,6 +1645,20 @@ public: gzpathCV(std::string const &conf) : componentDisabled(conf) {} }; +class colvar::aspathCV + : public colvar::componentDisabled +{ +public: + aspathCV(std::string const &conf) : componentDisabled(conf) {} +}; + +class colvar::azpathCV + : public colvar::componentDisabled +{ +public: + azpathCV(std::string const &conf) : componentDisabled(conf) {} +}; + #endif // C++11 checking // metrics functions for cvc implementations diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 97fb23b181..28d7677147 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -16,6 +16,8 @@ colvar::angle::angle(std::string const &conf) : cvc(conf) { function_type = "angle"; + init_as_angle(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); @@ -25,8 +27,6 @@ colvar::angle::angle(std::string const &conf) group3 = parse_group(conf, "group3"); init_total_force_params(conf); - - x.type(colvarvalue::type_scalar); } @@ -35,6 +35,8 @@ colvar::angle::angle(cvm::atom const &a1, cvm::atom const &a3) { function_type = "angle"; + init_as_angle(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); @@ -45,15 +47,6 @@ colvar::angle::angle(cvm::atom const &a1, register_atom_group(group1); register_atom_group(group2); register_atom_group(group3); - - x.type(colvarvalue::type_scalar); -} - - -colvar::angle::angle() -{ - function_type = "angle"; - x.type(colvarvalue::type_scalar); } @@ -148,13 +141,13 @@ colvar::dipole_angle::dipole_angle(std::string const &conf) : cvc(conf) { function_type = "dipole_angle"; + init_as_angle(); + group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); group3 = parse_group(conf, "group3"); init_total_force_params(conf); - - x.type(colvarvalue::type_scalar); } @@ -178,6 +171,7 @@ colvar::dipole_angle::dipole_angle(cvm::atom const &a1, colvar::dipole_angle::dipole_angle() { function_type = "dipole_angle"; + init_as_angle(); x.type(colvarvalue::type_scalar); } @@ -259,7 +253,7 @@ colvar::dihedral::dihedral(std::string const &conf) { function_type = "dihedral"; period = 360.0; - b_periodic = true; + enable(f_cvc_periodic); provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); @@ -285,7 +279,7 @@ colvar::dihedral::dihedral(cvm::atom const &a1, function_type = "dihedral"; period = 360.0; - b_periodic = true; + enable(f_cvc_periodic); provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); @@ -312,7 +306,7 @@ colvar::dihedral::dihedral() { function_type = "dihedral"; period = 360.0; - b_periodic = true; + enable(f_cvc_periodic); provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); x.type(colvarvalue::type_scalar); diff --git a/lib/colvars/colvarcomp_apath.cpp b/lib/colvars/colvarcomp_apath.cpp new file mode 100644 index 0000000000..591d8d0012 --- /dev/null +++ b/lib/colvars/colvarcomp_apath.cpp @@ -0,0 +1,192 @@ +#if (__cplusplus >= 201103L) + +#include +#include +#include +#include +#include + +#include "colvarmodule.h" +#include "colvarvalue.h" +#include "colvarparse.h" +#include "colvar.h" +#include "colvarcomp.h" + +colvar::aspathCV::aspathCV(std::string const &conf): CVBasedPath(conf) { + function_type = "aspathCV"; + cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); + std::vector p_weights(cv.size(), 1.0); + 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); + 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) { + if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { + use_explicit_gradients = false; + } + cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(weights[i_cv]) + std::string("\n")); + } +} + +void colvar::aspathCV::updateDistanceToReferenceFrames() { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + cv[i_cv]->calc_value(); + } + for (size_t i_frame = 0; i_frame < ref_cv.size(); ++i_frame) { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]); + colvarvalue current_cv_value(cv[i_cv]->value()); + if (current_cv_value.type() == colvarvalue::type_scalar) { + frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)), ref_cv_value.real_value); + } else { + frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value); + } + } + } +} + +void colvar::aspathCV::calc_value() { + computeValue(); + x = s; +} + +void colvar::aspathCV::calc_gradients() { + computeDerivatives(); + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + cv[i_cv]->calc_gradients(); + if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) && + !cv[i_cv]->is_enabled(f_cvc_scalable) && + !cv[i_cv]->is_enabled(f_cvc_scalable_com)) { + cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) { + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = dsdx[i_cv][j_elem] * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; + } + } + } + } + } +} + +void colvar::aspathCV::apply_force(colvarvalue const &force) { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) && + !cv[i_cv]->is_enabled(f_cvc_scalable) && + !cv[i_cv]->is_enabled(f_cvc_scalable_com) + ) { + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); + } + } else { + cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + colvarvalue cv_force = dsdx[i_cv] * force.real_value * factor_polynomial; + cv[i_cv]->apply_force(cv_force); + } + } +} + +colvar::aspathCV::~aspathCV() {} + +colvar::azpathCV::azpathCV(std::string const &conf): CVBasedPath(conf) { + function_type = "azpathCV"; + cvm::log(std::string("Total number of frames: ") + cvm::to_str(total_reference_frames) + std::string("\n")); + std::vector p_weights(cv.size(), 1.0); + 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); + 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) { + if (!cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { + use_explicit_gradients = false; + } + cvm::log(std::string("The weight of CV ") + cvm::to_str(i_cv) + std::string(" is ") + cvm::to_str(weights[i_cv]) + std::string("\n")); + } +} + +void colvar::azpathCV::updateDistanceToReferenceFrames() { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + cv[i_cv]->calc_value(); + } + for (size_t i_frame = 0; i_frame < ref_cv.size(); ++i_frame) { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]); + colvarvalue current_cv_value(cv[i_cv]->value()); + if (current_cv_value.type() == colvarvalue::type_scalar) { + frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)), ref_cv_value.real_value); + } else { + frame_element_distances[i_frame][i_cv] = 0.5 * cv[i_cv]->dist2_lgrad(cv[i_cv]->sup_coeff * current_cv_value, ref_cv_value); + } + } + } +} + +void colvar::azpathCV::calc_value() { + computeValue(); + x = z; +} + +void colvar::azpathCV::calc_gradients() { + computeDerivatives(); + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + cv[i_cv]->calc_gradients(); + if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) && + !cv[i_cv]->is_enabled(f_cvc_scalable) && + !cv[i_cv]->is_enabled(f_cvc_scalable_com)) { + cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + for (size_t j_elem = 0; j_elem < cv[i_cv]->value().size(); ++j_elem) { + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = dzdx[i_cv][j_elem] * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; + } + } + } + + } + } +} + +void colvar::azpathCV::apply_force(colvarvalue const &force) { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + if ( cv[i_cv]->is_enabled(f_cvc_explicit_gradient) && + !cv[i_cv]->is_enabled(f_cvc_scalable) && + !cv[i_cv]->is_enabled(f_cvc_scalable_com) + ) { + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); + } + } else { + cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + const colvarvalue cv_force = dzdx[i_cv] * force.real_value * factor_polynomial; + cv[i_cv]->apply_force(cv_force); + } + } +} + +colvar::azpathCV::~azpathCV() {} + +#endif diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index b6dd6535ef..f0fbabfc4f 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -97,6 +97,8 @@ colvar::coordnum::coordnum(std::string const &conf) function_type = "coordnum"; x.type(colvarvalue::type_scalar); + colvarproxy *proxy = cvm::main()->proxy; + group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); @@ -119,12 +121,12 @@ colvar::coordnum::coordnum(std::string const &conf) } bool const b_isotropic = get_keyval(conf, "cutoff", r0, - cvm::real(4.0 * cvm::unit_angstrom())); + cvm::real(4.0 * proxy->angstrom_value)); if (get_keyval(conf, "cutoff3", r0_vec, - cvm::rvector(4.0 * cvm::unit_angstrom(), - 4.0 * cvm::unit_angstrom(), - 4.0 * cvm::unit_angstrom()))) { + cvm::rvector(4.0 * proxy->angstrom_value, + 4.0 * proxy->angstrom_value, + 4.0 * proxy->angstrom_value))) { if (b_isotropic) { cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" " "at the same time.\n", @@ -174,14 +176,8 @@ colvar::coordnum::coordnum(std::string const &conf) } } -} - - -colvar::coordnum::coordnum() - : b_anisotropic(false), pairlist(NULL) -{ - function_type = "coordnum"; - x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, b_group2_center_only ? group1->size() : + group1->size() * group2->size()); } @@ -310,10 +306,13 @@ colvar::h_bond::h_bond(std::string const &conf) function_type = "h_bond"; x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, 1.0); - int a_num, d_num; - get_keyval(conf, "acceptor", a_num, -1); - get_keyval(conf, "donor", d_num, -1); + colvarproxy *proxy = cvm::main()->proxy; + + int a_num = -1, d_num = -1; + get_keyval(conf, "acceptor", a_num, a_num); + get_keyval(conf, "donor", d_num, a_num); if ( (a_num == -1) || (d_num == -1) ) { cvm::error("Error: either acceptor or donor undefined.\n"); @@ -326,7 +325,7 @@ colvar::h_bond::h_bond(std::string const &conf) atom_groups[0]->add_atom(acceptor); atom_groups[0]->add_atom(donor); - get_keyval(conf, "cutoff", r0, (3.3 * cvm::unit_angstrom())); + get_keyval(conf, "cutoff", r0, (3.3 * proxy->angstrom_value)); get_keyval(conf, "expNumer", en, 6); get_keyval(conf, "expDenom", ed, 8); @@ -352,6 +351,7 @@ colvar::h_bond::h_bond(cvm::atom const &acceptor, { function_type = "h_bond"; x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, 1.0); register_atom_group(new cvm::atom_group); atom_groups[0]->add_atom(acceptor); @@ -359,14 +359,6 @@ colvar::h_bond::h_bond(cvm::atom const &acceptor, } -colvar::h_bond::h_bond() - : cvc() -{ - function_type = "h_bond"; - x.type(colvarvalue::type_scalar); -} - - void colvar::h_bond::calc_value() { int const flags = coordnum::ef_null; @@ -406,9 +398,11 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf) function_type = "selfcoordnum"; x.type(colvarvalue::type_scalar); + colvarproxy *proxy = cvm::main()->proxy; + group1 = parse_group(conf, "group1"); - get_keyval(conf, "cutoff", r0, cvm::real(4.0 * cvm::unit_angstrom())); + get_keyval(conf, "cutoff", r0, cvm::real(4.0 * proxy->angstrom_value)); get_keyval(conf, "expNumer", en, 6); get_keyval(conf, "expDenom", ed, 12); @@ -437,14 +431,8 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf) } pairlist = new bool[(group1->size()-1) * (group1->size()-1)]; } -} - -colvar::selfcoordnum::selfcoordnum() - : pairlist(NULL) -{ - function_type = "selfcoordnum"; - x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, (group1->size()-1) * (group1->size()-1)); } @@ -553,6 +541,9 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf) { function_type = "groupcoordnum"; x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, 1.0); + + colvarproxy *proxy = cvm::main()->proxy; // group1 and group2 are already initialized by distance() if (group1->b_dummy || group2->b_dummy) { @@ -561,7 +552,7 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf) } bool const b_scale = get_keyval(conf, "cutoff", r0, - cvm::real(4.0 * cvm::unit_angstrom())); + cvm::real(4.0 * proxy->angstrom_value)); if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0, 4.0, 4.0), parse_silent)) { @@ -598,14 +589,6 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf) } -colvar::groupcoordnum::groupcoordnum() - : b_anisotropic(false) -{ - function_type = "groupcoordnum"; - x.type(colvarvalue::type_scalar); -} - - void colvar::groupcoordnum::calc_value() { // create fake atoms to hold the com coordinates diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 7cb4c30c31..c3c1ee55d5 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -19,6 +19,8 @@ colvar::distance::distance(std::string const &conf) : cvc(conf) { function_type = "distance"; + init_as_distance(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); @@ -27,8 +29,6 @@ colvar::distance::distance(std::string const &conf) group2 = parse_group(conf, "group2"); init_total_force_params(conf); - - x.type(colvarvalue::type_scalar); } @@ -36,10 +36,11 @@ colvar::distance::distance() : cvc() { function_type = "distance"; + init_as_distance(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); - x.type(colvarvalue::type_scalar); } @@ -177,7 +178,7 @@ colvar::distance_z::distance_z(std::string const &conf) // TODO detect PBC from MD engine (in simple cases) // and then update period in real time if (period != 0.0) - b_periodic = true; + enable(f_cvc_periodic); if ((wrap_center != 0.0) && (period == 0.0)) { cvm::error("Error: wrapAround was defined in a distanceZ component," @@ -317,7 +318,7 @@ cvm::real colvar::distance_z::dist2(colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; - if (b_periodic) { + if (is_enabled(f_cvc_periodic)) { cvm::real shift = cvm::floor(diff/period + 0.5); diff -= shift * period; } @@ -329,7 +330,7 @@ colvarvalue colvar::distance_z::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; - if (b_periodic) { + if (is_enabled(f_cvc_periodic)) { cvm::real shift = cvm::floor(diff/period + 0.5); diff -= shift * period; } @@ -341,7 +342,7 @@ colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1, colvarvalue const &x2) const { cvm::real diff = x1.real_value - x2.real_value; - if (b_periodic) { + if (is_enabled(f_cvc_periodic)) { cvm::real shift = cvm::floor(diff/period + 0.5); diff -= shift * period; } @@ -351,15 +352,13 @@ colvarvalue colvar::distance_z::dist2_rgrad(colvarvalue const &x1, void colvar::distance_z::wrap(colvarvalue &x_unwrapped) const { - if (!b_periodic) { + if (!is_enabled(f_cvc_periodic)) { // don't wrap if the period has not been set return; } - cvm::real shift = cvm::floor((x_unwrapped.real_value - wrap_center) / period + 0.5); x_unwrapped.real_value -= shift * period; - return; } @@ -368,10 +367,11 @@ colvar::distance_xy::distance_xy(std::string const &conf) : distance_z(conf) { function_type = "distance_xy"; + init_as_distance(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); - x.type(colvarvalue::type_scalar); } @@ -379,10 +379,11 @@ colvar::distance_xy::distance_xy() : distance_z() { function_type = "distance_xy"; + init_as_distance(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); enable(f_cvc_com_based); - x.type(colvarvalue::type_scalar); } @@ -557,6 +558,7 @@ colvar::distance_inv::distance_inv(std::string const &conf) : cvc(conf) { function_type = "distance_inv"; + init_as_distance(); group1 = parse_group(conf, "group1"); group2 = parse_group(conf, "group2"); @@ -585,16 +587,6 @@ colvar::distance_inv::distance_inv(std::string const &conf) "for distanceInv, because its value and gradients are computed " "simultaneously.\n"); } - - x.type(colvarvalue::type_scalar); -} - - -colvar::distance_inv::distance_inv() -{ - function_type = "distance_inv"; - exponent = 6; - x.type(colvarvalue::type_scalar); } @@ -807,6 +799,8 @@ colvar::gyration::gyration(std::string const &conf) : cvc(conf) { function_type = "gyration"; + init_as_distance(); + provide(f_cvc_inv_gradient); provide(f_cvc_Jacobian); atoms = parse_group(conf, "atoms"); @@ -818,17 +812,6 @@ colvar::gyration::gyration(std::string const &conf) 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)); } - - x.type(colvarvalue::type_scalar); -} - - -colvar::gyration::gyration() -{ - function_type = "gyration"; - provide(f_cvc_inv_gradient); - provide(f_cvc_Jacobian); - x.type(colvarvalue::type_scalar); } @@ -885,14 +868,7 @@ colvar::inertia::inertia(std::string const &conf) : gyration(conf) { function_type = "inertia"; - x.type(colvarvalue::type_scalar); -} - - -colvar::inertia::inertia() -{ - function_type = "inertia"; - x.type(colvarvalue::type_scalar); + init_as_distance(); } @@ -928,9 +904,10 @@ colvar::inertia_z::inertia_z(std::string const &conf) : inertia(conf) { function_type = "inertia_z"; + init_as_distance(); if (get_keyval(conf, "axis", axis, cvm::rvector(0.0, 0.0, 1.0))) { if (axis.norm2() == 0.0) { - cvm::error("Axis vector is zero!"); + cvm::error("Axis vector is zero!", INPUT_ERROR); return; } if (axis.norm2() != 1.0) { @@ -938,14 +915,6 @@ colvar::inertia_z::inertia_z(std::string const &conf) cvm::log("The normalized axis is: "+cvm::to_str(axis)+".\n"); } } - x.type(colvarvalue::type_scalar); -} - - -colvar::inertia_z::inertia_z() -{ - function_type = "inertia_z"; - x.type(colvarvalue::type_scalar); } @@ -982,9 +951,10 @@ simple_scalar_dist_functions(inertia) colvar::rmsd::rmsd(std::string const &conf) : cvc(conf) { - provide(f_cvc_inv_gradient); function_type = "rmsd"; - x.type(colvarvalue::type_scalar); + init_as_distance(); + + provide(f_cvc_inv_gradient); atoms = parse_group(conf, "atoms"); diff --git a/lib/colvars/colvarcomp_gpath.cpp b/lib/colvars/colvarcomp_gpath.cpp index 530aaf3284..d563403628 100644 --- a/lib/colvars/colvarcomp_gpath.cpp +++ b/lib/colvars/colvarcomp_gpath.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -119,9 +119,10 @@ colvar::CartesianBasedPath::~CartesianBasedPath() { (*it_comp_atoms) = nullptr; } } + atom_groups.clear(); } -void colvar::CartesianBasedPath::computeReferenceDistance(std::vector& result) { +void colvar::CartesianBasedPath::computeDistanceToReferenceFrames(std::vector& result) { for (size_t i_frame = 0; i_frame < reference_frames.size(); ++i_frame) { cvm::real frame_rmsd = 0.0; for (size_t i_atom = 0; i_atom < atoms->size(); ++i_atom) { @@ -147,13 +148,16 @@ colvar::gspath::gspath(std::string const &conf): CartesianBasedPath(conf) { } else { cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n")); } + if (total_reference_frames < 2) { + cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspath requires at least 2 frames.\n"); + } GeometricPathCV::GeometricPathBase::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame); cvm::log(std::string("Geometric pathCV(s) is initialized.\n")); cvm::log(std::string("Geometric pathCV(s) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n")); } -void colvar::gspath::updateReferenceDistances() { - computeReferenceDistance(frame_distances); +void colvar::gspath::updateDistanceToReferenceFrames() { + computeDistanceToReferenceFrames(frame_distances); } void colvar::gspath::prepareVectors() { @@ -290,14 +294,17 @@ colvar::gzpath::gzpath(std::string const &conf): CartesianBasedPath(conf) { if (b_use_z_square == true) { cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n")); } + if (total_reference_frames < 2) { + cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gzpath requires at least 2 frames.\n"); + } GeometricPathCV::GeometricPathBase::initialize(atoms->size(), cvm::atom_pos(), total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square); // Logging cvm::log(std::string("Geometric pathCV(z) is initialized.\n")); cvm::log(std::string("Geometric pathCV(z) loaded ") + cvm::to_str(reference_frames.size()) + std::string(" frames.\n")); } -void colvar::gzpath::updateReferenceDistances() { - computeReferenceDistance(frame_distances); +void colvar::gzpath::updateDistanceToReferenceFrames() { + computeDistanceToReferenceFrames(frame_distances); } void colvar::gzpath::prepareVectors() { @@ -538,6 +545,7 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { std::vector fields; split_string(line, token, fields); size_t num_value_required = 0; + cvm::log(std::string("Reading reference frame ") + cvm::to_str(total_reference_frames + 1) + std::string("\n")); for (size_t i_cv = 0; i_cv < tmp_cv.size(); ++i_cv) { const size_t value_size = tmp_cv[i_cv].size(); num_value_required += value_size; @@ -545,10 +553,9 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { if (num_value_required <= fields.size()) { size_t start_index = num_value_required - value_size; for (size_t i = start_index; i < num_value_required; ++i) { - tmp_cv[i_cv][i] = std::atof(fields[i].c_str()); - cvm::log(fields[i] + std::string(" ")); + tmp_cv[i_cv][i - start_index] = std::atof(fields[i].c_str()); + cvm::log(cvm::to_str(tmp_cv[i_cv][i - start_index])); } - cvm::log(std::string("\n")); } else { cvm::error("Error: incorrect format of path file.\n"); } @@ -558,6 +565,9 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { ++total_reference_frames; } } + if (total_reference_frames <= 1) { + cvm::error("Error: there is only 1 or 0 reference frame, which doesn't constitute a path.\n"); + } x.type(colvarvalue::type_scalar); use_explicit_gradients = true; for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { @@ -570,7 +580,7 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { } } -void colvar::CVBasedPath::computeReferenceDistance(std::vector& result) { +void colvar::CVBasedPath::computeDistanceToReferenceFrames(std::vector& result) { for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { cv[i_cv]->calc_value(); } @@ -593,6 +603,20 @@ void colvar::CVBasedPath::computeReferenceDistance(std::vector& resul } } +void colvar::CVBasedPath::computeDistanceBetweenReferenceFrames(std::vector& result) const { + if (ref_cv.size() < 2) return; + for (size_t i_frame = 1; i_frame < ref_cv.size(); ++i_frame) { + cvm::real dist_ij = 0.0; + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + colvarvalue ref_cv_value(ref_cv[i_frame][i_cv]); + colvarvalue prev_ref_cv_value(ref_cv[i_frame-1][i_cv]); + dist_ij += cv[i_cv]->dist2(ref_cv_value, prev_ref_cv_value); + } + dist_ij = cvm::sqrt(dist_ij); + result[i_frame-1] = dist_ij; + } +} + cvm::real colvar::CVBasedPath::getPolynomialFactorOfCVGradient(size_t i_cv) const { cvm::real factor_polynomial = 1.0; if (cv[i_cv]->value().type() == colvarvalue::type_scalar) { @@ -607,6 +631,7 @@ colvar::CVBasedPath::~CVBasedPath() { for (auto it = cv.begin(); it != cv.end(); ++it) { delete (*it); } + atom_groups.clear(); } colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) { @@ -625,6 +650,9 @@ colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) { } else { cvm::log(std::string("Geometric path s(σ) will use the neighbouring frame to compute s_(m+1)\n")); } + if (total_reference_frames < 2) { + cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gspathCV requires at least 2 frames.\n"); + } GeometricPathCV::GeometricPathBase::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame); x.type(colvarvalue::type_scalar); use_explicit_gradients = true; @@ -641,8 +669,8 @@ colvar::gspathCV::gspathCV(std::string const &conf): CVBasedPath(conf) { colvar::gspathCV::~gspathCV() {} -void colvar::gspathCV::updateReferenceDistances() { - computeReferenceDistance(frame_distances); +void colvar::gspathCV::updateDistanceToReferenceFrames() { + computeDistanceToReferenceFrames(frame_distances); } void colvar::gspathCV::prepareVectors() { @@ -766,6 +794,9 @@ colvar::gzpathCV::gzpathCV(std::string const &conf): CVBasedPath(conf) { if (b_use_z_square == true) { cvm::log(std::string("Geometric path z(σ) will use the square of distance from current frame to path compute z\n")); } + if (total_reference_frames < 2) { + cvm::error("Error: you have specified " + cvm::to_str(total_reference_frames) + " reference frames, but gzpathCV requires at least 2 frames.\n"); + } GeometricPathCV::GeometricPathBase::initialize(cv.size(), ref_cv[0], total_reference_frames, use_second_closest_frame, use_third_closest_frame, b_use_z_square); x.type(colvarvalue::type_scalar); use_explicit_gradients = true; @@ -783,8 +814,8 @@ colvar::gzpathCV::gzpathCV(std::string const &conf): CVBasedPath(conf) { colvar::gzpathCV::~gzpathCV() { } -void colvar::gzpathCV::updateReferenceDistances() { - computeReferenceDistance(frame_distances); +void colvar::gzpathCV::updateDistanceToReferenceFrames() { + computeDistanceToReferenceFrames(frame_distances); } void colvar::gzpathCV::prepareVectors() { diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index aa61cdf1dc..e0779b602e 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -30,6 +30,8 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); + colvarproxy *proxy = cvm::main()->proxy; + std::string segment_id; get_keyval(conf, "psfSegID", segment_id, std::string("MAIN")); @@ -91,7 +93,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) { cvm::real r0; size_t en, ed; - get_keyval(conf, "hBondCutoff", r0, (3.3 * cvm::unit_angstrom())); + get_keyval(conf, "hBondCutoff", r0, (3.3 * proxy->angstrom_value)); get_keyval(conf, "hBondExpNumer", en, 6); get_keyval(conf, "hBondExpDenom", ed, 8); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 3c8d0b4af6..c7a414b006 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -164,8 +164,8 @@ colvar::orientation_angle::orientation_angle(std::string const &conf) : orientation() { function_type = "orientation_angle"; + init_as_angle(); enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); init(conf); } @@ -176,15 +176,6 @@ int colvar::orientation_angle::init(std::string const &conf) } -colvar::orientation_angle::orientation_angle() - : orientation() -{ - function_type = "orientation_angle"; - enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); -} - - void colvar::orientation_angle::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -231,6 +222,7 @@ colvar::orientation_proj::orientation_proj(std::string const &conf) function_type = "orientation_proj"; enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); + init_scalar_boundaries(0.0, 1.0); init(conf); } @@ -241,15 +233,6 @@ int colvar::orientation_proj::init(std::string const &conf) } -colvar::orientation_proj::orientation_proj() - : orientation() -{ - function_type = "orientation_proj"; - enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); -} - - void colvar::orientation_proj::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -287,6 +270,7 @@ colvar::tilt::tilt(std::string const &conf) function_type = "tilt"; enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); + init_scalar_boundaries(-1.0, 1.0); init(conf); } @@ -307,15 +291,6 @@ int colvar::tilt::init(std::string const &conf) } -colvar::tilt::tilt() - : orientation() -{ - function_type = "tilt"; - enable(f_cvc_explicit_gradient); - x.type(colvarvalue::type_scalar); -} - - void colvar::tilt::calc_value() { atoms_cog = atoms->center_of_geometry(); @@ -358,7 +333,7 @@ colvar::spin_angle::spin_angle(std::string const &conf) { function_type = "spin_angle"; period = 360.0; - b_periodic = true; + enable(f_cvc_periodic); enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); init(conf); @@ -386,7 +361,7 @@ colvar::spin_angle::spin_angle() { function_type = "spin_angle"; period = 360.0; - b_periodic = true; + enable(f_cvc_periodic); enable(f_cvc_explicit_gradient); x.type(colvarvalue::type_scalar); } diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp index 276f2b39e7..8425982be4 100644 --- a/lib/colvars/colvardeps.cpp +++ b/lib/colvars/colvardeps.cpp @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvardeps.h b/lib/colvars/colvardeps.h index 198a24f330..47a91ab1ac 100644 --- a/lib/colvars/colvardeps.h +++ b/lib/colvars/colvardeps.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -219,7 +219,7 @@ public: /// Implements possible actions to be carried out /// when a given feature is enabled /// Base function does nothing, can be overloaded - virtual void do_feature_side_effects(int id) {} + virtual void do_feature_side_effects(int /* id */) {} // NOTE that all feature enums should start with f_*_active enum features_biases { @@ -296,6 +296,10 @@ public: f_cv_lower_boundary, /// \brief An upper boundary is defined f_cv_upper_boundary, + /// \brief The lower boundary is not defined from user's choice + f_cv_hard_lower_boundary, + /// \brief The upper boundary is not defined from user's choice + f_cv_hard_upper_boundary, /// \brief Provide a discretization of the values of the colvar to /// be used by the biases or in analysis (needs lower and upper /// boundary) @@ -310,6 +314,8 @@ public: f_cv_custom_function, /// \brief Colvar is periodic f_cv_periodic, + /// \brief The colvar has only one component + f_cv_single_cvc, /// \brief is scalar f_cv_scalar, f_cv_linear, @@ -321,20 +327,40 @@ public: }; enum features_cvc { + /// Computation of this CVC is enabled f_cvc_active, + /// This CVC computes a scalar value f_cvc_scalar, + /// Values of this CVC lie in a periodic interval + f_cvc_periodic, + /// This CVC provides a default value for the colvar's width + f_cvc_width, + /// This CVC provides a default value for the colvar's lower boundary + f_cvc_lower_boundary, + /// This CVC provides a default value for the colvar's upper boundary + f_cvc_upper_boundary, + /// CVC calculates atom gradients f_cvc_gradient, - /// \brief CVC calculates and stores explicit atom gradients + /// CVC calculates and stores explicit atom gradients f_cvc_explicit_gradient, + /// CVC calculates and stores inverse atom gradients (used for total force) f_cvc_inv_gradient, - /// \brief If enabled, calc_gradients() will call debug_gradients() for every group needed - f_cvc_debug_gradient, + /// CVC calculates the Jacobian term of the total-force expression f_cvc_Jacobian, - f_cvc_pbc_minimum_image, + /// The total force for this CVC will be computed from one group only f_cvc_one_site_total_force, + /// calc_gradients() will call debug_gradients() for every group needed + f_cvc_debug_gradient, + /// With PBCs, minimum-image convention will be used for distances + /// (does not affect the periodicity of CVC values, e.g. angles) + f_cvc_pbc_minimum_image, + /// This CVC is a function of centers of mass f_cvc_com_based, + /// This CVC can be computed in parallel f_cvc_scalable, + /// Centers-of-mass used in this CVC can be computed in parallel f_cvc_scalable_com, + /// Number of CVC features f_cvc_ntot }; diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index bc5ac3c57c..7c2fe3b33c 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 7792b7d0df..9d07ba0972 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -209,7 +209,8 @@ public: /// \brief "Almost copy-constructor": only copies configuration /// parameters from another grid, but doesn't reallocate stuff; /// setup() must be called after that; - colvar_grid(colvar_grid const &g) : nd(g.nd), + colvar_grid(colvar_grid const &g) : colvarparse(), + nd(g.nd), nx(g.nx), mult(g.mult), data(), @@ -284,8 +285,8 @@ public: } widths.push_back(cv[i]->width); - hard_lower_boundaries.push_back(cv[i]->hard_lower_boundary); - hard_upper_boundaries.push_back(cv[i]->hard_upper_boundary); + hard_lower_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_lower_boundary)); + hard_upper_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_upper_boundary)); periodic.push_back(cv[i]->periodic_boundaries()); // By default, get reported colvar value (for extended Lagrangian colvars) @@ -415,6 +416,20 @@ public: return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] ); } + /// \brief Report the fraction of bin beyond current_bin_scalar() + inline cvm::real current_bin_scalar_fraction(int const i) const + { + return value_to_bin_scalar_fraction(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i); + } + + /// \brief Use the lower boundary and the width to report the fraction of bin + /// beyond value_to_bin_scalar() that the provided value is in + inline cvm::real value_to_bin_scalar_fraction(colvarvalue const &value, const int i) const + { + cvm::real x = (value.real_value - lower_boundaries[i].real_value) / widths[i]; + return x - cvm::floor(x); + } + /// \brief Use the lower boundary and the width to report which bin /// the provided value is in and assign first or last bin if out of boundaries inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index baffc25c28..ecf0ae800d 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -32,6 +32,8 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) cv_traj_os = NULL; + xyz_reader_use_count = 0; + if (proxy == NULL) { proxy = proxy_in; // Pointer to the proxy object parse = new colvarparse(); // Parsing object for global options @@ -48,13 +50,21 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) cvm::log("Initializing the collective variables module, version "+ cvm::to_str(COLVARS_VERSION)+".\n"); cvm::log("Please cite Fiorin et al, Mol Phys 2013:\n " - "http://dx.doi.org/10.1080/00268976.2013.813594\n" + "https://dx.doi.org/10.1080/00268976.2013.813594\n" "in any publication based on this calculation.\n"); if (proxy->smp_enabled() == COLVARS_OK) { - cvm::log("SMP parallelism is available.\n"); + cvm::log("SMP parallelism is enabled; if needed, use \"smp off\" to override this.\n"); } +#if (__cplusplus >= 201103L) + cvm::log("This version was built with the C++11 standard or higher."); +#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" + "https://colvars.github.io/README-c++11.html"); +#endif + // set initial default values // "it_restart" will be set by the input state file, if any; @@ -260,6 +270,16 @@ int colvarmodule::parse_global_params(std::string const &conf) // TODO document and then echo this keyword parse->get_keyval(conf, "logLevel", log_level_, log_level_, colvarparse::parse_silent); + { + std::string units; + if (parse->get_keyval(conf, "units", units)) { + units = colvarparse::to_lower_cppstr(units); + int error_code = proxy->set_unit_system(units, (colvars.size() != 0)); + if (error_code != COLVARS_OK) { + return error_code; + } + } + } { std::string index_file_name; @@ -1160,8 +1180,7 @@ int colvarmodule::reset() } colvars.clear(); - index_groups.clear(); - index_group_names.clear(); + reset_index_groups(); proxy->reset(); @@ -1266,23 +1285,45 @@ std::istream & colvarmodule::read_restart(std::istream &is) // read global restart information std::string restart_conf; if (is >> colvarparse::read_block("configuration", restart_conf)) { + parse->get_keyval(restart_conf, "step", it_restart, static_cast(0), colvarparse::parse_restart); - it = it_restart; + it = it_restart; + std::string restart_version; + int restart_version_int = 0; parse->get_keyval(restart_conf, "version", restart_version, std::string(""), colvarparse::parse_restart); - if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) { - cvm::log("This state file was generated with version "+restart_version+"\n"); + if (restart_version.size()) { + if (restart_version != std::string(COLVARS_VERSION)) { + cvm::log("This state file was generated with version "+ + restart_version+"\n"); + } + restart_version_int = + proxy->get_version_from_string(restart_version.c_str()); } - if ((restart_version.size() == 0) || (restart_version.compare(std::string(COLVARS_VERSION)) < 0)) { + + if (restart_version_int < 20160810) { // check for total force change if (proxy->total_forces_enabled()) { warn_total_forces = true; } } + + std::string units_restart; + if (parse->get_keyval(restart_conf, "units", + units_restart, std::string(""), + colvarparse::parse_restart)) { + units_restart = colvarparse::to_lower_cppstr(units_restart); + if ((proxy->units.size() > 0) && (units_restart != proxy->units)) { + cvm::error("Error: the state file has units \""+units_restart+ + "\", but the current unit system is \""+proxy->units+ + "\".\n", INPUT_ERROR); + } + } + } is.clear(); parse->clear_keyword_registry(); @@ -1314,69 +1355,19 @@ std::istream & colvarmodule::read_restart(std::istream &is) if (warn_total_forces) { cvm::log(cvm::line_marker); - cvm::log("WARNING:\n"); - std::string const warning("### CHANGES IN THE DEFINITION OF SYSTEM FORCES (NOW TOTAL FORCES)\n\ -\n\ -Starting from the version 2016-08-10 of the Colvars module, \n\ -the role of system forces has been replaced by total forces.\n\ -\n\ -These include *all* forces acting on a collective variable, whether they\n\ -come from the force field potential or from external terms\n\ -(e.g. restraints), including forces applied by Colvars itself.\n\ -\n\ -In NAMD, forces applied by Colvars, IMD, SMD, TMD, symmetry\n\ -restraints and tclForces are now all counted in the total force.\n\ -\n\ -In LAMMPS, forces applied by Colvars itself are now counted in the total\n\ -force (all forces from other fixes were being counted already).\n\ -\n\ -\n\ -### WHEN IS THIS CHANGE RELEVANT\n\ -\n\ -This change affects results *only* when (1) outputSystemForce is\n\ -requested or (2) the ABF bias is used. All other usage cases are\n\ -*unaffected* (colvar restraints, metadynamics, etc).\n\ -\n\ -When system forces are reported (flag: outputSystemForce), their values\n\ -in the output may change, but the physical trajectory is never affected.\n\ -The physical results of ABF calculations may be affected in some cases.\n\ -\n\ -\n\ -### CHANGES TO ABF CALCULATIONS\n\ -\n\ -Compared to previous Colvars versions, the ABF method will now attempt\n\ -to cancel external forces (for example, boundary walls) and it may be\n\ -not possible to resume through a state file a simulation that was\n\ -performed with a previous version.\n\ -\n\ -There are three possible scenarios:\n\ -\n\ -1. No external forces are applied to the atoms used by ABF: results are\n\ -unchanged.\n\ -\n\ -2. Some of the atoms used by ABF experience external forces, but these\n\ -forces are not applied directly to the variables used by ABF\n\ -(e.g. another colvar that uses the same atoms, tclForces, etc): in this\n\ -case, we recommend beginning a new simulation.\n\ -\n\ -3. External forces are applied to one or more of the colvars used by\n\ -ABF, but no other forces are applied to their atoms: you may use the\n\ -subtractAppliedForce keyword inside the corresponding colvars to\n\ -continue the previous simulation.\n\n"); - cvm::log(warning); - cvm::log(cvm::line_marker); - + cvm::log("WARNING: The definition of system forces has changed. Please see:\n"); + cvm::log(" https://colvars.github.io/README-totalforce.html\n"); // update this ahead of time in this special case output_prefix() = proxy->input_prefix(); cvm::log("All output files will now be saved with the prefix \""+output_prefix()+".tmp.*\".\n"); - cvm::log(cvm::line_marker); cvm::log("Please review the important warning above. After that, you may rename:\n\ \""+output_prefix()+".tmp.colvars.state\"\n\ to:\n\ -\""+ proxy->input_prefix()+".colvars.state\"\n"); +\""+proxy->input_prefix()+".colvars.state\"\n\ +and load it to continue this simulation.\n"); output_prefix() = output_prefix()+".tmp"; write_restart_file(output_prefix()+".colvars.state"); - cvm::error("Exiting with error until issue is addressed.\n", FATAL_ERROR); + cvm::error("Exiting with error until issue is addressed.\n", INPUT_ERROR); } return is; @@ -1496,8 +1487,11 @@ std::ostream & colvarmodule::write_restart(std::ostream &os) << " step " << std::setw(it_width) << it << "\n" << " dt " << dt() << "\n" - << " version " << std::string(COLVARS_VERSION) << "\n" - << "}\n\n"; + << " version " << std::string(COLVARS_VERSION) << "\n"; + if (proxy->units.size() > 0) { + os << " units " << proxy->units << "\n"; + } + os << "}\n\n"; int error_code = COLVARS_OK; @@ -1718,32 +1712,61 @@ int cvm::read_index_file(char const *filename) while (is.good()) { char open, close; std::string group_name; + int index_of_group = -1; if ( (is >> open) && (open == '[') && (is >> group_name) && (is >> close) && (close == ']') ) { - for (std::list::iterator names_i = index_group_names.begin(); - names_i != index_group_names.end(); - names_i++) { - if (*names_i == group_name) { - cvm::error("Error: the group name \""+group_name+ - "\" appears in multiple index files.\n", - FILE_ERROR); + size_t i = 0; + for ( ; i < index_group_names.size(); i++) { + if (index_group_names[i] == group_name) { + // Found a group with the same name + index_of_group = i; } } - index_group_names.push_back(group_name); - index_groups.push_back(std::vector()); + if (index_of_group < 0) { + index_group_names.push_back(group_name); + index_groups.push_back(NULL); + index_of_group = index_groups.size()-1; + } } else { - cvm::error("Error: in parsing index file \""+ - std::string(filename)+"\".\n", - INPUT_ERROR); + return cvm::error("Error: in parsing index file \""+ + std::string(filename)+"\".\n", + INPUT_ERROR); } + std::vector *old_index_group = index_groups[index_of_group]; + std::vector *new_index_group = new std::vector(); + int atom_number = 1; size_t pos = is.tellg(); while ( (is >> atom_number) && (atom_number > 0) ) { - (index_groups.back()).push_back(atom_number); + new_index_group->push_back(atom_number); pos = is.tellg(); } + + if (old_index_group != NULL) { + bool equal = false; + if (new_index_group->size() == old_index_group->size()) { + if (std::equal(new_index_group->begin(), new_index_group->end(), + old_index_group->begin())) { + equal = true; + } + } + if (! equal) { + new_index_group->clear(); + delete new_index_group; + new_index_group = NULL; + return cvm::error("Error: the index group \""+group_name+ + "\" was redefined.\n", INPUT_ERROR); + } else { + old_index_group->clear(); + delete old_index_group; + old_index_group = NULL; + } + } + + index_groups[index_of_group] = new_index_group; + is.clear(); is.seekg(pos, std::ios::beg); std::string delim; @@ -1756,14 +1779,27 @@ int cvm::read_index_file(char const *filename) } } - cvm::log("The following index groups were read from the index file \""+ - std::string(filename)+"\":\n"); - std::list::iterator names_i = index_group_names.begin(); - std::list >::iterator lists_i = index_groups.begin(); - for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) { - cvm::log(" "+(*names_i)+" ("+cvm::to_str(lists_i->size())+" atoms).\n"); + cvm::log("The following index groups are currently defined:\n"); + size_t i = 0; + for ( ; i < index_group_names.size(); i++) { + cvm::log(" "+(index_group_names[i])+" ("+ + cvm::to_str((index_groups[i])->size())+" atoms)\n"); } - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + + return COLVARS_OK; +} + + +int colvarmodule::reset_index_groups() +{ + size_t i = 0; + for ( ; i < index_groups.size(); i++) { + delete index_groups[i]; + index_groups[i] = NULL; + } + index_group_names.clear(); + index_groups.clear(); + return COLVARS_OK; } @@ -1799,7 +1835,7 @@ int cvm::load_coords(char const *file_name, "for XYZ coordinate files.\n", INPUT_ERROR); } // For XYZ files, use internal parser - error_code |= cvm::load_coords_xyz(file_name, &sorted_pos, atoms); + error_code |= cvm::main()->load_coords_xyz(file_name, &sorted_pos, atoms); } else { // Otherwise, call proxy function for PDB error_code |= proxy->load_coords(file_name, @@ -1824,11 +1860,18 @@ int cvm::load_coords_xyz(char const *filename, unsigned int natoms; char symbol[256]; std::string line; + cvm::real x = 0.0, y = 0.0, z = 0.0; if ( ! (xyz_is >> natoms) ) { cvm::error("Error: cannot parse XYZ file " + std::string(filename) + ".\n", 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."); + } + // skip comment line cvm::getline(xyz_is, line); cvm::getline(xyz_is, line); @@ -1845,12 +1888,19 @@ int cvm::load_coords_xyz(char const *filename, next++; } xyz_is >> symbol; - xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; + 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); } } else { // Use all positions for ( ; pos_i != pos->end() ; pos_i++) { xyz_is >> symbol; - xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; + 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); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -1860,11 +1910,6 @@ int cvm::load_coords_xyz(char const *filename, // Wrappers to proxy functions: these may go in the future -cvm::real cvm::unit_angstrom() -{ - return proxy->unit_angstrom(); -} - cvm::real cvm::boltzmann() { @@ -1903,39 +1948,6 @@ cvm::real cvm::rand_gaussian(void) } - -bool cvm::replica_enabled() -{ - return proxy->replica_enabled(); -} - -int cvm::replica_index() -{ - return proxy->replica_index(); -} - -int cvm::replica_num() -{ - return proxy->replica_num(); -} - -void cvm::replica_comm_barrier() -{ - return proxy->replica_comm_barrier(); -} - -int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) -{ - return proxy->replica_comm_recv(msg_data,buf_len,src_rep); -} - -int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) -{ - return proxy->replica_comm_send(msg_data,msg_len,dest_rep); -} - - - template std::string _to_str(T const &x, size_t width, size_t prec) { diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 30620a2527..d3bc5e1b80 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -608,10 +608,6 @@ public: // proxy functions - /// \brief Value of the unit for atomic coordinates with respect to - /// angstroms (used by some variables for hard-coded default values) - static real unit_angstrom(); - /// \brief Boltmann constant static real boltzmann(); @@ -678,29 +674,23 @@ public: return 5; } - - // Replica exchange commands. - static bool replica_enabled(); - static int replica_index(); - static int replica_num(); - static void replica_comm_barrier(); - static int replica_comm_recv(char* msg_data, int buf_len, int src_rep); - static int replica_comm_send(char* msg_data, int msg_len, int dest_rep); - /// \brief Get the distance between two atomic positions with pbcs handled /// correctly static rvector position_distance(atom_pos const &pos1, atom_pos const &pos2); - /// \brief Names of groups from a Gromacs .ndx file to be read at startup - std::list index_group_names; + /// \brief Names of groups from one or more Gromacs .ndx files + std::vector index_group_names; - /// \brief Groups from a Gromacs .ndx file read at startup - std::list > index_groups; + /// \brief Groups from one or more Gromacs .ndx files + std::vector *> index_groups; /// \brief Read a Gromacs .ndx file int read_index_file(char const *filename); + /// Clear the index groups loaded so far + int reset_index_groups(); + /// \brief Select atom IDs from a file (usually PDB) \param filename name of /// the file \param atoms array into which atoms read from "filename" will be /// appended \param pdb_field (optional) if the file is a PDB and this @@ -726,11 +716,10 @@ public: std::string const &pdb_field, double pdb_field_value = 0.0); - /// \brief Load the coordinates for a group of atoms from an - /// XYZ file - static int load_coords_xyz(char const *filename, - std::vector *pos, - atom_group *atoms); + /// Load coordinates into an atom group from an XYZ file (assumes Angstroms) + int load_coords_xyz(char const *filename, + std::vector *pos, + atom_group *atoms); /// Frequency for collective variables trajectory output static size_t cv_traj_freq; @@ -771,6 +760,9 @@ private: /// Thread-specific depth std::vector depth_v; + /// Track how many times the XYZ reader has been used + int xyz_reader_use_count; + public: /// Get the current object depth in the hierarchy diff --git a/lib/colvars/colvarparams.cpp b/lib/colvars/colvarparams.cpp new file mode 100644 index 0000000000..e957a1841d --- /dev/null +++ b/lib/colvars/colvarparams.cpp @@ -0,0 +1,111 @@ +// -*- 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 +#include + +#include "colvarmodule.h" +#include "colvarvalue.h" +#include "colvarparams.h" + + + +colvarparams::colvarparams() +{} + + +colvarparams::~colvarparams() +{} + + +void colvarparams::register_param(std::string const ¶m_name, + void *param_ptr) +{ + param_map[param_name] = param_ptr; +} + + +void colvarparams::register_param_grad(std::string const ¶m_name, + colvarvalue *param_grad_ptr) +{ + param_grad_map[param_name] = param_grad_ptr; +} + + +int colvarparams::param_exists(std::string const ¶m_name) +{ + if (param_map.count(param_name) > 0) { + return COLVARS_OK; + } + return INPUT_ERROR; +} + + +std::vector colvarparams::get_param_names() +{ + std::vector result; + for (std::map::const_iterator elem = + param_map.begin(); elem != param_map.end(); elem++) { + result.push_back(elem->first); + } + return result; +} + + +std::vector colvarparams::get_param_grad_names() +{ + std::vector result; + for (std::map::const_iterator elem = + param_grad_map.begin(); elem != param_grad_map.end(); elem++) { + result.push_back(elem->first); + } + return result; +} + + +void const *colvarparams::get_param_ptr(std::string const ¶m_name) +{ + if (param_map.count(param_name) > 0) { + return param_map[param_name]; + } + cvm::error("Error: parameter \""+param_name+"\" not found.\n", INPUT_ERROR); + return NULL; +} + + +void const *colvarparams::get_param_grad_ptr(std::string const ¶m_name) +{ + if (param_grad_map.count(param_name) > 0) { + return param_grad_map[param_name]; + } + cvm::error("Error: gradient of parameter \""+param_name+"\" not found.\n", + INPUT_ERROR); + return NULL; +} + + +cvm::real colvarparams::get_param(std::string const ¶m_name) +{ + cvm::real const *ptr = + reinterpret_cast(get_param_ptr(param_name)); + return ptr != NULL ? *ptr : 0.0; +} + + +int colvarparams::set_param(std::string const ¶m_name, + void const *new_value) +{ + if (param_map.count(param_name) > 0) { + return cvm::error("Error: parameter \""+param_name+"\" cannot be " + "modified.\n", COLVARS_NOT_IMPLEMENTED); + } + return cvm::error("Error: parameter \""+param_name+"\" not found.\n", + INPUT_ERROR); +} diff --git a/lib/colvars/colvarparams.h b/lib/colvars/colvarparams.h new file mode 100644 index 0000000000..1d00300a8b --- /dev/null +++ b/lib/colvars/colvarparams.h @@ -0,0 +1,66 @@ +// -*- 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 COLVARPARAMS_H +#define COLVARPARAMS_H + +#include + +/// \file colvarparams.h Functions to handle scalar parameters used in objects + + +class colvarparams { + + public: + + /// Whether the parameter param_name exists + int param_exists(std::string const ¶m_name); + + /// Get a copy of the names of registered parameters + virtual std::vector get_param_names(); + + /// Get a copy of the names of registered parameter gradients + virtual std::vector get_param_grad_names(); + + /// Pointer to the parameter param_name + virtual void const *get_param_ptr(std::string const ¶m_name); + + /// Pointer to the gradient of parameter param_name + virtual void const *get_param_grad_ptr(std::string const ¶m_name); + + /// Value of the parameter param_name (must be a scalar) + virtual cvm::real get_param(std::string const ¶m_name); + + /// Set the named parameter to the given value + virtual int set_param(std::string const ¶m_name, void const *new_value); + + protected: + + /// Default constructor + colvarparams(); + + /// Default destructor + virtual ~colvarparams(); + + /// Pointers to relevant parameters that may be accessed by other objects + std::map param_map; + + /// Derivatives of the object with respect to internal parameters + std::map param_grad_map; + + /// Register the given parameter + void register_param(std::string const ¶m_name, void *param_ptr); + + /// Register the gradient of the given parameter + void register_param_grad(std::string const ¶m_name, + colvarvalue *param_grad_ptr); + +}; + +#endif diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 8041153494..77f03793a3 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -32,6 +32,41 @@ namespace { } +colvarparse::colvarparse() +{ + init(); +} + + +void colvarparse::init() +{ + config_string.clear(); + clear_keyword_registry(); +} + + +colvarparse::colvarparse(const std::string& conf) +{ + init(conf); +} + + +void colvarparse::init(std::string const &conf) +{ + if (! config_string.size()) { + init(); + config_string = conf; + } +} + + +colvarparse::~colvarparse() +{ + init(); +} + + + bool colvarparse::get_key_string_value(std::string const &conf, char const *key, std::string &data) { @@ -156,7 +191,7 @@ template<> int colvarparse::_get_keyval_scalar_value_(std::string const &key_str, std::string const &data, bool &value, - bool const &def_value) + bool const & /* def_value */) { if ( (data == std::string("on")) || (data == std::string("yes")) || @@ -176,8 +211,8 @@ int colvarparse::_get_keyval_scalar_value_(std::string const &key_str, template int colvarparse::_get_keyval_scalar_novalue_(std::string const &key_str, - TYPE &value, - Parse_Mode const &parse_mode) + TYPE & /* value */, + Parse_Mode const & /* parse_mode */) { return cvm::error("Error: improper or missing value " "for \""+key_str+"\".\n", INPUT_ERROR); diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index f4ead26601..f43f6527c8 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -12,10 +12,10 @@ #include #include -#include #include "colvarmodule.h" #include "colvarvalue.h" +#include "colvarparams.h" /// \file colvarparse.h Parsing functions for collective variables @@ -23,37 +23,24 @@ /// \brief Base class containing parsing functions; all objects which /// need to parse input inherit from this -class colvarparse { +class colvarparse : public colvarparams { public: /// Default constructor - inline colvarparse() - { - init(); - } + colvarparse(); /// Constructor that stores the object's config string - inline colvarparse(const std::string& conf) - { - init(conf); - } + colvarparse(const std::string& conf); /// Set the object ready to parse a new configuration string - inline void init() - { - config_string.clear(); - clear_keyword_registry(); - } + void init(); /// Set a new config string for this object - inline void init(std::string const &conf) - { - if (! config_string.size()) { - init(); - config_string = conf; - } - } + void init(std::string const &conf); + + /// Default destructor + virtual ~colvarparse(); /// Get the configuration string (includes comments) inline std::string const & get_config() const diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp index f077984556..63ae07a57c 100644 --- a/lib/colvars/colvarproxy.cpp +++ b/lib/colvars/colvarproxy.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -33,6 +33,7 @@ colvarproxy_system::colvarproxy_system() + : angstrom_value(0.) { reset_pbc_lattice(); } @@ -41,7 +42,7 @@ colvarproxy_system::colvarproxy_system() colvarproxy_system::~colvarproxy_system() {} -void colvarproxy_system::add_energy(cvm::real energy) {} +void colvarproxy_system::add_energy(cvm::real /* energy */) {} void colvarproxy_system::request_total_force(bool yesno) @@ -173,9 +174,9 @@ int colvarproxy_atoms::add_atom_slot(int atom_id) } -int colvarproxy_atoms::init_atom(cvm::residue_id const &residue, - std::string const &atom_name, - std::string const &segment_id) +int colvarproxy_atoms::init_atom(cvm::residue_id const & /* residue */, + std::string const & /* atom_name */, + std::string const & /* segment_id */) { cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n", COLVARS_NOT_IMPLEMENTED); @@ -204,9 +205,9 @@ void colvarproxy_atoms::clear_atom(int index) } -int colvarproxy_atoms::load_atoms(char const *filename, - cvm::atom_group &atoms, - std::string const &pdb_field, +int colvarproxy_atoms::load_atoms(char const * /* filename */, + cvm::atom_group & /* atoms */, + std::string const & /* pdb_field */, double) { return cvm::error("Error: loading atom identifiers from a file " @@ -215,10 +216,10 @@ int colvarproxy_atoms::load_atoms(char const *filename, } -int colvarproxy_atoms::load_coords(char const *filename, - std::vector &pos, - std::vector const &sorted_ids, - std::string const &pdb_field, +int colvarproxy_atoms::load_coords(char const * /* filename */, + std::vector & /* pos */, + std::vector const & /* sorted_ids */, + std::string const & /* pdb_field */, double) { return cvm::error("Error: loading atomic coordinates from a file " @@ -269,7 +270,7 @@ int colvarproxy_atom_groups::scalable_group_coms() } -int colvarproxy_atom_groups::init_atom_group(std::vector const &atoms_ids) +int colvarproxy_atom_groups::init_atom_group(std::vector const & /* atoms_ids */) { cvm::error("Error: initializing a group outside of the Colvars module " "is currently not supported.\n", @@ -456,51 +457,6 @@ int colvarproxy_smp::smp_unlock() -colvarproxy_replicas::colvarproxy_replicas() {} - - -colvarproxy_replicas::~colvarproxy_replicas() {} - - -bool colvarproxy_replicas::replica_enabled() -{ - return false; -} - - -int colvarproxy_replicas::replica_index() -{ - return 0; -} - - -int colvarproxy_replicas::replica_num() -{ - return 1; -} - - -void colvarproxy_replicas::replica_comm_barrier() {} - - -int colvarproxy_replicas::replica_comm_recv(char* msg_data, - int buf_len, - int src_rep) -{ - return COLVARS_NOT_IMPLEMENTED; -} - - -int colvarproxy_replicas::replica_comm_send(char* msg_data, - int msg_len, - int dest_rep) -{ - return COLVARS_NOT_IMPLEMENTED; -} - - - - colvarproxy_script::colvarproxy_script() { script = NULL; @@ -518,7 +474,7 @@ char const *colvarproxy_script::script_obj_to_str(unsigned char *obj) } -std::vector colvarproxy_script::script_obj_to_str_vector(unsigned char *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); @@ -532,19 +488,17 @@ int colvarproxy_script::run_force_callback() } -int colvarproxy_script::run_colvar_callback( - std::string const &name, - std::vector const &cvcs, - colvarvalue &value) +int colvarproxy_script::run_colvar_callback(std::string const & /* name */, + std::vector const & /* cvcs */, + colvarvalue & /* value */) { return COLVARS_NOT_IMPLEMENTED; } -int colvarproxy_script::run_colvar_gradient_callback( - std::string const &name, - std::vector const &cvcs, - std::vector > &gradient) +int colvarproxy_script::run_colvar_gradient_callback(std::string const & /* name */, + std::vector const & /* cvcs */, + std::vector > & /* gradient */) { return COLVARS_NOT_IMPLEMENTED; } @@ -599,9 +553,9 @@ int colvarproxy_tcl::tcl_run_force_callback() int colvarproxy_tcl::tcl_run_colvar_callback( - std::string const &name, - std::vector const &cvc_values, - colvarvalue &value) + std::string const & name, + std::vector const & cvc_values, + colvarvalue & value) { #if defined(COLVARS_TCL) @@ -636,9 +590,9 @@ int colvarproxy_tcl::tcl_run_colvar_callback( int colvarproxy_tcl::tcl_run_colvar_gradient_callback( - std::string const &name, - std::vector const &cvc_values, - std::vector > &gradient) + std::string const & name, + std::vector const & cvc_values, + std::vector > & gradient) { #if defined(COLVARS_TCL) @@ -776,7 +730,7 @@ int colvarproxy_io::close_output_stream(std::string const &output_name) } -int colvarproxy_io::backup_file(char const *filename) +int colvarproxy_io::backup_file(char const * /* filename */) { // TODO implement this using rename_file() return COLVARS_NOT_IMPLEMENTED; @@ -785,13 +739,33 @@ int colvarproxy_io::backup_file(char const *filename) int colvarproxy_io::remove_file(char const *filename) { + int error_code = COLVARS_OK; +#if defined(WIN32) && !defined(__CYGWIN__) + // Because the file may be open by other processes, rename it to filename.old + std::string const renamed_file(std::string(filename)+".old"); + // It may still be there from an interrupted run, so remove it to be safe + std::remove(renamed_file.c_str()); + int rename_exit_code = 0; + while ((rename_exit_code = std::rename(filename, + renamed_file.c_str())) != 0) { + if (errno == EINTR) continue; + error_code |= FILE_ERROR; + break; + } + // Ask to remove filename.old, but ignore any errors raised + std::remove(renamed_file.c_str()); +#else if (std::remove(filename)) { if (errno != ENOENT) { - return cvm::error("Error: in removing file \""+std::string(filename)+ - "\".\n.", - FILE_ERROR); + error_code |= FILE_ERROR; } } +#endif + if (error_code != COLVARS_OK) { + return cvm::error("Error: in removing file \""+std::string(filename)+ + "\".\n.", + error_code); + } return COLVARS_OK; } diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 21944106f1..7a43946c52 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -48,11 +48,45 @@ public: /// Destructor virtual ~colvarproxy_system(); - /// \brief Value of the unit for atomic coordinates with respect to - /// angstroms (used by some variables for hard-coded default values) - virtual cvm::real unit_angstrom() = 0; + /// \brief Name of the unit system used internally by Colvars (by default, that of the back-end). + /// Supported depending on the back-end: real (A, kcal/mol), metal (A, eV), electron (Bohr, Hartree), gromacs (nm, kJ/mol) + /// Note: calls to back-end PBC functions assume back-end length unit + /// We use different unit from back-end in VMD bc using PBC functions from colvarproxy base class + /// Colvars internal units are user specified, because the module exchanges info in unknown + /// composite dimensions with user input, while it only exchanges quantities of known + /// dimension with the back-end (length and forces) + std::string units; - /// \brief Boltzmann constant + /// \brief Request to set the units used internally by Colvars + virtual int set_unit_system(std::string const &units, bool check_only) = 0; + + /// \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 + /// * in VMD proxy, can only be changed when no variables are defined + /// as user-defined values in composite units must be compatible with that system + cvm::real angstrom_value; + + /// \brief Value of 1 Angstrom in the backend's unit for atomic coordinates + virtual cvm::real backend_angstrom_value() = 0; + + /// \brief Value of 1 kcal/mol in the internal Colvars unit for energy + cvm::real kcal_mol_value; + + /// \brief Convert a length from Angstrom to internal + inline cvm::real angstrom_to_internal(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.) { + // return l / backend_angstrom_value(); + // } + // return l * angstrom_value / backend_angstrom_value(); + // } + + /// \brief Boltzmann constant in internal Colvars units virtual cvm::real boltzmann() = 0; /// \brief Target temperature of the simulation (K units) @@ -212,7 +246,7 @@ public: } /// Read the current velocity of the given atom - inline cvm::rvector get_atom_velocity(int index) + inline cvm::rvector get_atom_velocity(int /* index */) { cvm::error("Error: reading the current velocity of an atom " "is not yet implemented.\n", @@ -355,7 +389,7 @@ public: } /// Read the current velocity of the given atom group - inline cvm::rvector get_atom_group_velocity(int index) + inline cvm::rvector get_atom_group_velocity(int /* index */) { cvm::error("Error: reading the current velocity of an atom group is not yet implemented.\n", COLVARS_NOT_IMPLEMENTED); @@ -446,15 +480,15 @@ public: virtual ~colvarproxy_replicas(); /// \brief Indicate if multi-replica support is available and active - virtual bool replica_enabled(); + virtual int replica_enabled(); /// \brief Index of this replica virtual int replica_index(); - /// \brief Total number of replica - virtual int replica_num(); + /// \brief Total number of replicas + virtual int num_replicas(); - /// \brief Synchronize replica + /// \brief Synchronize replica with others virtual void replica_comm_barrier(); /// \brief Receive data from other replica @@ -591,10 +625,10 @@ public: return backup_file(filename.c_str()); } - /// Remove the given file + /// Remove the given file (on Windows only, rename to filename.old) int remove_file(char const *filename); - /// Remove the given file + /// Remove the given file (on Windows only, rename to filename.old) inline int remove_file(std::string const &filename) { return remove_file(filename.c_str()); diff --git a/lib/colvars/colvarproxy_replicas.cpp b/lib/colvars/colvarproxy_replicas.cpp new file mode 100644 index 0000000000..1f336d3e44 --- /dev/null +++ b/lib/colvars/colvarproxy_replicas.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 "colvarmodule.h" +#include "colvarproxy.h" + + +colvarproxy_replicas::colvarproxy_replicas() {} + + +colvarproxy_replicas::~colvarproxy_replicas() {} + + +int colvarproxy_replicas::replica_enabled() +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_replicas::replica_index() +{ + return 0; +} + + +int colvarproxy_replicas::num_replicas() +{ + return 1; +} + + +void colvarproxy_replicas::replica_comm_barrier() {} + + +int colvarproxy_replicas::replica_comm_recv(char* /* msg_data */, + int /* buf_len */, + int /* src_rep */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_replicas::replica_comm_send(char* /* msg_data */, + int /* msg_len */, + int /* dest_rep */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index 77b3bd7e38..aa985bedd8 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -1,10 +1,3 @@ #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2019-08-05" -// 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. - +#define COLVARS_VERSION "2020-01-27" #endif diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index 1ad3283337..10b276c89e 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h index 341cb1f72c..69d52cbb51 100644 --- a/lib/colvars/colvarscript.h +++ b/lib/colvars/colvarscript.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -86,6 +86,7 @@ public: cv_printframe, cv_printframelabels, cv_frame, + cv_units, cv_colvar, cv_colvar_value, cv_colvar_update, @@ -269,9 +270,7 @@ extern "C" { "Clear the index groups loaded so far, allowing to replace them", 0, 0, { }, - cvm::main()->index_group_names.clear(); - cvm::main()->index_groups.clear(); - return COLVARS_OK; + return cvm::main()->reset_index_groups(); ) CVSCRIPT(cv_addenergy, @@ -292,6 +291,18 @@ extern "C" { return COLVARS_OK; ) + CVSCRIPT(cv_units, + "Get the current Colvars unit system", + 0, 1, + { }, + if (objc < 3) { + script->set_str_result(cvm::proxy->units); + return COLVARS_OK; + } else { + return cvm::proxy->set_unit_system(script->obj_to_str(objv[2]) , false); + } + ) + #ifndef COLVARSCRIPT_INIT_FN #ifdef __cplusplus } // extern "C" diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index b81a943eab..da00183323 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index 9973d92e98..e5f12154e6 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index accc5defec..7ab617bc44 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index 29b535a1a8..ca367dd43f 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp index 2c43a3fb26..2129461c9a 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.cpp +++ b/src/USER-COLVARS/colvarproxy_lammps.cpp @@ -2,7 +2,7 @@ // 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 +// 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. @@ -110,7 +110,7 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp, restart_output_prefix_str.erase(restart_output_prefix_str.rfind(".*"),2); // initialize multi-replica support, if available - if (replica_enabled()) { + if (replica_enabled() == COLVARS_OK) { MPI_Comm_rank(inter_comm, &inter_me); MPI_Comm_size(inter_comm, &inter_num); } @@ -131,6 +131,11 @@ void colvarproxy_lammps::init(const char *conf_file) cvm::to_str(COLVARPROXY_VERSION)+".\n"); my_angstrom = _lmp->force->angstrom; + // Front-end unit is the same as back-end + angstrom_value = my_angstrom; + + // my_kcal_mol = _lmp->force->qe2f / 23.060549; + // force->qe2f is 1eV expressed in LAMMPS' energy unit (1 if unit is eV, 23 if kcal/mol) my_boltzmann = _lmp->force->boltz; my_timestep = _lmp->update->dt * _lmp->force->femtosecond; @@ -325,6 +330,17 @@ void colvarproxy_lammps::fatal_error(std::string const &message) } +int colvarproxy_lammps::set_unit_system(std::string const &units_in, bool /*check_only*/) +{ + std::string lmp_units = _lmp->update->unit_style; + if (units_in != lmp_units) { + cvm::error("Error: Specified unit system for Colvars \"" + units_in + "\" is incompatible with LAMMPS internal units (" + lmp_units + ").\n"); + return COLVARS_ERROR; + } + return COLVARS_OK; +} + + int colvarproxy_lammps::backup_file(char const *filename) { if (std::string(filename).rfind(std::string(".colvars.state")) @@ -338,6 +354,24 @@ int colvarproxy_lammps::backup_file(char const *filename) // multi-replica support +int colvarproxy_lammps::replica_enabled() +{ + return (inter_comm != MPI_COMM_NULL) ? COLVARS_OK : COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_lammps::replica_index() +{ + return inter_me; +} + + +int colvarproxy_lammps::num_replicas() +{ + return inter_num; +} + + void colvarproxy_lammps::replica_comm_barrier() { MPI_Barrier(inter_comm); diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index 8ed03309a3..29ba5023ed 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -2,7 +2,7 @@ // 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 +// 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. @@ -48,13 +48,6 @@ class colvarproxy_lammps : public colvarproxy { bool total_force_requested; bool do_exit; - // std::vector colvars_atoms; - // std::vector colvars_atoms_ncopies; - // std::vector positions; - // std::vector total_forces; - // std::vector applied_forces; - // std::vector previous_applied_forces; - std::vector atoms_types; MPI_Comm inter_comm; // MPI comm with 1 root proc from each world @@ -98,10 +91,11 @@ class colvarproxy_lammps : public colvarproxy { // read additional config from string void add_config_string(const std::string &config); - // implementation of pure methods from base class - public: + // Request to set the units used internally by Colvars + int set_unit_system(std::string const &units_in, bool check_only); + + inline cvm::real backend_angstrom_value() { return my_angstrom; }; - inline cvm::real unit_angstrom() { return my_angstrom; }; inline cvm::real boltzmann() { return my_boltzmann; }; inline cvm::real temperature() { return t_target; }; inline cvm::real dt() { return my_timestep; }; // return _lmp->update->dt * _lmp->force->femtosecond; }; @@ -127,26 +121,12 @@ class colvarproxy_lammps : public colvarproxy { inline std::vector *modify_atom_types() { return &atoms_types; } - // implementation of optional methods from base class - public: + virtual int replica_enabled(); + virtual int replica_index(); + virtual int num_replicas(); - // Multi-replica support - // Indicate if multi-replica support is available and active - virtual bool replica_enabled() { return (inter_comm != MPI_COMM_NULL); } - - // Index of this replica - virtual int replica_index() { return inter_me; } - - // Total number of replica - virtual int replica_num() { return inter_num; } - - // Synchronize replica virtual void replica_comm_barrier(); - - // Receive data from other replica virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep); - - // Send data to other replica virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep); }; diff --git a/src/USER-COLVARS/colvarproxy_lammps_version.h b/src/USER-COLVARS/colvarproxy_lammps_version.h index d976f40510..7f52fa5e01 100644 --- a/src/USER-COLVARS/colvarproxy_lammps_version.h +++ b/src/USER-COLVARS/colvarproxy_lammps_version.h @@ -1,10 +1,3 @@ #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2019-08-01" -// 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. - +#define COLVARPROXY_VERSION "2019-12-04" #endif diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index 4a57e59fbd..9863e8c20f 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -2,7 +2,7 @@ // 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 +// 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. diff --git a/src/USER-COLVARS/fix_colvars.h b/src/USER-COLVARS/fix_colvars.h index 985457e746..2486446f35 100644 --- a/src/USER-COLVARS/fix_colvars.h +++ b/src/USER-COLVARS/fix_colvars.h @@ -2,7 +2,7 @@ // 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 +// 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.