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

Merge pull request #2356 from giacomofiorin/colvars-update 

Update Colvars to version 2020-07-07
This commit is contained in:
Axel Kohlmeyer 2020-09-16 23:37:02 -04:00 committed by GitHub
parent c14fd3131a 85d36ad6c4
commit 4e304177a1
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
50 changed files with 3249 additions and 1586 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
cmake/Modules/Packages/USER-COLVARS.cmake
View File

@ -2,6 +2,8 @@ set(COLVARS_SOURCE_DIR ${LAMMPS_LIB_SOURCE_DIR}/colvars)
file(GLOB COLVARS_SOURCES ${COLVARS_SOURCE_DIR}/[^.]*.cpp)
option(COLVARS_DEBUG "Debugging messages for Colvars (quite verbose)" OFF)
# Build Lepton by default
option(COLVARS_LEPTON "Build and link the Lepton library" ON)
@ -21,6 +23,11 @@ set_target_properties(colvars PROPERTIES OUTPUT_NAME lammps_colvars${LAMMPS_MACH
target_include_directories(colvars PUBLIC ${LAMMPS_LIB_SOURCE_DIR}/colvars)
target_link_libraries(lammps PRIVATE colvars)
if(COLVARS_DEBUG)
# Need to export the macro publicly to also affect the proxy
target_compile_definitions(colvars PUBLIC -DCOLVARS_DEBUG)
endif()
if(COLVARS_LEPTON)
target_link_libraries(lammps PRIVATE lepton)
target_compile_definitions(colvars PRIVATE -DLEPTON)

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

Binary file not shown.

2
examples/USER/colvars/out.colvars.state
View File

@ -1,7 +1,7 @@
configuration {
step 200
dt 2.000000e+00
version 2018-11-16
version 2020-07-07
}
colvar {

2
examples/USER/colvars/out.colvars.state.old
View File

@ -1,7 +1,7 @@
configuration {
step 100
dt 2.000000e+00
version 2018-11-16
version 2020-07-07
}
colvar {

2
examples/USER/colvars/out2.colvars.state
View File

@ -1,7 +1,7 @@
configuration {
step 300
dt 2.000000e+00
version 2018-11-16
version 2020-07-07
}
colvar {

2
examples/USER/colvars/peptide2.colvars.state
View File

@ -1,7 +1,7 @@
configuration {
step 100
dt 2.000000e+00
version 2018-11-16
version 2020-07-07
}
colvar {

26
lib/colvars/Makefile.common
View File

@ -38,6 +38,7 @@ COLVARS_SRCS = \
colvarcomp_gpath.cpp \
colvarcomp_protein.cpp \
colvarcomp_rotations.cpp \
colvarcomp_volmaps.cpp \
colvar.cpp \
colvardeps.cpp \
colvargrid.cpp \
@ -46,7 +47,12 @@ COLVARS_SRCS = \
colvarparse.cpp \
colvarproxy.cpp \
colvarproxy_replicas.cpp \
colvarproxy_tcl.cpp \
colvarproxy_volmaps.cpp \
colvarscript.cpp \
colvarscript_commands.cpp \
colvarscript_commands_bias.cpp \
colvarscript_commands_colvar.cpp \
colvartypes.cpp \
colvarvalue.cpp
@ -61,7 +67,7 @@ ifeq ($(COLVARS_LEPTON),no)
LEPTON_INCFLAGS =
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o)
else
LEPTON_INCFLAGS = -Ilepton/include -DLEPTON -DLEPTON_USE_STATIC_LIBRARIES
LEPTON_INCFLAGS = -Ilepton/include -DLEPTON
COLVARS_OBJS = $(COLVARS_SRCS:.cpp=.o) $(LEPTON_SRCS:.cpp=.o)
endif
@ -77,25 +83,9 @@ Makefile.deps: $(COLVARS_SRCS)
@echo > $@
@for src in $^ ; do \
obj=`basename $$src .cpp`.o ; \
$(CXX) -MM $(COLVARS_INCFLAGS) $(LEPTON_INCFLAGS) \
$(CXX) $(CXXFLAGS) -MM $(COLVARS_INCFLAGS) $(LEPTON_INCFLAGS) \
-MT '$$(COLVARS_OBJ_DIR)'$$obj $$src >> $@ ; \
done
include Makefile.deps
# Exceptions to pattern rule above for Lepton objects
lepton/src/CompiledExpression.o: lepton/src/CompiledExpression.cpp
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
lepton/src/ExpressionProgram.o: lepton/src/ExpressionProgram.cpp
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
lepton/src/ExpressionTreeNode.o: lepton/src/ExpressionTreeNode.cpp
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
lepton/src/Operation.o: lepton/src/Operation.cpp
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
lepton/src/ParsedExpression.o: lepton/src/ParsedExpression.cpp
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
lepton/src/Parser.o: lepton/src/Parser.cpp
$(CXX) $(CXXFLAGS) -Ilepton/include -DLEPTON_BUILDING_STATIC_LIBRARY -c -o $@ $<
include Makefile.lepton.deps # Hand-generated

133
lib/colvars/Makefile.deps
View File

@ -1,84 +1,103 @@
$(COLVARS_OBJ_DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h \
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \
colvarparse.h colvarparams.h colvaratoms.h colvardeps.h
colvarproxy_tcl.h colvarproxy_volmaps.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 colvarparams.h colvardeps.h colvarbias_abf.h colvarbias.h \
colvargrid.h colvar_UIestimator.h
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvar.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 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 colvarparams.h colvardeps.h colvargrid.h
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvarbias.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 colvarparams.h colvardeps.h \
colvarbias_histogram.h colvarbias.h colvargrid.h
colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.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 colvarparams.h colvardeps.h colvarbias_meta.h colvarbias.h \
colvargrid.h
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvar.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 \
colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvarbias_restraint.h colvarbias.h colvar.h colvarparse.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 colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h \
colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h
colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp colvarmodule.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 colvar_geometricpath.h
colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \
colvarmodule.h colvars_version.h colvarparse.h colvarvalue.h \
colvartypes.h colvarparams.h colvaratoms.h colvarproxy.h colvardeps.h \
colvar.h colvarcomp.h colvar_arithmeticpath.h colvar_geometricpath.h
colvartypes.h colvarparams.h colvaratoms.h colvarproxy.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvar.h \
colvarcomp.h colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \
colvars_version.h colvarvalue.h colvartypes.h colvar.h colvarparse.h \
colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \
colvar_arithmeticpath.h colvar_geometricpath.h
colvarproxy_tcl.h colvarproxy_volmaps.h colvar_arithmeticpath.h \
colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \
colvarmodule.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 \
colvar_geometricpath.h
colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp colvarmodule.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 colvar_geometricpath.h
colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \
colvarmodule.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 \
colvar_geometricpath.h
colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \
colvarmodule.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 \
colvar_geometricpath.h
colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvarcomp_volmaps.o: colvarcomp_volmaps.cpp \
colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \
colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \
colvaratoms.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h
$(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.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 colvar_geometricpath.h colvarscript.h \
colvarbias.h
colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_tcl.h \
colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h \
colvarscript.h colvarbias.h colvarscript_commands.h \
colvarscript_commands_colvar.h colvarscript_commands_bias.h
$(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvarmodule.h \
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvardeps.h \
colvarparse.h colvarparams.h
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \
colvarparams.h
$(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.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 colvar_geometricpath.h \
colvargrid.h
colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \
colvar_arithmeticpath.h colvar_geometricpath.h colvargrid.h
$(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \
colvars_version.h colvarparse.h colvarvalue.h colvartypes.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 \
colvar_geometricpath.h
colvarparams.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.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 \
colvarscript_commands.h colvarscript_commands_colvar.h \
colvarscript_commands_bias.h colvaratoms.h colvarcomp.h \
colvar_arithmeticpath.h colvar_geometricpath.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 \
@ -86,15 +105,43 @@ $(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.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 colvarparams.h \
colvardeps.h colvaratoms.h
colvarproxy_tcl.h colvarproxy_volmaps.h colvarscript.h colvarbias.h \
colvar.h colvarparse.h colvarparams.h colvardeps.h \
colvarscript_commands.h colvarscript_commands_colvar.h \
colvarscript_commands_bias.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 colvarparams.h colvardeps.h \
colvarproxy.h
colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h
$(COLVARS_OBJ_DIR)colvarproxy_tcl.o: colvarproxy_tcl.cpp colvarmodule.h \
colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvaratoms.h colvarparse.h \
colvarparams.h colvardeps.h
$(COLVARS_OBJ_DIR)colvarproxy_volmaps.o: colvarproxy_volmaps.cpp \
colvarmodule.h colvars_version.h colvarproxy_volmaps.h
$(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarproxy.h \
colvarmodule.h colvars_version.h colvartypes.h colvarvalue.h \
colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \
colvarparams.h colvarscript.h colvarbias.h colvar.h \
colvarscript_commands.h colvarscript_commands_colvar.h \
colvarscript_commands_bias.h
$(COLVARS_OBJ_DIR)colvarscript_commands.o: colvarscript_commands.cpp \
colvarproxy.h colvarmodule.h colvars_version.h colvartypes.h \
colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h \
colvarparse.h colvarparams.h colvarscript.h colvarbias.h colvar.h \
colvarscript_commands.h colvarscript_commands_colvar.h \
colvarscript_commands_bias.h
$(COLVARS_OBJ_DIR)colvarscript_commands_bias.o: \
colvarscript_commands_bias.cpp colvarproxy.h colvarmodule.h \
colvars_version.h colvartypes.h colvarvalue.h colvarproxy_tcl.h \
colvarproxy_volmaps.h colvardeps.h colvarparse.h colvarparams.h \
colvarscript.h colvarbias.h colvar.h colvarscript_commands.h \
colvarscript_commands_colvar.h colvarscript_commands_bias.h
$(COLVARS_OBJ_DIR)colvarscript_commands_colvar.o: \
colvarscript_commands_colvar.cpp colvarproxy.h colvarmodule.h \
colvars_version.h colvartypes.h colvarvalue.h colvarproxy_tcl.h \
colvarproxy_volmaps.h colvardeps.h colvarparse.h colvarparams.h \
colvarscript.h colvarbias.h colvar.h colvarscript_commands.h \
colvarscript_commands_colvar.h colvarscript_commands_bias.h
$(COLVARS_OBJ_DIR)colvartypes.o: colvartypes.cpp colvarmodule.h \
colvars_version.h colvartypes.h colvarparse.h colvarvalue.h \
colvarparams.h

103
lib/colvars/colvar.cpp
View File

@ -541,7 +541,7 @@ int colvar::init_grid_parameters(std::string const &conf)
cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
"consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
if (!get_keyval(conf, "lowerWall", lower_wall)) {
error_code != cvm::error("Error: the value of lowerWall must be set "
error_code |= cvm::error("Error: the value of lowerWall must be set "
"explicitly.\n", INPUT_ERROR);
}
lw_conf = std::string("\n\
@ -554,7 +554,7 @@ int colvar::init_grid_parameters(std::string const &conf)
cvm::log("Reading legacy options upperWall and upperWallConstant: "
"consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n");
if (!get_keyval(conf, "upperWall", upper_wall)) {
error_code != cvm::error("Error: the value of upperWall must be set "
error_code |= cvm::error("Error: the value of upperWall must be set "
"explicitly.\n", INPUT_ERROR);
}
uw_conf = std::string("\n\
@ -616,7 +616,7 @@ harmonicWalls {\n\
"are enabled).\n", INPUT_ERROR);
}
return COLVARS_OK;
return error_code;
}
@ -681,6 +681,9 @@ int colvar::init_extended_Lagrangian(std::string const &conf)
// Adjust Langevin sigma for slow time step if time_step_factor != 1
ext_sigma = cvm::sqrt(2.0 * cvm::boltzmann() * temp * ext_gamma * ext_mass / (cvm::dt() * cvm::real(time_step_factor)));
}
get_keyval_feature(this, conf, "reflectingLowerBoundary", f_cv_reflecting_lower_boundary, false);
get_keyval_feature(this, conf, "reflectingUpperBoundary", f_cv_reflecting_upper_boundary, false);
}
return COLVARS_OK;
@ -863,6 +866,8 @@ int colvar::init_components(std::string const &conf)
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");
error_code |= init_components_type<map_total>(conf, "total value of atomic map", "mapTotal");
if (!cvcs.size() || (error_code != COLVARS_OK)) {
cvm::error("Error: no valid components were provided "
"for this collective variable.\n",
@ -1040,85 +1045,93 @@ int colvar::init_dependencies() {
init_feature(f_cv_gradient, "gradient", f_type_dynamic);
require_feature_children(f_cv_gradient, f_cvc_gradient);
init_feature(f_cv_collect_gradient, "collect gradient", f_type_dynamic);
init_feature(f_cv_collect_gradient, "collect_gradient", f_type_dynamic);
require_feature_self(f_cv_collect_gradient, f_cv_gradient);
require_feature_self(f_cv_collect_gradient, f_cv_scalar);
// The following exlusion could be lifted by implementing the feature
exclude_feature_self(f_cv_collect_gradient, f_cv_scripted);
require_feature_children(f_cv_collect_gradient, f_cvc_explicit_gradient);
init_feature(f_cv_fdiff_velocity, "velocity from finite differences", f_type_dynamic);
init_feature(f_cv_fdiff_velocity, "velocity_from_finite_differences", f_type_dynamic);
// System force: either trivial (spring force); through extended Lagrangian, or calculated explicitly
init_feature(f_cv_total_force, "total force", f_type_dynamic);
init_feature(f_cv_total_force, "total_force", f_type_dynamic);
require_feature_alt(f_cv_total_force, f_cv_extended_Lagrangian, f_cv_total_force_calc);
// Deps for explicit total force calculation
init_feature(f_cv_total_force_calc, "total force calculation", f_type_dynamic);
init_feature(f_cv_total_force_calc, "total_force_calculation", f_type_dynamic);
require_feature_self(f_cv_total_force_calc, f_cv_scalar);
require_feature_self(f_cv_total_force_calc, f_cv_linear);
require_feature_children(f_cv_total_force_calc, f_cvc_inv_gradient);
require_feature_self(f_cv_total_force_calc, f_cv_Jacobian);
init_feature(f_cv_Jacobian, "Jacobian derivative", f_type_dynamic);
init_feature(f_cv_Jacobian, "Jacobian_derivative", f_type_dynamic);
require_feature_self(f_cv_Jacobian, f_cv_scalar);
require_feature_self(f_cv_Jacobian, f_cv_linear);
require_feature_children(f_cv_Jacobian, f_cvc_Jacobian);
init_feature(f_cv_hide_Jacobian, "hide Jacobian force", f_type_user);
init_feature(f_cv_hide_Jacobian, "hide_Jacobian_force", f_type_user);
require_feature_self(f_cv_hide_Jacobian, f_cv_Jacobian); // can only hide if calculated
init_feature(f_cv_extended_Lagrangian, "extended Lagrangian", f_type_user);
init_feature(f_cv_extended_Lagrangian, "extended_Lagrangian", f_type_user);
require_feature_self(f_cv_extended_Lagrangian, f_cv_scalar);
require_feature_self(f_cv_extended_Lagrangian, f_cv_gradient);
init_feature(f_cv_Langevin, "Langevin dynamics", f_type_user);
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_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);
init_feature(f_cv_output_energy, "output energy", f_type_user);
init_feature(f_cv_output_energy, "output_energy", f_type_user);
init_feature(f_cv_output_value, "output value", f_type_user);
init_feature(f_cv_output_value, "output_value", f_type_user);
init_feature(f_cv_output_velocity, "output velocity", f_type_user);
init_feature(f_cv_output_velocity, "output_velocity", f_type_user);
require_feature_self(f_cv_output_velocity, f_cv_fdiff_velocity);
init_feature(f_cv_output_applied_force, "output applied force", f_type_user);
init_feature(f_cv_output_applied_force, "output_applied_force", f_type_user);
init_feature(f_cv_output_total_force, "output total force", f_type_user);
init_feature(f_cv_output_total_force, "output_total_force", f_type_user);
require_feature_self(f_cv_output_total_force, f_cv_total_force);
init_feature(f_cv_subtract_applied_force, "subtract applied force from total force", f_type_user);
init_feature(f_cv_subtract_applied_force, "subtract_applied_force_from_total_force", f_type_user);
require_feature_self(f_cv_subtract_applied_force, f_cv_total_force);
init_feature(f_cv_lower_boundary, "lower boundary", f_type_user);
init_feature(f_cv_lower_boundary, "lower_boundary", f_type_user);
require_feature_self(f_cv_lower_boundary, f_cv_scalar);
init_feature(f_cv_upper_boundary, "upper boundary", f_type_user);
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);
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);
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_reflecting_lower_boundary, "reflecting_lower_boundary", f_type_user);
require_feature_self(f_cv_reflecting_lower_boundary, f_cv_lower_boundary);
require_feature_self(f_cv_reflecting_lower_boundary, f_cv_extended_Lagrangian);
init_feature(f_cv_reflecting_upper_boundary, "reflecting_upper_boundary", f_type_user);
require_feature_self(f_cv_reflecting_upper_boundary, f_cv_upper_boundary);
require_feature_self(f_cv_reflecting_upper_boundary, f_cv_extended_Lagrangian);
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);
init_feature(f_cv_runave, "running average", f_type_user);
init_feature(f_cv_runave, "running_average", f_type_user);
init_feature(f_cv_corrfunc, "correlation function", f_type_user);
init_feature(f_cv_corrfunc, "correlation_function", f_type_user);
init_feature(f_cv_scripted, "scripted", f_type_user);
init_feature(f_cv_custom_function, "custom function", f_type_user);
init_feature(f_cv_custom_function, "custom_function", f_type_user);
exclude_feature_self(f_cv_custom_function, f_cv_scripted);
init_feature(f_cv_periodic, "periodic", f_type_static);
@ -1129,7 +1142,7 @@ int colvar::init_dependencies() {
// because total forces are obtained from the previous time step,
// we cannot (currently) have colvar values and total forces for the same timestep
init_feature(f_cv_multiple_ts, "multiple timestep colvar", f_type_static);
init_feature(f_cv_multiple_ts, "multiple_timestep", f_type_static);
exclude_feature_self(f_cv_multiple_ts, f_cv_total_force_calc);
// check that everything is initialized
@ -1199,8 +1212,17 @@ colvar::~colvar()
(*ci)->remove_all_children();
delete *ci;
}
cvcs.clear();
// remove reference to this colvar from the CVM
while (biases.size() > 0) {
size_t const i = biases.size()-1;
cvm::log("Warning: before deleting colvar " + name
+ ", deleting related bias " + biases[i]->name);
delete biases[i];
}
biases.clear();
// remove reference to this colvar from the module
colvarmodule *cv = cvm::main();
for (std::vector<colvar *>::iterator cvi = cv->variables()->begin();
cvi != cv->variables()->end();
@ -1211,6 +1233,8 @@ colvar::~colvar()
}
}
cv->config_changed();
#ifdef LEPTON
for (std::vector<Lepton::CompiledExpression *>::iterator cei = value_evaluators.begin();
cei != value_evaluators.end();
@ -1599,6 +1623,15 @@ int colvar::calc_colvar_properties()
// just calculated from the cvcs
if ((cvm::step_relative() == 0 && !after_restart) || x_ext.type() == colvarvalue::type_notset) {
x_ext = x;
if (is_enabled(f_cv_reflecting_lower_boundary) && x_ext < lower_boundary) {
cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below.");
x_ext = lower_boundary;
}
if (is_enabled(f_cv_reflecting_upper_boundary) && x_ext > upper_boundary) {
cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above.");
x_ext = upper_boundary;
}
v_ext.reset(); // (already 0; added for clarity)
}
@ -1672,10 +1705,10 @@ cvm::real colvar::update_forces_energy()
cvm::log("Updating extended-Lagrangian degree of freedom.\n");
}
if (prev_timestep > -1) {
if (prev_timestep > -1L) {
// Keep track of slow timestep to integrate MTS colvars
// the colvar checks the interval after waking up twice
int n_timesteps = cvm::step_relative() - prev_timestep;
cvm::step_number n_timesteps = cvm::step_relative() - prev_timestep;
if (n_timesteps != 0 && n_timesteps != time_step_factor) {
cvm::error("Error: extended-Lagrangian " + description + " has timeStepFactor " +
cvm::to_str(time_step_factor) + ", but was activated after " + cvm::to_str(n_timesteps) +
@ -1737,6 +1770,14 @@ cvm::real colvar::update_forces_energy()
}
v_ext += (0.5 * dt) * f_ext / ext_mass;
x_ext += dt * v_ext;
cvm::real delta = 0; // Length of overshoot past either reflecting boundary
if ((is_enabled(f_cv_reflecting_lower_boundary) && (delta = x_ext - lower_boundary) < 0) ||
(is_enabled(f_cv_reflecting_upper_boundary) && (delta = x_ext - upper_boundary) > 0)) {
x_ext -= 2.0 * delta;
v_ext *= -1.0;
}
x_ext.apply_constraints();
this->wrap(x_ext);
} else {
@ -2082,7 +2123,7 @@ void colvar::wrap(colvarvalue &x_unwrapped) const
std::istream & colvar::read_state(std::istream &is)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
std::string conf;
if ( !(is >> colvarparse::read_block("colvar", &conf)) ) {
@ -2163,7 +2204,7 @@ std::istream & colvar::read_state(std::istream &is)
std::istream & colvar::read_traj(std::istream &is)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
if (is_enabled(f_cv_output_value)) {

3
lib/colvars/colvar.h
View File

@ -602,6 +602,9 @@ public:
class cartesian;
class orientation;
// components that do not handle any atoms directly
class map_total;
protected:
/// \brief Array of \link colvar::cvc \endlink objects

9
lib/colvars/colvar_UIestimator.h
View File

@ -116,7 +116,7 @@ namespace UIestimator {
int i;
for (i = 0; i < dimension; i++) {
temp[i] = round((round(y[i] / width[i] + EPSILON) - round(x[i] / width[i] + EPSILON)) + (y_size - 1) / 2 + EPSILON);
temp[i] = int(round((round(y[i] / width[i] + EPSILON) - round(x[i] / width[i] + EPSILON)) + (y_size - 1) / 2 + EPSILON));
}
int index = 0;
@ -305,12 +305,6 @@ namespace UIestimator {
int i;
if (step % output_freq == 0) {
calc_pmf();
write_files();
//write_interal_data();
}
for (i = 0; i < dimension; i++) {
// for dihedral RC, it is possible that x = 179 and y = -179, should correct it
// may have problem, need to fix
@ -381,6 +375,7 @@ namespace UIestimator {
bool written;
bool written_1D;
public:
// calculate gradients from the internal variables
void calc_pmf() {
int norm;

16
lib/colvars/colvaratoms.cpp
View File

@ -240,19 +240,19 @@ int cvm::atom_group::init_dependencies() {
}
init_feature(f_ag_active, "active", f_type_dynamic);
init_feature(f_ag_center, "translational fit", f_type_static);
init_feature(f_ag_rotate, "rotational fit", f_type_static);
init_feature(f_ag_fitting_group, "fitting group", f_type_static);
init_feature(f_ag_explicit_gradient, "explicit atom gradient", f_type_dynamic);
init_feature(f_ag_fit_gradients, "fit gradients", f_type_user);
init_feature(f_ag_center, "translational_fit", f_type_static);
init_feature(f_ag_rotate, "rotational_fit", f_type_static);
init_feature(f_ag_fitting_group, "fitting_group", f_type_static);
init_feature(f_ag_explicit_gradient, "explicit_atom_gradient", f_type_dynamic);
init_feature(f_ag_fit_gradients, "fit_gradients", f_type_user);
require_feature_self(f_ag_fit_gradients, f_ag_explicit_gradient);
init_feature(f_ag_atom_forces, "atomic forces", f_type_dynamic);
init_feature(f_ag_atom_forces, "atomic_forces", f_type_dynamic);
// parallel calculation implies that we have at least a scalable center of mass,
// but f_ag_scalable is kept as a separate feature to deal with future dependencies
init_feature(f_ag_scalable, "scalable group calculation", f_type_static);
init_feature(f_ag_scalable_com, "scalable group center of mass calculation", f_type_static);
init_feature(f_ag_scalable, "scalable_group", f_type_static);
init_feature(f_ag_scalable_com, "scalable_group_center_of_mass", f_type_static);
require_feature_self(f_ag_scalable, f_ag_scalable_com);
// check that everything is initialized

138
lib/colvars/colvarbias.cpp
View File

@ -7,6 +7,9 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <fstream>
#include <cstring>
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarvalue.h"
@ -23,8 +26,12 @@ colvarbias::colvarbias(char const *key)
init_dependencies();
rank = 1;
time_step_factor = 1;
has_data = false;
b_output_energy = false;
output_freq = cvm::restart_out_freq;
reset();
state_file_step = 0L;
matching_state = false;
@ -93,12 +100,18 @@ int colvarbias::init(std::string const &conf)
output_prefix = cvm::output_prefix();
get_keyval_feature(this, conf, "stepZeroData", f_cvb_step_zero_data, is_enabled(f_cvb_step_zero_data));
// Write energy to traj file?
get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy);
// Disabled by default in base class; default value can be overridden by derived class constructor
get_keyval_feature(this, conf, "bypassExtendedLagrangian", f_cvb_bypass_ext_lagrangian, is_enabled(f_cvb_bypass_ext_lagrangian), parse_silent);
// How often to write full output files?
get_keyval(conf, "outputFreq", output_freq, output_freq);
get_keyval(conf, "timeStepFactor", time_step_factor, 1);
// Disabled by default in base class; default value can be overridden by derived class constructor
get_keyval_feature(this, conf, "bypassExtendedLagrangian", f_cvb_bypass_ext_lagrangian, is_enabled(f_cvb_bypass_ext_lagrangian), parse_echo);
get_keyval(conf, "timeStepFactor", time_step_factor, time_step_factor);
if (time_step_factor < 1) {
cvm::error("Error: timeStepFactor must be 1 or greater.\n");
return COLVARS_ERROR;
@ -125,37 +138,40 @@ int colvarbias::init_dependencies() {
init_feature(f_cvb_awake, "awake", f_type_static);
require_feature_self(f_cvb_awake, f_cvb_active);
init_feature(f_cvb_apply_force, "apply force", f_type_user);
init_feature(f_cvb_step_zero_data, "step_zero_data", f_type_user);
init_feature(f_cvb_apply_force, "apply_force", f_type_user);
require_feature_children(f_cvb_apply_force, f_cv_gradient);
init_feature(f_cvb_bypass_ext_lagrangian, "bypass extended-Lagrangian coordinates", f_type_user);
init_feature(f_cvb_bypass_ext_lagrangian, "bypass_extended_Lagrangian_coordinates", f_type_user);
// The exclusion below prevents the inconsistency where biasing forces are applied onto
// the actual colvar, while total forces are measured on the extended coordinate
exclude_feature_self(f_cvb_bypass_ext_lagrangian, f_cvb_get_total_force);
init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic);
init_feature(f_cvb_get_total_force, "obtain_total_force", f_type_dynamic);
require_feature_children(f_cvb_get_total_force, f_cv_total_force);
init_feature(f_cvb_output_acc_work, "output accumulated work", f_type_user);
init_feature(f_cvb_output_acc_work, "output_accumulated_work", f_type_user);
require_feature_self(f_cvb_output_acc_work, f_cvb_apply_force);
init_feature(f_cvb_history_dependent, "history-dependent", f_type_static);
init_feature(f_cvb_history_dependent, "history_dependent", f_type_static);
init_feature(f_cvb_time_dependent, "time-dependent", f_type_static);
init_feature(f_cvb_time_dependent, "time_dependent", f_type_static);
init_feature(f_cvb_scalar_variables, "require scalar variables", f_type_static);
init_feature(f_cvb_scalar_variables, "require_scalar_variables", f_type_static);
require_feature_children(f_cvb_scalar_variables, f_cv_scalar);
init_feature(f_cvb_calc_pmf, "calculate a PMF", f_type_static);
init_feature(f_cvb_calc_pmf, "calculate_a_PMF", f_type_static);
init_feature(f_cvb_calc_ti_samples, "calculate TI samples", f_type_dynamic);
init_feature(f_cvb_calc_ti_samples, "calculate_TI_samples", f_type_dynamic);
require_feature_self(f_cvb_calc_ti_samples, f_cvb_get_total_force);
require_feature_children(f_cvb_calc_ti_samples, f_cv_grid);
init_feature(f_cvb_write_ti_samples, "write TI samples ", f_type_user);
init_feature(f_cvb_write_ti_samples, "write_TI_samples_", f_type_user);
require_feature_self(f_cvb_write_ti_samples, f_cvb_calc_ti_samples);
init_feature(f_cvb_write_ti_pmf, "write TI PMF", f_type_user);
init_feature(f_cvb_write_ti_pmf, "write_TI_PMF", f_type_user);
require_feature_self(f_cvb_write_ti_pmf, f_cvb_calc_ti_samples);
// check that everything is initialized
@ -237,6 +253,8 @@ int colvarbias::clear()
}
}
cv->config_changed();
return COLVARS_OK;
}
@ -302,6 +320,17 @@ int colvarbias::update()
}
bool colvarbias::can_accumulate_data()
{
colvarproxy *proxy = cvm::main()->proxy;
if (((cvm::step_relative() > 0) && !proxy->simulation_continuing()) ||
is_enabled(f_cvb_step_zero_data)) {
return true;
}
return false;
}
int colvarbias::calc_energy(std::vector<colvarvalue> const *)
{
bias_energy = 0.0;
@ -451,7 +480,7 @@ std::ostream & colvarbias::write_state(std::ostream &os)
std::istream & colvarbias::read_state(std::istream &is)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
std::string key, brace, conf;
if ( !(is >> key) || !(key == state_keyword || key == bias_type) ||
@ -501,12 +530,80 @@ std::istream & colvarbias::read_state(std::istream &is)
}
int colvarbias::write_state_prefix(std::string const &prefix)
{
std::string const filename =
cvm::state_file_prefix(prefix.c_str())+".colvars.state";
std::ostream *os = cvm::proxy->output_stream(filename.c_str());
int error_code = COLVARS_OK;
if (os != NULL) {
os->setf(std::ios::scientific, std::ios::floatfield);
error_code = write_state(*os).good() ? COLVARS_OK : FILE_ERROR;
} else {
error_code = FILE_ERROR;
}
cvm::proxy->close_output_stream(filename.c_str());
return error_code;
}
int colvarbias::write_state_string(std::string &output)
{
std::ostringstream os;
if (!write_state(os)) {
return cvm::error("Error: in writing state of bias \""+name+
"\" to buffer.\n", FILE_ERROR);
}
output = os.str();
return COLVARS_OK;
}
int colvarbias::read_state_prefix(std::string const &prefix)
{
std::string filename((prefix+std::string(".colvars.state")).c_str());
std::ifstream is(filename.c_str());
if (!is.good()) {
// try without the suffix
is.clear();
filename = prefix;
is.open(filename.c_str());
}
return read_state(is).good() ? COLVARS_OK :
cvm::error("Error: in reading state for \""+name+"\" from input file \""+
std::string(filename)+"\".\n", FILE_ERROR);
}
int colvarbias::read_state_string(char const *buffer)
{
if (buffer != NULL) {
size_t const buffer_size = strlen(buffer);
if (cvm::debug()) {
cvm::log("colvarbias::read_state_string() with argument:\n");
cvm::log(buffer);
}
if (buffer_size > 0) {
std::istringstream is;
is.rdbuf()->pubsetbuf(const_cast<char *>(buffer), buffer_size);
return read_state(is).good() ? COLVARS_OK :
cvm::error("Error: in reading state for \""+name+"\" from buffer.\n",
FILE_ERROR);
}
return COLVARS_OK;
}
return cvm::error("Error: NULL pointer for colvarbias::read_state_string()",
BUG_ERROR);
}
std::istream & colvarbias::read_state_data_key(std::istream &is, char const *key)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
std::string key_in;
if ( !(is >> key_in) ||
!(key_in == to_lower_cppstr(std::string(key))) ) {
!(to_lower_cppstr(key_in) == to_lower_cppstr(std::string(key))) ) {
cvm::error("Error: in reading restart configuration for "+
bias_type+" bias \""+this->name+"\" at position "+
cvm::to_str(static_cast<size_t>(is.tellg()))+
@ -546,7 +643,12 @@ std::ostream & colvarbias::write_traj(std::ostream &os)
colvarbias_ti::colvarbias_ti(char const *key)
: colvarbias(key)
{
colvarproxy *proxy = cvm::main()->proxy;
provide(f_cvb_calc_ti_samples);
if (!proxy->total_forces_same_step()) {
// Samples at step zero can not be collected
feature_states[f_cvb_step_zero_data].available = false;
}
ti_avg_forces = NULL;
ti_count = NULL;
}
@ -683,9 +785,11 @@ int colvarbias_ti::update_system_forces(std::vector<colvarvalue> const
(*subtract_forces)[i] : previous_colvar_forces[i]);
}
}
if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
ti_avg_forces->acc_value(ti_bin, ti_system_forces);
}
}
}
if (!proxy->total_forces_same_step()) {
// Set the index for use in the next iteration, when total forces come in

23
lib/colvars/colvarbias.h
View File

@ -57,6 +57,10 @@ public:
/// Some implementations may use calc_energy() and calc_forces()
virtual int update();
/// Returns true if the current step represent a valid increment, whose data
/// can be recorded (as opposed to e.g. a continuation step from a restart)
virtual bool can_accumulate_data();
/// 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
@ -144,14 +148,28 @@ public:
}
/// Read a keyword from the state data (typically a header)
/// \param Input stream
/// \param Keyword labeling the header block
std::istream & read_state_data_key(std::istream &is, char const *key);
/// Write the bias configuration to a restart file or other stream
/// Write the bias configuration to a state file or other stream
std::ostream & write_state(std::ostream &os);
/// Read the bias configuration from a restart file or other stream
std::istream & read_state(std::istream &is);
/// Write the bias state to a file with the given prefix
int write_state_prefix(std::string const &prefix);
/// Write the bias state to a string
int write_state_string(std::string &output);
/// Read the bias state from a file with this name or prefix
int read_state_prefix(std::string const &prefix);
/// Read the bias state from this string buffer
int read_state_string(char const *buffer);
/// Write a label to the trajectory file (comment line)
virtual std::ostream & write_traj_label(std::ostream &os);
@ -164,6 +182,9 @@ public:
return COLVARS_OK;
}
/// Frequency for writing output files
size_t output_freq;
/// Write any output files that this bias may have (e.g. PMF files)
virtual int write_output_files()
{

265
lib/colvars/colvarbias_abf.cpp
View File

@ -8,7 +8,6 @@
// Colvars repository at GitHub.
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias_abf.h"
@ -29,7 +28,13 @@ colvarbias_abf::colvarbias_abf(char const *key)
last_gradients(NULL),
last_samples(NULL)
{
colvarproxy *proxy = cvm::main()->proxy;
if (!proxy->total_forces_same_step()) {
// Samples at step zero can not be collected
feature_states[f_cvb_step_zero_data].available = false;
}
}
int colvarbias_abf::init(std::string const &conf)
{
@ -71,8 +76,19 @@ int colvarbias_abf::init(std::string const &conf)
// full_samples - min_samples >= 1 is guaranteed
get_keyval(conf, "inputPrefix", input_prefix, std::vector<std::string>());
get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
get_keyval(conf, "historyFreq", history_freq, 0);
if (history_freq != 0) {
if (output_freq == 0) {
cvm::error("Error: historyFreq must be a multiple of outputFreq.\n",
INPUT_ERROR);
} else {
if ((history_freq % output_freq) != 0) {
cvm::error("Error: historyFreq must be a multiple of outputFreq.\n",
INPUT_ERROR);
}
}
}
b_history_files = (history_freq > 0);
// shared ABF
@ -146,6 +162,7 @@ int colvarbias_abf::init(std::string const &conf)
for (i = 0; i < num_variables(); i++) {
if (max_force[i] < 0.0) {
cvm::error("Error: maxForce should be non-negative.");
return COLVARS_ERROR;
}
}
cap_force = true;
@ -184,25 +201,24 @@ int colvarbias_abf::init(std::string const &conf)
czar_gradients = new colvar_grid_gradient(colvars);
}
get_keyval(conf, "integrate", b_integrate, num_variables() <= 3); // Integrate for output if d<=3
if (b_integrate) {
// For now, we integrate on-the-fly iff the grid is < 3D
if ( num_variables() <= 3 ) {
if ( num_variables() > 3 ) {
cvm::error("Error: cannot integrate free energy in dimension > 3.\n");
return COLVARS_ERROR;
}
pmf = new integrate_potential(colvars, gradients);
if ( b_CZAR_estimator ) {
czar_pmf = new integrate_potential(colvars, czar_gradients);
}
get_keyval(conf, "integrate", b_integrate, true); // Integrate for output
if ( num_variables() > 1 ) {
// Projected ABF
get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0);
// Parameters for integrating initial (and final) gradient data
get_keyval(conf, "integrateInitMaxIterations", integrate_initial_iterations, 1e4);
get_keyval(conf, "integrateInitTol", integrate_initial_tol, 1e-6);
// for updating the integrated PMF on the fly
get_keyval(conf, "integrateMaxIterations", integrate_iterations, 100);
get_keyval(conf, "integrateTol", integrate_tol, 1e-4);
}
} else {
b_integrate = false;
get_keyval(conf, "integrateMaxIterations", integrate_iterations, 1e4, colvarparse::parse_silent);
get_keyval(conf, "integrateTol", integrate_tol, 1e-6, colvarparse::parse_silent);
// Projected ABF, updating the integrated PMF on the fly
get_keyval(conf, "pABFintegrateFreq", pabf_freq, 0, colvarparse::parse_silent);
get_keyval(conf, "pABFintegrateMaxIterations", pabf_integrate_iterations, 100, colvarparse::parse_silent);
get_keyval(conf, "pABFintegrateTol", pabf_integrate_tol, 1e-4, colvarparse::parse_silent);
}
// For shared ABF, we store a second set of grids.
@ -330,7 +346,7 @@ int colvarbias_abf::update()
force_bin = bin;
}
if (cvm::step_relative() > 0 || cvm::proxy->total_forces_same_step()) {
if (cvm::step_relative() > 0 || is_enabled(f_cvb_step_zero_data)) {
if (update_bias) {
// if (b_adiabatic_reweighting) {
@ -370,10 +386,10 @@ int colvarbias_abf::update()
if ( b_integrate ) {
if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
cvm::real err;
int iter = pmf->integrate(integrate_iterations, integrate_tol, err);
if ( iter == integrate_iterations ) {
cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol)
+ " in " + cvm::to_str(integrate_iterations)
int iter = pmf->integrate(pabf_integrate_iterations, pabf_integrate_tol, err);
if ( iter == pabf_integrate_iterations ) {
cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(pabf_integrate_tol)
+ " in " + cvm::to_str(pabf_integrate_iterations)
+ " steps. Residual error: " + cvm::to_str(err));
}
pmf->set_zero_minimum(); // TODO: do this only when necessary
@ -448,17 +464,6 @@ int colvarbias_abf::update()
output_prefix = cvm::output_prefix() + "." + this->name;
}
if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
if (cvm::debug()) cvm::log("ABF bias trying to write gradients and samples to disk");
write_gradients_samples(output_prefix);
}
if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
// file already exists iff cvm::step_relative() > 0
// otherwise, backup and replace
write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0));
}
if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) {
// Share gradients and samples for shared ABF.
replica_share();
@ -487,10 +492,10 @@ int colvarbias_abf::update()
eabf_UI.update(cvm::step_absolute(), x, y);
}
/// Add the bias energy for 1D ABF
bias_energy = calc_energy(NULL);
/// Compute the bias energy
int error_code = calc_energy(NULL);
return COLVARS_OK;
return error_code;
}
@ -570,83 +575,71 @@ int colvarbias_abf::replica_share() {
last_samples->copy_grid(*samples);
shared_last_step = cvm::step_absolute();
if (b_integrate) {
// Update divergence to account for newly shared gradients
pmf->set_div();
}
return COLVARS_OK;
}
void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool append)
{
std::string samples_out_name = prefix + ".count";
std::string gradients_out_name = prefix + ".grad";
std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
std::ostream *samples_os =
cvm::proxy->output_stream(samples_out_name, mode);
if (!samples_os) return;
samples->write_multicol(*samples_os);
cvm::proxy->close_output_stream(samples_out_name);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string samples_dx_out_name = prefix + ".count.dx";
std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
if (!samples_os) return;
samples->write_opendx(*samples_dx_os);
*samples_dx_os << std::endl;
cvm::proxy->close_output_stream(samples_dx_out_name);
template <class T> int colvarbias_abf::write_grid_to_file(T const *grid,
std::string const &filename,
bool close) {
std::ostream *os = cvm::proxy->output_stream(filename);
if (!os) {
return cvm::error("Error opening file " + filename + " for writing.\n", COLVARS_ERROR | FILE_ERROR);
}
grid->write_multicol(*os);
if (close) {
cvm::proxy->close_output_stream(filename);
} else {
// Insert empty line between frames in history files
*os << std::endl;
cvm::proxy->flush_output_stream(os);
}
std::ostream *gradients_os =
cvm::proxy->output_stream(gradients_out_name, mode);
if (!gradients_os) return;
gradients->write_multicol(*gradients_os);
cvm::proxy->close_output_stream(gradients_out_name);
// In dimension higher than 2, dx is easier to handle and visualize
// but we cannot write multiple frames in a dx file now
// (could be implemented as multiple dx files)
if (num_variables() > 2 && close) {
std::string dx = filename + ".dx";
std::ostream *dx_os = cvm::proxy->output_stream(dx);
if (!dx_os) {
return cvm::error("Error opening file " + dx + " for writing.\n", COLVARS_ERROR | FILE_ERROR);
}
grid->write_opendx(*dx_os);
// if (close) {
cvm::proxy->close_output_stream(dx);
// }
// else {
// // TODO, decide convention for multiple datasets in dx file
// *dx_os << std::endl;
// dx_os->flush();
// }
}
return COLVARS_OK;
}
void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool close)
{
write_grid_to_file<colvar_grid_count>(samples, prefix + ".count", close);
write_grid_to_file<colvar_grid_gradient>(gradients, prefix + ".grad", close);
if (b_integrate) {
// Do numerical integration (to high precision) and output a PMF
cvm::real err;
pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err);
pmf->integrate(integrate_iterations, integrate_tol, err);
pmf->set_zero_minimum();
std::string pmf_out_name = prefix + ".pmf";
std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
if (!pmf_os) return;
pmf->write_multicol(*pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string pmf_dx_out_name = prefix + ".pmf.dx";
std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
if (!pmf_dx_os) return;
pmf->write_opendx(*pmf_dx_os);
*pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(pmf_dx_out_name);
}
*pmf_os << std::endl;
cvm::proxy->close_output_stream(pmf_out_name);
write_grid_to_file<colvar_grid_scalar>(pmf, prefix + ".pmf", close);
}
if (b_CZAR_estimator) {
// Write eABF CZAR-related quantities
std::string z_samples_out_name = prefix + ".zcount";
std::ostream *z_samples_os =
cvm::proxy->output_stream(z_samples_out_name, mode);
if (!z_samples_os) return;
z_samples->write_multicol(*z_samples_os);
cvm::proxy->close_output_stream(z_samples_out_name);
write_grid_to_file<colvar_grid_count>(z_samples, prefix + ".zcount", close);
if (b_czar_window_file) {
std::string z_gradients_out_name = prefix + ".zgrad";
std::ostream *z_gradients_os =
cvm::proxy->output_stream(z_gradients_out_name, mode);
if (!z_gradients_os) return;
z_gradients->write_multicol(*z_gradients_os);
cvm::proxy->close_output_stream(z_gradients_out_name);
write_grid_to_file<colvar_grid_gradient>(z_gradients, prefix + ".zgrad", close);
}
// Calculate CZAR estimator of gradients
@ -657,39 +650,15 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
- cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n);
}
}
std::string czar_gradients_out_name = prefix + ".czar.grad";
std::ostream *czar_gradients_os =
cvm::proxy->output_stream(czar_gradients_out_name, mode);
if (!czar_gradients_os) return;
czar_gradients->write_multicol(*czar_gradients_os);
cvm::proxy->close_output_stream(czar_gradients_out_name);
write_grid_to_file<colvar_grid_gradient>(czar_gradients, prefix + ".czar.grad", close);
if (b_integrate) {
// Do numerical integration (to high precision) and output a PMF
cvm::real err;
czar_pmf->set_div();
czar_pmf->integrate(integrate_initial_iterations, integrate_initial_tol, err);
czar_pmf->integrate(integrate_iterations, integrate_tol, err);
czar_pmf->set_zero_minimum();
std::string czar_pmf_out_name = prefix + ".czar.pmf";
std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
if (!czar_pmf_os) return;
czar_pmf->write_multicol(*czar_pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
if (!czar_pmf_dx_os) return;
czar_pmf->write_opendx(*czar_pmf_dx_os);
*czar_pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
}
*czar_pmf_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_out_name);
write_grid_to_file<colvar_grid_scalar>(czar_pmf, prefix + ".czar.pmf", close);
}
}
return;
@ -836,25 +805,56 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
return is;
}
int colvarbias_abf::write_output_files()
{
write_gradients_samples(output_prefix);
if (cvm::debug()) {
cvm::log("ABF bias trying to write gradients and samples to disk");
}
if (shared_on && cvm::main()->proxy->replica_index() > 0) {
// No need to report the same data as replica 0, let it do the I/O job
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);
write_gradients_samples(output_prefix);
if (b_history_files) {
if ((cvm::step_absolute() % history_freq) == 0) {
write_gradients_samples(output_prefix + ".hist", false);
}
}
if (num_variables() != 1) return 0.0;
if (b_UI_estimator) {
eabf_UI.calc_pmf();
eabf_UI.write_files();
}
return COLVARS_OK;
}
int colvarbias_abf::calc_energy(std::vector<colvarvalue> const *values)
{
bias_energy = 0.0; // default value, overridden if a value can be calculated
if (num_variables() > 1 || values != NULL) {
// Use simple estimate: neglect effect of fullSamples,
// return value at center of bin
if (pmf != NULL) {
std::vector<int> const curr_bin = values ?
pmf->get_colvars_index(*values) :
pmf->get_colvars_index();
if (pmf->index_ok(curr_bin)) {
bias_energy = pmf->value(curr_bin);
}
}
return COLVARS_OK;
}
// Get the home bin.
int home0 = gradients->current_bin_scalar(0);
if (home0 < 0) return 0.0;
if (home0 < 0) return COLVARS_OK;
int gradient_len = (int)(gradients->number_of_points(0));
int home = (home0 < gradient_len) ? home0 : (gradient_len-1);
@ -886,5 +886,6 @@ int colvarbias_abf::calc_energy(std::vector<colvarvalue> const *values)
sum += fact*gradients->value(ix)/count*gradients->widths[0]*frac;
// The applied potential is the negative integral of force samples.
return -sum;
bias_energy = -sum;
return COLVARS_OK;
}

20
lib/colvars/colvarbias_abf.h
View File

@ -15,6 +15,7 @@
#include <sstream>
#include <iomanip>
#include "colvarproxy.h"
#include "colvarbias.h"
#include "colvargrid.h"
#include "colvar_UIestimator.h"
@ -56,8 +57,6 @@ private:
size_t full_samples;
/// Number of samples per bin before applying a scaled-down biasing force
size_t min_samples;
/// frequency for updating output files
int output_freq;
/// Write combined files with a history of all output data?
bool b_history_files;
/// Write CZAR output file for stratified eABF (.zgrad)
@ -74,13 +73,13 @@ private:
/// Frequency for updating pABF PMF (if zero, pABF is not used)
int pabf_freq;
/// Max number of CG iterations for integrating PMF at startup and for file output
int integrate_initial_iterations;
/// Tolerance for integrating PMF at startup and for file output
cvm::real integrate_initial_tol;
/// Max number of CG iterations for integrating PMF at on-the-fly pABF updates
int integrate_iterations;
/// Tolerance for integrating PMF at on-the-fly pABF updates
/// Tolerance for integrating PMF at startup and for file output
cvm::real integrate_tol;
/// Max number of CG iterations for integrating PMF at on-the-fly pABF updates
int pabf_integrate_iterations;
/// Tolerance for integrating PMF at on-the-fly pABF updates
cvm::real pabf_integrate_tol;
/// Cap the biasing force to be applied? (option maxForce)
bool cap_force;
@ -151,11 +150,16 @@ private:
virtual int bin_count(int bin_index);
/// Write human-readable FE gradients and sample count, and DX file in dim > 2
void write_gradients_samples(const std::string &prefix, bool append = false);
void write_gradients_samples(const std::string &prefix, bool close = true);
/// Read human-readable FE gradients and sample count (if not using restart)
void read_gradients_samples();
/// Template used in write_gradient_samples()
template <class T> int write_grid_to_file(T const *grid,
std::string const &name,
bool close);
virtual std::istream& read_state_data(std::istream&);
virtual std::ostream& write_state_data(std::ostream&);
virtual int write_output_files();

17
lib/colvars/colvarbias_histogram.cpp
View File

@ -30,9 +30,10 @@ int colvarbias_histogram::init(std::string const &conf)
size_t i;
get_keyval(conf, "outputFile", out_name, std::string(""));
get_keyval(conf, "outputFileDX", out_name_dx, std::string(""));
get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
get_keyval(conf, "outputFile", out_name, "");
// Write DX file by default only in dimension >= 3
std::string default_name_dx = this->num_variables() > 2 ? "" : "none";
get_keyval(conf, "outputFileDX", out_name_dx, default_name_dx);
/// with VMD, this may not be an error
// if ( output_freq == 0 ) {
@ -146,9 +147,11 @@ int colvarbias_histogram::update()
bin[i] = grid->current_bin_scalar(i);
}
if (can_accumulate_data()) {
if (grid->index_ok(bin)) {
grid->acc_value(bin, 1.0);
}
}
} else {
// update indices for vector/array values
size_t iv, i;
@ -163,10 +166,6 @@ int colvarbias_histogram::update()
}
}
if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
write_output_files();
}
error_code |= cvm::get_error();
return error_code;
}
@ -179,7 +178,7 @@ int colvarbias_histogram::write_output_files()
return COLVARS_OK;
}
if (out_name.size()) {
if (out_name.size() && out_name != "none") {
cvm::log("Writing the histogram file \""+out_name+"\".\n");
cvm::backup_file(out_name.c_str());
std::ostream *grid_os = cvm::proxy->output_stream(out_name);
@ -191,7 +190,7 @@ int colvarbias_histogram::write_output_files()
cvm::proxy->close_output_stream(out_name);
}
if (out_name_dx.size()) {
if (out_name_dx.size() && out_name_dx != "none") {
cvm::log("Writing the histogram file \""+out_name_dx+"\".\n");
cvm::backup_file(out_name_dx.c_str());
std::ostream *grid_os = cvm::proxy->output_stream(out_name_dx);

31
lib/colvars/colvarbias_meta.cpp
View File

@ -566,8 +566,8 @@ int colvarbias_meta::update_grid_params()
int colvarbias_meta::update_bias()
{
// add a new hill if the required time interval has passed
if ((cvm::step_absolute() % new_hill_freq) == 0 &&
is_enabled(f_cvb_history_dependent)) {
if (((cvm::step_absolute() % new_hill_freq) == 0) &&
can_accumulate_data() && is_enabled(f_cvb_history_dependent)) {
if (cvm::debug()) {
cvm::log("Metadynamics bias \""+this->name+"\""+
@ -883,7 +883,7 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
// TODO: improve it by looping over a small subgrid instead of the whole grid
std::vector<colvarvalue> colvar_values(num_variables());
std::vector<colvarvalue> new_colvar_values(num_variables());
std::vector<cvm::real> colvar_forces_scalar(num_variables());
std::vector<int> he_ix = he->new_index();
@ -902,17 +902,17 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
count++) {
size_t i;
for (i = 0; i < num_variables(); i++) {
colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
new_colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
}
// 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, &new_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, &new_colvar_values);
colvar_forces_scalar[i] = hills_forces_here[i].real_value;
}
hg->acc_force(hg_ix, &(colvar_forces_scalar.front()));
@ -935,16 +935,18 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first,
} else {
// TODO delete this (never used)
// simpler version, with just the energy
for ( ; (he->index_ok(he_ix)); ) {
for (size_t i = 0; i < num_variables(); i++) {
colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
new_colvar_values[i] = hills_energy->bin_to_value_scalar(he_ix[i], i);
}
hills_energy_here = 0.0;
calc_hills(h_first, h_last, hills_energy_here, &colvar_values);
calc_hills(h_first, h_last, hills_energy_here, &new_colvar_values);
he->acc_value(he_ix, hills_energy_here);
he->incr(he_ix);
@ -1240,7 +1242,7 @@ void colvarbias_meta::read_replica_files()
// test whether this is the end of the file
is.seekg(0, std::ios::end);
if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) {
if (is.tellg() > (replicas[ir])->replica_hills_file_pos + ((std::streampos) 1)) {
(replicas[ir])->update_status++;
} else {
(replicas[ir])->update_status = 0;
@ -1324,7 +1326,7 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
hills_energy_gradients = new colvar_grid_gradient(colvars);
}
size_t const hills_energy_pos = is.tellg();
std::streampos const hills_energy_pos = is.tellg();
std::string key;
if (!(is >> key)) {
if (hills_energy_backup != NULL) {
@ -1363,7 +1365,7 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is)
}
}
size_t const hills_energy_gradients_pos = is.tellg();
std::streampos const hills_energy_gradients_pos = is.tellg();
if (!(is >> key)) {
if (hills_energy_backup != NULL) {
delete hills_energy;
@ -1511,7 +1513,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is)
{
if (!is) return is; // do nothing if failbit is set
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
size_t i = 0;
std::string data;
@ -1958,11 +1960,12 @@ colvarbias_meta::hill::hill(cvm::step_number it_in,
colvarbias_meta::hill::hill(colvarbias_meta::hill const &h)
: sW(1.0),
: it(h.it),
hill_value(0.0),
sW(1.0),
W(h.W),
centers(h.centers),
sigmas(h.sigmas),
it(h.it),
replica(h.replica)
{}

2
lib/colvars/colvarbias_meta.h
View File

@ -257,7 +257,7 @@ protected:
std::string replica_hills_file;
/// Position within replica_hills_file (when reading it)
int replica_hills_file_pos;
std::streampos replica_hills_file_pos;
};

27
lib/colvars/colvarcomp.cpp
View File

@ -169,7 +169,7 @@ cvm::atom_group *colvar::cvc::parse_group(std::string const &conf,
group->check_keywords(group_conf, group_key);
if (cvm::get_error()) {
cvm::error("Error parsing definition for atom group \""+
std::string(group_key)+"\"\n.", INPUT_ERROR);
std::string(group_key)+"\".", INPUT_ERROR);
}
cvm::decrease_depth();
@ -200,40 +200,40 @@ int colvar::cvc::init_dependencies() {
init_feature(f_cvc_periodic, "periodic", f_type_static);
init_feature(f_cvc_width, "defined width", 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_lower_boundary, "defined_lower_boundary", f_type_static);
init_feature(f_cvc_upper_boundary, "defined upper 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);
init_feature(f_cvc_explicit_gradient, "explicit_gradient", f_type_static);
require_feature_children(f_cvc_explicit_gradient, f_ag_explicit_gradient);
init_feature(f_cvc_inv_gradient, "inverse gradient", f_type_dynamic);
init_feature(f_cvc_inv_gradient, "inverse_gradient", f_type_dynamic);
require_feature_self(f_cvc_inv_gradient, f_cvc_gradient);
init_feature(f_cvc_debug_gradient, "debug gradient", f_type_user);
init_feature(f_cvc_debug_gradient, "debug_gradient", f_type_user);
require_feature_self(f_cvc_debug_gradient, f_cvc_gradient);
require_feature_self(f_cvc_debug_gradient, f_cvc_explicit_gradient);
init_feature(f_cvc_Jacobian, "Jacobian derivative", f_type_dynamic);
init_feature(f_cvc_Jacobian, "Jacobian_derivative", f_type_dynamic);
require_feature_self(f_cvc_Jacobian, f_cvc_inv_gradient);
// 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);
init_feature(f_cvc_one_site_total_force, "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_com_based, "function_of_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_pbc_minimum_image, "use_minimum-image_with_PBCs", f_type_user);
init_feature(f_cvc_scalable, "scalable calculation", f_type_static);
init_feature(f_cvc_scalable, "scalable_calculation", f_type_static);
require_feature_self(f_cvc_scalable, f_cvc_scalable_com);
init_feature(f_cvc_scalable_com, "scalable calculation of centers of mass", f_type_static);
init_feature(f_cvc_scalable_com, "scalable_calculation_of_centers_of_mass", f_type_static);
require_feature_self(f_cvc_scalable_com, f_cvc_com_based);
@ -626,6 +626,7 @@ void colvar::cvc::wrap(colvarvalue & /* x_unwrapped */) const
}
// Static members
std::vector<colvardeps::feature *> colvar::cvc::cvc_features;

35
lib/colvars/colvarcomp.h
View File

@ -1358,8 +1358,14 @@ protected:
cvm::atom_group *atoms;
/// Reference coordinates (for RMSD calculation only)
/// Includes sets with symmetry permutations (n_permutations * n_atoms)
std::vector<cvm::atom_pos> ref_pos;
/// Number of permutations of symmetry-related atoms
size_t n_permutations;
/// Index of the permutation yielding the smallest RMSD (0 for identity)
size_t best_perm_index;
public:
/// Constructor
@ -1661,6 +1667,33 @@ public:
#endif // C++11 checking
// \brief Colvar component: total value of a scalar map
// (usually implemented as a grid by the simulation engine)
class colvar::map_total
: public colvar::cvc
{
public:
map_total();
map_total(std::string const &conf);
virtual ~map_total() {}
virtual int init(std::string const &conf);
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
protected:
/// Identifier of the map object (as used by the simulation engine)
std::string map_name;
/// Index of the map objet in the proxy arrays
int volmap_index;
};
// metrics functions for cvc implementations
// simple definitions of the distance functions; these are useful only
@ -1691,6 +1724,6 @@ public:
{ \
return this->dist2_lgrad(x2, x1); \
} \
\
#endif

58
lib/colvars/colvarcomp_distances.cpp
View File

@ -7,6 +7,8 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <algorithm>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
@ -1045,6 +1047,44 @@ colvar::rmsd::rmsd(std::string const &conf)
// this is only required for ABF, but we do both groups here for better caching
atoms->rot.request_group2_gradients(atoms->size());
}
std::string perm_conf;
size_t pos = 0; // current position in config string
n_permutations = 1;
while (key_lookup(conf, "atomPermutation", &perm_conf, &pos)) {
std::vector<size_t> perm;
if (perm_conf.size()) {
std::istringstream is(perm_conf);
size_t index;
while (is >> index) {
std::vector<int> const &ids = atoms->ids();
size_t const ia = std::find(ids.begin(), ids.end(), index-1) - ids.begin();
if (ia == atoms->size()) {
cvm::error("Error: atom id " + cvm::to_str(index) +
" is not a member of group \"atoms\".");
return;
}
if (std::find(perm.begin(), perm.end(), ia) != perm.end()) {
cvm::error("Error: atom id " + cvm::to_str(index) +
" is mentioned more than once in atomPermutation list.");
return;
}
perm.push_back(ia);
}
if (perm.size() != atoms->size()) {
cvm::error("Error: symmetry permutation in input contains " + cvm::to_str(perm.size()) +
" indices, but group \"atoms\" contains " + cvm::to_str(atoms->size()) + " atoms.");
return;
}
cvm::log("atomPermutation = " + cvm::to_str(perm));
n_permutations++;
// Record a copy of reference positions in new order
for (size_t ia = 0; ia < atoms->size(); ia++) {
ref_pos.push_back(ref_pos[perm[ia]]);
}
}
}
}
@ -1056,6 +1096,20 @@ void colvar::rmsd::calc_value()
for (size_t ia = 0; ia < atoms->size(); ia++) {
x.real_value += ((*atoms)[ia].pos - ref_pos[ia]).norm2();
}
best_perm_index = 0;
// Compute sum of squares for each symmetry permutation of atoms, keep the smallest
size_t ref_pos_index = atoms->size();
for (size_t ip = 1; ip < n_permutations; ip++) {
cvm::real value = 0.0;
for (size_t ia = 0; ia < atoms->size(); ia++) {
value += ((*atoms)[ia].pos - ref_pos[ref_pos_index++]).norm2();
}
if (value < x.real_value) {
x.real_value = value;
best_perm_index = ip;
}
}
x.real_value /= cvm::real(atoms->size()); // MSD
x.real_value = cvm::sqrt(x.real_value);
}
@ -1067,8 +1121,10 @@ void colvar::rmsd::calc_gradients()
0.5 / (x.real_value * cvm::real(atoms->size())) :
0.0;
// Use the appropriate symmetry permutation of reference positions to calculate gradients
size_t const start = atoms->size() * best_perm_index;
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[ia]));
(*atoms)[ia].grad = (drmsddx2 * 2.0 * ((*atoms)[ia].pos - ref_pos[start + ia]));
}
}

60
lib/colvars/colvarcomp_volmaps.cpp Normal file
View File

@ -0,0 +1,60 @@
// -*- 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 "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
colvar::map_total::map_total()
: cvc(), volmap_index(-1)
{
function_type = "map_total";
x.type(colvarvalue::type_scalar);
}
colvar::map_total::map_total(std::string const &conf)
: cvc(), volmap_index(-1)
{
function_type = "map_total";
x.type(colvarvalue::type_scalar);
map_total::init(conf);
}
int colvar::map_total::init(std::string const &conf)
{
int error_code = cvc::init(conf);
get_keyval(conf, "mapName", map_name, map_name);
volmap_index = (cvm::proxy)->init_volmap(map_name);
error_code |= volmap_index > 0 ? COLVARS_OK : INPUT_ERROR;
return error_code;
}
void colvar::map_total::calc_value()
{
x.real_value = (cvm::proxy)->get_volmap_value(volmap_index);
}
void colvar::map_total::calc_gradients()
{
// Atomic coordinates are not available here
}
void colvar::map_total::apply_force(colvarvalue const &force)
{
(cvm::proxy)->apply_volmap_force(volmap_index, force.real_value);
}

4
lib/colvars/colvardeps.cpp
View File

@ -115,7 +115,9 @@ bool colvardeps::get_keyval_feature(colvarparse *cvp,
}
bool value;
bool const found = cvp->get_keyval(conf, key, value, def_value, parse_mode);
if (value) enable(feature_id);
// If the default value is on, this function should be able to disable the feature!
set_enabled(feature_id, value);
return found;
}

6
lib/colvars/colvardeps.h
View File

@ -227,6 +227,8 @@ public:
f_cvb_active,
/// \brief Bias is awake (active on its own accord) this timestep
f_cvb_awake,
/// Accumulates data starting from step 0 of a simulation run
f_cvb_step_zero_data,
/// \brief will apply forces
f_cvb_apply_force,
/// \brief force this bias to act on actual value for extended-Lagrangian coordinates
@ -302,6 +304,10 @@ public:
f_cv_hard_lower_boundary,
/// \brief The upper boundary is not defined from user's choice
f_cv_hard_upper_boundary,
/// \brief Reflecting lower boundary condition
f_cv_reflecting_lower_boundary,
/// \brief Reflecting upper boundary condition
f_cv_reflecting_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)

47
lib/colvars/colvargrid.h
View File

@ -239,19 +239,21 @@ public:
}
/// \brief Constructor from a vector of colvars
/// \param add_extra_bin requests that non-periodic dimensions are extended
/// by 1 bin to accommodate the integral (PMF) of another gridded quantity (gradient)
colvar_grid(std::vector<colvar *> const &colvars,
T const &t = T(),
size_t mult_i = 1,
bool margin = false)
bool add_extra_bin = false)
: has_data(false)
{
this->init_from_colvars(colvars, t, mult_i, margin);
this->init_from_colvars(colvars, t, mult_i, add_extra_bin);
}
int init_from_colvars(std::vector<colvar *> const &colvars,
T const &t = T(),
size_t mult_i = 1,
bool margin = false)
bool add_extra_bin = false)
{
if (cvm::debug()) {
cvm::log("Reading grid configuration from collective variables.\n");
@ -299,7 +301,7 @@ public:
use_actual_value[i-1] = true;
}
if (margin) {
if (add_extra_bin) {
if (periodic[i]) {
// Shift the grid by half the bin width (values at edges instead of center of bins)
lower_boundaries.push_back(cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
@ -699,7 +701,7 @@ public:
}
if (scale_factor != 1.0)
for (size_t i = 0; i < data.size(); i++) {
data[i] += scale_factor * other_grid.data[i];
data[i] += static_cast<T>(scale_factor * other_grid.data[i]);
}
else
// skip multiplication if possible
@ -712,7 +714,7 @@ public:
/// \brief Return the value suitable for output purposes (so that it
/// may be rescaled or manipulated without changing it permanently)
virtual inline T value_output(std::vector<int> const &ix,
size_t const &imult = 0)
size_t const &imult = 0) const
{
return value(ix, imult);
}
@ -820,6 +822,8 @@ public:
std::vector<int> old_nx = nx;
std::vector<colvarvalue> old_lb = lower_boundaries;
std::vector<colvarvalue> old_ub = upper_boundaries;
std::vector<cvm::real> old_w = widths;
{
size_t nd_in = 0;
@ -858,12 +862,15 @@ public:
if (! periodic.size()) periodic.assign(nd, false);
if (! widths.size()) widths.assign(nd, 1.0);
cvm::real eps = 1.e-10;
bool new_params = false;
if (old_nx.size()) {
for (size_t i = 0; i < nd; i++) {
if ( (old_nx[i] != nx[i]) ||
(cvm::sqrt(cv[i]->dist2(old_lb[i],
lower_boundaries[i])) > 1.0E-10) ) {
if (old_nx[i] != nx[i] ||
cvm::sqrt(cv[i]->dist2(old_lb[i], lower_boundaries[i])) > eps ||
cvm::sqrt(cv[i]->dist2(old_ub[i], upper_boundaries[i])) > eps ||
cvm::fabs(old_w[i] - widths[i]) > eps) {
new_params = true;
}
}
@ -928,7 +935,7 @@ public:
/// \brief Read grid entry in restart file
std::istream & read_restart(std::istream &is)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
std::string key, conf;
if ((is >> key) && (key == std::string("grid_parameters"))) {
is.seekg(start_pos, std::ios::beg);
@ -955,7 +962,7 @@ public:
/// represented in memory
/// \param buf_size Number of values per line
std::ostream & write_raw(std::ostream &os,
size_t const buf_size = 3)
size_t const buf_size = 3) const
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
@ -981,7 +988,7 @@ public:
/// \brief Read data written by colvar_grid::write_raw()
std::istream & read_raw(std::istream &is)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
for (size_t imult = 0; imult < mult; imult++) {
@ -1004,7 +1011,7 @@ public:
/// \brief Write the grid in a format which is both human readable
/// and suitable for visualization e.g. with gnuplot
void write_multicol(std::ostream &os)
void write_multicol(std::ostream &os) const
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
@ -1145,7 +1152,7 @@ public:
/// \brief Write the grid data without labels, as they are
/// represented in memory
std::ostream & write_opendx(std::ostream &os)
std::ostream & write_opendx(std::ostream &os) const
{
// write the header
os << "object 1 class gridpositions counts";
@ -1208,7 +1215,7 @@ public:
/// Constructor from a vector of colvars
colvar_grid_count(std::vector<colvar *> &colvars,
size_t const &def_count = 0,
bool margin = false);
bool add_extra_bin = false);
/// Increment the counter at given position
inline void incr_count(std::vector<int> const &ix)
@ -1360,7 +1367,7 @@ public:
/// Constructor from a vector of colvars
colvar_grid_scalar(std::vector<colvar *> &colvars,
bool margin = 0);
bool add_extra_bin = false);
/// Accumulate the value
inline void acc_value(std::vector<int> const &ix,
@ -1419,9 +1426,9 @@ public:
// 100 101 110 111 000 001 010 011
grad[0] = 0.25 * ((p[4] + p[5] + p[6] + p[7]) - (p[0] + p[1] + p[2] + p[3])) / widths[0];
// 010 011 110 111 000 001 100 101
grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[0];
grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[1];
// 001 011 101 111 000 010 100 110
grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[0];
grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[2];
} else {
cvm::error("Finite differences available in dimension 2 and 3 only.");
}
@ -1430,7 +1437,7 @@ public:
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual cvm::real value_output(std::vector<int> const &ix,
size_t const &imult = 0)
size_t const &imult = 0) const
{
if (imult > 0) {
cvm::error("Error: trying to access a component "
@ -1574,7 +1581,7 @@ public:
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual inline cvm::real value_output(std::vector<int> const &ix,
size_t const &imult = 0)
size_t const &imult = 0) const
{
if (samples)
return (samples->value(ix) > 0) ?

283
lib/colvars/colvarmodule.cpp
View File

@ -50,7 +50,7 @@ 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 "
"https://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) {
@ -80,7 +80,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
colvarmodule::rotation::crossing_threshold = 1.0e-02;
cv_traj_freq = 100;
restart_out_freq = proxy->restart_frequency();
restart_out_freq = proxy->default_restart_frequency();
// by default overwrite the existing trajectory file
cv_traj_append = false;
@ -91,7 +91,7 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
colvarmodule * colvarmodule::main()
{
return proxy->colvars;
return proxy ? proxy->colvars : NULL;
}
@ -245,9 +245,6 @@ int colvarmodule::parse_config(std::string &conf)
// Update any necessary proxy data
proxy->setup();
// configuration might have changed, better redo the labels
cv_traj_write_labels = true;
return get_error();
}
@ -265,6 +262,12 @@ int colvarmodule::append_new_config(std::string const &new_conf)
}
void colvarmodule::config_changed()
{
cv_traj_write_labels = true;
}
int colvarmodule::parse_global_params(std::string const &conf)
{
// TODO document and then echo this keyword
@ -332,8 +335,13 @@ int colvarmodule::parse_global_params(std::string const &conf)
scripting_after_biases, scripting_after_biases);
if (use_scripted_forces && !proxy->force_script_defined) {
if (proxy->simulation_running()) {
return cvm::error("User script for scripted colvar forces not found.",
INPUT_ERROR);
} else {
// Not necessary if we are not applying biases in a real simulation (eg. VMD)
cvm::log("Warning: User script for scripted colvar forces not found.");
}
}
return cvm::get_error();
@ -369,6 +377,11 @@ int colvarmodule::parse_colvars(std::string const &conf)
colvar_conf = "";
}
if (pos > 0) {
// One or more new variables were added
config_changed();
}
if (!colvars.size()) {
cvm::log("Warning: no collective variables defined.\n");
}
@ -419,6 +432,10 @@ int colvarmodule::parse_biases_type(std::string const &conf,
}
bias_conf = "";
}
if (conf_saved_pos > 0) {
// One or more new biases were added
config_changed();
}
return COLVARS_OK;
}
@ -656,6 +673,7 @@ std::string colvarmodule::read_colvar(std::string const &name)
return ss.str();
}
cvm::real colvarmodule::energy_difference(std::string const &bias_name,
std::string const &conf)
{
@ -672,75 +690,6 @@ cvm::real colvarmodule::energy_difference(std::string const &bias_name,
return energy_diff;
}
int colvarmodule::bias_current_bin(std::string const &bias_name)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
ret = b->current_bin();
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::bias_bin_num(std::string const &bias_name)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
ret = b->bin_num();
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::bias_bin_count(std::string const &bias_name, size_t bin_index)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
ret = b->bin_count(bin_index);
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::bias_share(std::string const &bias_name)
{
cvm::increase_depth();
int ret;
colvarbias *b = bias_by_name(bias_name);
if (b != NULL) {
b->replica_share();
ret = COLVARS_OK;
} else {
cvm::error("Error: bias not found.\n");
ret = COLVARS_ERROR;
}
cvm::decrease_depth();
return ret;
}
int colvarmodule::calc()
{
@ -753,11 +702,6 @@ int colvarmodule::calc()
}
error_code |= calc_colvars();
// set biasing forces to zero before biases are calculated and summed over
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end(); cvi++) {
(*cvi)->reset_bias_force();
}
error_code |= calc_biases();
error_code |= update_colvar_forces();
@ -768,17 +712,45 @@ int colvarmodule::calc()
error_code |= write_traj_files();
}
// write restart files, if needed
// write restart files and similar data
if (restart_out_freq && (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
if (restart_out_name.size()) {
// Write restart file, if different from main output
error_code |= write_restart_file(restart_out_name);
} else {
error_code |= write_restart_file(output_prefix()+".colvars.state");
}
write_output_files();
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
// TODO remove this when corrFunc becomes a bias
error_code |= (*cvi)->write_output_files();
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
error_code |= (*bi)->write_state_to_replicas();
}
cvm::decrease_depth();
}
// Write output files for biases, at the specified frequency for each
cvm::increase_depth();
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ((*bi)->output_freq > 0) {
if ((cvm::step_relative() > 0) &&
((cvm::step_absolute() % (*bi)->output_freq) == 0) ) {
error_code |= (*bi)->write_output_files();
}
}
}
cvm::decrease_depth();
error_code |= end_of_step();
@ -885,6 +857,12 @@ int colvarmodule::calc_biases()
if (cvm::debug() && num_biases())
cvm::log("Updating collective variable biases.\n");
// set biasing forces to zero before biases are calculated and summed over
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end(); cvi++) {
(*cvi)->reset_bias_force();
}
std::vector<colvarbias *>::iterator bi;
int error_code = COLVARS_OK;
@ -1031,12 +1009,32 @@ int colvarmodule::write_restart_file(std::string const &out_name)
return cvm::error("Error: in writing restart file.\n", FILE_ERROR);
}
proxy->close_output_stream(out_name);
if (cv_traj_os != NULL) {
// Take the opportunity to flush colvars.traj
proxy->flush_output_stream(cv_traj_os);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::write_restart_string(std::string &output)
{
cvm::log("Saving state to output buffer.\n");
std::ostringstream os;
if (!write_restart(os)) {
return cvm::error("Error: in writing restart to buffer.\n", FILE_ERROR);
}
output = os.str();
return COLVARS_OK;
}
int colvarmodule::write_traj_files()
{
if (cvm::debug()) {
cvm::log("colvarmodule::write_traj_files()\n");
}
if (cv_traj_os == NULL) {
if (open_traj_file(cv_traj_name) != COLVARS_OK) {
return cvm::get_error();
@ -1057,15 +1055,12 @@ int colvarmodule::write_traj_files()
write_traj(*cv_traj_os);
}
if (restart_out_freq && (cv_traj_os != NULL)) {
// flush the trajectory file if we are at the restart frequency
if ( (cvm::step_relative() > 0) &&
if (restart_out_freq && (cv_traj_os != NULL) &&
((cvm::step_absolute() % restart_out_freq) == 0)) {
cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+
cv_traj_name+"\".\n");
proxy->flush_output_stream(cv_traj_os);
}
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@ -1182,53 +1177,64 @@ int colvarmodule::reset()
reset_index_groups();
proxy->flush_output_streams();
proxy->reset();
if (cv_traj_os != NULL) {
// Do not close traj file here, as we might not be done with it yet.
proxy->flush_output_stream(cv_traj_os);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::setup_input()
{
std::string restart_in_name("");
// read the restart configuration, if available
if (proxy->input_prefix().size()) {
// read the restart file
restart_in_name = proxy->input_prefix();
// Read a state file
std::string restart_in_name(proxy->input_prefix()+
std::string(".colvars.state"));
std::ifstream input_is(restart_in_name.c_str());
if (!input_is.good()) {
// try by adding the suffix
// Try without the suffix
input_is.clear();
restart_in_name = restart_in_name+std::string(".colvars.state");
restart_in_name = proxy->input_prefix();
input_is.open(restart_in_name.c_str());
}
// Now that the file has been opened, clear this for the next round
proxy->input_prefix().clear();
if (!input_is.good()) {
cvm::error("Error: in opening input file \""+
return cvm::error("Error: in opening input state file \""+
std::string(restart_in_name)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
} else {
cvm::log(cvm::line_marker);
cvm::log("Restarting from file \""+restart_in_name+"\".\n");
cvm::log("Loading state from file \""+restart_in_name+"\".\n");
read_restart(input_is);
if (cvm::get_error() != COLVARS_OK) {
return COLVARS_ERROR;
} else {
proxy->input_prefix().clear();
}
cvm::log(cvm::line_marker);
return cvm::get_error();
}
}
if (proxy->input_buffer() != NULL) {
// Read a string buffer
char const *buffer = proxy->input_buffer();
size_t const buffer_size = strlen(proxy->input_buffer());
// Clear proxy pointer for the next round
proxy->input_buffer() = NULL;
if (buffer_size > 0) {
std::istringstream input_is;
// Replace the buffer of input_is; work around the lack of const in
// pubsetbuf's prototype (which also needs to support output streams)
input_is.rdbuf()->pubsetbuf(const_cast<char *>(buffer), buffer_size);
cvm::log(cvm::line_marker);
cvm::log("Loading state from input buffer.\n");
read_restart(input_is);
cvm::log(cvm::line_marker);
return cvm::get_error();
}
}
return COLVARS_OK;
}
int colvarmodule::setup_output()
@ -1277,6 +1283,20 @@ int colvarmodule::setup_output()
}
std::string colvarmodule::state_file_prefix(char const *filename)
{
std::string const filename_str(filename);
std::string const prefix =
filename_str.substr(0, filename_str.find(".colvars.state"));
if (prefix.size() == 0) {
cvm::error("Error: invalid filename/prefix value \""+filename_str+"\".",
INPUT_ERROR);
}
return prefix;
}
std::istream & colvarmodule::read_restart(std::istream &is)
{
bool warn_total_forces = false;
@ -1340,7 +1360,7 @@ std::istream & colvarmodule::read_restart(std::istream &is)
std::istream & colvarmodule::read_objects_state(std::istream &is)
{
size_t pos = 0;
std::streampos pos = 0;
std::string word;
while (is.good()) {
@ -1365,7 +1385,7 @@ std::istream & colvarmodule::read_objects_state(std::istream &is)
"collective variable \""+(*cvi)->name+"\".\n",
INPUT_ERROR);
}
if (static_cast<size_t>(is.tellg()) > pos) break; // found it
if (is.tellg() > pos) break; // found it
}
cvm::decrease_depth();
@ -1386,13 +1406,13 @@ std::istream & colvarmodule::read_objects_state(std::istream &is)
(*bi)->name+"\".\n",
INPUT_ERROR);
}
if (static_cast<size_t>(is.tellg()) > pos) break; // found it
if (is.tellg() > pos) break; // found it
}
cvm::decrease_depth();
}
}
if (static_cast<size_t>(is.tellg()) == pos) {
if (is.tellg() == pos) {
// This block has not been read by any object: discard it and move on
// to the next one
is >> colvarparse::read_block(word, NULL);
@ -1438,34 +1458,21 @@ int colvarmodule::backup_file(char const *filename)
int colvarmodule::write_output_files()
{
int error_code = COLVARS_OK;
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
error_code |= (*cvi)->write_output_files();
}
cvm::decrease_depth();
cvm::increase_depth();
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
// Only write output files if they have not already been written this time step
if ((*bi)->output_freq == 0 || (cvm::step_absolute() % (*bi)->output_freq) != 0) {
error_code |= (*bi)->write_output_files();
}
error_code |= (*bi)->write_state_to_replicas();
}
cvm::decrease_depth();
if (cv_traj_os != NULL) {
// do not close, there may be another run command
proxy->flush_output_stream(cv_traj_os);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int colvarmodule::read_traj(char const *traj_filename,
long traj_read_begin,
long traj_read_end)
@ -1581,12 +1588,10 @@ int colvarmodule::open_traj_file(std::string const &file_name)
// (re)open trajectory file
if (cv_traj_append) {
cvm::log("Appending to colvar trajectory file \""+file_name+
"\".\n");
cvm::log("Appending to trajectory file \""+file_name+"\".\n");
cv_traj_os = (cvm::proxy)->output_stream(file_name, std::ios::app);
} else {
cvm::log("Writing to colvar trajectory file \""+file_name+
"\".\n");
cvm::log("Opening trajectory file \""+file_name+"\".\n");
proxy->backup_file(file_name.c_str());
cv_traj_os = (cvm::proxy)->output_stream(file_name);
}
@ -1603,6 +1608,7 @@ int colvarmodule::open_traj_file(std::string const &file_name)
int colvarmodule::close_traj_file()
{
if (cv_traj_os != NULL) {
cvm::log("Closing trajectory file \""+cv_traj_name+"\".\n");
proxy->close_output_stream(cv_traj_name);
cv_traj_os = NULL;
}
@ -1668,7 +1674,7 @@ std::ostream & colvarmodule::write_traj(std::ostream &os)
}
void cvm::log(std::string const &message, int min_log_level)
void colvarmodule::log(std::string const &message, int min_log_level)
{
if (cvm::log_level() < min_log_level) return;
// allow logging when the module is not fully initialized
@ -1681,13 +1687,13 @@ void cvm::log(std::string const &message, int min_log_level)
}
void cvm::increase_depth()
void colvarmodule::increase_depth()
{
(depth())++;
}
void cvm::decrease_depth()
void colvarmodule::decrease_depth()
{
if (depth() > 0) {
(depth())--;
@ -1695,7 +1701,7 @@ void cvm::decrease_depth()
}
size_t & cvm::depth()
size_t & colvarmodule::depth()
{
// NOTE: do not call log() or error() here, to avoid recursion
colvarmodule *cv = cvm::main();
@ -1726,16 +1732,19 @@ void colvarmodule::set_error_bits(int code)
proxy->smp_unlock();
}
bool colvarmodule::get_error_bit(int code)
{
return bool(errorCode & code);
}
void colvarmodule::clear_error()
{
proxy->smp_lock();
errorCode = COLVARS_OK;
proxy->smp_unlock();
proxy->clear_error_msgs();
}
@ -1749,9 +1758,7 @@ int colvarmodule::error(std::string const &message, int code)
int colvarmodule::fatal_error(std::string const &message)
{
set_error_bits(FATAL_ERROR);
proxy->fatal_error(message);
return get_error();
return error(message, FATAL_ERROR);
}
@ -1793,7 +1800,7 @@ int cvm::read_index_file(char const *filename)
std::vector<int> *new_index_group = new std::vector<int>();
int atom_number = 1;
size_t pos = is.tellg();
std::streampos pos = is.tellg();
while ( (is >> atom_number) && (atom_number > 0) ) {
new_index_group->push_back(atom_number);
pos = is.tellg();

19
lib/colvars/colvarmodule.h
View File

@ -221,7 +221,6 @@ public:
static void clear_error();
/// Current step number
static step_number it;
/// Starting step number for this run
@ -371,6 +370,9 @@ public:
/// anything that triggers another call
int append_new_config(std::string const &conf);
/// Signals to the module object that the configuration has changed
void config_changed();
private:
/// Configuration string read so far by the module (includes comments)
@ -441,6 +443,9 @@ public:
/// Write the output restart file
std::ostream & write_restart(std::ostream &os);
/// Strips .colvars.state from filename and checks that it is not empty
static std::string state_file_prefix(char const *filename);
/// Open a trajectory file if requested (and leave it open)
int open_traj_file(std::string const &file_name);
/// Close it (note: currently unused)
@ -459,6 +464,9 @@ public:
/// Backup a file before writing it
static int backup_file(char const *filename);
/// Write the state into a string
int write_restart_string(std::string &output);
/// Look up a bias by name; returns NULL if not found
static colvarbias * bias_by_name(std::string const &name);
@ -479,15 +487,6 @@ public:
/// currently works for harmonic (force constant and/or centers)
real energy_difference(std::string const &bias_name, std::string const &conf);
/// Give the total number of bins for a given bias.
int bias_bin_num(std::string const &bias_name);
/// Calculate the bin index for a given bias.
int bias_current_bin(std::string const &bias_name);
//// Give the count at a given bin index.
int bias_bin_count(std::string const &bias_name, size_t bin_index);
//// Share among replicas.
int bias_share(std::string const &bias_name);
/// Main worker function
int calc();

2
lib/colvars/colvarparse.cpp
View File

@ -853,7 +853,7 @@ colvarparse::read_block::~read_block()
std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
{
size_t start_pos = is.tellg();
std::streampos start_pos = is.tellg();
std::string read_key, next;
if ( !(is >> read_key) || !(read_key == rb.key) ||

384
lib/colvars/colvarproxy.cpp
View File

@ -20,11 +20,6 @@
#include <omp.h>
#endif
#if defined(NAMD_TCL) || defined(VMDTCL)
#define COLVARS_TCL
#include <tcl.h>
#endif
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarscript.h"
@ -33,8 +28,9 @@
colvarproxy_system::colvarproxy_system()
: angstrom_value(0.)
{
angstrom_value = 0.0;
total_force_requested = false;
reset_pbc_lattice();
}
@ -67,7 +63,7 @@ bool colvarproxy_system::total_forces_same_step() const
inline int round_to_integer(cvm::real x)
{
return cvm::floor(x+0.5);
return int(cvm::floor(x+0.5));
}
@ -135,6 +131,14 @@ cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1,
}
int colvarproxy_system::get_molid(int &)
{
cvm::error("Error: only VMD allows the use of multiple \"molecules\", "
"i.e. multiple molecular systems.", COLVARS_NOT_IMPLEMENTED);
return -1;
}
colvarproxy_atoms::colvarproxy_atoms()
{
@ -505,144 +509,12 @@ int colvarproxy_script::run_colvar_gradient_callback(std::string const & /* name
colvarproxy_tcl::colvarproxy_tcl()
colvarproxy_io::colvarproxy_io()
{
tcl_interp_ = NULL;
input_buffer_ = NULL;
}
colvarproxy_tcl::~colvarproxy_tcl()
{
}
void colvarproxy_tcl::init_tcl_pointers()
{
cvm::error("Error: Tcl support is currently unavailable "
"outside NAMD or VMD.\n", COLVARS_NOT_IMPLEMENTED);
}
char const *colvarproxy_tcl::tcl_obj_to_str(unsigned char *obj)
{
#if defined(COLVARS_TCL)
return Tcl_GetString(reinterpret_cast<Tcl_Obj *>(obj));
#else
return NULL;
#endif
}
int colvarproxy_tcl::tcl_run_force_callback()
{
#if defined(COLVARS_TCL)
Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
std::string cmd = std::string("calc_colvar_forces ")
+ cvm::to_str(cvm::step_absolute());
int err = Tcl_Eval(tcl_interp, cmd.c_str());
if (err != TCL_OK) {
cvm::log(std::string("Error while executing calc_colvar_forces:\n"));
cvm::error(Tcl_GetStringResult(tcl_interp));
return COLVARS_ERROR;
}
return cvm::get_error();
#else
return COLVARS_NOT_IMPLEMENTED;
#endif
}
int colvarproxy_tcl::tcl_run_colvar_callback(
std::string const & name,
std::vector<const colvarvalue *> const & cvc_values,
colvarvalue & value)
{
#if defined(COLVARS_TCL)
Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
size_t i;
std::string cmd = std::string("calc_") + name;
for (i = 0; i < cvc_values.size(); i++) {
cmd += std::string(" {") + (*(cvc_values[i])).to_simple_string() +
std::string("}");
}
int err = Tcl_Eval(tcl_interp, cmd.c_str());
const char *result = Tcl_GetStringResult(tcl_interp);
if (err != TCL_OK) {
return cvm::error(std::string("Error while executing ")
+ cmd + std::string(":\n") +
std::string(Tcl_GetStringResult(tcl_interp)), COLVARS_ERROR);
}
std::istringstream is(result);
if (value.from_simple_string(is.str()) != COLVARS_OK) {
cvm::log("Error parsing colvar value from script:");
cvm::error(result);
return COLVARS_ERROR;
}
return cvm::get_error();
#else
return COLVARS_NOT_IMPLEMENTED;
#endif
}
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)
{
#if defined(COLVARS_TCL)
Tcl_Interp *const tcl_interp = reinterpret_cast<Tcl_Interp *>(tcl_interp_);
size_t i;
std::string cmd = std::string("calc_") + name + "_gradient";
for (i = 0; i < cvc_values.size(); i++) {
cmd += std::string(" {") + (*(cvc_values[i])).to_simple_string() +
std::string("}");
}
int err = Tcl_Eval(tcl_interp, cmd.c_str());
if (err != TCL_OK) {
return cvm::error(std::string("Error while executing ")
+ cmd + std::string(":\n") +
std::string(Tcl_GetStringResult(tcl_interp)), COLVARS_ERROR);
}
Tcl_Obj **list;
int n;
Tcl_ListObjGetElements(tcl_interp, Tcl_GetObjResult(tcl_interp),
&n, &list);
if (n != int(gradient.size())) {
cvm::error("Error parsing list of gradient values from script: found "
+ cvm::to_str(n) + " values instead of " +
cvm::to_str(gradient.size()));
return COLVARS_ERROR;
}
for (i = 0; i < gradient.size(); i++) {
std::istringstream is(Tcl_GetString(list[i]));
if (gradient[i].from_simple_string(is.str()) != COLVARS_OK) {
cvm::log("Gradient matrix size: " + cvm::to_str(gradient[i].size()));
cvm::log("Gradient string: " + cvm::to_str(Tcl_GetString(list[i])));
cvm::error("Error parsing gradient value from script", COLVARS_ERROR);
return COLVARS_ERROR;
}
}
return cvm::get_error();
#else
return COLVARS_NOT_IMPLEMENTED;
#endif
}
colvarproxy_io::colvarproxy_io() {}
colvarproxy_io::~colvarproxy_io() {}
@ -658,78 +530,6 @@ int colvarproxy_io::set_frame(long int)
}
std::ostream * colvarproxy_io::output_stream(std::string const &output_name,
std::ios_base::openmode mode)
{
if (cvm::debug()) {
cvm::log("Using colvarproxy::output_stream()\n");
}
std::ostream *os = get_output_stream(output_name);
if (os != NULL) return os;
if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
backup_file(output_name);
}
std::ofstream *osf = new std::ofstream(output_name.c_str(), mode);
if (!osf->is_open()) {
cvm::error("Error: cannot write to file/channel \""+output_name+"\".\n",
FILE_ERROR);
return NULL;
}
output_stream_names.push_back(output_name);
output_files.push_back(osf);
return osf;
}
std::ostream *colvarproxy_io::get_output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
return *osi;
}
}
return NULL;
}
int colvarproxy_io::flush_output_stream(std::ostream *os)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osi == os) {
((std::ofstream *) (*osi))->flush();
return COLVARS_OK;
}
}
return cvm::error("Error: trying to flush an output file/channel "
"that wasn't open.\n", BUG_ERROR);
}
int colvarproxy_io::close_output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
((std::ofstream *) (*osi))->close();
delete *osi;
output_files.erase(osi);
output_stream_names.erase(osni);
return COLVARS_OK;
}
}
return cvm::error("Error: trying to close an output file/channel "
"that wasn't open.\n", BUG_ERROR);
}
int colvarproxy_io::backup_file(char const * /* filename */)
{
// TODO implement this using rename_file()
@ -796,11 +596,33 @@ colvarproxy::colvarproxy()
{
colvars = NULL;
b_simulation_running = true;
b_simulation_continuing = false;
b_delete_requested = false;
}
colvarproxy::~colvarproxy() {}
colvarproxy::~colvarproxy()
{
close_files();
}
int colvarproxy::close_files()
{
if (smp_enabled() == COLVARS_OK && smp_thread_id() > 0) {
// Nothing to do on non-master threads
return COLVARS_OK;
}
std::list<std::string>::iterator osni = output_stream_names.begin();
std::list<std::ostream *>::iterator osi = output_files.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
((std::ofstream *) (*osi))->close();
delete *osi;
}
output_files.clear();
output_stream_names.clear();
return COLVARS_OK;
}
int colvarproxy::reset()
@ -838,9 +660,37 @@ int colvarproxy::update_output()
}
size_t colvarproxy::restart_frequency()
int colvarproxy::post_run()
{
return 0;
int error_code = COLVARS_OK;
if (colvars->output_prefix().size()) {
error_code |= colvars->write_restart_file(cvm::output_prefix()+".colvars.state");
error_code |= colvars->write_output_files();
}
error_code |= flush_output_streams();
return error_code;
}
void colvarproxy::add_error_msg(std::string const &message)
{
std::istringstream is(message);
std::string line;
while (std::getline(is, line)) {
error_output += line+"\n";
}
}
void colvarproxy::clear_error_msgs()
{
error_output.clear();
}
std::string const & colvarproxy::get_error_msgs()
{
return error_output;
}
@ -852,3 +702,105 @@ int colvarproxy::get_version_from_string(char const *version_string)
is >> newint;
return newint;
}
void colvarproxy::smp_stream_error()
{
cvm::error("Error: trying to access an output stream from a "
"multi-threaded region (bug). For a quick workaround, use "
"\"smp off\" in the Colvars config.\n", BUG_ERROR);
}
std::ostream * colvarproxy::output_stream(std::string const &output_name,
std::ios_base::openmode mode)
{
if (cvm::debug()) {
cvm::log("Using colvarproxy::output_stream()\n");
}
std::ostream *os = get_output_stream(output_name);
if (os != NULL) return os;
if (!(mode & (std::ios_base::app | std::ios_base::ate))) {
backup_file(output_name);
}
std::ofstream *osf = new std::ofstream(output_name.c_str(), mode);
if (!osf->is_open()) {
cvm::error("Error: cannot write to file/channel \""+output_name+"\".\n",
FILE_ERROR);
return NULL;
}
output_stream_names.push_back(output_name);
output_files.push_back(osf);
return osf;
}
std::ostream *colvarproxy::get_output_stream(std::string const &output_name)
{
if (smp_enabled() == COLVARS_OK) {
if (smp_thread_id() > 0) smp_stream_error();
}
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
return *osi;
}
}
return NULL;
}
int colvarproxy::flush_output_stream(std::ostream *os)
{
if (smp_enabled() == COLVARS_OK) {
if (smp_thread_id() > 0) smp_stream_error();
}
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osi == os) {
((std::ofstream *) (*osi))->flush();
return COLVARS_OK;
}
}
return cvm::error("Error: trying to flush an output file/channel "
"that wasn't open.\n", BUG_ERROR);
}
int colvarproxy::flush_output_streams()
{
if (smp_enabled() == COLVARS_OK && smp_thread_id() > 0)
return COLVARS_OK;
std::list<std::ostream *>::iterator osi = output_files.begin();
for ( ; osi != output_files.end(); osi++) {
((std::ofstream *) (*osi))->flush();
}
return COLVARS_OK;
}
int colvarproxy::close_output_stream(std::string const &output_name)
{
if (smp_enabled() == COLVARS_OK) {
if (smp_thread_id() > 0) smp_stream_error();
}
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
((std::ofstream *) (*osi))->close();
delete *osi;
output_files.erase(osi);
output_stream_names.erase(osni);
return COLVARS_OK;
}
}
return cvm::error("Error: trying to close an output file/channel "
"that wasn't open.\n", BUG_ERROR);
}

169
lib/colvars/colvarproxy.h
View File

@ -16,7 +16,8 @@
#include "colvarmodule.h"
#include "colvartypes.h"
#include "colvarvalue.h"
#include "colvarproxy_tcl.h"
#include "colvarproxy_volmaps.h"
/// \file colvarproxy.h
/// \brief Colvars proxy classes
@ -29,7 +30,7 @@
///
/// To interface to a new MD engine, the simplest solution is to derive a new
/// class from \link colvarproxy \endlink. Currently implemented are: \link
/// colvarproxy_lammps \endlink, \link colvarproxy_namd \endlink, \link
/// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link
/// colvarproxy_vmd \endlink.
@ -121,8 +122,15 @@ public:
/// Are total forces from the current step available?
virtual bool total_forces_same_step() const;
/// Get the molecule ID when called in VMD; raise error otherwise
/// \param molid Set this argument equal to the current VMD molid
virtual int get_molid(int &molid);
protected:
/// Whether the total forces have been requested
bool total_force_requested;
/// \brief Type of boundary conditions
///
/// Orthogonal and triclinic cells are made available to objects.
@ -542,46 +550,6 @@ public:
};
/// Methods for using Tcl within Colvars
class colvarproxy_tcl {
public:
/// Constructor
colvarproxy_tcl();
/// Destructor
virtual ~colvarproxy_tcl();
/// Is Tcl available? (trigger initialization if needed)
int tcl_available();
/// Tcl implementation of script_obj_to_str()
char const *tcl_obj_to_str(unsigned char *obj);
/// Run a user-defined colvar forces script
int tcl_run_force_callback();
int tcl_run_colvar_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value);
int tcl_run_colvar_gradient_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient);
protected:
/// Pointer to Tcl interpreter object
void *tcl_interp_;
/// Set Tcl pointers
virtual void init_tcl_pointers();
};
/// Methods for data input/output
class colvarproxy_io {
@ -601,21 +569,6 @@ public:
// Returns error code
virtual int set_frame(long int);
/// \brief Returns a reference to the given output channel;
/// if this is not open already, then open it
virtual std::ostream *output_stream(std::string const &output_name,
std::ios_base::openmode mode =
std::ios_base::out);
/// Returns a reference to output_name if it exists, NULL otherwise
virtual std::ostream *get_output_stream(std::string const &output_name);
/// \brief Flushes the given output channel
virtual int flush_output_stream(std::ostream *os);
/// \brief Closes the given output channel
virtual int close_output_stream(std::string const &output_name);
/// \brief Rename the given file, before overwriting it
virtual int backup_file(char const *filename);
@ -644,30 +597,49 @@ public:
return rename_file(filename.c_str(), newfilename.c_str());
}
/// \brief Prefix of the input state file
/// Prefix of the input state file to be read next
inline std::string & input_prefix()
{
return input_prefix_str;
}
/// \brief Prefix to be used for output restart files
inline std::string & restart_output_prefix()
{
return restart_output_prefix_str;
}
/// \brief Prefix to be used for output files (final system
/// configuration)
/// Default prefix to be used for all output files (final configuration)
inline std::string & output_prefix()
{
return output_prefix_str;
}
/// Prefix of the restart (checkpoint) file to be written next
inline std::string & restart_output_prefix()
{
return restart_output_prefix_str;
}
/// Default restart frequency (as set by the simulation engine)
inline int default_restart_frequency() const
{
return restart_frequency_engine;
}
/// Buffer from which the input state information may be read
inline char const * & input_buffer()
{
return input_buffer_;
}
protected:
/// \brief Prefix to be used for input files (restarts, not
/// configuration)
std::string input_prefix_str, output_prefix_str, restart_output_prefix_str;
/// Prefix of the input state file to be read next
std::string input_prefix_str;
/// Default prefix to be used for all output files (final configuration)
std::string output_prefix_str;
/// Prefix of the restart (checkpoint) file to be written next
std::string restart_output_prefix_str;
/// How often the simulation engine will write its own restart
int restart_frequency_engine;
/// \brief Currently opened output files: by default, these are ofstream objects.
/// Allows redefinition to implement different output mechanisms
@ -675,6 +647,8 @@ protected:
/// \brief Identifiers for output_stream objects: by default, these are the names of the files
std::list<std::string> output_stream_names;
/// Buffer from which the input state information may be read
char const *input_buffer_;
};
@ -687,6 +661,7 @@ class colvarproxy
: public colvarproxy_system,
public colvarproxy_atoms,
public colvarproxy_atom_groups,
public colvarproxy_volmaps,
public colvarproxy_smp,
public colvarproxy_replicas,
public colvarproxy_script,
@ -717,6 +692,9 @@ public:
/// \brief Reset proxy state, e.g. requested atoms
virtual int reset();
/// Close any open files to prevent data loss
int close_files();
/// (Re)initialize required member data after construction
virtual int setup();
@ -731,14 +709,17 @@ public:
/// Print a message to the main log
virtual void log(std::string const &message) = 0;
/// Print a message to the main log and let the rest of the program handle the error
/// Print a message to the main log and/or let the host code know about it
virtual void error(std::string const &message) = 0;
/// Print a message to the main log and exit with error code
virtual void fatal_error(std::string const &message) = 0;
/// Record error message (used by VMD to collect them after a script call)
void add_error_msg(std::string const &message);
/// \brief Restarts will be written each time this number of steps has passed
virtual size_t restart_frequency();
/// Retrieve accumulated error messages
std::string const & get_error_msgs();
/// As the name says
void clear_error_msgs();
/// Whether a simulation is running (warn against irrecovarable errors)
inline bool simulation_running() const
@ -746,6 +727,17 @@ public:
return b_simulation_running;
}
/// Is the current step a repetition of a step just executed?
/// This is set to true when the step 0 of a new "run" command is being
/// executed, regardless of whether a state file has been loaded.
inline bool simulation_continuing() const
{
return b_simulation_continuing;
}
/// Called at the end of a simulation segment (i.e. "run" command)
int post_run();
/// Convert a version string "YYYY-MM-DD" into an integer
int get_version_from_string(char const *version_string);
@ -755,17 +747,46 @@ public:
return version_int;
}
/// \brief Returns a reference to the given output channel;
/// if this is not open already, then open it
virtual std::ostream *output_stream(std::string const &output_name,
std::ios_base::openmode mode =
std::ios_base::out);
/// Returns a reference to output_name if it exists, NULL otherwise
virtual std::ostream *get_output_stream(std::string const &output_name);
/// \brief Flushes the given output channel
virtual int flush_output_stream(std::ostream *os);
/// \brief Flushes all output channels
virtual int flush_output_streams();
/// \brief Closes the given output channel
virtual int close_output_stream(std::string const &output_name);
protected:
/// Collected error messages
std::string error_output;
/// Whether a simulation is running (warn against irrecovarable errors)
bool b_simulation_running;
/// Is the current step a repetition of a step just executed?
/// This is set to true when the step 0 of a new "run" command is being
/// executed, regardless of whether a state file has been loaded.
bool b_simulation_continuing;
/// Whether the entire module should be deallocated by the host engine
bool b_delete_requested;
/// Integer representing the version string (allows comparisons)
int version_int;
/// Raise when the output stream functions are used on threads other than 0
void smp_stream_error();
};

160
lib/colvars/colvarproxy_tcl.cpp Normal file
View File

@ -0,0 +1,160 @@
// -*- 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>
#if defined(NAMD_TCL) || defined(VMDTCL)
#define COLVARS_TCL
#include <tcl.h>
#endif
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarproxy_tcl.h"
#include "colvaratoms.h"
colvarproxy_tcl::colvarproxy_tcl()
{
tcl_interp_ = NULL;
}
colvarproxy_tcl::~colvarproxy_tcl()
{
}
void colvarproxy_tcl::init_tcl_pointers()
{
cvm::error("Error: Tcl support is not available in this build.\n",
COLVARS_NOT_IMPLEMENTED);
}
char const *colvarproxy_tcl::tcl_get_str(void *obj)
{
#if defined(COLVARS_TCL)
return Tcl_GetString(reinterpret_cast<Tcl_Obj *>(obj));
#else
return NULL;
#endif
}
int colvarproxy_tcl::tcl_run_force_callback()
{
#if defined(COLVARS_TCL)
Tcl_Interp *const tcl_interp =
reinterpret_cast<Tcl_Interp *>(get_tcl_interp());
std::string cmd = std::string("calc_colvar_forces ")
+ cvm::to_str(cvm::step_absolute());
int err = Tcl_Eval(tcl_interp, cmd.c_str());
if (err != TCL_OK) {
cvm::log(std::string("Error while executing calc_colvar_forces:\n"));
cvm::error(Tcl_GetStringResult(tcl_interp));
return COLVARS_ERROR;
}
return cvm::get_error();
#else
return COLVARS_NOT_IMPLEMENTED;
#endif
}
int colvarproxy_tcl::tcl_run_colvar_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvc_values,
colvarvalue &value)
{
#if defined(COLVARS_TCL)
Tcl_Interp *const tcl_interp =
reinterpret_cast<Tcl_Interp *>(get_tcl_interp());
size_t i;
std::string cmd = std::string("calc_") + name;
for (i = 0; i < cvc_values.size(); i++) {
cmd += std::string(" {") + (*(cvc_values[i])).to_simple_string() +
std::string("}");
}
int err = Tcl_Eval(tcl_interp, cmd.c_str());
const char *result = Tcl_GetStringResult(tcl_interp);
if (err != TCL_OK) {
return cvm::error(std::string("Error while executing ")
+ cmd + std::string(":\n") +
std::string(Tcl_GetStringResult(tcl_interp)),
COLVARS_ERROR);
}
std::istringstream is(result);
if (value.from_simple_string(is.str()) != COLVARS_OK) {
cvm::log("Error parsing colvar value from script:");
cvm::error(result);
return COLVARS_ERROR;
}
return cvm::get_error();
#else
return COLVARS_NOT_IMPLEMENTED;
#endif
}
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)
{
#if defined(COLVARS_TCL)
Tcl_Interp *const tcl_interp =
reinterpret_cast<Tcl_Interp *>(get_tcl_interp());
size_t i;
std::string cmd = std::string("calc_") + name + "_gradient";
for (i = 0; i < cvc_values.size(); i++) {
cmd += std::string(" {") + (*(cvc_values[i])).to_simple_string() +
std::string("}");
}
int err = Tcl_Eval(tcl_interp, cmd.c_str());
if (err != TCL_OK) {
return cvm::error(std::string("Error while executing ")
+ cmd + std::string(":\n") +
std::string(Tcl_GetStringResult(tcl_interp)),
COLVARS_ERROR);
}
Tcl_Obj **list;
int n;
Tcl_ListObjGetElements(tcl_interp, Tcl_GetObjResult(tcl_interp),
&n, &list);
if (n != int(gradient.size())) {
cvm::error("Error parsing list of gradient values from script: found "
+ cvm::to_str(n) + " values instead of " +
cvm::to_str(gradient.size()));
return COLVARS_ERROR;
}
for (i = 0; i < gradient.size(); i++) {
std::istringstream is(Tcl_GetString(list[i]));
if (gradient[i].from_simple_string(is.str()) != COLVARS_OK) {
cvm::log("Gradient matrix size: " + cvm::to_str(gradient[i].size()));
cvm::log("Gradient string: " + cvm::to_str(Tcl_GetString(list[i])));
cvm::error("Error parsing gradient value from script", COLVARS_ERROR);
return COLVARS_ERROR;
}
}
return cvm::get_error();
#else
return COLVARS_NOT_IMPLEMENTED;
#endif
}

64
lib/colvars/colvarproxy_tcl.h Normal file
View File

@ -0,0 +1,64 @@
// -*- 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 COLVARPROXY_TCL_H
#define COLVARPROXY_TCL_H
#include <vector>
/// Methods for using Tcl within Colvars
class colvarproxy_tcl {
public:
/// Constructor
colvarproxy_tcl();
/// Destructor
virtual ~colvarproxy_tcl();
/// Is Tcl available? (trigger initialization if needed)
int tcl_available();
/// Tcl implementation of script_obj_to_str()
char const *tcl_get_str(void *obj);
/// Tcl implementation of run_force_callback()
int tcl_run_force_callback();
/// Tcl implementation of run_colvar_callback()
int tcl_run_colvar_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value);
/// Tcl implementation of run_colvar_gradient_callback()
int tcl_run_colvar_gradient_callback(
std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient);
/// Get a pointer to the Tcl interpreter
inline void *get_tcl_interp()
{
return tcl_interp_;
}
protected:
/// Pointer to Tcl interpreter object
void *tcl_interp_;
/// Set Tcl pointers
virtual void init_tcl_pointers();
};
#endif

81
lib/colvars/colvarproxy_volmaps.cpp Normal file
View File

@ -0,0 +1,81 @@
// -*- 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_volmaps.h"
colvarproxy_volmaps::colvarproxy_volmaps() {}
colvarproxy_volmaps::~colvarproxy_volmaps() {}
int colvarproxy_volmaps::volmaps_available()
{
return COLVARS_NOT_IMPLEMENTED;
}
int colvarproxy_volmaps::reset()
{
for (size_t i = 0; i < volmaps_ids.size(); i++) {
clear_volmap(i);
}
volmaps_ids.clear();
volmaps_ncopies.clear();
volmaps_values.clear();
volmaps_new_colvar_forces.clear();
return COLVARS_OK;
}
int colvarproxy_volmaps::add_volmap_slot(int volmap_id)
{
volmaps_ids.push_back(volmap_id);
volmaps_ncopies.push_back(1);
volmaps_values.push_back(0.0);
volmaps_new_colvar_forces.push_back(0.0);
return (volmaps_ids.size() - 1);
}
int colvarproxy_volmaps::init_volmap(int volmap_id)
{
return cvm::error("Error: access to volumetric maps is unavailable "
"in this build.\n",
COLVARS_NOT_IMPLEMENTED);
}
int colvarproxy_volmaps::init_volmap(const char *volmap_name)
{
return cvm::error("Error: access to volumetric maps is unavailable "
"in this build.\n",
COLVARS_NOT_IMPLEMENTED);
}
int colvarproxy_volmaps::init_volmap(const std::string &volmap_name)
{
return init_volmap(volmap_name.c_str());
}
void colvarproxy_volmaps::clear_volmap(int index)
{
if (((size_t) index) >= volmaps_ids.size()) {
cvm::error("Error: trying to unrequest a volumetric map that was not "
"previously requested.\n", INPUT_ERROR);
}
if (volmaps_ncopies[index] > 0) {
volmaps_ncopies[index] -= 1;
}
}

76
lib/colvars/colvarproxy_volmaps.h Normal file
View File

@ -0,0 +1,76 @@
// -*- c++ -*-
#ifndef COLVARPROXY_VOLMAPS_H
#define COLVARPROXY_VOLMAPS_H
/// \brief Container of grid-based objects
class colvarproxy_volmaps {
public:
/// Contructor
colvarproxy_volmaps();
/// Destructor
virtual ~colvarproxy_volmaps();
/// Clear volumetric map data
int reset();
/// \brief Whether this implementation has capability to use volumetric maps
virtual int volmaps_available();
/// Create a slot for a volumetric map not requested yet
int add_volmap_slot(int volmap_id);
/// Request and prepare this volumetric map for use by Colvars
virtual int init_volmap(int volmap_id);
/// Request and prepare this volumetric map for use by Colvars
virtual int init_volmap(char const *volmap_name);
/// Request and prepare this volumetric map for use by Colvars
int init_volmap(std::string const &volmap_name);
/// \brief Used by the CVC destructors
virtual void clear_volmap(int index);
/// Get the numeric ID of the given volumetric map (for the MD program)
inline int get_volmap_id(int index) const
{
return volmaps_ids[index];
}
/// Read the current value of the volumetric map
inline cvm::real get_volmap_value(int index) const
{
return volmaps_values[index];
}
/// Request that this force is applied to the given volumetric map
inline void apply_volmap_force(int index, cvm::real const &new_force)
{
volmaps_new_colvar_forces[index] += new_force;
}
protected:
/// \brief Array of numeric IDs of volumetric maps
std::vector<int> volmaps_ids;
/// \brief Keep track of how many times each vol map is used by a
/// separate colvar object
std::vector<size_t> volmaps_ncopies;
/// \brief Current values of the vol maps
std::vector<cvm::real> volmaps_values;
/// \brief Forces applied from colvars, to be communicated to the MD
/// integrator
std::vector<cvm::real> volmaps_new_colvar_forces;
};
#endif

2
lib/colvars/colvars_version.h
View File

@ -1,3 +1,3 @@
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2020-02-25"
#define COLVARS_VERSION "2020-07-07"
#endif

999
lib/colvars/colvarscript.cpp
View File

File diff suppressed because it is too large Load Diff

470
lib/colvars/colvarscript.h
View File

@ -8,7 +8,7 @@
// Colvars repository at GitHub.
#ifndef COLVARSCRIPT_H
//#define COLVARSCRIPT_H // Delay definition until later
#define COLVARSCRIPT_H
#include <string>
#include <vector>
@ -29,7 +29,7 @@ class colvarscript {
private:
colvarproxy *proxy;
colvarproxy *proxy_;
colvarmodule *colvars;
inline colvarscript() {} // no-argument construction forbidden
@ -39,7 +39,8 @@ public:
friend class colvarproxy;
colvarscript(colvarproxy *p);
inline ~colvarscript() {}
~colvarscript();
/// If an error is caught by the proxy through fatal_error(), this is set to
/// COLVARSCRIPT_ERROR
@ -49,112 +50,196 @@ public:
/// error message
std::string result;
/// Run script command with given positional arguments (objects)
/// Run a script command with space-separated positional arguments (objects)
int run(int objc, unsigned char *const objv[]);
/// Set the return value of the script command to the given string
inline void set_str_result(std::string const &s)
/// Get the string result of the current scripting call
inline std::string const &str_result() const
{
result = s;
return result;
}
/// Build and return a short help
std::string help_string(void) const;
/// Modify the string result of the current scripting call
inline std::string &modify_str_result()
{
return result;
}
/// Set the return value to the given string
int set_result_str(std::string const &s);
/// Clear the string result
int clear_str_result();
/// Add the given string to the error message of the script interface
void add_error_msg(std::string const &s);
/// Commands available
enum command {
#define CVSCRIPT_ENUM_COMM(COMM) COMM,
#undef CVSCRIPT
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_ENUM_COMM(COMM)
#ifdef COLVARSCRIPT_COMMANDS_H
#undef COLVARSCRIPT_COMMANDS_H
#endif
#include "colvarscript_commands.h"
#undef COLVARSCRIPT_COMMANDS_H
#undef CVSCRIPT
#undef CVSCRIPT_ENUM_COMM
cv_n_commands
};
/// Type of object handling a script command
enum Object_type {
use_module,
use_colvar,
use_bias
};
/// Return the prefix of the individual command for each object function
std::string get_cmd_prefix(Object_type t);
/// Get a pointer to the i-th argument of the command (NULL if not given)
template<Object_type T>
unsigned char *get_cmd_arg(int iarg, int objc, unsigned char *const objv[]);
/// Instantiation of get_cmd_arg<> for module-level commands
unsigned char *get_module_cmd_arg(int iarg, int objc,
unsigned char *const objv[]);
/// Instantiation of get_cmd_arg<> for colvar-level commands
unsigned char *get_colvar_cmd_arg(int iarg, int objc,
unsigned char *const objv[]);
/// Instantiation of get_cmd_arg<> for bias-level commands
unsigned char *get_bias_cmd_arg(int iarg, int objc,
unsigned char *const objv[]);
/// Check the argument count of the command
template<Object_type T>
int check_cmd_nargs(char const *cmd, int objc,
int n_args_min, int n_args_max);
/// Instantiation of check_cmd_nargs<> for module-level commands
int check_module_cmd_nargs(char const *cmd, int objc,
int n_args_min, int n_args_max);
/// Instantiation of check_cmd_nargs<> for colvar-level commands
int check_colvar_cmd_nargs(char const *cmd, int objc,
int n_args_min, int n_args_max);
/// Instantiation of get_cmd_arg<> for bias-level commands
int check_bias_cmd_nargs(char const *cmd, int objc,
int n_args_min, int n_args_max);
/// Number of positional arguments to shift for each object type
template<colvarscript::Object_type T>
int cmd_arg_shift();
/// Use scripting language to get the string representation of an object
inline char const *obj_to_str(unsigned char *const obj)
{
return cvm::proxy->script_obj_to_str(obj);
return (obj == NULL ? NULL : proxy_->script_obj_to_str(obj));
}
enum command {
cv_help,
cv_version,
cv_config,
cv_getconfig,
cv_configfile,
cv_reset,
cv_resetindexgroups,
cv_delete,
cv_list,
cv_list_biases,
cv_load,
cv_save,
cv_update,
cv_addenergy,
cv_getenergy,
cv_printframe,
cv_printframelabels,
cv_frame,
cv_units,
cv_colvar,
cv_colvar_value,
cv_colvar_update,
cv_colvar_type,
cv_colvar_delete,
cv_colvar_addforce,
cv_colvar_getappliedforce,
cv_colvar_gettotalforce,
cv_colvar_cvcflags,
cv_colvar_getconfig,
cv_colvar_get,
cv_colvar_set,
cv_bias,
cv_bias_energy,
cv_bias_update,
cv_bias_delete,
cv_bias_getconfig,
cv_bias_get,
cv_bias_set,
cv_n_commands
};
/// Get names of all commands
inline char const **get_command_names() const
{
return cmd_names;
}
/// Get help string for a command (does not specify how it is launched)
/// \param cmd Name of the command's function (e.g. "cv_units")
std::string get_command_help(char const *cmd);
/// Get summary of command line syntax for all commands of a given context
/// \param t One of use_module, use_colvar or use_bias
std::string get_cmdline_help_summary(Object_type t);
/// Get a description of how the command should be used in a command line
/// \param t One of use_module, use_colvar or use_bias
/// \param c Value of the \link command \endlink enum
std::string get_command_cmdline_syntax(Object_type t, command c);
/// Get the command line syntax following by the help string
/// \param t One of use_module, use_colvar or use_bias
/// \param cmd Name of the subcommand (e.g. "units")
std::string get_command_cmdline_help(Object_type t, std::string const &cmd);
/// Set error code for unsupported script operation
int unsupported_op();
/// Pointer to the Colvars main object
inline colvarmodule *module()
{
return this->colvars;
}
/// Pointer to the colvarproxy object (interface with host engine)
inline colvarproxy *proxy()
{
return this->proxy_;
}
private:
/// Set up all script API functions
int init_commands();
/// Set up a single script API function
int init_command(colvarscript::command const &comm,
char const *name, char const *help,
int n_args_min, int n_args_max, char const *arghelp,
int (*fn)(void *, int, unsigned char * const *));
/// Execute a script command
inline int exec_command(command c,
void *pobj,
int objc, unsigned char * const *objv)
{
return (*(comm_fns[c]))(pobj, objc, objv);
return (*(cmd_fns[c]))(pobj, objc, objv);
}
/// Get help for a command (TODO reformat for each language?)
inline std::string command_help(colvarscript::command c) const
{
return comm_help[c];
}
/// Clear all object results
inline void clear_results()
{
result.clear();
}
private:
/// Run subcommands on colvar
int proc_colvar(colvar *cv, int argc, unsigned char *const argv[]);
/// Run subcommands on bias
int proc_bias(colvarbias *b, int argc, unsigned char *const argv[]);
public: // TODO this function will be removed soon
/// Run subcommands on base colvardeps object (colvar, bias, ...)
int proc_features(colvardeps *obj,
int argc, unsigned char *const argv[]);
private: // TODO
/// Internal identifiers of command strings
std::map<std::string, command> comm_str_map;
std::map<std::string, command> cmd_str_map;
/// Inverse of cmd_str_map (to be exported outside this class)
char const **cmd_names;
/// Help strings for each command
std::vector<std::string> comm_help;
std::vector<std::string> cmd_help;
/// Number of arguments for each command
std::vector<size_t> comm_n_args;
/// Minimum number of arguments for each command
std::vector<size_t> cmd_n_args_min;
/// Arguments for each command
std::vector< std::vector<std::string> > comm_args;
/// Maximum number of arguments for each command
std::vector<size_t> cmd_n_args_max;
/// Help strings for each command argument
std::vector< std::vector<std::string> > cmd_arghelp;
/// Implementations of each command
std::vector<int (*)(void *, int, unsigned char * const *)> comm_fns;
std::vector<int (*)(void *, int, unsigned char * const *)> cmd_fns;
/// Get a pointer to the implementation of the given command
inline int (*get_cmd_fn(std::string const &cmd_key))(void *,
int,
unsigned char * const *)
{
if (cmd_str_map.count(cmd_key) > 0) {
return cmd_fns[cmd_str_map[cmd_key]];
}
return NULL;
}
};
@ -165,12 +250,14 @@ inline static colvarscript *colvarscript_obj()
return cvm::main()->proxy->script;
}
/// Get a pointer to the colvar object pointed to by pobj
inline static colvar *colvar_obj(void *pobj)
{
return reinterpret_cast<colvar *>(pobj);
}
/// Get a pointer to the colvarbias object pointed to by pobj
inline static colvarbias *colvarbias_obj(void *pobj)
{
@ -178,137 +265,124 @@ inline static colvarbias *colvarbias_obj(void *pobj)
}
#define CVSCRIPT_COMM_FNAME(COMM) cvscript_ ## COMM
#define CVSCRIPT_COMM_PROTO(COMM) \
int CVSCRIPT_COMM_FNAME(COMM)(void *, int, unsigned char *const *);
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_PROTO(COMM)
#undef COLVARSCRIPT_H
#endif // #ifndef COLVARSCRIPT_H
#ifdef COLVARSCRIPT_CPP
#define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
extern "C" int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \
int objc, \
unsigned char *const objv[]) \
{ \
colvarscript *script = colvarscript_obj(); \
script->clear_results(); \
if (objc < 2+N_ARGS_MIN) /* "cv" and "COMM" are 1st and 2nd */ { \
script->set_str_result("Missing arguments\n" + \
script->command_help(colvarscript::COMM)); \
return COLVARSCRIPT_ERROR; \
} \
if (objc > 2+N_ARGS_MAX) { \
script->set_str_result("Too many arguments\n" + \
script->command_help(colvarscript::COMM)); \
return COLVARSCRIPT_ERROR; \
} \
FN_BODY; \
template<colvarscript::Object_type T>
unsigned char *colvarscript::get_cmd_arg(int iarg,
int objc,
unsigned char *const objv[])
{
int const shift = cmd_arg_shift<T>();
return (shift+iarg < objc) ? objv[shift+iarg] : NULL;
}
#undef CVSCRIPT
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)
#endif // #ifdef COLVARSCRIPT_CPP
#ifdef COLVARSCRIPT_INIT_FN
#define CVSCRIPT_COMM_INIT(COMM,HELP,ARGS) { \
comm_str_map[#COMM] = COMM; \
comm_help[COMM] = HELP; \
comm_fns[COMM] = &(CVSCRIPT_COMM_FNAME(COMM)); \
inline unsigned char *colvarscript::get_module_cmd_arg(int iarg, int objc,
unsigned char *const objv[])
{
return get_cmd_arg<use_module>(iarg, objc, objv);
}
#undef CVSCRIPT
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_INIT(COMM,HELP,ARGS)
#endif
#if !defined(COLVARSCRIPT_H) || defined(COLVARSCRIPT_INIT_FN)
#define COLVARSCRIPT_H
#ifndef COLVARSCRIPT_INIT_FN
#ifdef __cplusplus
extern "C" {
#endif
#endif
// Add optional arguments for command-specific help?
CVSCRIPT(cv_help,
"Print the help message",
0, 0,
{},
script->set_str_result(script->help_string());
return COLVARS_OK;
)
CVSCRIPT(cv_config,
"Read configuration from the given string",
1, 1,
{ "conf (str) - Configuration string" },
std::string const conf(script->obj_to_str(objv[2]));
if (cvm::main()->read_config_string(conf) == COLVARS_OK) {
return COLVARS_OK;
inline unsigned char *colvarscript::get_colvar_cmd_arg(int iarg, int objc,
unsigned char *const objv[])
{
return get_cmd_arg<use_colvar>(iarg, objc, objv);
}
script->set_str_result("Error parsing configuration string");
inline unsigned char *colvarscript::get_bias_cmd_arg(int iarg, int objc,
unsigned char *const objv[])
{
return get_cmd_arg<use_bias>(iarg, objc, objv);
}
template<colvarscript::Object_type T>
int colvarscript::check_cmd_nargs(char const *cmd,
int objc,
int n_args_min,
int n_args_max)
{
int const shift = cmd_arg_shift<T>();
if (objc < shift+n_args_min) {
add_error_msg("Missing arguments for script function \""+std::string(cmd)+
"\":\n"+get_command_help(cmd));
return COLVARSCRIPT_ERROR;
)
CVSCRIPT(cv_getconfig,
"Get the module's configuration string read so far",
0, 0,
{ },
script->set_str_result(cvm::main()->get_config());
return COLVARS_OK;
)
CVSCRIPT(cv_resetindexgroups,
"Clear the index groups loaded so far, allowing to replace them",
0, 0,
{ },
return cvm::main()->reset_index_groups();
)
CVSCRIPT(cv_addenergy,
"Add an energy to the MD engine",
1, 1,
{ "E (float) - Amount of energy to add" },
cvm::main()->total_bias_energy +=
strtod(script->obj_to_str(objv[2]), NULL);
return COLVARS_OK;
)
CVSCRIPT(cv_getenergy,
"Get the current Colvars energy",
1, 1,
{ "E (float) - Store the energy in this variable" },
double *energy = reinterpret_cast<double *>(objv[2]);
*energy = cvm::main()->total_bias_energy;
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);
}
)
if (objc > shift+n_args_max) {
add_error_msg("Too many arguments for script function \""+std::string(cmd)+
"\":\n"+get_command_help(cmd));
return COLVARSCRIPT_ERROR;
}
return COLVARSCRIPT_OK;
}
#ifndef COLVARSCRIPT_INIT_FN
#ifdef __cplusplus
} // extern "C"
#endif
#endif
#undef CVSCRIPT
inline int colvarscript::check_module_cmd_nargs(char const *cmd,
int objc,
int n_args_min,
int n_args_max)
{
return check_cmd_nargs<use_module>(cmd, objc, n_args_min, n_args_max);
}
inline int colvarscript::check_colvar_cmd_nargs(char const *cmd,
int objc,
int n_args_min,
int n_args_max)
{
return check_cmd_nargs<use_colvar>(cmd, objc, n_args_min, n_args_max);
}
inline int colvarscript::check_bias_cmd_nargs(char const *cmd,
int objc,
int n_args_min,
int n_args_max)
{
return check_cmd_nargs<use_bias>(cmd, objc, n_args_min, n_args_max);
}
template<colvarscript::Object_type T>
int colvarscript::cmd_arg_shift()
{
int shift = 0;
if (T == use_module) {
// "cv" and "COMMAND" are 1st and 2nd argument, and shift is equal to 2
shift = 2;
} else if (T == use_colvar) {
// Same as above with additional arguments "colvar" and "NAME"
shift = 4;
} else if (T == use_bias) {
shift = 4;
}
return shift;
}
extern "C" {
#if defined(COLVARS_TCL)
/// Run the script API via Tcl command-line interface
/// \param clientData Not used
/// \param my_interp Pointer to Tcl_Interp object (read from Colvars if NULL)
/// \param objc Number of Tcl command parameters
/// \param objv Array of command parameters
/// \return Result of the script command
int tcl_run_colvarscript_command(ClientData clientData,
Tcl_Interp *interp_in,
int objc, Tcl_Obj *const objv[]);
#endif // #if defined(COLVARS_TCL)
/// Generic wrapper for string-based scripting
int run_colvarscript_command(int objc, unsigned char *const objv[]);
/// Get the string result of a script call
const char * get_colvarscript_result();
}
#endif // #ifndef COLVARSCRIPT_H

65
lib/colvars/colvarscript_commands.cpp Normal file
View File

@ -0,0 +1,65 @@
// -*- 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 <vector>
#include <cstdlib>
#include <string.h>
#include "colvarproxy.h"
#include "colvardeps.h"
#include "colvarscript.h"
#include "colvarscript_commands.h"
extern "C"
int cvscript_n_commands()
{
return static_cast<int>(colvarscript::cv_n_commands);
}
extern "C"
char const **cvscript_command_names()
{
colvarscript *script = colvarscript_obj();
return script->get_command_names();
}
// Instantiate the body of all script commands
#define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \
int objc, unsigned char *const objv[]) \
{ \
if (cvm::debug()) { \
cvm::log("Executing script function \""+std::string(#COMM)+"\""); \
} \
colvarscript *script = colvarscript_obj(); \
script->clear_str_result(); \
if (script->check_module_cmd_nargs(#COMM, \
objc, N_ARGS_MIN, N_ARGS_MAX) != \
COLVARSCRIPT_OK) { \
return COLVARSCRIPT_ERROR; \
} \
FN_BODY; \
}
#undef CVSCRIPT
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)
// Skips the colvar- and bias- specific commands
#define COLVARSCRIPT_COMMANDS_GLOBAL
#undef COLVARSCRIPT_COMMANDS_H
#include "colvarscript_commands.h"
#undef CVSCRIPT_COMM_FN
#undef CVSCRIPT

407
lib/colvars/colvarscript_commands.h Normal file
View File

@ -0,0 +1,407 @@
// -*- 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 COLVARSCRIPT_COMMANDS_H
#define COLVARSCRIPT_COMMANDS_H
// The following is a complete definition of the scripting API.
// The CVSCRIPT macro is used in four distinct contexts:
// 1) Expand to the functions' prototypes (when included generically)
// 2) List colvarscript::command entries (when included in colvarscript.h)
// 3) Implement colvarscript::init() (when included in colvarscript.cpp)
// 4) Define the functions' bodies (when included in colvarscript_commands.cpp)
// Each command is created by an instance of the CVSCRIPT macro
// The arguments of the CVSCRIPT macro are:
// COMM = the id of the command (must be a member of colvarscript::command)
// HELP = a one-line description (C string literal) for the command
// N_ARGS_MIN = the lowest number of arguments allowed
// N_ARGS_MAX = the highest number of arguments allowed
// ARGS = multi-line string literal describing each parameter; each line
// follows the format "name : type - description"
// FN_BODY = the implementation of the function; this should be a thin wrapper
// over existing functions; the "script" pointer to the colvarscript
// object is already set by the CVSCRIPT_COMM_FN macro; see also the
// functions in colvarscript_commands.h.
#ifndef CVSCRIPT_COMM_FNAME
#define CVSCRIPT_COMM_FNAME(COMM) cvscript_ ## COMM
#endif
// If CVSCRIPT is not defined, this file yields the function prototypes
#ifndef CVSCRIPT
#ifdef __cplusplus
#define CVSCRIPT_COMM_PROTO(COMM) \
extern "C" int CVSCRIPT_COMM_FNAME(COMM)(void *, \
int, unsigned char *const *);
#else
#define CVSCRIPT_COMM_PROTO(COMM) \
int CVSCRIPT_COMM_FNAME(COMM)(void *, int, unsigned char *const *);
#endif
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_PROTO(COMM)
// Utility functions used to query the command database
extern "C" {
/// Get the number of colvarscript commands
int cvscript_n_commands();
/// Get the names of all commands (array of strings)
char const ** cvscript_command_names();
}
#endif
CVSCRIPT(cv_addenergy,
"Add an energy to the MD engine (no effect in VMD)",
1, 1,
"E : float - Amount of energy to add",
char const *Earg =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
cvm::main()->total_bias_energy += strtod(Earg, NULL);
return cvm::get_error(); // TODO Make this multi-language
)
CVSCRIPT(cv_bias,
"Prefix for bias-specific commands",
0, 0,
"",
// This cannot be executed from a command line
return COLVARS_OK;
)
CVSCRIPT(cv_colvar,
"Prefix for colvar-specific commands",
0, 0,
"",
// This cannot be executed from a command line
return COLVARS_OK;
)
CVSCRIPT(cv_config,
"Read configuration from the given string",
1, 1,
"conf : string - Configuration string",
char const *conf_str =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
std::string const conf(conf_str);
if (cvm::main()->read_config_string(conf) == COLVARS_OK) {
return COLVARS_OK;
}
script->add_error_msg("Error parsing configuration string");
return COLVARSCRIPT_ERROR;
)
CVSCRIPT(cv_configfile,
"Read configuration from a file",
1, 1,
"conf_file : string - Path to configuration file",
char const *conf_file_name =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
if (script->module()->read_config_file(conf_file_name) == COLVARS_OK) {
return COLVARS_OK;
} else {
script->add_error_msg("Error parsing configuration file");
return COLVARSCRIPT_ERROR;
}
)
CVSCRIPT(cv_delete,
"Delete this Colvars module instance (VMD only)",
0, 0,
"",
return script->proxy()->request_deletion();
)
CVSCRIPT(cv_frame,
"Get or set current frame number (VMD only)",
0, 1,
"frame : integer - Frame number",
char const *arg =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
if (arg == NULL) {
long int f = -1;
if (script->proxy()->get_frame(f) == COLVARS_OK) {
script->set_result_str(cvm::to_str(f));
return COLVARS_OK;
} else {
script->add_error_msg("Frame number is not available");
return COLVARSCRIPT_ERROR;
}
} else {
int const f = strtol(const_cast<char *>(arg), NULL, 10);
int error_code = script->proxy()->set_frame(f);
if (error_code == COLVARS_NO_SUCH_FRAME) {
script->add_error_msg("Invalid frame number: \""+std::string(arg)+
"\"\n");
}
return error_code;
}
return COLVARS_OK;
)
CVSCRIPT(cv_getconfig,
"Get the module's configuration string read so far",
0, 0,
"",
script->set_result_str(cvm::main()->get_config());
return COLVARS_OK;
)
CVSCRIPT(cv_getenergy,
"Get the current Colvars energy",
0, 0,
"",
script->set_result_str(cvm::to_str(cvm::main()->total_bias_energy));
return COLVARS_OK;
)
CVSCRIPT(cv_help,
"Get the help string of the Colvars scripting interface",
0, 1,
"command : string - Get the help string of this specific command",
unsigned char *const cmdobj =
script->get_module_cmd_arg(0, objc, objv);
if (cmdobj) {
std::string const cmdstr(script->obj_to_str(cmdobj));
if (cmdstr.size()) {
if (cmdstr == std::string("colvar")) {
script->set_result_str(script->get_cmdline_help_summary(colvarscript::use_colvar));
} else if (cmdstr == std::string("bias")) {
script->set_result_str(script->get_cmdline_help_summary(colvarscript::use_bias));
} else {
script->set_result_str(script->get_command_cmdline_help(colvarscript::use_module,
cmdstr));
}
return cvm::get_error();
} else {
return COLVARSCRIPT_ERROR;
}
} else {
script->set_result_str(script->get_cmdline_help_summary(colvarscript::use_module));
return COLVARS_OK;
}
)
CVSCRIPT(cv_list,
"Return a list of all variables or biases",
0, 1,
"param : string - \"colvars\" or \"biases\"; default is \"colvars\"",
std::string res;
unsigned char *const kwarg = script->get_module_cmd_arg(0, objc, objv);
std::string const kwstr = kwarg ? script->obj_to_str(kwarg) :
std::string("colvars");
if (kwstr == "colvars") {
for (std::vector<colvar *>::iterator cvi = script->module()->variables()->begin();
cvi != script->module()->variables()->end();
++cvi) {
res += (cvi == script->module()->variables()->begin() ? "" : " ") + (*cvi)->name;
}
script->set_result_str(res);
return COLVARS_OK;
} else if (kwstr == "biases") {
for (std::vector<colvarbias *>::iterator bi = script->module()->biases.begin();
bi != script->module()->biases.end();
++bi) {
res += (bi == script->module()->biases.begin() ? "" : " ") + (*bi)->name;
}
script->set_result_str(res);
return COLVARS_OK;
} else {
script->add_error_msg("Wrong arguments to command \"list\"\n");
return COLVARSCRIPT_ERROR;
}
)
CVSCRIPT(cv_listcommands,
"Get the list of script functions, prefixed with \"cv_\", \"colvar_\" or \"bias_\"",
0, 0,
"",
int const n_commands = cvscript_n_commands();
char const **command_names = cvscript_command_names();
std::string result;
for (int i = 0; i < n_commands; i++) {
if (i > 0) result.append(1, ' ');
result.append(std::string(command_names[i]));
}
script->set_result_str(result);
return COLVARS_OK;
)
CVSCRIPT(cv_load,
"Load data from a state file into all matching colvars and biases",
1, 1,
"prefix : string - Path to existing state file or input prefix",
char const *arg =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
script->proxy()->input_prefix() = cvm::state_file_prefix(arg);
if (script->module()->setup_input() == COLVARS_OK) {
return COLVARS_OK;
} else {
script->add_error_msg("Error loading state file");
return COLVARSCRIPT_ERROR;
}
)
CVSCRIPT(cv_loadfromstring,
"Load state data from a string into all matching colvars and biases",
1, 1,
"buffer : string - String buffer containing the state information",
char const *arg =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
script->proxy()->input_buffer() = arg;
if (script->module()->setup_input() == COLVARS_OK) {
return COLVARS_OK;
} else {
script->add_error_msg("Error loading state string");
return COLVARSCRIPT_ERROR;
}
)
CVSCRIPT(cv_molid,
"Get or set the molecule ID on which Colvars is defined (VMD only)",
0, 1,
"molid : integer - Molecule ID; -1 means undefined",
char const *arg =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
if (arg == NULL) {
int molid = -1;
script->proxy()->get_molid(molid);
script->set_result_str(cvm::to_str(molid));
return COLVARS_OK;
} else {
script->add_error_msg("Error: To change the molecule ID in VMD, use cv delete first.");
return COLVARS_NOT_IMPLEMENTED;
}
)
CVSCRIPT(cv_printframe,
"Return the values that would be written to colvars.traj",
0, 0,
"",
std::ostringstream os;
script->module()->write_traj(os);
script->set_result_str(os.str());
return COLVARS_OK;
)
CVSCRIPT(cv_printframelabels,
"Return the labels that would be written to colvars.traj",
0, 0,
"",
std::ostringstream os;
script->module()->write_traj_label(os);
script->set_result_str(os.str());
return COLVARS_OK;
)
CVSCRIPT(cv_reset,
"Delete all internal configuration",
0, 0,
"",
return script->module()->reset();
)
CVSCRIPT(cv_resetindexgroups,
"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;
)
CVSCRIPT(cv_save,
"Change the prefix of all output files and save them",
1, 1,
"prefix : string - Output prefix with trailing \".colvars.state\" gets removed)",
std::string const prefix =
cvm::state_file_prefix(script->obj_to_str(script->get_module_cmd_arg(0, objc, objv)));
script->proxy()->output_prefix() = prefix;
int error_code = COLVARS_OK;
error_code |= script->module()->setup_output();
error_code |= script->module()->write_restart_file(prefix+
".colvars.state");
error_code |= script->module()->write_output_files();
return error_code;
)
CVSCRIPT(cv_savetostring,
"Write the Colvars state to a string and return it",
0, 0,
"",
return script->module()->write_restart_string(script->modify_str_result());
)
CVSCRIPT(cv_units,
"Get or set the current Colvars unit system",
0, 1,
"units : string - The new unit system",
char const *argstr =
script->obj_to_str(script->get_module_cmd_arg(0, objc, objv));
if (argstr) {
return cvm::proxy->set_unit_system(argstr, false);
} else {
script->set_result_str(cvm::proxy->units);
return COLVARS_OK;
}
)
CVSCRIPT(cv_update,
"Recalculate colvars and biases",
0, 0,
"",
int error_code = script->proxy()->update_input();
if (error_code) {
script->add_error_msg("Error updating the Colvars module (input)");
return error_code;
}
error_code |= script->module()->calc();
if (error_code) {
script->add_error_msg("Error updating the Colvars module (calc)");
return error_code;
}
error_code |= script->proxy()->update_output();
if (error_code) {
script->add_error_msg("Error updating the Colvars module (output)");
}
return error_code;
)
CVSCRIPT(cv_version,
"Get the Colvars Module version number",
0, 0,
"",
script->set_result_str(COLVARS_VERSION);
return COLVARS_OK;
)
// This guard allows compiling colvar and bias function bodies in their
// respecitve files instead of colvarscript_commands.o
#ifndef COLVARSCRIPT_COMMANDS_GLOBAL
#include "colvarscript_commands_colvar.h"
#include "colvarscript_commands_bias.h"
#endif
#endif // #ifndef COLVARSCRIPT_COMMANDS_H

49
lib/colvars/colvarscript_commands_bias.cpp Normal file
View File

@ -0,0 +1,49 @@
// -*- 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 <vector>
#include <cstdlib>
#include <stdlib.h>
#include <string.h>
#include "colvarproxy.h"
#include "colvardeps.h"
#include "colvarscript.h"
#include "colvarscript_commands.h"
// Instantiate the body of all bias-specific script commands
#define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \
int objc, unsigned char *const objv[]) \
{ \
if (cvm::debug()) { \
cvm::log("Executing script function \""+std::string(#COMM)+"\""); \
} \
colvarscript *script = colvarscript_obj(); \
script->clear_str_result(); \
if (script->check_bias_cmd_nargs(#COMM, \
objc, N_ARGS_MIN, N_ARGS_MAX) != \
COLVARSCRIPT_OK) { \
return COLVARSCRIPT_ERROR; \
} \
colvarbias *this_bias = colvarbias_obj(pobj); \
FN_BODY; \
}
#undef CVSCRIPT
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)
#include "colvarscript_commands_bias.h"
#undef CVSCRIPT_COMM_FN
#undef CVSCRIPT

173
lib/colvars/colvarscript_commands_bias.h Normal file
View File

@ -0,0 +1,173 @@
// -*- 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.
CVSCRIPT(bias_bin,
"Get the current grid bin index (1D ABF only for now)",
0, 0,
"",
script->set_result_str(cvm::to_str(this_bias->current_bin()));
return COLVARS_OK;
)
CVSCRIPT(bias_bincount,
"Get the number of samples at the given grid bin (1D ABF only for now)",
0, 1,
"index : integer - Grid index; defaults to current bin",
int index = this_bias->current_bin();
char const *indexarg =
script->obj_to_str(script->get_bias_cmd_arg(0, objc, objv));
if (indexarg) {
std::string const param(indexarg);
if (!(std::istringstream(param) >> index)) {
script->add_error_msg("bincount: error parsing bin index");
return COLVARSCRIPT_ERROR;
}
}
script->set_result_str(cvm::to_str(this_bias->bin_count(index)));
return COLVARS_OK;
)
CVSCRIPT(bias_binnum,
"Get the total number of grid points of this bias (1D ABF only for now)",
0, 0,
"",
int r = this_bias->bin_num();
if (r < 0) {
script->add_error_msg("Error: calling bin_num() for bias " +
this_bias->name);
return COLVARSCRIPT_ERROR;
}
script->set_result_str(cvm::to_str(r));
return COLVARS_OK;
)
CVSCRIPT(bias_delete,
"Delete this bias",
0, 0,
"",
delete this_bias;
return COLVARS_OK;
)
CVSCRIPT(bias_energy,
"Get the current energy of this bias",
0, 0,
"",
script->set_result_str(cvm::to_str(this_bias->get_energy()));
return COLVARS_OK;
)
CVSCRIPT(bias_get,
"Get the value of the given feature for this bias",
1, 1,
"feature : string - Name of the feature",
return script->proc_features(this_bias, objc, objv);
)
CVSCRIPT(bias_getconfig,
"Return the configuration string of this bias",
0, 0,
"",
script->set_result_str(this_bias->get_config());
return COLVARS_OK;
)
CVSCRIPT(bias_help,
"Get a help summary or the help string of one bias subcommand",
0, 1,
"command : string - Get the help string of this specific command",
unsigned char *const cmdobj =
script->get_colvar_cmd_arg(0, objc, objv);
if (this_bias) {
}
if (cmdobj) {
std::string const cmdstr(script->obj_to_str(cmdobj));
if (cmdstr.size()) {
script->set_result_str(script->get_command_cmdline_help(colvarscript::use_bias,
cmdstr));
return COLVARS_OK;
} else {
return COLVARSCRIPT_ERROR;
}
} else {
script->set_result_str(script->get_cmdline_help_summary(colvarscript::use_bias));
return COLVARS_OK;
}
)
CVSCRIPT(bias_load,
"Load data into this bias",
1, 1,
"prefix : string - Read from a file with this name or prefix",
char const *arg =
script->obj_to_str(script->get_bias_cmd_arg(0, objc, objv));
return this_bias->read_state_prefix(std::string(arg));
)
CVSCRIPT(bias_loadfromstring,
"Load state data into this bias from a string",
1, 1,
"buffer : string - String buffer containing the state information",
char const *buffer = script->obj_to_str(script->get_bias_cmd_arg(0, objc, objv));
return this_bias->read_state_string(buffer);
)
CVSCRIPT(bias_save,
"Save data from this bias into a file with the given prefix",
1, 1,
"prefix : string - Prefix for the state file of this bias",
std::string const prefix =
cvm::state_file_prefix(script->obj_to_str(script->get_bias_cmd_arg(0, objc, objv)));
return this_bias->write_state_prefix(prefix);
)
CVSCRIPT(bias_savetostring,
"Save data from this bias into a string and return it",
0, 0,
"",
return this_bias->write_state_string(script->modify_str_result());
)
CVSCRIPT(bias_set,
"Set the given feature of this bias to a new value",
2, 2,
"feature : string - Name of the feature\n"
"value : string - String representation of the new feature value",
return script->proc_features(this_bias, objc, objv);
)
CVSCRIPT(bias_share,
"Share bias information with other replicas (multiple-walker scheme)",
0, 0,
"",
if (this_bias->replica_share() != COLVARS_OK) {
script->add_error_msg("Error: calling replica_share() for bias " +
this_bias->name);
return COLVARSCRIPT_ERROR;
}
return COLVARS_OK;
)
CVSCRIPT(bias_state,
"Print a string representation of the feature state of this bias",
0, 0,
"",
this_bias->print_state();
return COLVARS_OK;
)
CVSCRIPT(bias_update,
"Recompute this bias and return its up-to-date energy",
0, 0,
"",
this_bias->update();
script->set_result_str(cvm::to_str(this_bias->get_energy()));
return COLVARS_OK;
)

49
lib/colvars/colvarscript_commands_colvar.cpp Normal file
View File

@ -0,0 +1,49 @@
// -*- 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 <vector>
#include <cstdlib>
#include <stdlib.h>
#include <string.h>
#include "colvarproxy.h"
#include "colvardeps.h"
#include "colvarscript.h"
#include "colvarscript_commands.h"
// Instantiate the body of all colvar-specific script commands
#define CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
int CVSCRIPT_COMM_FNAME(COMM)(void *pobj, \
int objc, unsigned char *const objv[]) \
{ \
if (cvm::debug()) { \
cvm::log("Executing script function \""+std::string(#COMM)+"\""); \
} \
colvarscript *script = colvarscript_obj(); \
script->clear_str_result(); \
if (script->check_colvar_cmd_nargs(#COMM, \
objc, N_ARGS_MIN, N_ARGS_MAX) != \
COLVARSCRIPT_OK) { \
return COLVARSCRIPT_ERROR; \
} \
colvar *this_colvar = colvar_obj(pobj); \
FN_BODY; \
}
#undef CVSCRIPT
#define CVSCRIPT(COMM,HELP,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY) \
CVSCRIPT_COMM_FN(COMM,N_ARGS_MIN,N_ARGS_MAX,ARGS,FN_BODY)
#include "colvarscript_commands_colvar.h"
#undef CVSCRIPT_COMM_FN
#undef CVSCRIPT

239
lib/colvars/colvarscript_commands_colvar.h Normal file
View File

@ -0,0 +1,239 @@
// -*- 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.
CVSCRIPT(colvar_addforce,
"Apply the given force onto this colvar and return the same",
1, 1,
"force : float or array - Applied force; must match colvar dimensionality",
std::string const f_str(script->obj_to_str(script->get_colvar_cmd_arg(0, objc, objv)));
std::istringstream is(f_str);
is.width(cvm::cv_width);
is.precision(cvm::cv_prec);
colvarvalue force(this_colvar->value());
force.is_derivative();
if (force.from_simple_string(is.str()) != COLVARS_OK) {
script->add_error_msg("addforce : error parsing force value");
return COLVARSCRIPT_ERROR;
}
this_colvar->add_bias_force(force);
script->set_result_str(force.to_simple_string());
return COLVARS_OK;
)
CVSCRIPT(colvar_cvcflags,
"Enable or disable individual components by setting their active flags",
1, 1,
"flags : integer array - Zero/nonzero value disables/enables the CVC",
std::string const flags_str(script->obj_to_str(script->get_colvar_cmd_arg(0, objc, objv)));
std::istringstream is(flags_str);
std::vector<bool> flags;
int flag;
while (is >> flag) {
flags.push_back(flag != 0);
}
int res = this_colvar->set_cvc_flags(flags);
if (res != COLVARS_OK) {
script->add_error_msg("Error setting CVC flags");
return COLVARSCRIPT_ERROR;
}
script->set_result_str("0");
return COLVARS_OK;
)
CVSCRIPT(colvar_delete,
"Delete this colvar, along with all biases that depend on it",
0, 0,
"",
delete this_colvar;
return COLVARS_OK;
)
CVSCRIPT(colvar_get,
"Get the value of the given feature for this colvar",
1, 1,
"feature : string - Name of the feature",
return script->proc_features(this_colvar, objc, objv);
)
CVSCRIPT(colvar_getappliedforce,
"Return the total of the forces applied to this colvar",
0, 0,
"",
script->set_result_str((this_colvar->applied_force()).to_simple_string());
return COLVARS_OK;
)
CVSCRIPT(colvar_getatomgroups,
"Return the atom indices used by this colvar as a list of lists",
0, 0,
"",
std::string result;
std::vector<std::vector<int> > lists = this_colvar->get_atom_lists();
std::vector<std::vector<int> >::iterator li = lists.begin();
for ( ; li != lists.end(); ++li) {
result += "{";
std::vector<int>::iterator lj = (*li).begin();
for ( ; lj != (*li).end(); ++lj) {
result += cvm::to_str(*lj);
result += " ";
}
result += "} ";
}
script->set_result_str(result);
return COLVARS_OK;
)
CVSCRIPT(colvar_getatomids,
"Return the list of atom indices used by this colvar",
0, 0,
"",
std::string result;
std::vector<int>::iterator li = this_colvar->atom_ids.begin();
for ( ; li != this_colvar->atom_ids.end(); ++li) {
result += cvm::to_str(*li);
result += " ";
}
script->set_result_str(result);
return COLVARS_OK;
)
CVSCRIPT(colvar_getconfig,
"Return the configuration string of this colvar",
0, 0,
"",
script->set_result_str(this_colvar->get_config());
return COLVARS_OK;
)
CVSCRIPT(colvar_getgradients,
"Return the atomic gradients of this colvar",
0, 0,
"",
std::string result;
std::vector<cvm::rvector>::iterator li =
this_colvar->atomic_gradients.begin();
for ( ; li != this_colvar->atomic_gradients.end(); ++li) {
result += "{";
int j;
for (j = 0; j < 3; ++j) {
result += cvm::to_str((*li)[j]);
result += " ";
}
result += "} ";
}
script->set_result_str(result);
return COLVARS_OK;
)
CVSCRIPT(colvar_gettotalforce,
"Return the sum of internal and external forces to this colvar",
0, 0,
"",
script->set_result_str((this_colvar->total_force()).to_simple_string());
return COLVARS_OK;
)
CVSCRIPT(colvar_help,
"Get a help summary or the help string of one colvar subcommand",
0, 1,
"command : string - Get the help string of this specific command",
unsigned char *const cmdobj =
script->get_colvar_cmd_arg(0, objc, objv);
if (this_colvar) {
}
if (cmdobj) {
std::string const cmdstr(script->obj_to_str(cmdobj));
if (cmdstr.size()) {
script->set_result_str(script->get_command_cmdline_help(colvarscript::use_colvar,
cmdstr));
return cvm::get_error();
} else {
return COLVARSCRIPT_ERROR;
}
} else {
script->set_result_str(script->get_cmdline_help_summary(colvarscript::use_colvar));
return COLVARS_OK;
}
)
CVSCRIPT(colvar_modifycvcs,
"Modify configuration of individual components by passing string arguments",
1, 1,
"confs : sequence of strings - New configurations; empty strings are skipped",
std::vector<std::string> const confs(script->proxy()->script_obj_to_str_vector(script->get_colvar_cmd_arg(0, objc, objv)));
cvm::increase_depth();
int res = this_colvar->update_cvc_config(confs);
cvm::decrease_depth();
if (res != COLVARS_OK) {
script->add_error_msg("Error setting CVC flags");
return COLVARSCRIPT_ERROR;
}
script->set_result_str("0");
return COLVARS_OK;
)
CVSCRIPT(colvar_run_ave,
"Get the current running average of the value of this colvar",
0, 0,
"",
script->set_result_str(this_colvar->run_ave().to_simple_string());
return COLVARS_OK;
)
CVSCRIPT(colvar_set,
"Set the given feature of this colvar to a new value",
2, 2,
"feature : string - Name of the feature\n"
"value : string - String representation of the new feature value",
return script->proc_features(this_colvar, objc, objv);
)
CVSCRIPT(colvar_state,
"Print a string representation of the feature state of this colvar",
0, 0,
"",
this_colvar->print_state();
return COLVARS_OK;
)
CVSCRIPT(colvar_type,
"Get the type description of this colvar",
0, 0,
"",
script->set_result_str(this_colvar->value().type_desc(this_colvar->value().value_type));
return COLVARS_OK;
)
CVSCRIPT(colvar_update,
"Recompute this colvar and return its up-to-date value",
0, 0,
"",
this_colvar->calc();
this_colvar->update_forces_energy();
script->set_result_str((this_colvar->value()).to_simple_string());
return COLVARS_OK;
)
CVSCRIPT(colvar_value,
"Get the current value of this colvar",
0, 0,
"",
script->set_result_str(this_colvar->value().to_simple_string());
return COLVARS_OK;
)
CVSCRIPT(colvar_width,
"Get the width of this colvar",
0, 0,
"",
script->set_result_str(cvm::to_str(this_colvar->width, 0,
cvm::cv_prec));
return COLVARS_OK;
)

4
lib/colvars/colvartypes.cpp
View File

@ -74,7 +74,7 @@ std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v)
std::istream & operator >> (std::istream &is, colvarmodule::rvector &v)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
char sep;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> v.x) || !(is >> sep) || !(sep == ',') ||
@ -130,7 +130,7 @@ std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q)
std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
{
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
std::string euler("");

2
lib/colvars/colvartypes.h
View File

@ -270,7 +270,7 @@ public:
cvm::vector1d<T> &v)
{
if (v.size() == 0) return is;
size_t const start_pos = is.tellg();
std::streampos const start_pos = is.tellg();
char sep;
if ( !(is >> sep) || !(sep == '(') ) {
is.clear();

32
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -71,11 +71,9 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
_random = new LAMMPS_NS::RanPark(lmp,seed);
first_timestep=true;
total_force_requested=false;
previous_step=-1;
t_target=temp;
do_exit=false;
restart_every=0;
// User-scripted forces are not available in LAMMPS
force_script_defined = false;
@ -94,15 +92,17 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
// check if it is possible to save output configuration
if ((!output_prefix_str.size()) && (!restart_output_prefix_str.size())) {
fatal_error("Error: neither the final output state file or "
error("Error: neither the final output state file or "
"the output restart file could be defined, exiting.\n");
}
// try to extract a restart prefix from a potential restart command.
LAMMPS_NS::Output *outp = _lmp->output;
if ((outp->restart_every_single > 0) && (outp->restart1 != 0)) {
restart_frequency_engine = outp->restart_every_single;
restart_output_prefix_str = std::string(outp->restart1);
} else if ((outp->restart_every_double > 0) && (outp->restart2a != 0)) {
restart_frequency_engine = outp->restart_every_double;
restart_output_prefix_str = std::string(outp->restart2a);
}
// trim off unwanted stuff from the restart prefix
@ -207,12 +207,18 @@ double colvarproxy_lammps::compute()
first_timestep = false;
} else {
// Use the time step number from LAMMPS Update object
if ( _lmp->update->ntimestep - previous_step == 1 )
if ( _lmp->update->ntimestep - previous_step == 1 ) {
colvars->it++;
// Other cases could mean:
b_simulation_continuing = false;
} else {
// Cases covered by this condition:
// - run 0
// - beginning of a new run statement
// then the internal counter should not be incremented
// The internal counter is not incremented, and the objects are made
// aware of this via the following flag
b_simulation_continuing = true;
}
}
previous_step = _lmp->update->ntimestep;
@ -275,13 +281,6 @@ void colvarproxy_lammps::serialize_status(std::string &rst)
rst = os.str();
}
void colvarproxy_lammps::write_output_files()
{
// TODO skip output if undefined
colvars->write_restart_file(cvm::output_prefix()+".colvars.state");
colvars->write_output_files();
}
// set status from string
bool colvarproxy_lammps::deserialize_status(std::string &rst)
{
@ -322,13 +321,6 @@ void colvarproxy_lammps::log(std::string const &message)
void colvarproxy_lammps::error(std::string const &message)
{
// In LAMMPS, all errors are fatal
fatal_error(message);
}
void colvarproxy_lammps::fatal_error(std::string const &message)
{
log(message);
_lmp->error->one(FLERR,

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

@ -41,11 +41,9 @@ class colvarproxy_lammps : public colvarproxy {
// state of LAMMPS properties
double t_target, my_timestep, my_boltzmann, my_angstrom;
double bias_energy;
int restart_every;
int previous_step;
bool first_timestep;
bool total_force_requested;
bool do_exit;
std::vector<int> atoms_types;
@ -59,7 +57,7 @@ class colvarproxy_lammps : public colvarproxy {
const char *, const int, const double, MPI_Comm);
virtual ~colvarproxy_lammps();
void init(const char*);
int setup();
virtual int setup();
// disable default and copy constructor
private:
@ -82,9 +80,6 @@ class colvarproxy_lammps : public colvarproxy {
// set status from string
bool deserialize_status(std::string &);
// Write files expected from Colvars (called by post_run())
void write_output_files();
// read additional config from file
int add_config_file(char const *config_filename);
@ -103,14 +98,11 @@ class colvarproxy_lammps : public colvarproxy {
inline cvm::real temperature() { return t_target; };
inline cvm::real dt() { return my_timestep; }; // return _lmp->update->dt * _lmp->force->femtosecond; };
inline size_t restart_frequency() { return restart_every; };
void add_energy(cvm::real energy) { bias_energy += energy; };
void request_total_force(bool yesno) { total_force_requested = yesno; };
void log(std::string const &message);
void error(std::string const &message);
void fatal_error(std::string const &message);
cvm::rvector position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2) const;

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

@ -1,3 +1,3 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2019-12-04"
#define COLVARPROXY_VERSION "2020-04-07"
#endif

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

@ -499,7 +499,7 @@ int FixColvars::modify_param(int narg, char **arg)
if (me == 0) {
if (! proxy)
error->one(FLERR,"Cannot use fix_modify before initialization");
return proxy->add_config_file(arg[1]) == COLVARS_OK ? 0 : 2;
return proxy->add_config_file(arg[1]) == COLVARS_OK ? 2 : 0;
}
return 2;
} else if (strcmp(arg[0],"config") == 0) {
@ -508,7 +508,7 @@ int FixColvars::modify_param(int narg, char **arg)
if (! proxy)
error->one(FLERR,"Cannot use fix_modify before initialization");
std::string const conf(arg[1]);
return proxy->add_config_string(conf) == COLVARS_OK ? 0 : 2;
return proxy->add_config_string(conf) == COLVARS_OK ? 2 : 0;
}
return 2;
} else if (strcmp(arg[0],"load") == 0) {
@ -516,7 +516,7 @@ int FixColvars::modify_param(int narg, char **arg)
if (me == 0) {
if (! proxy)
error->one(FLERR,"Cannot use fix_modify before initialization");
return proxy->read_state_file(arg[1]) == COLVARS_OK ? 0 : 2;
return proxy->read_state_file(arg[1]) == COLVARS_OK ? 2 : 0;
}
return 2;
}
@ -984,7 +984,7 @@ void FixColvars::restart(char *buf)
void FixColvars::post_run()
{
if (me == 0) {
proxy->write_output_files();
proxy->post_run();
}
}