jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update Colvars library to version 2020-01-27 

One new feature (arithmetic path variables) and several smaller enhancements
or bugfixes.

    Fix error check when loading an old state file
    https://github.com/Colvars/colvars/pull/317

    Get default values for grid boundaries when available
    https://github.com/Colvars/colvars/pull/310

    Allow redefining existing index groups (warn when this happens)
    https://github.com/Colvars/colvars/pull/302

    Simplified replica communication syntax in metadynamics
    https://github.com/Colvars/colvars/pull/301

    Obtain the bias_energy for ABF biases
    https://github.com/Colvars/colvars/pull/294

    Fix reading path file with vector CVCs
    https://github.com/Colvars/colvars/pull/288

    Fix segfault at deleting CVBasedPath
    https://github.com/Colvars/colvars/pull/286

    Parrinello's (arithmetic) pathCV in CV space
    https://github.com/Colvars/colvars/pull/274

    Fix race condition when starting a bundle of walkers
    https://github.com/Colvars/colvars/pull/279
This commit is contained in:
Giacomo Fiorin 2020-01-27 13:39:58 -05:00
parent 5eef3b1828
commit 9427fc50a5
56 changed files with 1876 additions and 1017 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/src/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

3
lib/colvars/Makefile.common
View File

@ -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

237
lib/colvars/Makefile.deps
View File

@ -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

193
lib/colvars/colvar.cpp
View File

@ -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<Lepton::ParsedExpression> 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<colvarvalue const *>(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<colvarvalue const *>(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<typename def_class_name> 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<typename def_class_name> 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<distance>(conf, "distance", "distance");
error_code |= init_components_type<distance_vec>(conf, "distance vector", "distanceVec");
@ -793,9 +854,11 @@ int colvar::init_components(std::string const &conf)
error_code |= init_components_type<eigenvector>(conf, "eigenvector", "eigenvector");
error_code |= init_components_type<gspath>(conf, "geometrical path collective variables (s)", "gspath");
error_code |= init_components_type<gzpath>(conf, "geometrical path collective variables (z)", "gzpath");
error_code |= init_components_type<linearCombination>(conf, "linear combination of other collective variables", "subColvar");
error_code |= init_components_type<linearCombination>(conf, "linear combination of other collective variables", "linearCombination");
error_code |= init_components_type<gspathCV>(conf, "geometrical path collective variables (s) for other CVs", "gspathCV");
error_code |= init_components_type<gzpathCV>(conf, "geometrical path collective variables (z) for other CVs", "gzpathCV");
error_code |= init_components_type<aspathCV>(conf, "arithmetic path collective variables (s) for other CVs", "aspathCV");
error_code |= init_components_type<azpathCV>(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<std::string> const &confs)
}
int colvar::cvc_param_exists(std::string const &param_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 &param_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 &param_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 &param_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 &param_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

73
lib/colvars/colvar.h
View File

@ -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<bool> 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<std::string> const &confs);
/// Whether this named parameter exists (in the first and only component)
int cvc_param_exists(std::string const &param_name);
/// Get the value of the named parameter (from the first and only component)
cvm::real get_cvc_param(std::string const &param_name);
/// Get a pointer to the named parameter (from the first and only component)
void const *get_cvc_param_ptr(std::string const &param_name);
/// Pointer to the derivative of the variable with respect to param_name
colvarvalue const *get_cvc_param_grad(std::string const &param_name);
/// Set the named parameter in the first and only component to the given value
int set_cvc_param(std::string const &param_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<int> atom_ids;
@ -628,34 +654,35 @@ public:
/// For scalar variables only!
std::vector<cvm::rvector> 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<std::vector<int> > 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;

2
lib/colvars/colvar_UIestimator.h
View File

@ -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.

129
lib/colvars/colvar_arithmeticpath.h Normal file
View File

@ -0,0 +1,129 @@
#ifndef ARITHMETICPATHCV_H
#define ARITHMETICPATHCV_H
#include "colvarmodule.h"
#include <vector>
#include <cmath>
#include <iostream>
#include <limits>
namespace ArithmeticPathCV {
using std::vector;
enum path_sz {S, Z};
template <typename element_type, typename scalar_type, path_sz path_type>
class ArithmeticPathBase {
public:
ArithmeticPathBase() {}
virtual ~ArithmeticPathBase() {}
virtual void initialize(size_t p_num_elements, size_t p_total_frames, double p_lambda, const vector<element_type>& p_element, const vector<double>& p_weights);
virtual void updateDistanceToReferenceFrames() = 0;
virtual void computeValue();
virtual void computeDerivatives();
virtual void compute();
protected:
scalar_type lambda;
vector<scalar_type> weights;
size_t num_elements;
size_t total_frames;
vector< vector<element_type> > frame_element_distances;
scalar_type s;
scalar_type z;
vector<element_type> dsdx;
vector<element_type> dzdx;
private:
// intermediate variables
vector<scalar_type> s_numerator_frame;
vector<scalar_type> s_denominator_frame;
scalar_type numerator_s;
scalar_type denominator_s;
scalar_type normalization_factor;
};
template <typename element_type, typename scalar_type, path_sz path_type>
void ArithmeticPathBase<element_type, scalar_type, path_type>::initialize(size_t p_num_elements, size_t p_total_frames, double p_lambda, const vector<element_type>& p_element, const vector<double>& 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<scalar_type>(total_frames - 1);
}
template <typename element_type, typename scalar_type, path_sz path_type>
void ArithmeticPathBase<element_type, scalar_type, path_type>::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<scalar_type>(i_frame) * exponent_tmp;
denominator_s += exponent_tmp;
s_numerator_frame[i_frame] = static_cast<scalar_type>(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 <typename element_type, typename scalar_type, path_sz path_type>
void ArithmeticPathBase<element_type, scalar_type, path_type>::compute() {
computeValue();
computeDerivatives();
}
template <typename element_type, typename scalar_type, path_sz path_type>
void ArithmeticPathBase<element_type, scalar_type, path_type>::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<scalar_type>::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

32
lib/colvars/colvar_geometricpath.h
View File

@ -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 <vector>
#include <cmath>
#include <iostream>
@ -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<element_type>& 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<element_type, scalar_type, path_type>::initialize(size_t
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::initialize(size_t vector_size, const std::vector<element_type>& elements, size_t total_frames, bool p_use_second_closest_frame, bool p_use_third_closest_frame, bool p_use_z_square) {
void GeometricPathBase<element_type, scalar_type, path_type>::initialize(size_t /* vector_size */, const std::vector<element_type>& 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<element_type, scalar_type, path_type>::initialize(size_t
m = static_cast<scalar_type>(1.0);
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::prepareVectors() {
std::cout << "Warning: you should not call the prepareVectors() in base class!\n";
std::cout << std::flush;
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::updateReferenceDistances() {
std::cout << "Warning: you should not call the updateReferenceDistances() in base class!\n";
std::cout << std::flush;
}
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::compute() {
computeValue();
@ -182,7 +172,7 @@ void GeometricPathBase<element_type, scalar_type, path_type>::determineClosestFr
// sigma(z) is on the right side of the closest frame
sign = -1;
}
if (std::abs(static_cast<long>(frame_index[0]) - static_cast<long>(frame_index[1])) > 1) {
if (cvm::fabs(static_cast<long>(frame_index[0]) - static_cast<long>(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<element_type, scalar_type, path_type>::determineClosestFr
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::computeValue() {
updateReferenceDistances();
updateDistanceToReferenceFrames();
determineClosestFrames();
prepareVectors();
v1v1 = scalar_type();
@ -218,14 +208,14 @@ void GeometricPathBase<element_type, scalar_type, path_type>::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<element_type, scalar_type, path_type>::computeValue() {
template <typename element_type, typename scalar_type, path_sz path_type>
void GeometricPathBase<element_type, scalar_type, path_type>::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

45
lib/colvars/colvaratoms.cpp
View File

@ -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<std::string> const &index_group_names =
cvm::main()->index_group_names;
std::vector<std::vector<int> *> const &index_groups =
cvm::main()->index_groups;
std::list<std::string>::iterator names_i = cv->index_group_names.begin();
std::list<std::vector<int> >::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<int> 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;
}

2
lib/colvars/colvaratoms.h
View File

@ -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.

47
lib/colvars/colvarbias.cpp
View File

@ -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<colvarvalue> const *)
{
bias_energy = 0.0;
return COLVARS_OK;
}
int colvarbias::calc_forces(std::vector<colvarvalue> 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<size_t>(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;
}

21
lib/colvars/colvarbias.h
View File

@ -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<colvarvalue> const &values =
std::vector<colvarvalue>(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<colvarvalue> 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<colvarvalue> const *values);
/// Send forces to the collective variables
virtual void communicate_forces();

84
lib/colvars/colvarbias_abf.cpp
View File

@ -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<colvarvalue> 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<int> 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<int> 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;
}

6
lib/colvars/colvarbias_abf.h
View File

@ -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<colvarvalue> const *values);
};
#endif

10
lib/colvars/colvarbias_alb.cpp
View File

@ -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<cvm::real> (colvar_centers[i])) - 1) * ssd[i] / (fmax(update_calls,2) - 1);
<< -2.0 * (means[i] / (static_cast<cvm::real>(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;
}

2
lib/colvars/colvarbias_alb.h
View File

@ -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.

4
lib/colvars/colvarbias_histogram.cpp
View File

@ -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();

2
lib/colvars/colvarbias_histogram.h
View File

@ -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.

289
lib/colvars/colvarbias_meta.cpp
View File

@ -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<int> 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<colvarvalue> const &values)
int colvarbias_meta::calc_energy(std::vector<colvarvalue> const *values)
{
size_t ir = 0;
@ -586,8 +619,8 @@ int colvarbias_meta::calc_energy(std::vector<colvarvalue> const &values)
replicas[ir]->bias_energy = 0.0;
}
std::vector<int> const curr_bin = values.size() ?
hills_energy->get_colvars_index(values) :
std::vector<int> 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<colvarvalue> 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<colvarvalue> const &values)
}
int colvarbias_meta::calc_forces(std::vector<colvarvalue> const &values)
int colvarbias_meta::calc_forces(std::vector<colvarvalue> 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<colvarvalue> const &values)
}
}
std::vector<int> const curr_bin = values.size() ?
hills_energy->get_colvars_index(values) :
std::vector<int> 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<colvarvalue> const &values)
void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last,
cvm::real &energy,
std::vector<colvarvalue> const &colvar_values)
std::vector<colvarvalue> const *values)
{
size_t i = 0;
std::vector<colvarvalue> 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<colvarvalue> &forces,
std::vector<colvarvalue> const &values)
std::vector<colvarvalue> 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);
}

13
lib/colvars/colvarbias_meta.h
View File

@ -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<colvarvalue> const &values =
std::vector<colvarvalue>(0));
virtual int calc_forces(std::vector<colvarvalue> const &values =
std::vector<colvarvalue>(0));
virtual int calc_energy(std::vector<colvarvalue> const *values);
virtual int calc_forces(std::vector<colvarvalue> 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<colvarvalue> const &values = std::vector<colvarvalue>(0));
std::vector<colvarvalue> 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<colvarvalue> &forces,
std::vector<colvarvalue> const &values = std::vector<colvarvalue>(0));
std::vector<colvarvalue> const *values);
/// Height of new hills

4
lib/colvars/colvarbias_restraint.cpp
View File

@ -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;

12
lib/colvars/colvarbias_restraint.h
View File

@ -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);

117
lib/colvars/colvarcomp.cpp
View File

@ -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<void *>(&sup_coeff));
register_param("componentExp", reinterpret_cast<void *>(&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<void *>(&period));
register_param("wrapAround", reinterpret_cast<void *>(&wrap_center));
get_keyval_feature(dynamic_cast<colvarparse *>(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<void *>(&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<void *>(&lower_boundary));
register_param("upperBoundary", reinterpret_cast<void *>(&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 &param_name)
{
colvarvalue const *ptr =
reinterpret_cast<colvarvalue const *>(get_param_grad_ptr(param_name));
return ptr != NULL ? ptr : NULL;
}
int colvar::cvc::set_param(std::string const &param_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<cvm::real const *>(new_value));
}
if (param_name.compare("componentExp") == 0) {
sup_np = *(reinterpret_cast<int const *>(new_value));
}
if (is_enabled(f_cvc_periodic)) {
if (param_name.compare("period") == 0) {
period = *(reinterpret_cast<cvm::real const *>(new_value));
}
if (param_name.compare("wrapAround") == 0) {
wrap_center = *(reinterpret_cast<cvm::real const *>(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;
}

110
lib/colvars/colvarcomp.h
View File

@ -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 <functional>
#endif // C++11 checking
#endif
#include <map>
@ -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<cvm::atom_group *> 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 &param_name);
/// Set the named parameter to the given value
virtual int set_param(std::string const &param_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<cvm::real>& result);
virtual void computeDistanceToReferenceFrames(std::vector<cvm::real>& 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<cvm::real>& result);
virtual void computeDistanceToReferenceFrames(std::vector<cvm::real>& result);
/// Helper function to determine the distance between reference frames
virtual void computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& 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<colvarvalue, cvm::real, GeometricPathCV::path_sz::S>
{
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<colvarvalue, cvm::real, GeometricPathCV::path_sz::Z>
{
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<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::S>
{
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<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::Z>
{
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

28
lib/colvars/colvarcomp_angles.cpp
View File

@ -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);

192
lib/colvars/colvarcomp_apath.cpp Normal file
View File

@ -0,0 +1,192 @@
#if (__cplusplus >= 201103L)
#include <numeric>
#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <limits>
#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<cvm::real> p_weights(cv.size(), 1.0);
get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
std::vector<cvm::real> 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<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::S>::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<cvm::real> p_weights(cv.size(), 1.0);
get_keyval(conf, "weights", p_weights, std::vector<cvm::real>(cv.size(), 1.0));
x.type(colvarvalue::type_scalar);
use_explicit_gradients = true;
std::vector<cvm::real> 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<colvarvalue, cvm::real, ArithmeticPathCV::path_sz::Z>::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

67
lib/colvars/colvarcomp_coordnums.cpp
View File

@ -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

76
lib/colvars/colvarcomp_distances.cpp
View File

@ -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");

59
lib/colvars/colvarcomp_gpath.cpp
View File

@ -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<cvm::real>& result) {
void colvar::CartesianBasedPath::computeDistanceToReferenceFrames(std::vector<cvm::real>& 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<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::S>::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<cvm::atom_pos, cvm::real, GeometricPathCV::path_sz::Z>::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<std::string> 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<cvm::real>& result) {
void colvar::CVBasedPath::computeDistanceToReferenceFrames(std::vector<cvm::real>& 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<cvm::real>& resul
}
}
void colvar::CVBasedPath::computeDistanceBetweenReferenceFrames(std::vector<cvm::real>& 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<colvarvalue, cvm::real, GeometricPathCV::path_sz::S>::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<colvarvalue, cvm::real, GeometricPathCV::path_sz::Z>::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() {

6
lib/colvars/colvarcomp_protein.cpp
View File

@ -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);

37
lib/colvars/colvarcomp_rotations.cpp
View File

@ -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);
}

2
lib/colvars/colvardeps.cpp
View File

@ -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
lib/colvars/colvardeps.h
View File

@ -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
};

2
lib/colvars/colvargrid.cpp
View File

@ -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
lib/colvars/colvargrid.h
View File

@ -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<T> const &g) : nd(g.nd),
colvar_grid(colvar_grid<T> 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

266
lib/colvars/colvarmodule.cpp
View File

@ -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<step_number>(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<std::string>::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<int>());
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<int> *old_index_group = index_groups[index_of_group];
std::vector<int> *new_index_group = new std::vector<int>();
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<std::string>::iterator names_i = index_group_names.begin();
std::list<std::vector<int> >::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<typename T> std::string _to_str(T const &x,
size_t width, size_t prec)
{

38
lib/colvars/colvarmodule.h
View File

@ -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<std::string> index_group_names;
/// \brief Names of groups from one or more Gromacs .ndx files
std::vector<std::string> index_group_names;
/// \brief Groups from a Gromacs .ndx file read at startup
std::list<std::vector<int> > index_groups;
/// \brief Groups from one or more Gromacs .ndx files
std::vector<std::vector<int> *> 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<rvector> *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<rvector> *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<size_t> 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

111
lib/colvars/colvarparams.cpp Normal file
View File

@ -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 <sstream>
#include <iostream>
#include <algorithm>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparams.h"
colvarparams::colvarparams()
{}
colvarparams::~colvarparams()
{}
void colvarparams::register_param(std::string const &param_name,
void *param_ptr)
{
param_map[param_name] = param_ptr;
}
void colvarparams::register_param_grad(std::string const &param_name,
colvarvalue *param_grad_ptr)
{
param_grad_map[param_name] = param_grad_ptr;
}
int colvarparams::param_exists(std::string const &param_name)
{
if (param_map.count(param_name) > 0) {
return COLVARS_OK;
}
return INPUT_ERROR;
}
std::vector<std::string> colvarparams::get_param_names()
{
std::vector<std::string> result;
for (std::map<std::string, void const *>::const_iterator elem =
param_map.begin(); elem != param_map.end(); elem++) {
result.push_back(elem->first);
}
return result;
}
std::vector<std::string> colvarparams::get_param_grad_names()
{
std::vector<std::string> result;
for (std::map<std::string, colvarvalue const *>::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 &param_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 &param_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 &param_name)
{
cvm::real const *ptr =
reinterpret_cast<cvm::real const *>(get_param_ptr(param_name));
return ptr != NULL ? *ptr : 0.0;
}
int colvarparams::set_param(std::string const &param_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);
}

66
lib/colvars/colvarparams.h Normal file
View File

@ -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 <map>
/// \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 &param_name);
/// Get a copy of the names of registered parameters
virtual std::vector<std::string> get_param_names();
/// Get a copy of the names of registered parameter gradients
virtual std::vector<std::string> get_param_grad_names();
/// Pointer to the parameter param_name
virtual void const *get_param_ptr(std::string const &param_name);
/// Pointer to the gradient of parameter param_name
virtual void const *get_param_grad_ptr(std::string const &param_name);
/// Value of the parameter param_name (must be a scalar)
virtual cvm::real get_param(std::string const &param_name);
/// Set the named parameter to the given value
virtual int set_param(std::string const &param_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<std::string, void const *> param_map;
/// Derivatives of the object with respect to internal parameters
std::map<std::string, colvarvalue const *> param_grad_map;
/// Register the given parameter
void register_param(std::string const &param_name, void *param_ptr);
/// Register the gradient of the given parameter
void register_param_grad(std::string const &param_name,
colvarvalue *param_grad_ptr);
};
#endif

43
lib/colvars/colvarparse.cpp
View File

@ -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<typename TYPE>
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);

33
lib/colvars/colvarparse.h
View File

@ -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 <cstring>
#include <string>
#include <map>
#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

128
lib/colvars/colvarproxy.cpp
View File

@ -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<cvm::atom_pos> &pos,
std::vector<int> const &sorted_ids,
std::string const &pdb_field,
int colvarproxy_atoms::load_coords(char const * /* filename */,
std::vector<cvm::atom_pos> & /* pos */,
std::vector<int> 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<int> const &atoms_ids)
int colvarproxy_atom_groups::init_atom_group(std::vector<int> 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<std::string> colvarproxy_script::script_obj_to_str_vector(unsigned char *obj)
std::vector<std::string> 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 colvarvalue *> const &cvcs,
colvarvalue &value)
int colvarproxy_script::run_colvar_callback(std::string const & /* name */,
std::vector<const colvarvalue *> const & /* cvcs */,
colvarvalue & /* value */)
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_script::run_colvar_gradient_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient)
int colvarproxy_script::run_colvar_gradient_callback(std::string const & /* name */,
std::vector<const colvarvalue *> const & /* cvcs */,
std::vector<cvm::matrix2d<cvm::real> > & /* 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 colvarvalue *> const &cvc_values,
colvarvalue &value)
std::string const & name,
std::vector<const colvarvalue *> 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 colvarvalue *> const &cvc_values,
std::vector<cvm::matrix2d<cvm::real> > &gradient)
std::string const & name,
std::vector<const colvarvalue *> const & cvc_values,
std::vector<cvm::matrix2d<cvm::real> > & 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;
}

60
lib/colvars/colvarproxy.h
View File

@ -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());

56
lib/colvars/colvarproxy_replicas.cpp Normal file
View File

@ -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;
}

9
lib/colvars/colvars_version.h
View File

@ -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

2
lib/colvars/colvarscript.cpp
View File

@ -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
lib/colvars/colvarscript.h
View File

@ -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"

2
lib/colvars/colvartypes.cpp
View File

@ -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.

2
lib/colvars/colvartypes.h
View File

@ -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.

2
lib/colvars/colvarvalue.cpp
View File

@ -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.

2
lib/colvars/colvarvalue.h
View File

@ -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
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -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);

36
src/USER-COLVARS/colvarproxy_lammps.h
View File

@ -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<int> colvars_atoms;
// std::vector<size_t> colvars_atoms_ncopies;
// std::vector<struct commdata> positions;
// std::vector<struct commdata> total_forces;
// std::vector<struct commdata> applied_forces;
// std::vector<struct commdata> previous_applied_forces;
std::vector<int> 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<int> *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);
};

9
src/USER-COLVARS/colvarproxy_lammps_version.h
View File

@ -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

2
src/USER-COLVARS/fix_colvars.cpp
View File

@ -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.

2
src/USER-COLVARS/fix_colvars.h
View File

@ -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.