diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index 59d12f9253..4f85537b35 100644 Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ diff --git a/lib/colvars/Makefile.common b/lib/colvars/Makefile.common index 356a7f4a91..9e7be12958 100644 --- a/lib/colvars/Makefile.common +++ b/lib/colvars/Makefile.common @@ -50,7 +50,9 @@ COLVARS_SRCS = \ colvarparams.cpp \ colvarparse.cpp \ colvarproxy.cpp \ + colvarproxy_io.cpp \ colvarproxy_replicas.cpp \ + colvarproxy_system.cpp \ colvarproxy_tcl.cpp \ colvarproxy_volmaps.cpp \ colvarscript.cpp \ diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps index 6619653af0..6bf5171710 100644 --- a/lib/colvars/Makefile.deps +++ b/lib/colvars/Makefile.deps @@ -1,371 +1,140 @@ $(COLVARS_OBJ_DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvarparse.h colvarparams.h \ - colvaratoms.h colvardeps.h + colvarproxy_io.h colvarproxy_system.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 colvar.h colvarvalue.h colvartypes.h colvarparse.h \ - colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarbias_abf.h colvarproxy.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvarbias.h colvargrid.h \ - colvar_UIestimator.h + colvarparams.h colvardeps.h colvarbias_abf.h colvarproxy.h \ + colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarbias_alb.h + colvarparse.h colvarparams.h colvardeps.h colvarbias_alb.h $(COLVARS_OBJ_DIR)colvarbias.o: colvarbias.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvarbias.h colvar.h \ - colvarparse.h colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvargrid.h + colvarproxy_io.h colvarproxy_system.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 colvarproxy_tcl.h colvarproxy_volmaps.h colvar.h \ - colvarparse.h colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarbias_histogram.h colvarbias.h \ - colvargrid.h + colvarvalue.h colvarproxy_io.h colvarproxy_system.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_histogram_reweight_amd.o: \ colvarbias_histogram_reweight_amd.cpp \ colvarbias_histogram_reweight_amd.h colvarbias_histogram.h colvarbias.h \ colvar.h colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvargrid.h colvarproxy.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h + colvarparse.h colvarparams.h colvardeps.h colvargrid.h colvarproxy.h \ + colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h $(COLVARS_OBJ_DIR)colvarbias_meta.o: colvarbias_meta.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvar.h colvarparse.h \ - colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarbias_meta.h colvarbias.h \ - colvargrid.h + colvarproxy_io.h colvarproxy_system.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 colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvarbias_restraint.h colvarbias.h colvar.h colvarparse.h \ - colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h + colvarvalue.h colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvarbias_restraint.h colvarbias.h colvar.h \ + colvarparse.h colvarparams.h colvardeps.h $(COLVARS_OBJ_DIR)colvarcomp_alchlambda.o: colvarcomp_alchlambda.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvar.h colvardeps.h \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvarproxy_io.h colvarproxy_system.h \ + colvarproxy_tcl.h colvarproxy_volmaps.h colvar_arithmeticpath.h \ + colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \ colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \ - colvarparse.h colvarparams.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvarproxy_io.h colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvarproxy_io.h colvarproxy_system.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 \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvar.h \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvar_arithmeticpath.h \ - colvar_geometricpath.h + colvarproxy_io.h colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \ + colvarproxy_io.h colvarproxy_system.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 \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvarproxy_io.h colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_neuralnetwork.o: \ colvarcomp_neuralnetwork.cpp colvarmodule.h colvars_version.h \ colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \ - colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ + colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_io.h \ + colvarproxy_system.h colvarproxy_tcl.h colvarproxy_volmaps.h \ colvar_arithmeticpath.h colvar_geometricpath.h \ colvar_neuralnetworkcompute.h $(COLVARS_OBJ_DIR)colvarcomp_combination.o: colvarcomp_combination.cpp \ colvarcomp.h colvarmodule.h colvars_version.h colvar.h colvarvalue.h \ - colvartypes.h colvarparse.h colvarparams.h colvardeps.h \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvaratoms.h colvarproxy.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvar_arithmeticpath.h \ - colvar_geometricpath.h + colvartypes.h colvarparse.h colvarparams.h colvardeps.h colvaratoms.h \ + colvarproxy.h colvarproxy_io.h colvarproxy_system.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 \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvarproxy_io.h colvarproxy_system.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 \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvarproxy_io.h colvarproxy_system.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 \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h + colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ + colvaratoms.h colvarproxy.h colvarproxy_io.h colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ + colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h colvarproxy_io.h \ + colvarproxy_system.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 \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \ - colvarparams.h + colvarproxy_io.h colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar_arithmeticpath.h colvar_geometricpath.h colvargrid.h + colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvar_arithmeticpath.h colvar_geometricpath.h \ + colvargrid.h colvargrid_def.h $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \ colvars_version.h colvarparse.h colvarvalue.h colvartypes.h \ - colvarparams.h colvarproxy.h colvarproxy_tcl.h colvarproxy_volmaps.h \ - colvar.h colvardeps.h ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarbias.h colvarbias_abf.h \ - colvargrid.h colvar_UIestimator.h colvarbias_alb.h \ - colvarbias_histogram.h colvarbias_histogram_reweight_amd.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 colvarmodule_refs.h + colvarparams.h colvarproxy.h colvarproxy_io.h colvarproxy_system.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_histogram_reweight_amd.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 colvarmodule_refs.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 \ @@ -373,97 +142,53 @@ $(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 \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvarscript.h colvarbias.h \ - colvar.h colvarparse.h colvarparams.h colvardeps.h \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarscript_commands.h \ + colvarproxy_io.h colvarproxy_system.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 colvarmodule_utils.h +$(COLVARS_OBJ_DIR)colvarproxy_io.o: colvarproxy_io.cpp colvarmodule.h \ + colvars_version.h colvarproxy_io.h $(COLVARS_OBJ_DIR)colvarproxy_replicas.o: colvarproxy_replicas.cpp \ colvarmodule.h colvars_version.h colvarproxy.h colvartypes.h \ - colvarvalue.h colvarproxy_tcl.h colvarproxy_volmaps.h + colvarvalue.h colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h +$(COLVARS_OBJ_DIR)colvarproxy_system.o: colvarproxy_system.cpp \ + colvarmodule.h colvars_version.h colvartypes.h colvarproxy_system.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 + colvarproxy_io.h colvarproxy_system.h colvarproxy_tcl.h \ + colvarproxy_volmaps.h colvaratoms.h colvarparse.h colvarparams.h \ + colvardeps.h $(COLVARS_OBJ_DIR)colvarproxy_volmaps.o: colvarproxy_volmaps.cpp \ colvarmodule.h colvars_version.h colvarproxy_volmaps.h \ colvarmodule_utils.h $(COLVARS_OBJ_DIR)colvarscript.o: colvarscript.cpp colvarproxy.h \ colvarmodule.h colvars_version.h colvartypes.h colvarvalue.h \ - colvarproxy_tcl.h colvarproxy_volmaps.h colvardeps.h colvarparse.h \ - colvarparams.h colvarscript.h colvarbias.h colvar.h \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarscript_commands.h \ + colvarproxy_io.h colvarproxy_system.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 \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarscript_commands.h \ + colvarvalue.h colvarproxy_io.h colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarscript_commands.h \ - colvarscript_commands_colvar.h colvarscript_commands_bias.h + colvars_version.h colvartypes.h colvarvalue.h colvarproxy_io.h \ + colvarproxy_system.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 ../lepton/include/Lepton.h \ - ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarscript_commands.h \ - colvarscript_commands_colvar.h colvarscript_commands_bias.h + colvars_version.h colvartypes.h colvarvalue.h colvarproxy_io.h \ + colvarproxy_system.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 ../../src/math_eigen_impl.h @@ -471,15 +196,6 @@ $(COLVARS_OBJ_DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h $(COLVARS_OBJ_DIR)colvar_neuralnetworkcompute.o: \ colvar_neuralnetworkcompute.cpp colvar_neuralnetworkcompute.h \ - ../lepton/include/Lepton.h ../lepton/include/lepton/CompiledExpression.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/windowsIncludes.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/ExpressionProgram.h \ - ../lepton/include/lepton/ExpressionTreeNode.h \ - ../lepton/include/lepton/Operation.h \ - ../lepton/include/lepton/CustomFunction.h \ - ../lepton/include/lepton/Exception.h \ - ../lepton/include/lepton/ParsedExpression.h \ - ../lepton/include/lepton/Parser.h colvarparse.h colvarmodule.h \ - colvars_version.h colvarvalue.h colvartypes.h colvarparams.h + colvarparse.h colvarmodule.h colvars_version.h colvarvalue.h \ + colvartypes.h colvarparams.h colvarproxy.h colvarproxy_io.h \ + colvarproxy_system.h colvarproxy_tcl.h colvarproxy_volmaps.h diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 85b8443d72..700d3752ac 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include #include "colvarmodule.h" #include "colvarvalue.h" @@ -24,8 +26,6 @@ std::mapproxy; get_keyval_feature(this, conf, "extendedLagrangian", f_cv_extended_Lagrangian, false); if (is_enabled(f_cv_extended_Lagrangian)) { @@ -646,7 +646,8 @@ int colvar::init_extended_Lagrangian(std::string const &conf) x_ext.type(colvarvalue::type_notset); v_ext.type(value()); fr.type(value()); - const bool temp_provided = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); + const bool temp_provided = get_keyval(conf, "extendedTemp", temp, + proxy->target_temperature()); if (is_enabled(f_cv_external)) { // In the case of an "external" coordinate, there is no coupling potential: // only the fictitious mass is meaningful @@ -669,14 +670,14 @@ int colvar::init_extended_Lagrangian(std::string const &conf) cvm::error("Error: \"extendedFluctuation\" must be positive.\n", COLVARS_INPUT_ERROR); return COLVARS_INPUT_ERROR; } - ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); + ext_force_k = proxy->boltzmann() * temp / (tolerance * tolerance); cvm::log("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " [E]/U^2\n"); get_keyval(conf, "extendedTimeConstant", extended_period, 200.0); if (extended_period <= 0.0) { cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", COLVARS_INPUT_ERROR); } - ext_mass = (cvm::boltzmann() * temp * extended_period * extended_period) + ext_mass = (proxy->boltzmann() * temp * extended_period * extended_period) / (4.0 * PI * PI * tolerance * tolerance); cvm::log("Computed fictitious mass: " + cvm::to_str(ext_mass) + " [E]/(U/fs)^2 (U: colvar unit)\n"); } @@ -697,7 +698,7 @@ int colvar::init_extended_Lagrangian(std::string const &conf) enable(f_cv_Langevin); ext_gamma *= 1.0e-3; // correct as long as input is required in ps-1 and cvm::dt() is in fs // 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))); + ext_sigma = cvm::sqrt(2.0 * proxy->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); @@ -1664,6 +1665,7 @@ int colvar::calc_cvc_Jacobians(int first_cvc, size_t num_cvcs) int colvar::collect_cvc_Jacobians() { + colvarproxy *proxy = cvm::main()->proxy; if (is_enabled(f_cv_Jacobian)) { fj.reset(); for (size_t i = 0; i < cvcs.size(); i++) { @@ -1676,7 +1678,7 @@ int colvar::collect_cvc_Jacobians() // linear combination is assumed fj += (cvcs[i])->Jacobian_derivative() * (cvcs[i])->sup_coeff / active_cvc_square_norm; } - fj *= cvm::boltzmann() * cvm::temperature(); + fj *= proxy->boltzmann() * proxy->target_temperature(); } return COLVARS_OK; @@ -2411,8 +2413,8 @@ std::ostream & colvar::write_state(std::ostream &os) { os << "}\n\n"; - if (runave_os) { - cvm::main()->proxy->flush_output_stream(runave_os); + if (runave_outfile.size() > 0) { + cvm::main()->proxy->flush_output_stream(runave_outfile); } return os; @@ -2536,9 +2538,13 @@ int colvar::write_output_files() } cvm::log("Writing correlation function to file \""+acf_outfile+"\".\n"); cvm::backup_file(acf_outfile.c_str()); - std::ostream *acf_os = cvm::proxy->output_stream(acf_outfile); - if (!acf_os) return cvm::get_error(); - error_code |= write_acf(*acf_os); + std::ostream &acf_os = cvm::proxy->output_stream(acf_outfile, + "colvar ACF file"); + if (!acf_os) { + error_code |= COLVARS_FILE_ERROR; + } else { + error_code |= write_acf(acf_os); + } cvm::proxy->close_output_stream(acf_outfile); } } @@ -2807,6 +2813,7 @@ int colvar::write_acf(std::ostream &os) int colvar::calc_runave() { int error_code = COLVARS_OK; + colvarproxy *proxy = cvm::main()->proxy; if (x_history.empty()) { @@ -2831,22 +2838,22 @@ int colvar::calc_runave() if ((*x_history_p).size() >= runave_length-1) { - if (runave_os == NULL) { - if (runave_outfile.size() == 0) { - runave_outfile = std::string(cvm::output_prefix()+"."+ - this->name+".runave.traj"); - } + if (runave_outfile.size() == 0) { + runave_outfile = std::string(cvm::output_prefix()+"."+ + this->name+".runave.traj"); + } + if (! proxy->output_stream_exists(runave_outfile)) { size_t const this_cv_width = x.output_width(cvm::cv_width); - cvm::proxy->backup_file(runave_outfile); - runave_os = cvm::proxy->output_stream(runave_outfile); - runave_os->setf(std::ios::scientific, std::ios::floatfield); - *runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2) - << " " - << cvm::wrap_string("running average", this_cv_width) - << " " - << cvm::wrap_string("running stddev", this_cv_width) - << "\n"; + std::ostream &runave_os = proxy->output_stream(runave_outfile, + "colvar running average"); + runave_os.setf(std::ios::scientific, std::ios::floatfield); + runave_os << "# " << cvm::wrap_string("step", cvm::it_width-2) + << " " + << cvm::wrap_string("running average", this_cv_width) + << " " + << cvm::wrap_string("running stddev", this_cv_width) + << "\n"; } runave = x; @@ -2866,12 +2873,17 @@ int colvar::calc_runave() } runave_variance *= 1.0 / cvm::real(runave_length-1); - *runave_os << std::setw(cvm::it_width) << cvm::step_relative() - << " " - << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) - << runave << " " - << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) - << cvm::sqrt(runave_variance) << "\n"; + if (runave_outfile.size() > 0) { + std::ostream &runave_os = proxy->output_stream(runave_outfile); + runave_os << std::setw(cvm::it_width) << cvm::step_relative() + << " " + << std::setprecision(cvm::cv_prec) + << std::setw(cvm::cv_width) + << runave << " " + << std::setprecision(cvm::cv_prec) + << std::setw(cvm::cv_width) + << cvm::sqrt(runave_variance) << "\n"; + } } history_add_value(runave_length, *x_history_p, x); diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 7b2863a3b5..9af26dedd3 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -551,8 +551,6 @@ protected: size_t runave_stride; /// Name of the file to write the running average std::string runave_outfile; - /// File to write the running average - std::ostream *runave_os; /// Current value of the running average colvarvalue runave; /// Current value of the square deviation from the running average diff --git a/lib/colvars/colvar_UIestimator.h b/lib/colvars/colvar_UIestimator.h index ca1ecc8c43..30a90a5799 100644 --- a/lib/colvars/colvar_UIestimator.h +++ b/lib/colvars/colvar_UIestimator.h @@ -1,4 +1,4 @@ -// -*- c++ -*- +// -*- Mode:c++; c-basic-offset: 4; -*- // This file is part of the Collective Variables module (Colvars). // The original version of Colvars and its updates are located at: @@ -21,6 +21,7 @@ // only for colvar module! // when integrated into other code, just remove this line and "...cvm::backup_file(...)" #include "colvarmodule.h" +#include "colvarproxy.h" namespace UIestimator { const int Y_SIZE = 21; // defines the range of extended CV with respect to a given CV @@ -380,6 +381,8 @@ namespace UIestimator { public: // calculate gradients from the internal variables void calc_pmf() { + colvarproxy *proxy = cvm::main()->proxy; + int norm; int i, j, k; @@ -455,7 +458,7 @@ namespace UIestimator { std::vector grad_temp(dimension, 0); for (k = 0; k < dimension; k++) { diff_av[k] /= (norm > 0 ? norm : 1); - av[k] = cvm::boltzmann() * temperature * av[k] / (norm > 0 ? norm : 1); + av[k] = proxy->boltzmann() * temperature * av[k] / (norm > 0 ? norm : 1); grad_temp[k] = av[k] - krestr[k] * diff_av[k]; } grad.set_value(loop_flag_x, grad_temp); @@ -515,15 +518,16 @@ namespace UIestimator { // only for colvars module! if (written_1D) cvm::backup_file(pmf_filename.c_str()); - std::ostream* ofile_pmf = cvm::proxy->output_stream(pmf_filename.c_str()); + std::ostream &ofile_pmf = cvm::proxy->output_stream(pmf_filename, + "PMF file"); std::vector position(1, 0); for (double i = lowerboundary[0]; i < upperboundary[0] + EPSILON; i += width[0]) { - *ofile_pmf << i << " "; + ofile_pmf << i << " "; position[0] = i + EPSILON; - *ofile_pmf << oneD_pmf.get_value(position) << std::endl; + ofile_pmf << oneD_pmf.get_value(position) << std::endl; } - cvm::proxy->close_output_stream(pmf_filename.c_str()); + cvm::proxy->close_output_stream(pmf_filename); written_1D = true; } @@ -541,7 +545,8 @@ namespace UIestimator { void write_interal_data() { std::string internal_filename = output_filename + ".UI.internal"; - std::ostream* ofile_internal = cvm::proxy->output_stream(internal_filename.c_str()); + std::ostream &ofile_internal = cvm::proxy->output_stream(internal_filename, + "UI internal file"); std::vector loop_flag(dimension, 0); for (int i = 0; i < dimension; i++) { @@ -551,11 +556,11 @@ namespace UIestimator { int n = 0; while (n >= 0) { for (int j = 0; j < dimension; j++) { - *ofile_internal << loop_flag[j] + 0.5 * width[j] << " "; + ofile_internal << loop_flag[j] + 0.5 * width[j] << " "; } for (int k = 0; k < dimension; k++) { - *ofile_internal << grad.get_value(loop_flag)[k] << " "; + ofile_internal << grad.get_value(loop_flag)[k] << " "; } std::vector ii(dimension,0); @@ -563,10 +568,10 @@ namespace UIestimator { for (double j = loop_flag[1] - 10; j< loop_flag[1] + 10 + EPSILON; j+=width[1]) { ii[0] = i; ii[1] = j; - *ofile_internal << i <<" "<output_stream(grad_filename.c_str()); - std::ostream* ofile_hist = cvm::proxy->output_stream(hist_filename.c_str(), std::ios::app); - std::ostream* ofile_count = cvm::proxy->output_stream(count_filename.c_str()); + std::ostream &ofile = cvm::proxy->output_stream(grad_filename, + "gradient file"); + std::ostream &ofile_hist = cvm::proxy->output_stream(hist_filename, + "gradient history file"); + std::ostream &ofile_count = cvm::proxy->output_stream(count_filename, + "count file"); - writehead(*ofile); - writehead(*ofile_hist); - writehead(*ofile_count); + writehead(ofile); + writehead(ofile_hist); + writehead(ofile_count); if (dimension == 1) { calc_1D_pmf(); @@ -617,19 +625,19 @@ namespace UIestimator { i = 0; while (i >= 0) { for (j = 0; j < dimension; j++) { - *ofile << loop_flag[j] + 0.5 * width[j] << " "; - *ofile_hist << loop_flag[j] + 0.5 * width[j] << " "; - *ofile_count << loop_flag[j] + 0.5 * width[j] << " "; + ofile << loop_flag[j] + 0.5 * width[j] << " "; + ofile_hist << loop_flag[j] + 0.5 * width[j] << " "; + ofile_count << loop_flag[j] + 0.5 * width[j] << " "; } if (restart == false) { for (j = 0; j < dimension; j++) { - *ofile << grad.get_value(loop_flag)[j] << " "; - *ofile_hist << grad.get_value(loop_flag)[j] << " "; + ofile << grad.get_value(loop_flag)[j] << " "; + ofile_hist << grad.get_value(loop_flag)[j] << " "; } - *ofile << std::endl; - *ofile_hist << std::endl; - *ofile_count << count.get_value(loop_flag) << " " < upperboundary[i] - width[i] + EPSILON) { loop_flag[i] = lowerboundary[i]; i--; - *ofile << std::endl; - *ofile_hist << std::endl; - *ofile_count << std::endl; + ofile << std::endl; + ofile_hist << std::endl; + ofile_count << std::endl; } else break; } } cvm::proxy->close_output_stream(grad_filename.c_str()); - cvm::proxy->close_output_stream(hist_filename.c_str()); + // cvm::proxy->close_output_stream(hist_filename.c_str()); cvm::proxy->close_output_stream(count_filename.c_str()); written = true; @@ -677,6 +685,7 @@ namespace UIestimator { int dimension_temp; int i, j, k, l, m; + colvarproxy *proxy = cvm::main()->proxy; std::vector loop_bin_size(dimension, 0); std::vector position_temp(dimension, 0); std::vector grad_temp(dimension, 0); @@ -687,8 +696,14 @@ namespace UIestimator { std::string count_filename = filename[i] + ".UI.count"; std::string grad_filename = filename[i] + ".UI.grad"; - std::ifstream count_file(count_filename.c_str(), std::ios::in); - std::ifstream grad_file(grad_filename.c_str(), std::ios::in); + std::istream &count_file = + proxy->input_stream(count_filename, "count filename"); + std::istream &grad_file = + proxy->input_stream(grad_filename, "gradient filename"); + + if (!count_file || !grad_file) { + return; + } count_file >> sharp >> dimension_temp; grad_file >> sharp >> dimension_temp; @@ -724,8 +739,8 @@ namespace UIestimator { input_count.increase_value(position_temp, count_temp); } - count_file.close(); - grad_file.close(); + proxy->close_input_stream(count_filename); + proxy->close_input_stream(grad_filename); } } }; diff --git a/lib/colvars/colvar_arithmeticpath.h b/lib/colvars/colvar_arithmeticpath.h index cf613af389..bea86a1014 100644 --- a/lib/colvars/colvar_arithmeticpath.h +++ b/lib/colvars/colvar_arithmeticpath.h @@ -5,7 +5,6 @@ #include #include -#include #include #include diff --git a/lib/colvars/colvar_geometricpath.h b/lib/colvars/colvar_geometricpath.h index 7f8acfd233..9ff78261f2 100644 --- a/lib/colvars/colvar_geometricpath.h +++ b/lib/colvars/colvar_geometricpath.h @@ -12,10 +12,8 @@ #include #include -#include #include #include -#include namespace GeometricPathCV { diff --git a/lib/colvars/colvar_neuralnetworkcompute.cpp b/lib/colvars/colvar_neuralnetworkcompute.cpp index d0f9633652..a1ad717946 100644 --- a/lib/colvars/colvar_neuralnetworkcompute.cpp +++ b/lib/colvars/colvar_neuralnetworkcompute.cpp @@ -1,6 +1,19 @@ +// -*- Mode:c++; c-basic-offset: 4; -*- + +// This file is part of the Collective Variables module (Colvars). +// The original version of Colvars and its updates are located at: +// https://github.com/Colvars/colvars +// Please update all Colvars source files before making any changes. +// If you wish to distribute your changes, please submit them to the +// Colvars repository at GitHub. + +#include +#include + #if (__cplusplus >= 201103L) #include "colvar_neuralnetworkcompute.h" #include "colvarparse.h" +#include "colvarproxy.h" namespace neuralnetworkCV { std::map, std::function>> activation_function_map @@ -124,12 +137,10 @@ void denseLayer::readFromFile(const std::string& weights_file, const std::string m_weights.clear(); m_biases.clear(); std::string line; - std::ifstream ifs_weights(weights_file.c_str()); - if (!ifs_weights) { - throw std::runtime_error("Cannot open file " + weights_file); - } + colvarproxy *proxy = cvm::main()->proxy; + auto &ifs_weights = proxy->input_stream(weights_file, "weights file"); while (std::getline(ifs_weights, line)) { - if (ifs_weights.bad()) { + if (!ifs_weights) { throw std::runtime_error("I/O error while reading " + weights_file); } std::vector splitted_data; @@ -146,13 +157,12 @@ void denseLayer::readFromFile(const std::string& weights_file, const std::string m_weights.push_back(weights_tmp); } } + proxy->close_input_stream(weights_file); + // parse biases file - std::ifstream ifs_biases(biases_file.c_str()); - if (!ifs_biases) { - throw std::runtime_error("Cannot open file " + biases_file); - } + auto &ifs_biases = proxy->input_stream(biases_file, "biases file"); while (std::getline(ifs_biases, line)) { - if (ifs_biases.bad()) { + if (!ifs_biases) { throw std::runtime_error("I/O error while reading " + biases_file); } std::vector splitted_data; @@ -167,6 +177,8 @@ void denseLayer::readFromFile(const std::string& weights_file, const std::string m_biases.push_back(bias); } } + proxy->close_input_stream(biases_file); + m_input_size = m_weights[0].size(); m_output_size = m_weights.size(); } diff --git a/lib/colvars/colvar_neuralnetworkcompute.h b/lib/colvars/colvar_neuralnetworkcompute.h index a48f8c7f14..5a56887431 100644 --- a/lib/colvars/colvar_neuralnetworkcompute.h +++ b/lib/colvars/colvar_neuralnetworkcompute.h @@ -1,3 +1,12 @@ +// -*- 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. + #if (__cplusplus >= 201103L) #ifndef NEURALNETWORKCOMPUTE_H #define NEURALNETWORKCOMPUTE_H @@ -5,7 +14,6 @@ #include #include #include -#include #include #include #include diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index f950bf5965..e31041ea6b 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -10,6 +10,8 @@ #include #include #include +#include +#include #include "colvarmodule.h" #include "colvarproxy.h" @@ -29,12 +31,8 @@ cvm::atom::atom() cvm::atom::atom(int atom_number) { - colvarproxy *p = cvm::proxy; + colvarproxy *p = cvm::main()->proxy; index = p->init_atom(atom_number); - if (cvm::debug()) { - cvm::log("The index of this atom in the colvarproxy arrays is "+ - cvm::to_str(index)+".\n"); - } id = p->get_atom_id(index); update_mass(); update_charge(); @@ -46,12 +44,8 @@ cvm::atom::atom(cvm::residue_id const &residue, std::string const &atom_name, std::string const &segment_id) { - colvarproxy *p = cvm::proxy; + colvarproxy *p = cvm::main()->proxy; index = p->init_atom(residue, atom_name, segment_id); - if (cvm::debug()) { - cvm::log("The index of this atom in the colvarproxy_namd arrays is "+ - cvm::to_str(index)+".\n"); - } id = p->get_atom_id(index); update_mass(); update_charge(); @@ -62,7 +56,9 @@ cvm::atom::atom(cvm::residue_id const &residue, cvm::atom::atom(atom const &a) : index(a.index) { - id = (cvm::proxy)->get_atom_id(index); + colvarproxy *p = cvm::main()->proxy; + id = p->get_atom_id(index); + p->increase_refcount(index); update_mass(); update_charge(); reset_data(); @@ -72,7 +68,7 @@ cvm::atom::atom(atom const &a) cvm::atom::~atom() { if (index >= 0) { - (cvm::proxy)->clear_atom(index); + (cvm::main()->proxy)->clear_atom(index); } } @@ -80,7 +76,7 @@ cvm::atom::~atom() cvm::atom & cvm::atom::operator = (cvm::atom const &a) { index = a.index; - id = (cvm::proxy)->get_atom_id(index); + id = (cvm::main()->proxy)->get_atom_id(index); update_mass(); update_charge(); reset_data(); @@ -113,7 +109,7 @@ cvm::atom_group::atom_group(std::vector const &atoms_in) cvm::atom_group::~atom_group() { if (is_enabled(f_ag_scalable) && !b_dummy) { - (cvm::proxy)->clear_atom_group(index); + (cvm::main()->proxy)->clear_atom_group(index); index = -1; } @@ -330,7 +326,7 @@ void cvm::atom_group::update_total_mass() } if (is_enabled(f_ag_scalable)) { - total_mass = (cvm::proxy)->get_atom_group_mass(index); + total_mass = (cvm::main()->proxy)->get_atom_group_mass(index); } else { total_mass = 0.0; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { @@ -351,7 +347,7 @@ void cvm::atom_group::update_total_charge() } if (is_enabled(f_ag_scalable)) { - total_charge = (cvm::proxy)->get_atom_group_charge(index); + total_charge = (cvm::main()->proxy)->get_atom_group_charge(index); } else { total_charge = 0.0; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index bb2f2310f7..18969fc531 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -7,7 +7,7 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -#include +#include #include #include "colvarmodule.h" @@ -42,6 +42,8 @@ colvarbias::colvarbias(char const *key) int colvarbias::init(std::string const &conf) { + int error_code = COLVARS_OK; + name = bias_type + cvm::to_str(rank); colvarparse::set_string(conf); @@ -59,9 +61,9 @@ int colvarbias::init(std::string const &conf) if (bias_with_name != NULL) { if ((bias_with_name->rank != this->rank) || (bias_with_name->bias_type != this->bias_type)) { - cvm::error("Error: this bias cannot have the same name, \""+this->name+ - "\", as another bias.\n", COLVARS_INPUT_ERROR); - return COLVARS_INPUT_ERROR; + error_code |= cvm::error("Error: this bias cannot have the same name, \""+ + this->name+"\", as another bias.\n", + COLVARS_INPUT_ERROR); } } description = "bias " + name; @@ -71,9 +73,9 @@ int colvarbias::init(std::string const &conf) std::vector colvar_names; if (get_keyval(conf, "colvars", colvar_names)) { if (num_variables()) { - cvm::error("Error: cannot redefine the colvars that a bias was already defined on.\n", - COLVARS_INPUT_ERROR); - return COLVARS_INPUT_ERROR; + error_code |= cvm::error("Error: cannot redefine the colvars that " + "a bias was already defined on.\n", + COLVARS_INPUT_ERROR); } for (i = 0; i < colvar_names.size(); i++) { add_colvar(colvar_names[i]); @@ -82,8 +84,8 @@ int colvarbias::init(std::string const &conf) } if (!num_variables()) { - cvm::error("Error: no collective variables specified.\n", COLVARS_INPUT_ERROR); - return COLVARS_INPUT_ERROR; + error_code |= cvm::error("Error: no collective variables specified.\n", + COLVARS_INPUT_ERROR); } } else { @@ -112,8 +114,8 @@ int colvarbias::init(std::string const &conf) 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; + error_code |= cvm::error("Error: timeStepFactor must be 1 or greater.\n", + COLVARS_INPUT_ERROR); } // Use the scaling factors from a grid? @@ -124,25 +126,17 @@ int colvarbias::init(std::string const &conf) std::string biasing_force_scaling_factors_in_filename; get_keyval(conf, "scaledBiasingForceFactorsGrid", biasing_force_scaling_factors_in_filename, std::string()); - std::ifstream is; - cvm::log("Reading scaling factors for the forces of bias " + - name + " from " + biasing_force_scaling_factors_in_filename); - is.open(biasing_force_scaling_factors_in_filename.c_str()); - if (!is.is_open()) { - cvm::error("Error opening the grid file " + - biasing_force_scaling_factors_in_filename + " for reading"); - } biasing_force_scaling_factors = new colvar_grid_scalar(colvars); - biasing_force_scaling_factors->read_multicol(is, true); + error_code |= biasing_force_scaling_factors->read_multicol(biasing_force_scaling_factors_in_filename, + "grid file"); biasing_force_scaling_factors_bin.assign(num_variables(), 0); - is.close(); } // Now that children are defined, we can solve dependencies enable(f_cvb_active); if (cvm::debug()) print_state(); - return COLVARS_OK; + return error_code; } @@ -377,10 +371,11 @@ int colvarbias::calc_forces(std::vector const *) } -void colvarbias::communicate_forces() +int colvarbias::communicate_forces() { + int error_code = COLVARS_OK; if (! is_enabled(f_cvb_apply_force)) { - return; + return error_code; } cvm::real biasing_force_factor = 1.0; size_t i = 0; @@ -409,6 +404,7 @@ void colvarbias::communicate_forces() } previous_colvar_forces[i] = colvar_forces[i]; } + return error_code; } @@ -573,11 +569,11 @@ 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()); + std::ostream &os = cvm::proxy->output_stream(filename.c_str(), "bias state file"); int error_code = COLVARS_OK; - if (os != NULL) { - os->setf(std::ios::scientific, std::ios::floatfield); - error_code = write_state(*os).good() ? COLVARS_OK : COLVARS_FILE_ERROR; + if (os) { + os.setf(std::ios::scientific, std::ios::floatfield); + error_code = write_state(os) ? COLVARS_OK : COLVARS_FILE_ERROR; } else { error_code = COLVARS_FILE_ERROR; } @@ -600,17 +596,19 @@ int colvarbias::write_state_string(std::string &output) 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(); + std::string filename(prefix+std::string(".colvars.state")); + std::istream *is = &(cvm::main()->proxy->input_stream(filename, + "bias state file", + false)); + if (!*is) { filename = prefix; - is.open(filename.c_str()); + is = &(cvm::main()->proxy->input_stream(filename, "bias state file")); } - return read_state(is).good() ? COLVARS_OK : - cvm::error("Error: in reading state for \""+name+"\" from input file \""+ - std::string(filename)+"\".\n", COLVARS_FILE_ERROR); + + if (read_state(*is)) { + return cvm::main()->proxy->close_input_stream(filename); + } + return COLVARS_FILE_ERROR; } @@ -899,6 +897,8 @@ std::istream & colvarbias_ti::read_state_data(std::istream &is) int colvarbias_ti::write_output_files() { + int error_code = COLVARS_OK; + if (!has_data) { // nothing to write return COLVARS_OK; @@ -906,38 +906,30 @@ int colvarbias_ti::write_output_files() std::string const ti_output_prefix = cvm::output_prefix()+"."+this->name; - std::ostream *os = NULL; - if (is_enabled(f_cvb_write_ti_samples)) { std::string const ti_count_file_name(ti_output_prefix+".ti.count"); - os = cvm::proxy->output_stream(ti_count_file_name); - if (os) { - ti_count->write_multicol(*os); - cvm::proxy->close_output_stream(ti_count_file_name); - } + error_code |= ti_count->write_multicol(ti_count_file_name, "TI count file"); std::string const ti_grad_file_name(ti_output_prefix+".ti.force"); - os = cvm::proxy->output_stream(ti_grad_file_name); - if (os) { - ti_avg_forces->write_multicol(*os); - cvm::proxy->close_output_stream(ti_grad_file_name); - } + error_code |= ti_avg_forces->write_multicol(ti_grad_file_name, "TI gradient file"); } if (is_enabled(f_cvb_write_ti_pmf)) { std::string const pmf_file_name(ti_output_prefix+".ti.pmf"); cvm::log("Writing TI PMF to file \""+pmf_file_name+"\".\n"); - os = cvm::proxy->output_stream(pmf_file_name); + std::ostream &os = cvm::proxy->output_stream(pmf_file_name, "TI PMF"); if (os) { // get the FE gradient ti_avg_forces->multiply_constant(-1.0); - ti_avg_forces->write_1D_integral(*os); + ti_avg_forces->write_1D_integral(os); ti_avg_forces->multiply_constant(-1.0); cvm::proxy->close_output_stream(pmf_file_name); + } else { + error_code |= COLVARS_FILE_ERROR; } } - return COLVARS_OK; + return error_code; } diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index ecb3e1eff1..e865902be6 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -73,7 +73,7 @@ public: virtual int calc_forces(std::vector const *values); /// Send forces to the collective variables - void communicate_forces(); + int communicate_forces(); /// Carry out operations needed before next step is run virtual int end_of_step(); diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 5cb5e7af24..7d6b7b7fef 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -7,6 +7,8 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include + #include "colvarmodule.h" #include "colvar.h" #include "colvarbias_abf.h" @@ -38,17 +40,17 @@ colvarbias_abf::colvarbias_abf(char const *key) int colvarbias_abf::init(std::string const &conf) { + colvarproxy *proxy = cvm::main()->proxy; + colvarbias::init(conf); cvm::main()->cite_feature("ABF colvar bias implementation"); - colvarproxy *proxy = cvm::main()->proxy; - enable(f_cvb_scalar_variables); enable(f_cvb_calc_pmf); - // TODO relax this in case of VMD plugin - if (cvm::temperature() == 0.0) + if ((proxy->target_temperature() == 0.0) && proxy->simulation_running()) { cvm::log("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n"); + } // ************* parsing general ABF options *********************** @@ -277,7 +279,7 @@ int colvarbias_abf::init(std::string const &conf) cvm::restart_out_freq, UI_restart, // whether restart from a .count and a .grad file input_prefix, // the prefixes of input files - cvm::temperature()); + proxy->target_temperature()); } } @@ -599,17 +601,17 @@ int colvarbias_abf::replica_share() { template int colvarbias_abf::write_grid_to_file(T const *grid, std::string const &filename, bool close) { - std::ostream *os = cvm::proxy->output_stream(filename); + std::ostream &os = cvm::proxy->output_stream(filename); if (!os) { return cvm::error("Error opening file " + filename + " for writing.\n", COLVARS_ERROR | COLVARS_FILE_ERROR); } - grid->write_multicol(*os); + 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); + os << std::endl; + cvm::proxy->flush_output_stream(filename); } // In dimension higher than 2, dx is easier to handle and visualize @@ -617,11 +619,11 @@ template int colvarbias_abf::write_grid_to_file(T const *grid, // (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); + std::ostream &dx_os = cvm::proxy->output_stream(dx); if (!dx_os) { return cvm::error("Error opening file " + dx + " for writing.\n", COLVARS_ERROR | COLVARS_FILE_ERROR); } - grid->write_opendx(*dx_os); + grid->write_opendx(dx_os); // if (close) { cvm::proxy->close_output_stream(dx); // } @@ -637,6 +639,8 @@ template int colvarbias_abf::write_grid_to_file(T const *grid, void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool close) { + colvarproxy *proxy = cvm::main()->proxy; + write_grid_to_file(samples, prefix + ".count", close); write_grid_to_file(gradients, prefix + ".grad", close); @@ -660,7 +664,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool clo czar_gradients->index_ok(ix); czar_gradients->incr(ix)) { for (size_t n = 0; n < czar_gradients->multiplicity(); n++) { czar_gradients->set_value(ix, z_gradients->value_output(ix, n) - - cvm::temperature() * cvm::boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n); + - proxy->target_temperature() * proxy->boltzmann() * z_samples->log_gradient_finite_diff(ix, n), n); } } write_grid_to_file(czar_gradients, prefix + ".czar.grad", close); @@ -698,8 +702,10 @@ int colvarbias_abf::bin_count(int bin_index) { } -void colvarbias_abf::read_gradients_samples() +int colvarbias_abf::read_gradients_samples() { + int error_code = COLVARS_OK; + std::string samples_in_name, gradients_in_name, z_samples_in_name, z_gradients_in_name; for ( size_t i = 0; i < input_prefix.size(); i++ ) { @@ -708,43 +714,30 @@ void colvarbias_abf::read_gradients_samples() z_samples_in_name = input_prefix[i] + ".zcount"; z_gradients_in_name = input_prefix[i] + ".zgrad"; // For user-provided files, the per-bias naming scheme may not apply + cvm::log("Reading sample count from " + samples_in_name + + " and gradient from " + gradients_in_name); - std::ifstream is; + error_code |= samples->read_multicol(samples_in_name, + "ABF samples file", + true); - cvm::log("Reading sample count from " + samples_in_name + " and gradient from " + gradients_in_name); - is.open(samples_in_name.c_str()); - if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading"); - samples->read_multicol(is, true); - is.close(); - is.clear(); - - is.open(gradients_in_name.c_str()); - if (!is.is_open()) { - cvm::error("Error opening ABF gradient file " + - gradients_in_name + " for reading", COLVARS_INPUT_ERROR); - } else { - gradients->read_multicol(is, true); - is.close(); - } + error_code |= gradients->read_multicol(gradients_in_name, + "ABF gradient file", + true); if (b_CZAR_estimator) { // Read eABF z-averaged data for CZAR cvm::log("Reading z-histogram from " + z_samples_in_name + " and z-gradient from " + z_gradients_in_name); - - is.clear(); - is.open(z_samples_in_name.c_str()); - if (!is.is_open()) cvm::error("Error opening eABF z-histogram file " + z_samples_in_name + " for reading"); - z_samples->read_multicol(is, true); - is.close(); - is.clear(); - - is.open(z_gradients_in_name.c_str()); - if (!is.is_open()) cvm::error("Error opening eABF z-gradient file " + z_gradients_in_name + " for reading"); - z_gradients->read_multicol(is, true); - is.close(); + error_code |= z_samples->read_multicol(z_samples_in_name, + "eABF z-histogram file", + true); + error_code |= z_gradients->read_multicol(z_gradients_in_name, + "eABF z-gradient file", + true); } } - return; + + return error_code; } diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index 3a3120a058..f5d5bd267f 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -155,7 +155,7 @@ private: 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(); + int read_gradients_samples(); /// Template used in write_gradient_samples() template int write_grid_to_file(T const *grid, diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 5f39b0d93a..b432659bf4 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -7,9 +7,12 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include +#include #include #include "colvarmodule.h" +#include "colvarproxy.h" #include "colvarbias.h" #include "colvarbias_alb.h" @@ -36,6 +39,7 @@ colvarbias_alb::colvarbias_alb(char const *key) int colvarbias_alb::init(std::string const &conf) { + colvarproxy *proxy = cvm::main()->proxy; colvarbias::init(conf); cvm::main()->cite_feature("ALB colvar bias implementation"); @@ -110,10 +114,12 @@ int colvarbias_alb::init(std::string const &conf) if (!get_keyval(conf, "forceRange", max_coupling_range, max_coupling_range)) { //set to default for (i = 0; i < num_variables(); i++) { - if (cvm::temperature() > 0) - max_coupling_range[i] = 3 * cvm::temperature() * cvm::boltzmann(); - else - max_coupling_range[i] = 3 * cvm::boltzmann(); + if (proxy->target_temperature() > 0.0) { + max_coupling_range[i] = 3 * proxy->target_temperature() * + proxy->boltzmann(); + } else { + max_coupling_range[i] = 3 * proxy->boltzmann(); + } } } @@ -139,6 +145,7 @@ colvarbias_alb::~colvarbias_alb() int colvarbias_alb::update() { + colvarproxy *proxy = cvm::main()->proxy; bias_energy = 0.0; update_calls++; @@ -214,10 +221,11 @@ int colvarbias_alb::update() temp = 2. * (means[i] / (static_cast (colvar_centers[i])) - 1) * ssd[i] / (update_calls - 1); - if (cvm::temperature() > 0) - step_size = temp / (cvm::temperature() * cvm::boltzmann()); - else - step_size = temp / cvm::boltzmann(); + if (proxy->target_temperature() > 0.0) { + step_size = temp / (proxy->target_temperature() * proxy->boltzmann()); + } else { + step_size = temp / proxy->boltzmann(); + } means[i] = 0; ssd[i] = 0; diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index 2719489560..84f1a5bdee 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -179,31 +179,19 @@ int colvarbias_histogram::write_output_files() return COLVARS_OK; } + int error_code = COLVARS_OK; + 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); - if (!grid_os) { - return cvm::error("Error opening histogram file "+out_name+ - " for writing.\n", COLVARS_FILE_ERROR); - } - grid->write_multicol(*grid_os); - cvm::proxy->close_output_stream(out_name); + error_code |= grid->write_multicol(out_name, "histogram output file"); } 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); - if (!grid_os) { - return cvm::error("Error opening histogram file "+out_name_dx+ - " for writing.\n", COLVARS_FILE_ERROR); - } - grid->write_opendx(*grid_os); - cvm::proxy->close_output_stream(out_name_dx); + error_code |= grid->write_opendx(out_name_dx, "histogram DX output file"); } - return COLVARS_OK; + return error_code; } diff --git a/lib/colvars/colvarbias_histogram_reweight_amd.cpp b/lib/colvars/colvarbias_histogram_reweight_amd.cpp index 97c278abcb..85f1bde35c 100644 --- a/lib/colvars/colvarbias_histogram_reweight_amd.cpp +++ b/lib/colvars/colvarbias_histogram_reweight_amd.cpp @@ -81,6 +81,7 @@ int colvarbias_reweightaMD::init(std::string const &conf) { } int colvarbias_reweightaMD::update() { + colvarproxy *proxy = cvm::main()->proxy; int error_code = COLVARS_OK; if (cvm::step_relative() >= start_after_steps) { // update base class @@ -110,7 +111,7 @@ int colvarbias_reweightaMD::update() { grid->acc_value(previous_bin, reweighting_factor); if (b_use_cumulant_expansion) { const cvm::real dV = cvm::logn(reweighting_factor) * - cvm::temperature() * cvm::boltzmann(); + proxy->target_temperature() * proxy->boltzmann(); grid_dV->acc_value(previous_bin, dV); grid_dV_square->acc_value(previous_bin, dV * dV); } @@ -129,7 +130,7 @@ int colvarbias_reweightaMD::update() { grid->acc_value(previous_bin, reweighting_factor); if (b_use_cumulant_expansion) { const cvm::real dV = cvm::logn(reweighting_factor) * - cvm::temperature() * cvm::boltzmann(); + proxy->target_temperature() * proxy->boltzmann(); grid_dV->acc_value(previous_bin, dV); grid_dV_square->acc_value(previous_bin, dV * dV); } @@ -138,10 +139,6 @@ int colvarbias_reweightaMD::update() { } previous_bin.assign(num_variables(), 0); - if (output_freq && (cvm::step_absolute() % output_freq) == 0) { - write_output_files(); - } - error_code |= cvm::get_error(); } return error_code; @@ -160,7 +157,7 @@ int colvarbias_reweightaMD::write_output_files() { (cvm::step_absolute() % history_freq) == 0; if (write_history) { error_code |= write_exponential_reweighted_pmf( - out_name_pmf + ".hist", (cvm::step_relative() > 0)); + out_name_pmf + ".hist", true); error_code |= write_count(out_count_prefix + ".hist", (cvm::step_relative() > 0)); } @@ -170,7 +167,7 @@ int colvarbias_reweightaMD::write_output_files() { error_code |= write_cumulant_expansion_pmf(out_name_cumulant_pmf); if (write_history) { error_code |= write_cumulant_expansion_pmf( - out_name_cumulant_pmf + ".hist", (cvm::step_relative() > 0)); + out_name_cumulant_pmf + ".hist", true); } } error_code |= cvm::get_error(); @@ -178,17 +175,13 @@ int colvarbias_reweightaMD::write_output_files() { } int colvarbias_reweightaMD::write_exponential_reweighted_pmf( - const std::string& p_output_prefix, bool append) { + const std::string& p_output_prefix, bool keep_open) { const std::string output_pmf = p_output_prefix + ".pmf"; + cvm::log("Writing the accelerated MD PMF file \"" + output_pmf + "\".\n"); - if (!append) { - cvm::backup_file(output_pmf.c_str()); - } - const std::ios::openmode mode = (append ? std::ios::app : std::ios::out); - std::ostream *pmf_grid_os = cvm::proxy->output_stream(output_pmf, mode); + std::ostream &pmf_grid_os = cvm::proxy->output_stream(output_pmf, "PMF file"); if (!pmf_grid_os) { - return cvm::error("Error opening PMF file " + output_pmf + - " for writing.\n", COLVARS_FILE_ERROR); + return COLVARS_FILE_ERROR; } pmf_grid_exp_avg->copy_grid(*grid); // compute the average @@ -200,19 +193,18 @@ int colvarbias_reweightaMD::write_exponential_reweighted_pmf( } } hist_to_pmf(pmf_grid_exp_avg, grid_count); - pmf_grid_exp_avg->write_multicol(*pmf_grid_os); - cvm::proxy->close_output_stream(output_pmf); + pmf_grid_exp_avg->write_multicol(pmf_grid_os); + if (!keep_open) { + cvm::proxy->close_output_stream(output_pmf); + } + if (b_write_gradients) { const std::string output_grad = p_output_prefix + ".grad"; cvm::log("Writing the accelerated MD gradients file \"" + output_grad + "\".\n"); - if (!append) { - cvm::backup_file(output_grad.c_str()); - } - std::ostream *grad_grid_os = cvm::proxy->output_stream(output_grad, mode); + std::ostream &grad_grid_os = cvm::proxy->output_stream(output_grad, "gradient file"); if (!grad_grid_os) { - return cvm::error("Error opening grad file " + output_grad + - " for writing.\n", COLVARS_FILE_ERROR); + return COLVARS_FILE_ERROR; } for (std::vector ix = grad_grid_exp_avg->new_index(); grad_grid_exp_avg->index_ok(ix); grad_grid_exp_avg->incr(ix)) { @@ -221,38 +213,37 @@ int colvarbias_reweightaMD::write_exponential_reweighted_pmf( ix, pmf_grid_exp_avg->gradient_finite_diff(ix, n), n); } } - grad_grid_exp_avg->write_multicol(*grad_grid_os); - cvm::proxy->close_output_stream(output_grad); + grad_grid_exp_avg->write_multicol(grad_grid_os); + if (!keep_open) { + cvm::proxy->close_output_stream(output_grad); + } } + return COLVARS_OK; } int colvarbias_reweightaMD::write_cumulant_expansion_pmf( - const std::string& p_output_prefix, bool append) { + const std::string& p_output_prefix, bool keep_open) { const std::string output_pmf = p_output_prefix + ".pmf"; cvm::log("Writing the accelerated MD PMF file using cumulant expansion: \"" + output_pmf + "\".\n"); - if (!append) cvm::backup_file(output_pmf.c_str()); - const std::ios::openmode mode = (append ? std::ios::app : std::ios::out); - std::ostream *pmf_grid_cumulant_os = cvm::proxy->output_stream(output_pmf, mode); + std::ostream &pmf_grid_cumulant_os = cvm::proxy->output_stream(output_pmf, "PMF file"); if (!pmf_grid_cumulant_os) { - return cvm::error("Error opening PMF file " + output_pmf + - " for writing.\n", COLVARS_FILE_ERROR); + return COLVARS_FILE_ERROR; } compute_cumulant_expansion_factor(grid_dV, grid_dV_square, grid_count, pmf_grid_cumulant); hist_to_pmf(pmf_grid_cumulant, grid_count); - pmf_grid_cumulant->write_multicol(*pmf_grid_cumulant_os); - cvm::proxy->close_output_stream(output_pmf); + pmf_grid_cumulant->write_multicol(pmf_grid_cumulant_os); + if (!keep_open) { + cvm::proxy->close_output_stream(output_pmf); + } + if (b_write_gradients) { const std::string output_grad = p_output_prefix + ".grad"; cvm::log("Writing the accelerated MD gradients file \"" + output_grad + "\".\n"); - if (!append) { - cvm::backup_file(output_grad.c_str()); - } - std::ostream *grad_grid_os = cvm::proxy->output_stream(output_grad, mode); + std::ostream &grad_grid_os = cvm::proxy->output_stream(output_grad, "grad file"); if (!grad_grid_os) { - return cvm::error("Error opening grad file " + output_grad + - " for writing.\n", COLVARS_FILE_ERROR); + return COLVARS_FILE_ERROR; } for (std::vector ix = grad_grid_cumulant->new_index(); grad_grid_cumulant->index_ok(ix); grad_grid_cumulant->incr(ix)) { @@ -261,31 +252,33 @@ int colvarbias_reweightaMD::write_cumulant_expansion_pmf( ix, pmf_grid_cumulant->gradient_finite_diff(ix, n), n); } } - grad_grid_cumulant->write_multicol(*grad_grid_os); + grad_grid_cumulant->write_multicol(grad_grid_os); cvm::proxy->close_output_stream(output_grad); } return COLVARS_OK; } -int colvarbias_reweightaMD::write_count(const std::string& p_output_prefix, bool append) { +int colvarbias_reweightaMD::write_count(const std::string& p_output_prefix, bool keep_open) { const std::string output_name = p_output_prefix + ".count"; cvm::log("Writing the accelerated MD count file \""+output_name+"\".\n"); - if (!append) cvm::backup_file(output_name.c_str()); - const std::ios::openmode mode = (append ? std::ios::app : std::ios::out); - std::ostream *grid_count_os = cvm::proxy->output_stream(output_name, mode); + std::ostream &grid_count_os = cvm::proxy->output_stream(output_name, "count file"); if (!grid_count_os) { - return cvm::error("Error opening count file "+output_name+ - " for writing.\n", COLVARS_FILE_ERROR); + return COLVARS_FILE_ERROR; + } + grid_count->write_multicol(grid_count_os); + if (!keep_open) { + cvm::proxy->close_output_stream(output_name); } - grid_count->write_multicol(*grid_count_os); - cvm::proxy->close_output_stream(output_name); return COLVARS_OK; } void colvarbias_reweightaMD::hist_to_pmf( - colvar_grid_scalar* hist, const colvar_grid_scalar* hist_count) const { + colvar_grid_scalar* hist, + const colvar_grid_scalar* hist_count) const +{ + colvarproxy *proxy = cvm::main()->proxy; if (hist->raw_data_num() == 0) return; - const cvm::real kbt = cvm::boltzmann() * cvm::temperature(); + const cvm::real kbt = proxy->boltzmann() * proxy->target_temperature(); bool first_min_element = true; bool first_max_element = true; cvm::real min_element = 0; @@ -329,12 +322,15 @@ void colvarbias_reweightaMD::hist_to_pmf( } } + void colvarbias_reweightaMD::compute_cumulant_expansion_factor( const colvar_grid_scalar* hist_dV, const colvar_grid_scalar* hist_dV_square, const colvar_grid_scalar* hist_count, - colvar_grid_scalar* cumulant_expansion_factor) const { - const cvm::real beta = 1.0 / (cvm::boltzmann() * cvm::temperature()); + colvar_grid_scalar* cumulant_expansion_factor) const +{ + colvarproxy *proxy = cvm::main()->proxy; + const cvm::real beta = 1.0 / (proxy->boltzmann() * proxy->target_temperature()); size_t i = 0; for (i = 0; i < hist_dV->raw_data_num(); ++i) { const cvm::real count = hist_count->value(i); diff --git a/lib/colvars/colvarbias_histogram_reweight_amd.h b/lib/colvars/colvarbias_histogram_reweight_amd.h index 8d74b5315d..f126738305 100644 --- a/lib/colvars/colvarbias_histogram_reweight_amd.h +++ b/lib/colvars/colvarbias_histogram_reweight_amd.h @@ -51,21 +51,21 @@ public: /// @brief output the PMF by the exponential average estimator /// @param[in] p_output_prefix the prefix of the output file - /// @param[in] append append the output to a .hist file if true + /// @param[in] keep_open Allow writing the history of the PMF virtual int write_exponential_reweighted_pmf( - const std::string& p_output_prefix, bool append = false); + const std::string& p_output_prefix, bool keep_open = false); /// @brief output the PMF by the cumulant expansion estimator /// @param[in] p_output_prefix the prefix of the output file - /// @param[in] append append the output to a .hist file if true + /// @param[in] keep_open Allow writing the history of the expansion virtual int write_cumulant_expansion_pmf( - const std::string& p_output_prefix, bool append = false); + const std::string& p_output_prefix, bool keep_open = false); /// @brief output the biased sampling /// @param[in] p_output_prefix the prefix of the output file - /// @param[in] append append the output to a .hist file if true + /// @param[in] keep_open Allow writing the history of the samples virtual int write_count( - const std::string& p_output_prefix, bool append = false); + const std::string& p_output_prefix, bool keep_open = false); protected: /// Current accelMD factor is the from previous frame std::vector previous_bin; diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 0b31479276..90497bf150 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -36,7 +36,6 @@ colvarbias_meta::colvarbias_meta(char const *key) : colvarbias(key), colvarbias_ti(key) { new_hills_begin = hills.end(); - hills_traj_os = NULL; hill_weight = 0.0; hill_width = 0.0; @@ -268,28 +267,31 @@ int colvarbias_meta::init_well_tempered_params(std::string const &conf) int colvarbias_meta::init_ebmeta_params(std::string const &conf) { + int error_code = COLVARS_OK; // for ebmeta target_dist = NULL; get_keyval(conf, "ebMeta", ebmeta, false); if(ebmeta){ cvm::main()->cite_feature("Ensemble-biased metadynamics (ebMetaD)"); if (use_grids && expand_grids) { - cvm::error("Error: expandBoundaries is not supported with " - "ebMeta please allocate wide enough boundaries for " - "each colvar ahead of time and set targetdistfile " - "accordingly.\n", COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: expandBoundaries is not supported with " + "ebMeta; please allocate wide enough boundaries " + "for each colvar ahead of time and set " + "targetDistFile accordingly.\n", + COLVARS_INPUT_ERROR); } target_dist = new colvar_grid_scalar(); - target_dist->init_from_colvars(colvars); + error_code |= target_dist->init_from_colvars(colvars); std::string target_dist_file; get_keyval(conf, "targetDistFile", target_dist_file); - std::ifstream targetdiststream(target_dist_file.c_str()); - target_dist->read_multicol(targetdiststream); + error_code |= target_dist->read_multicol(target_dist_file, + "ebMeta target histogram"); cvm::real min_val = target_dist->minimum_value(); cvm::real max_val = target_dist->maximum_value(); - if(min_val<0){ - cvm::error("Error: Target distribution of EBMetaD " - "has negative values!.\n", COLVARS_INPUT_ERROR); + if (min_val < 0.0) { + error_code |= cvm::error("Error: Target distribution of EBMetaD " + "has negative values!.\n", + COLVARS_INPUT_ERROR); } cvm::real target_dist_min_val; get_keyval(conf, "targetDistMinVal", target_dist_min_val, 1/1000000.0); @@ -301,17 +303,19 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf) cvm::log("NOTE: targetDistMinVal is set to zero, the minimum value of the target \n"); cvm::log(" distribution will be set as the minimum positive value.\n"); cvm::real min_pos_val = target_dist->minimum_pos_value(); - if(min_pos_val<=0){ - cvm::error("Error: Target distribution of EBMetaD has negative " - "or zero minimum positive value!.\n", COLVARS_INPUT_ERROR); + if (min_pos_val <= 0.0){ + error_code |= cvm::error("Error: Target distribution of EBMetaD has " + "negative or zero minimum positive value.\n", + COLVARS_INPUT_ERROR); } - if(min_val==0){ + if (min_val == 0.0){ cvm::log("WARNING: Target distribution has zero values.\n"); cvm::log("Zeros will be converted to the minimum positive value.\n"); target_dist->remove_small_values(min_pos_val); } } else { - cvm::error("Error: targetDistMinVal must be a value between 0 and 1!.\n", COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: targetDistMinVal must be a value " + "between 0 and 1.\n", COLVARS_INPUT_ERROR); } } // normalize target distribution and multiply by effective volume = exp(differential entropy) @@ -321,23 +325,18 @@ int colvarbias_meta::init_ebmeta_params(std::string const &conf) get_keyval(conf, "ebMetaEquilSteps", ebmeta_equil_steps, ebmeta_equil_steps); } - return COLVARS_OK; + return error_code; } colvarbias_meta::~colvarbias_meta() { colvarbias_meta::clear_state_data(); - colvarproxy *proxy = cvm::proxy; + colvarproxy *proxy = cvm::main()->proxy; - if (proxy->get_output_stream(replica_hills_file)) { - proxy->close_output_stream(replica_hills_file); - } + proxy->close_output_stream(replica_hills_file); - if (hills_traj_os) { - proxy->close_output_stream(hills_traj_file_name()); - hills_traj_os = NULL; - } + proxy->close_output_stream(hills_traj_file_name()); if (target_dist) { delete target_dist; @@ -392,9 +391,12 @@ colvarbias_meta::add_hill(colvarbias_meta::hill const &h) } // output to trajectory (if specified) - if (hills_traj_os) { - *hills_traj_os << (hills.back()).output_traj(); - cvm::proxy->flush_output_stream(hills_traj_os); + if (b_hills_traj) { + // Open trajectory file or access the one already open + std::ostream &hills_traj_os = + cvm::proxy->output_stream(hills_traj_file_name()); + hills_traj_os << (hills.back()).output_traj(); + cvm::proxy->flush_output_stream(hills_traj_file_name()); } has_data = true; @@ -424,12 +426,14 @@ colvarbias_meta::delete_hill(hill_iter &h) } } - if (hills_traj_os) { + if (b_hills_traj) { // output to the trajectory - *hills_traj_os << "# DELETED this hill: " - << (hills.back()).output_traj() - << "\n"; - cvm::proxy->flush_output_stream(hills_traj_os); + std::ostream &hills_traj_os = + cvm::proxy->output_stream(hills_traj_file_name()); + hills_traj_os << "# DELETED this hill: " + << (hills.back()).output_traj() + << "\n"; + cvm::proxy->flush_output_stream(hills_traj_file_name()); } return hills.erase(h); @@ -573,6 +577,7 @@ int colvarbias_meta::update_grid_params() int colvarbias_meta::update_bias() { + colvarproxy *proxy = cvm::main()->proxy; // add a new hill if the required time interval has passed if (((cvm::step_absolute() % new_hill_freq) == 0) && can_accumulate_data() && is_enabled(f_cvb_history_dependent)) { @@ -603,7 +608,7 @@ int colvarbias_meta::update_bias() } else { calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here, NULL); } - hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); + hills_scale *= cvm::exp(-1.0*hills_energy_sum_here/(bias_temperature*proxy->boltzmann())); } switch (comm) { @@ -618,10 +623,10 @@ int colvarbias_meta::update_bias() case multiple_replicas: add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale, colvar_values, colvar_sigmas, replica_id)); - std::ostream *replica_hills_os = - cvm::proxy->get_output_stream(replica_hills_file); + std::ostream &replica_hills_os = + cvm::proxy->output_stream(replica_hills_file); if (replica_hills_os) { - *replica_hills_os << hills.back(); + replica_hills_os << hills.back(); } else { return cvm::error("Error: in metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ @@ -979,25 +984,24 @@ void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first int colvarbias_meta::replica_share() { + int error_code = COLVARS_OK; colvarproxy *proxy = cvm::proxy; // sync with the other replicas (if needed) if (comm == multiple_replicas) { // reread the replicas registry - update_replicas_registry(); + error_code |= update_replicas_registry(); // empty the output buffer - std::ostream *replica_hills_os = - proxy->get_output_stream(replica_hills_file); - if (replica_hills_os) { - proxy->flush_output_stream(replica_hills_os); - } - read_replica_files(); + error_code |= proxy->flush_output_stream(replica_hills_file); + error_code |= read_replica_files(); } - return COLVARS_OK; + return error_code; } -void colvarbias_meta::update_replicas_registry() +int colvarbias_meta::update_replicas_registry() { + int error_code = COLVARS_OK; + if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ ": updating the list of replicas, currently containing "+ @@ -1012,8 +1016,9 @@ void colvarbias_meta::update_replicas_registry() while (colvarparse::getline_nocomments(reg_file, line)) replicas_registry.append(line+"\n"); } else { - cvm::error("Error: failed to open file \""+replicas_registry_file+ - "\" for reading.\n", COLVARS_FILE_ERROR); + error_code |= cvm::error("Error: failed to open file \""+ + replicas_registry_file+"\" for reading.\n", + COLVARS_FILE_ERROR); } } @@ -1081,8 +1086,8 @@ void colvarbias_meta::update_replicas_registry() } } } else { - cvm::error("Error: cannot read the replicas registry file \""+ - replicas_registry+"\".\n", COLVARS_FILE_ERROR); + error_code |= cvm::error("Error: cannot read the replicas registry file \""+ + replicas_registry+"\".\n", COLVARS_FILE_ERROR); } // now (re)read the list file of each replica @@ -1122,10 +1127,12 @@ void colvarbias_meta::update_replicas_registry() if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+ cvm::to_str(replicas.size())+" elements.\n"); + + return error_code; } -void colvarbias_meta::read_replica_files() +int colvarbias_meta::read_replica_files() { // Note: we start from the 2nd replica. for (size_t ir = 1; ir < replicas.size(); ir++) { @@ -1250,6 +1257,7 @@ void colvarbias_meta::read_replica_files() " steps. Ensure that it is still running.\n"); } } + return COLVARS_OK; } @@ -1413,6 +1421,10 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) size_t const old_hills_size = hills.size(); hill_iter old_hills_end = hills.end(); hill_iter old_hills_off_grid_end = hills_off_grid.end(); + if (cvm::debug()) { + cvm::log("Before reading hills from the state file, there are "+ + cvm::to_str(hills.size())+" hills in memory.\n"); + } // Read any hills following the grid data (if any) while (read_hill(is)) { @@ -1431,6 +1443,10 @@ std::istream & colvarbias_meta::read_state_data(std::istream& is) if (existing_hills) { hills.erase(hills.begin(), old_hills_end); hills_off_grid.erase(hills_off_grid.begin(), old_hills_off_grid_end); + if (cvm::debug()) { + cvm::log("After pruning the old hills, there are now "+ + cvm::to_str(hills.size())+" hills in memory.\n"); + } } if (rebin_grids) { @@ -1614,6 +1630,8 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) int colvarbias_meta::setup_output() { + int error_code = COLVARS_OK; + output_prefix = cvm::output_prefix(); if (cvm::main()->num_biases_feature(colvardeps::f_cvb_calc_pmf) > 1) { // if this is not the only free energy integrator, append @@ -1626,8 +1644,11 @@ int colvarbias_meta::setup_output() // TODO: one may want to specify the path manually for intricated filesystems? char *pwd = new char[3001]; - if (GETCWD(pwd, 3000) == NULL) - cvm::error("Error: cannot get the path of the current working directory.\n"); + if (GETCWD(pwd, 3000) == NULL) { + return cvm::error("Error: cannot get the path of the current working directory.\n", + COLVARS_BUG_ERROR); + } + replica_list_file = (std::string(pwd)+std::string(PATHSEP)+ this->name+"."+replica_id+".files.txt"); @@ -1680,38 +1701,35 @@ int colvarbias_meta::setup_output() // if we're running without grids, use a growing list of "hills" files // otherwise, just one state file and one "hills" file as buffer - std::ostream *list_os = - cvm::proxy->output_stream(replica_list_file, - (use_grids ? std::ios_base::trunc : - std::ios_base::app)); - if (!list_os) { - return cvm::get_error(); + std::ostream &list_os = cvm::proxy->output_stream(replica_list_file); + if (list_os) { + list_os << "stateFile " << replica_state_file << "\n"; + list_os << "hillsFile " << replica_hills_file << "\n"; + cvm::proxy->close_output_stream(replica_list_file); + } else { + error_code |= COLVARS_FILE_ERROR; } - *list_os << "stateFile " << replica_state_file << "\n"; - *list_os << "hillsFile " << replica_hills_file << "\n"; - cvm::proxy->close_output_stream(replica_list_file); // finally, add a new record for this replica to the registry if (! registered_replica) { - std::ostream *reg_os = - cvm::proxy->output_stream(replicas_registry_file, - std::ios::app); + std::ofstream reg_os(replicas_registry_file.c_str(), std::ios::app); if (!reg_os) { return cvm::get_error(); } - *reg_os << replica_id << " " << replica_list_file << "\n"; + reg_os << replica_id << " " << replica_list_file << "\n"; cvm::proxy->close_output_stream(replicas_registry_file); } } if (b_hills_traj) { + std::ostream &hills_traj_os = + cvm::proxy->output_stream(hills_traj_file_name()); if (!hills_traj_os) { - hills_traj_os = cvm::proxy->output_stream(hills_traj_file_name()); - if (!hills_traj_os) return cvm::get_error(); + error_code |= COLVARS_FILE_ERROR; } } - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return error_code; } @@ -1804,6 +1822,7 @@ int colvarbias_meta::write_output_files() void colvarbias_meta::write_pmf() { + colvarproxy *proxy = cvm::main()->proxy; // allocate a new grid to store the pmf colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy); pmf->setup(); @@ -1820,7 +1839,7 @@ void colvarbias_meta::write_pmf() cvm::real target_val=target_dist->value(i); if (target_val>0) { pmf_val=pmf->value(i); - pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val); + pmf_val=pmf_val + proxy->target_temperature() * proxy->boltzmann() * cvm::logn(target_val); } pmf->set_value(i,pmf_val); } @@ -1830,7 +1849,7 @@ void colvarbias_meta::write_pmf() pmf->add_constant(-1.0 * max); pmf->multiply_constant(-1.0); if (well_tempered) { - cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature; + cvm::real const well_temper_scale = (bias_temperature + proxy->target_temperature()) / bias_temperature; pmf->multiply_constant(well_temper_scale); } { @@ -1839,10 +1858,7 @@ void colvarbias_meta::write_pmf() (dump_fes_save ? "."+cvm::to_str(cvm::step_absolute()) : "") + ".pmf"); - cvm::proxy->backup_file(fes_file_name); - std::ostream *fes_os = cvm::proxy->output_stream(fes_file_name); - pmf->write_multicol(*fes_os); - cvm::proxy->close_output_stream(fes_file_name); + pmf->write_multicol(fes_file_name, "PMF file"); } } @@ -1861,7 +1877,7 @@ void colvarbias_meta::write_pmf() cvm::real target_val=target_dist->value(i); if (target_val>0) { pmf_val=pmf->value(i); - pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val); + pmf_val=pmf_val + proxy->target_temperature() * proxy->boltzmann() * cvm::logn(target_val); } pmf->set_value(i,pmf_val); } @@ -1871,17 +1887,14 @@ void colvarbias_meta::write_pmf() pmf->add_constant(-1.0 * max); pmf->multiply_constant(-1.0); if (well_tempered) { - cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature; + cvm::real const well_temper_scale = (bias_temperature + proxy->target_temperature()) / bias_temperature; pmf->multiply_constant(well_temper_scale); } std::string const fes_file_name(this->output_prefix + (dump_fes_save ? "."+cvm::to_str(cvm::step_absolute()) : "") + ".pmf"); - cvm::proxy->backup_file(fes_file_name); - std::ostream *fes_os = cvm::proxy->output_stream(fes_file_name); - pmf->write_multicol(*fes_os); - cvm::proxy->close_output_stream(fes_file_name); + pmf->write_multicol(fes_file_name, "partial PMF file"); } delete pmf; @@ -1902,9 +1915,9 @@ int colvarbias_meta::write_replica_state_file() // Write to temporary state file std::string const tmp_state_file(replica_state_file+".tmp"); error_code |= proxy->remove_file(tmp_state_file); - std::ostream *rep_state_os = cvm::proxy->output_stream(tmp_state_file); + std::ostream &rep_state_os = cvm::proxy->output_stream(tmp_state_file); if (rep_state_os) { - if (!write_state(*rep_state_os)) { + if (!write_state(rep_state_os)) { error_code |= cvm::error("Error: in writing to temporary file \""+ tmp_state_file+"\".\n", COLVARS_FILE_ERROR); } @@ -1921,13 +1934,13 @@ int colvarbias_meta::reopen_replica_buffer_file() { int error_code = COLVARS_OK; colvarproxy *proxy = cvm::proxy; - if (proxy->get_output_stream(replica_hills_file) != NULL) { + if (proxy->output_stream(replica_hills_file)) { error_code |= proxy->close_output_stream(replica_hills_file); } error_code |= proxy->remove_file(replica_hills_file); - std::ostream *replica_hills_os = proxy->output_stream(replica_hills_file); + std::ostream &replica_hills_os = proxy->output_stream(replica_hills_file); if (replica_hills_os) { - replica_hills_os->setf(std::ios::scientific, std::ios::floatfield); + replica_hills_os.setf(std::ios::scientific, std::ios::floatfield); } else { error_code |= COLVARS_FILE_ERROR; } diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 565534f22b..bac2ff74d4 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -12,8 +12,7 @@ #include #include -#include -#include +#include #include "colvarbias.h" #include "colvargrid.h" @@ -84,9 +83,7 @@ protected: size_t new_hill_freq; /// Write the hill logfile - bool b_hills_traj; - /// Logfile of hill management (creation and deletion) - std::ostream *hills_traj_os; + bool b_hills_traj; /// Name of the hill logfile std::string const hills_traj_file_name() const; @@ -214,10 +211,10 @@ protected: std::string replica_file_name; /// \brief Read the existing replicas on registry - virtual void update_replicas_registry(); + virtual int update_replicas_registry(); /// \brief Read new data from replicas' files - virtual void read_replica_files(); + virtual int read_replica_files(); /// Write full state information to be read by other replicas virtual int write_replica_state_file(); diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 0b374db9a8..a7963c4f8f 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -7,6 +7,10 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include +#include +#include + #include "colvarmodule.h" #include "colvarproxy.h" #include "colvarvalue.h" @@ -196,6 +200,8 @@ int colvarbias_restraint_moving::init(std::string const &conf) if (b_chg_centers || b_chg_force_k) { + first_step = cvm::step_absolute(); + get_keyval(conf, "targetNumSteps", target_nsteps, target_nsteps); if (!target_nsteps) { cvm::error("Error: targetNumSteps must be non-zero.\n", COLVARS_INPUT_ERROR); @@ -335,13 +341,17 @@ int colvarbias_restraint_centers_moving::update_centers(cvm::real lambda) int colvarbias_restraint_centers_moving::update() { + if (!cvm::main()->proxy->simulation_running()) { + return COLVARS_OK; + } + if (b_chg_centers) { if (target_nstages) { // Staged update if (stage <= target_nstages) { if ((cvm::step_relative() > 0) && - ((cvm::step_absolute() % target_nsteps) == 1)) { + (((cvm::step_absolute() - first_step) % target_nsteps) == 1)) { cvm::real const lambda = cvm::real(stage)/cvm::real(target_nstages); update_centers(lambda); @@ -358,9 +368,9 @@ int colvarbias_restraint_centers_moving::update() } } else { // Continuous update - if (cvm::step_absolute() <= target_nsteps) { + if (cvm::step_absolute() - first_step <= target_nsteps) { cvm::real const lambda = - cvm::real(cvm::step_absolute())/cvm::real(target_nsteps); + cvm::real(cvm::step_absolute() - first_step)/cvm::real(target_nsteps); update_centers(lambda); } else { for (size_t i = 0; i < num_variables(); i++) { @@ -389,10 +399,13 @@ int colvarbias_restraint_centers_moving::update() int colvarbias_restraint_centers_moving::update_acc_work() { + if (!cvm::main()->proxy->simulation_running()) { + return COLVARS_OK; + } if (b_chg_centers) { if (is_enabled(f_cvb_output_acc_work)) { if ((cvm::step_relative() > 0) && - (cvm::step_absolute() <= target_nsteps)) { + (cvm::step_absolute() - first_step <= target_nsteps)) { for (size_t i = 0; i < num_variables(); i++) { // project forces on the calculated increments at this step acc_work += colvar_forces[i] * centers_incr[i]; @@ -496,10 +509,11 @@ colvarbias_restraint_k_moving::colvarbias_restraint_k_moving(char const *key) colvarbias_restraint_moving(key) { b_chg_force_k = false; + b_decoupling = false; target_equil_steps = 0; target_force_k = -1.0; starting_force_k = -1.0; - force_k_exp = 1.0; + lambda_exp = 1.0; restraint_FE = 0.0; force_k_incr = 0.0; } @@ -509,18 +523,29 @@ int colvarbias_restraint_k_moving::init(std::string const &conf) { colvarbias_restraint_k::init(conf); + get_keyval(conf, "decoupling", b_decoupling, b_decoupling); + if (b_decoupling) { + starting_force_k = 0.0; + target_force_k = force_k; + b_chg_force_k = true; + } + if (get_keyval(conf, "targetForceConstant", target_force_k, target_force_k)) { + if (b_decoupling) { + cvm::error("Error: targetForceConstant may not be specified together with decoupling.\n", COLVARS_INPUT_ERROR); + return COLVARS_ERROR; + } starting_force_k = force_k; b_chg_force_k = true; } - if (b_chg_force_k) { - // parse moving restraint options - colvarbias_restraint_moving::init(conf); - } else { + if (!b_chg_force_k) { return COLVARS_OK; } + // parse moving restraint options + colvarbias_restraint_moving::init(conf); + get_keyval(conf, "targetEquilSteps", target_equil_steps, target_equil_steps); if (get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule) && @@ -534,12 +559,13 @@ int colvarbias_restraint_k_moving::init(std::string const &conf) target_nstages = lambda_schedule.size() - 1; } - if (get_keyval(conf, "targetForceExponent", force_k_exp, force_k_exp)) { - if (! b_chg_force_k) - cvm::log("Warning: not changing force constant: targetForceExponent will be ignored\n"); + if ((get_keyval(conf, "targetForceExponent", lambda_exp, lambda_exp, parse_deprecated) + || get_keyval(conf, "lambdaExponent", lambda_exp, lambda_exp)) + && !b_chg_force_k) { + cvm::error("Error: cannot set lambdaExponent unless a changing force constant is active.\n", COLVARS_INPUT_ERROR); } - if (force_k_exp < 1.0) { - cvm::log("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n"); + if (lambda_exp < 1.0) { + cvm::log("Warning: for all practical purposes, lambdaExponent should be 1.0 or greater.\n"); } return COLVARS_OK; @@ -548,24 +574,27 @@ int colvarbias_restraint_k_moving::init(std::string const &conf) int colvarbias_restraint_k_moving::update() { + if (!cvm::main()->proxy->simulation_running()) { + return COLVARS_OK; + } if (b_chg_force_k) { cvm::real lambda; if (target_nstages) { - if (cvm::step_absolute() == 0) { + if (cvm::step_absolute() == first_step) { // Setup first stage of staged variable force constant calculation if (lambda_schedule.size()) { lambda = lambda_schedule[0]; } else { - lambda = 0.0; + lambda = (b_decoupling ? 1.0 : 0.0); } force_k = starting_force_k + (target_force_k - starting_force_k) - * cvm::pow(lambda, force_k_exp); + * cvm::pow(lambda, lambda_exp); cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda) - + ", k = " + cvm::to_str(force_k)); + + ", k = " + cvm::to_str(force_k)+"\n"); } // TI calculation: estimate free energy derivative @@ -574,9 +603,10 @@ int colvarbias_restraint_k_moving::update() lambda = lambda_schedule[stage]; } else { lambda = cvm::real(stage) / cvm::real(target_nstages); + if (b_decoupling) lambda = 1.0 - lambda; } - if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) { + if (target_equil_steps == 0 || (cvm::step_absolute() - first_step) % target_nsteps >= target_equil_steps) { // Start averaging after equilibration period, if requested // Derivative of energy with respect to force_k @@ -584,17 +614,17 @@ int colvarbias_restraint_k_moving::update() for (size_t i = 0; i < num_variables(); i++) { dU_dk += d_restraint_potential_dk(i); } - restraint_FE += force_k_exp * cvm::pow(lambda, force_k_exp - 1.0) + restraint_FE += lambda_exp * cvm::pow(lambda, lambda_exp - 1.0) * (target_force_k - starting_force_k) * dU_dk; } // Finish current stage... - if (cvm::step_absolute() % target_nsteps == 0 && - cvm::step_absolute() > 0) { + if ((cvm::step_absolute() - first_step) % target_nsteps == 0 && + cvm::step_absolute() > first_step) { cvm::log("Restraint " + this->name + " Lambda= " + cvm::to_str(lambda) + " dA/dLambda= " - + cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps))); + + cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps))+"\n"); // ...and move on to the next one if (stage < target_nstages) { @@ -605,23 +635,24 @@ int colvarbias_restraint_k_moving::update() lambda = lambda_schedule[stage]; } else { lambda = cvm::real(stage) / cvm::real(target_nstages); + if (b_decoupling) lambda = 1.0 - lambda; } force_k = starting_force_k + (target_force_k - starting_force_k) - * cvm::pow(lambda, force_k_exp); + * cvm::pow(lambda, lambda_exp); cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda) - + ", k = " + cvm::to_str(force_k)); + + ", k = " + cvm::to_str(force_k)+"\n"); } } - } else if (cvm::step_absolute() <= target_nsteps) { - + } else if (cvm::step_absolute() - first_step <= target_nsteps) { // update force constant (slow growth) - lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps); + lambda = cvm::real(cvm::step_absolute() - first_step) / cvm::real(target_nsteps); + if (b_decoupling) lambda = 1.0 - lambda; cvm::real const force_k_old = force_k; force_k = starting_force_k + (target_force_k - starting_force_k) - * cvm::pow(lambda, force_k_exp); + * cvm::pow(lambda, lambda_exp); force_k_incr = force_k - force_k_old; } } @@ -632,6 +663,9 @@ int colvarbias_restraint_k_moving::update() int colvarbias_restraint_k_moving::update_acc_work() { + if (!cvm::main()->proxy->simulation_running()) { + return COLVARS_OK; + } if (b_chg_force_k) { if (is_enabled(f_cvb_output_acc_work)) { if (cvm::step_relative() > 0) { @@ -980,7 +1014,7 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf) } // Initialize starting value of the force constant (in case it's changing) - starting_force_k = force_k; + starting_force_k = (b_decoupling ? 0.0 : force_k); if (lower_walls.size() > 0) { for (i = 0; i < num_variables(); i++) { @@ -1302,6 +1336,8 @@ colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key) int colvarbias_restraint_histogram::init(std::string const &conf) { + int error_code = COLVARS_OK; + colvarbias::init(conf); enable(f_cvb_apply_force); @@ -1312,18 +1348,19 @@ int colvarbias_restraint_histogram::init(std::string const &conf) get_keyval(conf, "width", width, width); if (width <= 0.0) { - cvm::error("Error: \"width\" must be positive.\n", COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: \"width\" must be positive.\n", + COLVARS_INPUT_ERROR); } get_keyval(conf, "gaussianWidth", gaussian_width, 2.0 * width, colvarparse::parse_silent); get_keyval(conf, "gaussianSigma", gaussian_width, 2.0 * width); if (lower_boundary >= upper_boundary) { - cvm::error("Error: the upper boundary, "+ - cvm::to_str(upper_boundary)+ - ", is not higher than the lower boundary, "+ - cvm::to_str(lower_boundary)+".\n", - COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: the upper boundary, "+ + cvm::to_str(upper_boundary)+ + ", is not higher than the lower boundary, "+ + cvm::to_str(lower_boundary)+".\n", + COLVARS_INPUT_ERROR); } cvm::real const nbins = (upper_boundary - lower_boundary) / width; @@ -1347,22 +1384,30 @@ int colvarbias_restraint_histogram::init(std::string const &conf) get_keyval(conf, "refHistogramFile", ref_p_file, std::string("")); if (ref_p_file.size()) { if (inline_ref_p) { - cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n", - COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n", + COLVARS_INPUT_ERROR); } else { - std::ifstream is(ref_p_file.c_str()); + + std::istream &is = + cvm::main()->proxy->input_stream(ref_p_file, + "reference histogram file"); + std::string data_s = ""; std::string line; while (getline_nocomments(is, line)) { data_s.append(line+"\n"); } if (data_s.size() == 0) { - cvm::error("Error: file \""+ref_p_file+"\" empty or unreadable.\n", COLVARS_FILE_ERROR); + error_code |= cvm::error("Error: file \""+ref_p_file+ + "\" empty or unreadable.\n", + COLVARS_FILE_ERROR); } - is.close(); + error_code |= cvm::main()->proxy->close_input_stream(ref_p_file); + cvm::vector1d data; if (data.from_simple_string(data_s) != 0) { - cvm::error("Error: could not read histogram from file \""+ref_p_file+"\".\n"); + error_code |= cvm::error("Error: could not read histogram from file \""+ + ref_p_file+"\".\n"); } if (data.size() == 2*ref_p.size()) { // file contains both x and p(x) @@ -1373,11 +1418,13 @@ int colvarbias_restraint_histogram::init(std::string const &conf) } else if (data.size() == ref_p.size()) { ref_p = data; } else { - cvm::error("Error: file \""+ref_p_file+"\" contains a histogram of different length.\n", - COLVARS_INPUT_ERROR); + error_code |= cvm::error("Error: file \""+ref_p_file+ + "\" contains a histogram of different length.\n", + COLVARS_INPUT_ERROR); } } } + cvm::real const ref_integral = ref_p.sum() * width; if (cvm::fabs(ref_integral - 1.0) > 1.0e-03) { cvm::log("Reference distribution not normalized, normalizing to unity.\n"); @@ -1387,7 +1434,7 @@ int colvarbias_restraint_histogram::init(std::string const &conf) get_keyval(conf, "writeHistogram", b_write_histogram, false); get_keyval(conf, "forceConstant", force_k, 1.0); - return COLVARS_OK; + return error_code; } @@ -1492,27 +1539,29 @@ int colvarbias_restraint_histogram::update() int colvarbias_restraint_histogram::write_output_files() { if (b_write_histogram) { + colvarproxy *proxy = cvm::main()->proxy; std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat"); - std::ostream *os = cvm::proxy->output_stream(file_name); - *os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width) - << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3) - << ")\n"; + std::ostream &os = proxy->output_stream(file_name, + "histogram output file"); + os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width) + << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3) + << ")\n"; - os->setf(std::ios::fixed, std::ios::floatfield); + os.setf(std::ios::fixed, std::ios::floatfield); size_t igrid; for (igrid = 0; igrid < p.size(); igrid++) { cvm::real const x_grid = (lower_boundary + (igrid+1)*width); - *os << " " - << std::setprecision(cvm::cv_prec) - << std::setw(cvm::cv_width) - << x_grid - << " " - << std::setprecision(cvm::cv_prec) - << std::setw(cvm::cv_width) - << p[igrid] << "\n"; + os << " " + << std::setprecision(cvm::cv_prec) + << std::setw(cvm::cv_width) + << x_grid + << " " + << std::setprecision(cvm::cv_prec) + << std::setw(cvm::cv_width) + << p[igrid] << "\n"; } - cvm::proxy->close_output_stream(file_name); + proxy->close_output_stream(file_name); } return COLVARS_OK; } diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index b2c9055dfe..45a96d14f8 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -123,6 +123,9 @@ protected: /// \brief Changing force constant? bool b_chg_force_k; + /// \brief Perform decoupling of the restraint? + bool b_decoupling; + /// \brief Number of stages over which to perform the change /// If zero, perform a continuous change int target_nstages; @@ -137,6 +140,9 @@ protected: /// or restraint centers cvm::step_number target_nsteps; + /// \brief Timestep at which the restraint starts moving + cvm::step_number first_step; + /// \brief Accumulated work (computed when outputAccumulatedWork == true) cvm::real acc_work; }; @@ -207,7 +213,7 @@ protected: cvm::real starting_force_k; /// \brief Exponent for varying the force constant - cvm::real force_k_exp; + cvm::real lambda_exp; /// \brief Intermediate quantity to compute the restraint free energy /// (in TI, would be the accumulating FE derivative) diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 9c93de4b49..87d3e27a8f 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -616,7 +616,7 @@ public: dipole_magnitude (std::string const &conf); dipole_magnitude (cvm::atom const &a1); dipole_magnitude(); - virtual inline ~dipole_magnitude() {} + virtual ~dipole_magnitude() {} virtual void calc_value(); virtual void calc_gradients(); //virtual void calc_force_invgrads(); diff --git a/lib/colvars/colvarcomp_combination.cpp b/lib/colvars/colvarcomp_combination.cpp index 8aea869983..64745b6472 100644 --- a/lib/colvars/colvarcomp_combination.cpp +++ b/lib/colvars/colvarcomp_combination.cpp @@ -141,13 +141,15 @@ void colvar::linearCombination::apply_force(colvarvalue const &force) { colvar::customColvar::customColvar(std::string const &conf): linearCombination(conf) { use_custom_function = false; // code swipe from colvar::init_custom_function -#ifdef LEPTON std::string expr_in, expr; + size_t pos = 0; // current position in config string +#ifdef LEPTON std::vector pexprs; Lepton::ParsedExpression pexpr; double *ref; - size_t pos = 0; // current position in config string +#endif if (key_lookup(conf, "customFunction", &expr_in, &pos)) { +#ifdef LEPTON use_custom_function = true; cvm::log("This colvar uses a custom function.\n"); do { @@ -208,11 +210,15 @@ colvar::customColvar::customColvar(std::string const &conf): linearCombination(c } else { x.type(colvarvalue::type_scalar); } - } else { - cvm::log(std::string{"Warning: no customFunction specified.\n"}); - cvm::log(std::string{"Warning: use linear combination instead.\n"}); - } +#else + cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" + "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", + COLVARS_INPUT_ERROR); #endif + } else { + cvm::log("Warning: no customFunction specified.\n"); + cvm::log("Warning: use linear combination instead.\n"); + } } colvar::customColvar::~customColvar() { @@ -227,96 +233,111 @@ colvar::customColvar::~customColvar() { } void colvar::customColvar::calc_value() { -#ifdef LEPTON - for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { - cv[i_cv]->calc_value(); - } - x.reset(); - size_t l = 0; - for (size_t i = 0; i < x.size(); ++i) { - for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { - const colvarvalue& current_cv_value = cv[i_cv]->value(); - for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { - if (current_cv_value.type() == colvarvalue::type_scalar) { - *(value_eval_var_refs[l++]) = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)); - } else { - *(value_eval_var_refs[l++]) = cv[i_cv]->sup_coeff * current_cv_value[j_elem]; - } - } - } - x[i] = value_evaluators[i]->evaluate(); - } -#endif if (!use_custom_function) { colvar::linearCombination::calc_value(); + } else { +#ifdef LEPTON + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + cv[i_cv]->calc_value(); + } + x.reset(); + size_t l = 0; + for (size_t i = 0; i < x.size(); ++i) { + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + const colvarvalue& current_cv_value = cv[i_cv]->value(); + for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { + if (current_cv_value.type() == colvarvalue::type_scalar) { + *(value_eval_var_refs[l++]) = cv[i_cv]->sup_coeff * (cvm::pow(current_cv_value.real_value, cv[i_cv]->sup_np)); + } else { + *(value_eval_var_refs[l++]) = cv[i_cv]->sup_coeff * current_cv_value[j_elem]; + } + } + } + x[i] = value_evaluators[i]->evaluate(); + } +#else + cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" + "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", + COLVARS_INPUT_ERROR); +#endif } } void colvar::customColvar::calc_gradients() { + if (!use_custom_function) { + colvar::linearCombination::calc_gradients(); + } else { #ifdef LEPTON - size_t r = 0; // index in the vector of variable references - size_t e = 0; // index of the gradient evaluator - for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { // for each CV - cv[i_cv]->calc_gradients(); - if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { - const colvarvalue& current_cv_value = cv[i_cv]->value(); - const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); - for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { // for each element in this CV - for (size_t c = 0; c < x.size(); ++c) { // for each custom function expression - for (size_t k = 0; k < cv.size(); ++k) { // this is required since we need to feed all CV values to this expression - const cvm::real factor_polynomial_k = getPolynomialFactorOfCVGradient(k); - for (size_t l = 0; l < cv[k]->value().size(); ++l) { - *(grad_eval_var_refs[r++]) = factor_polynomial_k * cv[k]->value()[l]; + size_t r = 0; // index in the vector of variable references + size_t e = 0; // index of the gradient evaluator + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { // for each CV + cv[i_cv]->calc_gradients(); + if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { + const colvarvalue& current_cv_value = cv[i_cv]->value(); + const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { // for each element in this CV + for (size_t c = 0; c < x.size(); ++c) { // for each custom function expression + for (size_t k = 0; k < cv.size(); ++k) { // this is required since we need to feed all CV values to this expression + const cvm::real factor_polynomial_k = getPolynomialFactorOfCVGradient(k); + for (size_t l = 0; l < cv[k]->value().size(); ++l) { + *(grad_eval_var_refs[r++]) = factor_polynomial_k * cv[k]->value()[l]; + } } - } - const double expr_grad = gradient_evaluators[e++]->evaluate(); - for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { - for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { - (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = expr_grad * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; + const double expr_grad = gradient_evaluators[e++]->evaluate(); + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + for (size_t l_atom = 0; l_atom < (cv[i_cv]->atom_groups)[k_ag]->size(); ++l_atom) { + (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad = expr_grad * factor_polynomial * (*(cv[i_cv]->atom_groups)[k_ag])[l_atom].grad; + } } } } } } - } +#else + cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" + "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", + COLVARS_INPUT_ERROR); #endif - if (!use_custom_function) { - colvar::linearCombination::calc_gradients(); } } void colvar::customColvar::apply_force(colvarvalue const &force) { -#ifdef LEPTON - size_t r = 0; // index in the vector of variable references - size_t e = 0; // index of the gradient evaluator - for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { - // If this CV us explicit gradients, then atomic gradients is already calculated - // We can apply the force to atom groups directly - if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { - for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { - (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); - } - } else { - const colvarvalue& current_cv_value = cv[i_cv]->value(); - colvarvalue cv_force(current_cv_value.type()); - const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); - for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { - for (size_t c = 0; c < x.size(); ++c) { - for (size_t k = 0; k < cv.size(); ++k) { - const cvm::real factor_polynomial_k = getPolynomialFactorOfCVGradient(k); - for (size_t l = 0; l < cv[k]->value().size(); ++l) { - *(grad_eval_var_refs[r++]) = factor_polynomial_k * cv[k]->value()[l]; - } - } - cv_force[j_elem] += factor_polynomial * gradient_evaluators[e++]->evaluate() * force.real_value; - } - } - cv[i_cv]->apply_force(cv_force); - } - } -#endif if (!use_custom_function) { colvar::linearCombination::apply_force(force); + } else { +#ifdef LEPTON + size_t r = 0; // index in the vector of variable references + size_t e = 0; // index of the gradient evaluator + for (size_t i_cv = 0; i_cv < cv.size(); ++i_cv) { + // If this CV us explicit gradients, then atomic gradients is already calculated + // We can apply the force to atom groups directly + if (cv[i_cv]->is_enabled(f_cvc_explicit_gradient)) { + for (size_t k_ag = 0 ; k_ag < cv[i_cv]->atom_groups.size(); ++k_ag) { + (cv[i_cv]->atom_groups)[k_ag]->apply_colvar_force(force.real_value); + } + } else { + const colvarvalue& current_cv_value = cv[i_cv]->value(); + colvarvalue cv_force(current_cv_value.type()); + const cvm::real factor_polynomial = getPolynomialFactorOfCVGradient(i_cv); + for (size_t j_elem = 0; j_elem < current_cv_value.size(); ++j_elem) { + for (size_t c = 0; c < x.size(); ++c) { + for (size_t k = 0; k < cv.size(); ++k) { + const cvm::real factor_polynomial_k = getPolynomialFactorOfCVGradient(k); + for (size_t l = 0; l < cv[k]->value().size(); ++l) { + *(grad_eval_var_refs[r++]) = factor_polynomial_k * cv[k]->value()[l]; + } + } + cv_force[j_elem] += factor_polynomial * gradient_evaluators[e++]->evaluate() * force.real_value; + } + } + cv[i_cv]->apply_force(cv_force); + } + } +#else + cvm::error("customFunction requires the Lepton library, but it is not enabled during compilation.\n" + "Please refer to the Compilation Notes section of the Colvars manual for more information.\n", + COLVARS_INPUT_ERROR); +#endif } } diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index 80ec3177a9..3d618ff805 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -121,12 +121,12 @@ colvar::coordnum::coordnum(std::string const &conf) } bool const b_isotropic = get_keyval(conf, "cutoff", r0, - cvm::real(4.0 * proxy->angstrom_value)); + cvm::real(proxy->angstrom_to_internal(4.0))); if (get_keyval(conf, "cutoff3", r0_vec, - cvm::rvector(4.0 * proxy->angstrom_value, - 4.0 * proxy->angstrom_value, - 4.0 * proxy->angstrom_value))) { + cvm::rvector(proxy->angstrom_to_internal(4.0), + proxy->angstrom_to_internal(4.0), + proxy->angstrom_to_internal(4.0)))) { if (b_isotropic) { cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" " "at the same time.\n", @@ -328,7 +328,7 @@ colvar::h_bond::h_bond(std::string const &conf) atom_groups[0]->add_atom(acceptor); atom_groups[0]->add_atom(donor); - get_keyval(conf, "cutoff", r0, (3.3 * proxy->angstrom_value)); + get_keyval(conf, "cutoff", r0, proxy->angstrom_to_internal(3.3)); get_keyval(conf, "expNumer", en, 6); get_keyval(conf, "expDenom", ed, 8); @@ -405,7 +405,7 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf) group1 = parse_group(conf, "group1"); - get_keyval(conf, "cutoff", r0, cvm::real(4.0 * proxy->angstrom_value)); + get_keyval(conf, "cutoff", r0, cvm::real(proxy->angstrom_to_internal(4.0))); get_keyval(conf, "expNumer", en, 6); get_keyval(conf, "expDenom", ed, 12); @@ -556,7 +556,7 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf) } bool const b_scale = get_keyval(conf, "cutoff", r0, - cvm::real(4.0 * proxy->angstrom_value)); + cvm::real(proxy->angstrom_to_internal(4.0))); if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0, 4.0, 4.0), parse_silent)) { diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index fa27e5345e..d96cb33482 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -1055,6 +1055,7 @@ colvar::rmsd::rmsd(std::string const &conf) n_permutations = 1; while (key_lookup(conf, "atomPermutation", &perm_conf, &pos)) { + cvm::main()->cite_feature("Symmetry-adapted RMSD"); std::vector perm; if (perm_conf.size()) { std::istringstream is(perm_conf); diff --git a/lib/colvars/colvarcomp_gpath.cpp b/lib/colvars/colvarcomp_gpath.cpp index bff432fdbe..c73ab4f6c1 100644 --- a/lib/colvars/colvarcomp_gpath.cpp +++ b/lib/colvars/colvarcomp_gpath.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include "colvarmodule.h" #include "colvarvalue.h" @@ -426,9 +427,8 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { std::string path_filename; get_keyval(conf, "pathFile", path_filename); cvm::log(std::string("Reading path file: ") + path_filename + std::string("\n")); - std::ifstream ifs_path(path_filename); - if (!ifs_path.is_open()) { - cvm::error("Error: failed to open path file.\n"); + auto &ifs_path = cvm::main()->proxy->input_stream(path_filename); + if (!ifs_path) { return; } std::string line; @@ -459,6 +459,7 @@ colvar::CVBasedPath::CVBasedPath(std::string const &conf): cvc(conf) { ++total_reference_frames; } } + cvm::main()->proxy->close_input_stream(path_filename); if (total_reference_frames <= 1) { cvm::error("Error: there is only 1 or 0 reference frame, which doesn't constitute a path.\n"); return; diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index b92cb533b7..80ef9b855c 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -86,7 +86,7 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) { cvm::real r0; size_t en, ed; - get_keyval(conf, "hBondCutoff", r0, (3.3 * proxy->angstrom_value)); + get_keyval(conf, "hBondCutoff", r0, proxy->angstrom_to_internal(3.3)); get_keyval(conf, "hBondExpNumer", en, 6); get_keyval(conf, "hBondExpDenom", ed, 8); @@ -333,10 +333,11 @@ colvar::dihedPC::dihedPC(std::string const &conf) return; } - std::ifstream vecFile; - vecFile.open(vecFileName.c_str()); - if (!vecFile.good()) { - cvm::error("Error opening dihedral PCA vector file " + vecFileName + " for reading"); + std::istream &vecFile = + cvm::main()->proxy->input_stream(vecFileName, + "dihedral PCA vector file"); + if (!vecFile) { + return; } // TODO: adapt to different formats by setting this flag @@ -375,7 +376,7 @@ colvar::dihedPC::dihedPC(std::string const &conf) } } */ - vecFile.close(); + cvm::main()->proxy->close_input_stream(vecFileName); } else { get_keyval(conf, "vector", coeffs, coeffs); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 2e0d25ab81..f469174a62 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -67,15 +67,17 @@ int colvar::orientation::init(std::string const &conf) } - cvm::log("Centering the reference coordinates: it is " - "assumed that each atom is the closest " - "periodic image to the center of geometry.\n"); cvm::rvector ref_cog(0.0, 0.0, 0.0); size_t i; for (i = 0; i < ref_pos.size(); i++) { ref_cog += ref_pos[i]; } ref_cog /= cvm::real(ref_pos.size()); + cvm::log("Centering the reference coordinates on the origin by subtracting " + "the center of geometry at "+ + cvm::to_str(-1.0 * ref_cog)+"; it is " + "assumed that each atom is the closest " + "periodic image to the center of geometry.\n"); for (i = 0; i < ref_pos.size(); i++) { ref_pos[i] -= ref_cog; } diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index 96eb7cb5e1..2fd385ccec 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -7,14 +7,18 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include +#include + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" #include "colvar.h" #include "colvarcomp.h" #include "colvargrid.h" +#include "colvargrid_def.h" + -#include colvar_grid_count::colvar_grid_count() : colvar_grid() @@ -33,6 +37,42 @@ colvar_grid_count::colvar_grid_count(std::vector &colvars, : colvar_grid(colvars, def_count, 1, margin) {} +std::istream & colvar_grid_count::read_multicol(std::istream &is, bool add) +{ + return colvar_grid::read_multicol(is, add); +} + +int colvar_grid_count::read_multicol(std::string const &filename, + std::string description, + bool add) +{ + return colvar_grid::read_multicol(filename, description, add); +} + +std::ostream & colvar_grid_count::write_multicol(std::ostream &os) const +{ + return colvar_grid::write_multicol(os); +} + +int colvar_grid_count::write_multicol(std::string const &filename, + std::string description) const +{ + return colvar_grid::write_multicol(filename, description); +} + +std::ostream & colvar_grid_count::write_opendx(std::ostream &os) const +{ + return colvar_grid::write_opendx(os); +} + +int colvar_grid_count::write_opendx(std::string const &filename, + std::string description) const +{ + return colvar_grid::write_opendx(filename, description); +} + + + colvar_grid_scalar::colvar_grid_scalar() : colvar_grid(), samples(NULL) {} @@ -56,6 +96,41 @@ colvar_grid_scalar::~colvar_grid_scalar() { } +std::istream & colvar_grid_scalar::read_multicol(std::istream &is, bool add) +{ + return colvar_grid::read_multicol(is, add); +} + +int colvar_grid_scalar::read_multicol(std::string const &filename, + std::string description, + bool add) +{ + return colvar_grid::read_multicol(filename, description, add); +} + +std::ostream & colvar_grid_scalar::write_multicol(std::ostream &os) const +{ + return colvar_grid::write_multicol(os); +} + +int colvar_grid_scalar::write_multicol(std::string const &filename, + std::string description) const +{ + return colvar_grid::write_multicol(filename, description); +} + +std::ostream & colvar_grid_scalar::write_opendx(std::ostream &os) const +{ + return colvar_grid::write_opendx(os); +} + +int colvar_grid_scalar::write_opendx(std::string const &filename, + std::string description) const +{ + return colvar_grid::write_opendx(filename, description); +} + + cvm::real colvar_grid_scalar::maximum_value() const { cvm::real max = data[0]; @@ -145,10 +220,9 @@ colvar_grid_gradient::colvar_grid_gradient(std::string &filename) samples(NULL), weights(NULL) { - std::ifstream is; - is.open(filename.c_str()); - if (!is.is_open()) { - cvm::error("Error opening multicol gradient file " + filename + " for reading.\n"); + std::istream &is = cvm::main()->proxy->input_stream(filename, + "gradient file"); + if (!is) { return; } @@ -203,7 +277,42 @@ colvar_grid_gradient::colvar_grid_gradient(std::string &filename) is.clear(); is.seekg(0); read_multicol(is); - is.close(); + cvm::main()->proxy->close_input_stream(filename); +} + + +std::istream & colvar_grid_gradient::read_multicol(std::istream &is, bool add) +{ + return colvar_grid::read_multicol(is, add); +} + +int colvar_grid_gradient::read_multicol(std::string const &filename, + std::string description, + bool add) +{ + return colvar_grid::read_multicol(filename, description, add); +} + +std::ostream & colvar_grid_gradient::write_multicol(std::ostream &os) const +{ + return colvar_grid::write_multicol(os); +} + +int colvar_grid_gradient::write_multicol(std::string const &filename, + std::string description) const +{ + return colvar_grid::write_multicol(filename, description); +} + +std::ostream & colvar_grid_gradient::write_opendx(std::ostream &os) const +{ + return colvar_grid::write_opendx(os); +} + +int colvar_grid_gradient::write_opendx(std::string const &filename, + std::string description) const +{ + return colvar_grid::write_opendx(filename, description); } diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index f34c5eccab..e2fc1e0fe2 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -732,8 +732,8 @@ 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 const &ix, - size_t const &imult = 0) const + virtual T value_output(std::vector const &ix, + size_t const &imult = 0) const { return value(ix, imult); } @@ -741,10 +741,10 @@ public: /// \brief Get the value from a formatted output and transform it /// into the internal representation (the two may be different, /// e.g. when using colvar_grid_count) - virtual inline void value_input(std::vector const &ix, - T const &t, - size_t const &imult = 0, - bool add = false) + virtual void value_input(std::vector const &ix, + T const &t, + size_t const &imult = 0, + bool add = false) { if ( add ) data[address(ix) + imult] += t; @@ -803,7 +803,7 @@ public: } } - /// \brief Write the grid parameters (number of colvars, boundaries, width and number of points) + /// Write the grid parameters (number of colvars, boundaries, width and number of points) std::ostream & write_params(std::ostream &os) { size_t i; @@ -1028,188 +1028,27 @@ public: return is; } - /// \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) const - { - std::streamsize const w = os.width(); - std::streamsize const p = os.precision(); + /// Read a grid written by write_multicol(), incrementing if add is true + std::istream & read_multicol(std::istream &is, bool add = false); - // Data in the header: nColvars, then for each - // xiMin, dXi, nPoints, periodic + /// Read a grid written by write_multicol(), incrementing if add is true + int read_multicol(std::string const &filename, + std::string description = "grid file", + bool add = false); - os << std::setw(2) << "# " << nd << "\n"; - for (size_t i = 0; i < nd; i++) { - os << "# " - << std::setw(10) << lower_boundaries[i] - << std::setw(10) << widths[i] - << std::setw(10) << nx[i] << " " - << periodic[i] << "\n"; - } + /// Write grid in a format which is both human-readable and gnuplot-friendly + std::ostream & write_multicol(std::ostream &os) const; + /// Write grid in a format which is both human-readable and gnuplot-friendly + int write_multicol(std::string const &filename, + std::string description = "grid file") const; - for (std::vector ix = new_index(); index_ok(ix); incr(ix) ) { + /// Write the grid data without labels, as they are represented in memory + std::ostream & write_opendx(std::ostream &os) const; - if (ix.back() == 0) { - // if the last index is 0, add a new line to mark the new record - os << "\n"; - } - - for (size_t i = 0; i < nd; i++) { - os << " " - << std::setw(w) << std::setprecision(p) - << bin_to_value_scalar(ix[i], i); - } - os << " "; - for (size_t imult = 0; imult < mult; imult++) { - os << " " - << std::setw(w) << std::setprecision(p) - << value_output(ix, imult); - } - os << "\n"; - } - } - - /// \brief Read a grid written by colvar_grid::write_multicol() - /// Adding data if add is true, replacing if false - std::istream & read_multicol(std::istream &is, bool add = false) - { - // Data in the header: nColvars, then for each - // xiMin, dXi, nPoints, periodic flag - - std::string hash; - cvm::real lower, width, x; - size_t n, periodic_flag; - bool remap; - std::vector new_value; - std::vector nx_read; - std::vector bin; - - if ( cv.size() > 0 && cv.size() != nd ) { - cvm::error("Cannot read grid file: number of variables in file differs from number referenced by grid.\n"); - return is; - } - - if ( !(is >> hash) || (hash != "#") ) { - cvm::error("Error reading grid at position "+ - cvm::to_str(static_cast(is.tellg()))+ - " in stream(read \"" + hash + "\")\n"); - return is; - } - - is >> n; - if ( n != nd ) { - cvm::error("Error reading grid: wrong number of collective variables.\n"); - return is; - } - - nx_read.resize(n); - bin.resize(n); - new_value.resize(mult); - - if (this->has_parent_data && add) { - new_data.resize(data.size()); - } - - remap = false; - for (size_t i = 0; i < nd; i++ ) { - if ( !(is >> hash) || (hash != "#") ) { - cvm::error("Error reading grid at position "+ - cvm::to_str(static_cast(is.tellg()))+ - " in stream(read \"" + hash + "\")\n"); - return is; - } - - is >> lower >> width >> nx_read[i] >> periodic_flag; - - - if ( (cvm::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) || - (cvm::fabs(width - widths[i] ) > 1.0e-10) || - (nx_read[i] != nx[i]) ) { - cvm::log("Warning: reading from different grid definition (colvar " - + cvm::to_str(i+1) + "); remapping data on new grid.\n"); - remap = true; - } - } - - if ( remap ) { - // re-grid data - while (is.good()) { - bool end_of_file = false; - - for (size_t i = 0; i < nd; i++ ) { - if ( !(is >> x) ) end_of_file = true; - bin[i] = value_to_bin_scalar(x, i); - } - if (end_of_file) break; - - for (size_t imult = 0; imult < mult; imult++) { - is >> new_value[imult]; - } - - if ( index_ok(bin) ) { - for (size_t imult = 0; imult < mult; imult++) { - value_input(bin, new_value[imult], imult, add); - } - } - } - } else { - // do not re-grid the data but assume the same grid is used - for (std::vector ix = new_index(); index_ok(ix); incr(ix) ) { - for (size_t i = 0; i < nd; i++ ) { - is >> x; - } - for (size_t imult = 0; imult < mult; imult++) { - is >> new_value[imult]; - value_input(ix, new_value[imult], imult, add); - } - } - } - has_data = true; - return is; - } - - /// \brief Write the grid data without labels, as they are - /// represented in memory - std::ostream & write_opendx(std::ostream &os) const - { - // write the header - os << "object 1 class gridpositions counts"; - size_t icv; - for (icv = 0; icv < num_variables(); icv++) { - os << " " << number_of_points(icv); - } - os << "\n"; - - os << "origin"; - for (icv = 0; icv < num_variables(); icv++) { - os << " " << (lower_boundaries[icv].real_value + 0.5 * widths[icv]); - } - os << "\n"; - - for (icv = 0; icv < num_variables(); icv++) { - os << "delta"; - for (size_t icv2 = 0; icv2 < num_variables(); icv2++) { - if (icv == icv2) os << " " << widths[icv]; - else os << " " << 0.0; - } - os << "\n"; - } - - os << "object 2 class gridconnections counts"; - for (icv = 0; icv < num_variables(); icv++) { - os << " " << number_of_points(icv); - } - os << "\n"; - - os << "object 3 class array type double rank 0 items " - << number_of_points() << " data follows\n"; - - write_raw(os); - - os << "object \"collective variables scalar field\" class field\n"; - return os; - } + /// Write the grid data without labels, as they are represented in memory + int write_opendx(std::string const &filename, + std::string description = "grid file") const; }; @@ -1224,7 +1063,7 @@ public: colvar_grid_count(); /// Destructor - virtual inline ~colvar_grid_count() + virtual ~colvar_grid_count() {} /// Constructor @@ -1249,13 +1088,33 @@ public: return new_data[address(ix) + imult]; } - /// \brief Get the value from a formatted output and transform it - /// into the internal representation (it may have been rescaled or - /// manipulated) - virtual inline void value_input(std::vector const &ix, - size_t const &t, - size_t const &imult = 0, - bool add = false) + /// Read a grid written by write_multicol(), incrementin if data is true + std::istream & read_multicol(std::istream &is, bool add = false); + + /// Read a grid written by write_multicol(), incrementin if data is true + int read_multicol(std::string const &filename, + std::string description = "grid file", + bool add = false); + + /// Write grid in a format which is both human-readable and gnuplot-friendly + std::ostream & write_multicol(std::ostream &os) const; + + /// Write grid in a format which is both human-readable and gnuplot-friendly + int write_multicol(std::string const &filename, + std::string description = "grid file") const; + + /// Write the grid data without labels, as they are represented in memory + std::ostream & write_opendx(std::ostream &os) const; + + /// Write the grid data without labels, as they are represented in memory + int write_opendx(std::string const &filename, + std::string description = "grid file") const; + + /// Enter or add a value, but also handle parent grid + virtual void value_input(std::vector const &ix, + size_t const &t, + size_t const &imult = 0, + bool add = false) { (void) imult; if (add) { @@ -1322,7 +1181,7 @@ public: /// \brief Return the gradient of discrete count from finite differences /// on the *same* grid for dimension n inline cvm::real gradient_finite_diff(const std::vector &ix0, - int n = 0) + int n = 0) { cvm::real A0, A1, A2; std::vector ix = ix0; @@ -1380,7 +1239,7 @@ public: colvar_grid_scalar(colvar_grid_scalar const &g); /// Destructor - ~colvar_grid_scalar(); + virtual ~colvar_grid_scalar(); /// Constructor from specific sizes arrays colvar_grid_scalar(std::vector const &nx_i); @@ -1402,6 +1261,28 @@ public: has_data = true; } + /// Read a grid written by write_multicol(), incrementin if data is true + std::istream & read_multicol(std::istream &is, bool add = false); + + /// Read a grid written by write_multicol(), incrementin if data is true + int read_multicol(std::string const &filename, + std::string description = "grid file", + bool add = false); + + /// Write grid in a format which is both human-readable and gnuplot-friendly + std::ostream & write_multicol(std::ostream &os) const; + + /// Write grid in a format which is both human-readable and gnuplot-friendly + int write_multicol(std::string const &filename, + std::string description = "grid file") const; + + /// Write the grid data without labels, as they are represented in memory + std::ostream & write_opendx(std::ostream &os) const; + + /// Write the grid data without labels, as they are represented in memory + int write_opendx(std::string const &filename, + std::string description = "grid file") const; + /// \brief Return the gradient of the scalar field from finite differences /// Input coordinates are those of gradient grid, shifted wrt scalar grid /// Should not be called on edges of scalar grid, provided the latter has margins @@ -1458,7 +1339,7 @@ public: /// \brief Return the gradient of discrete count from finite differences /// on the *same* grid for dimension n inline cvm::real gradient_finite_diff(const std::vector &ix0, - int n = 0) + int n = 0) { cvm::real A0, A1, A2; std::vector ix = ix0; @@ -1517,9 +1398,7 @@ public: } } - /// \brief Get the value from a formatted output and transform it - /// into the internal representation (it may have been rescaled or - /// manipulated) + /// Enter or add value but also deal with count grid virtual void value_input(std::vector const &ix, cvm::real const &new_value, size_t const &imult = 0, @@ -1580,7 +1459,7 @@ public: colvar_grid_gradient(); /// Destructor - virtual inline ~colvar_grid_gradient() + virtual ~colvar_grid_gradient() {} /// Constructor from specific sizes arrays @@ -1592,6 +1471,28 @@ public: /// Constructor from a multicol file colvar_grid_gradient(std::string &filename); + /// Read a grid written by write_multicol(), incrementin if data is true + virtual std::istream & read_multicol(std::istream &is, bool add = false); + + /// Read a grid written by write_multicol(), incrementin if data is true + virtual int read_multicol(std::string const &filename, + std::string description = "grid file", + bool add = false); + + /// Write grid in a format which is both human-readable and gnuplot-friendly + virtual std::ostream & write_multicol(std::ostream &os) const; + + /// Write grid in a format which is both human-readable and gnuplot-friendly + virtual int write_multicol(std::string const &filename, + std::string description = "grid file") const; + + /// Write the grid data without labels, as they are represented in memory + virtual std::ostream & write_opendx(std::ostream &os) const; + + /// Write the grid data without labels, as they are represented in memory + virtual int write_opendx(std::string const &filename, + std::string description = "grid file") const; + /// \brief Get a vector with the binned value(s) indexed by ix, normalized if applicable inline void vector_value(std::vector const &ix, std::vector &v) const { @@ -1647,8 +1548,8 @@ 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 const &ix, - size_t const &imult = 0) const + virtual cvm::real value_output(std::vector const &ix, + size_t const &imult = 0) const { if (samples) return (samples->value(ix) > 0) ? @@ -1661,10 +1562,10 @@ public: /// \brief Get the value from a formatted output and transform it /// into the internal representation (it may have been rescaled or /// manipulated) - virtual inline void value_input(std::vector const &ix, - cvm::real const &new_value, - size_t const &imult = 0, - bool add = false) + virtual void value_input(std::vector const &ix, + cvm::real const &new_value, + size_t const &imult = 0, + bool add = false) { if (add) { if (samples) diff --git a/lib/colvars/colvargrid_def.h b/lib/colvars/colvargrid_def.h new file mode 100644 index 0000000000..f2245f3d81 --- /dev/null +++ b/lib/colvars/colvargrid_def.h @@ -0,0 +1,261 @@ +// -*- 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. + +/// \file Definition of the more complex members of colvar_grid<> template + +#ifndef COLVARGRID_DEF_H +#define COLVARGRID_DEF_H + +#include +#include + +#include "colvarmodule.h" +#include "colvarproxy.h" +#include "colvar.h" +#include "colvargrid.h" + + +template +std::istream & colvar_grid::read_multicol(std::istream &is, bool add) +{ + // Data in the header: nColvars, then for each + // xiMin, dXi, nPoints, periodic flag + + std::string hash; + cvm::real lower, width, x; + size_t n, periodic_flag; + bool remap; + std::vector new_value; + std::vector nx_read; + std::vector bin; + + if ( cv.size() > 0 && cv.size() != nd ) { + cvm::error("Cannot read grid file: number of variables in file differs from number referenced by grid.\n"); + return is; + } + + if ( !(is >> hash) || (hash != "#") ) { + cvm::error("Error reading grid at position "+ + cvm::to_str(static_cast(is.tellg()))+ + " in stream(read \"" + hash + "\")\n", COLVARS_INPUT_ERROR); + return is; + } + + is >> n; + if ( n != nd ) { + cvm::error("Error reading grid: wrong number of collective variables.\n"); + return is; + } + + nx_read.resize(n); + bin.resize(n); + new_value.resize(mult); + + if (this->has_parent_data && add) { + new_data.resize(data.size()); + } + + remap = false; + for (size_t i = 0; i < nd; i++ ) { + if ( !(is >> hash) || (hash != "#") ) { + cvm::error("Error reading grid at position "+ + cvm::to_str(static_cast(is.tellg()))+ + " in stream(read \"" + hash + "\")\n"); + return is; + } + + is >> lower >> width >> nx_read[i] >> periodic_flag; + + + if ( (cvm::fabs(lower - lower_boundaries[i].real_value) > 1.0e-10) || + (cvm::fabs(width - widths[i] ) > 1.0e-10) || + (nx_read[i] != nx[i]) ) { + cvm::log("Warning: reading from different grid definition (colvar " + + cvm::to_str(i+1) + "); remapping data on new grid.\n"); + remap = true; + } + } + + if ( remap ) { + // re-grid data + while (is.good()) { + bool end_of_file = false; + + for (size_t i = 0; i < nd; i++ ) { + if ( !(is >> x) ) end_of_file = true; + bin[i] = value_to_bin_scalar(x, i); + } + if (end_of_file) break; + + for (size_t imult = 0; imult < mult; imult++) { + is >> new_value[imult]; + } + + if ( index_ok(bin) ) { + for (size_t imult = 0; imult < mult; imult++) { + value_input(bin, new_value[imult], imult, add); + } + } + } + } else { + // do not re-grid the data but assume the same grid is used + for (std::vector ix = new_index(); index_ok(ix); incr(ix) ) { + for (size_t i = 0; i < nd; i++ ) { + is >> x; + } + for (size_t imult = 0; imult < mult; imult++) { + is >> new_value[imult]; + value_input(ix, new_value[imult], imult, add); + } + } + } + has_data = true; + return is; +} + + +template +int colvar_grid::read_multicol(std::string const &filename, + std::string description, + bool add) +{ + std::istream &is = cvm::main()->proxy->input_stream(filename, description); + if (!is) { + return COLVARS_FILE_ERROR; + } + if (colvar_grid::read_multicol(is, add)) { + cvm::main()->proxy->close_input_stream(filename); + return COLVARS_OK; + } + return COLVARS_FILE_ERROR; +} + + +template +std::ostream & colvar_grid::write_multicol(std::ostream &os) const +{ + // Save the output formats + std::ios_base::fmtflags prev_flags(os.flags()); + + // Data in the header: nColvars, then for each + // xiMin, dXi, nPoints, periodic + + os << std::setw(2) << "# " << nd << "\n"; + // Write the floating numbers in full precision + os.setf(std::ios::scientific, std::ios::floatfield); + for (size_t i = 0; i < nd; i++) { + os << "# " + << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << lower_boundaries[i] << " " + << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) << widths[i] << " " + << std::setw(10) << nx[i] << " " + << periodic[i] << "\n"; + } + + for (std::vector ix = new_index(); index_ok(ix); incr(ix) ) { + + if (ix.back() == 0) { + // if the last index is 0, add a new line to mark the new record + os << "\n"; + } + + for (size_t i = 0; i < nd; i++) { + os << " " + << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) + << bin_to_value_scalar(ix[i], i); + } + os << " "; + for (size_t imult = 0; imult < mult; imult++) { + os << " " + << std::setw(cvm::cv_width) << std::setprecision(cvm::cv_prec) + << value_output(ix, imult); + } + os << "\n"; + } + + // Restore the output formats + os.flags(prev_flags); + + return os; +} + + +template +int colvar_grid::write_multicol(std::string const &filename, + std::string description) const +{ + int error_code = COLVARS_OK; + std::ostream &os = cvm::main()->proxy->output_stream(filename, description); + if (!os) { + return COLVARS_FILE_ERROR; + } + error_code |= colvar_grid::write_multicol(os) ? COLVARS_OK : + COLVARS_FILE_ERROR; + cvm::main()->proxy->close_output_stream(filename); + return error_code; +} + + +template +std::ostream & colvar_grid::write_opendx(std::ostream &os) const +{ + // write the header + os << "object 1 class gridpositions counts"; + size_t icv; + for (icv = 0; icv < num_variables(); icv++) { + os << " " << number_of_points(icv); + } + os << "\n"; + + os << "origin"; + for (icv = 0; icv < num_variables(); icv++) { + os << " " << (lower_boundaries[icv].real_value + 0.5 * widths[icv]); + } + os << "\n"; + + for (icv = 0; icv < num_variables(); icv++) { + os << "delta"; + for (size_t icv2 = 0; icv2 < num_variables(); icv2++) { + if (icv == icv2) os << " " << widths[icv]; + else os << " " << 0.0; + } + os << "\n"; + } + + os << "object 2 class gridconnections counts"; + for (icv = 0; icv < num_variables(); icv++) { + os << " " << number_of_points(icv); + } + os << "\n"; + + os << "object 3 class array type double rank 0 items " + << number_of_points() << " data follows\n"; + + write_raw(os); + + os << "object \"collective variables scalar field\" class field\n"; + return os; +} + + +template +int colvar_grid::write_opendx(std::string const &filename, + std::string description) const +{ + int error_code = COLVARS_OK; + std::ostream &os = cvm::main()->proxy->output_stream(filename, description); + if (!os) { + return COLVARS_FILE_ERROR; + } + error_code |= colvar_grid::write_opendx(os) ? COLVARS_OK : + COLVARS_FILE_ERROR; + cvm::main()->proxy->close_output_stream(filename); + return error_code; +} + +#endif diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 6ad945202e..68cd402e1c 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -7,7 +7,10 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. +#include +#include #include +#include #include #include #include @@ -71,8 +74,6 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) depth_s = 0; log_level_ = 10; - cv_traj_os = NULL; - xyz_reader_use_count = 0; num_biases_types_used_ = @@ -132,8 +133,6 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) cv_traj_freq = 100; restart_out_freq = proxy->default_restart_frequency(); - // by default overwrite the existing trajectory file - cv_traj_append = false; cv_traj_write_labels = true; // Removes the need for proxy specializations to create this @@ -190,12 +189,12 @@ int colvarmodule::read_config_file(char const *config_filename) std::string(config_filename)+"\":\n"); // open the configfile - config_s.open(config_filename); - if (!config_s.is_open()) { - cvm::error("Error: in opening configuration file \""+ - std::string(config_filename)+"\".\n", - COLVARS_FILE_ERROR); - return COLVARS_ERROR; + std::istream &config_s = proxy->input_stream(config_filename, + "configuration file/string"); + if (!config_s) { + return cvm::error("Error: in opening configuration file \""+ + std::string(config_filename)+"\".\n", + COLVARS_FILE_ERROR); } // read the config file into a string @@ -206,7 +205,7 @@ int colvarmodule::read_config_file(char const *config_filename) if (line.find_first_not_of(colvarparse::white_space) != std::string::npos) conf.append(line+"\n"); } - config_s.close(); + proxy->close_input_stream(config_filename); return parse_config(conf); } @@ -299,9 +298,6 @@ int colvarmodule::parse_config(std::string &conf) cvm::log("Collective variables module (re)initialized.\n"); cvm::log(cvm::line_marker); - // Update any necessary proxy data - proxy->setup(); - if (source_Tcl_script.size() > 0) { run_tcl_script(source_Tcl_script); } @@ -385,10 +381,6 @@ int colvarmodule::parse_global_params(std::string const &conf) parse->get_keyval(conf, "colvarsRestartFrequency", restart_out_freq, restart_out_freq); - // Deprecate append flag - parse->get_keyval(conf, "colvarsTrajAppend", - cv_traj_append, cv_traj_append, colvarparse::parse_silent); - parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, use_scripted_forces); @@ -433,6 +425,7 @@ int colvarmodule::parse_colvars(std::string const &conf) cvm::log("Error while constructing colvar number " + cvm::to_str(colvars.size()) + " : deleting."); delete colvars.back(); // the colvar destructor updates the colvars array + cvm::decrease_depth(); return COLVARS_ERROR; } cvm::decrease_depth(); @@ -855,11 +848,13 @@ int colvarmodule::calc_colvars() // so they can activate colvars as needed std::vector::iterator bi; for (bi = biases.begin(); bi != biases.end(); bi++) { - int tsf = (*bi)->get_time_step_factor(); - if (tsf > 0 && (step_absolute() % tsf == 0)) { - (*bi)->enable(colvardeps::f_cvb_awake); - } else { - (*bi)->disable(colvardeps::f_cvb_awake); + int const tsf = (*bi)->get_time_step_factor(); + if (tsf > 1) { + if (step_absolute() % tsf == 0) { + (*bi)->enable(colvardeps::f_cvb_awake); + } else { + (*bi)->disable(colvardeps::f_cvb_awake); + } } } @@ -870,12 +865,14 @@ int colvarmodule::calc_colvars() variables_active()->clear(); variables_active()->reserve(variables()->size()); for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) { - // Wake up or put to sleep variables + // Wake up or put to sleep variables with MTS int tsf = (*cvi)->get_time_step_factor(); - if (tsf > 0 && (step_absolute() % tsf == 0)) { - (*cvi)->enable(colvardeps::f_cv_awake); - } else { - (*cvi)->disable(colvardeps::f_cv_awake); + if (tsf > 1) { + if (step_absolute() % tsf == 0) { + (*cvi)->enable(colvardeps::f_cv_awake); + } else { + (*cvi)->disable(colvardeps::f_cv_awake); + } } if ((*cvi)->is_enabled()) { @@ -1013,10 +1010,7 @@ int colvarmodule::update_colvar_forces() cvm::log("Collecting forces from all biases.\n"); cvm::increase_depth(); for (bi = biases_active()->begin(); bi != biases_active()->end(); bi++) { - (*bi)->communicate_forces(); - if (cvm::get_error()) { - return COLVARS_ERROR; - } + error_code |= (*bi)->communicate_forces(); } cvm::decrease_depth(); @@ -1040,9 +1034,6 @@ int colvarmodule::update_colvar_forces() for (cvi = variables()->begin(); cvi != variables()->end(); cvi++) { // Inactive colvars will only reset their forces and return 0 energy total_colvar_energy += (*cvi)->update_forces_energy(); - if (cvm::get_error()) { - return COLVARS_ERROR; - } } cvm::decrease_depth(); if (cvm::debug()) @@ -1064,7 +1055,7 @@ int colvarmodule::update_colvar_forces() } cvm::decrease_depth(); - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return error_code; } @@ -1089,17 +1080,15 @@ int colvarmodule::calc_scripted_forces() int colvarmodule::write_restart_file(std::string const &out_name) { cvm::log("Saving collective variables state to \""+out_name+"\".\n"); - proxy->backup_file(out_name); - std::ostream *restart_out_os = proxy->output_stream(out_name); - if (!restart_out_os) return cvm::get_error(); - if (!write_restart(*restart_out_os)) { + std::ostream &restart_out_os = proxy->output_stream(out_name, "state file"); + if (!restart_out_os) return COLVARS_FILE_ERROR; + if (!write_restart(restart_out_os)) { return cvm::error("Error: in writing restart file.\n", COLVARS_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); - } + + // Take the opportunity to flush colvars.traj + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -1118,38 +1107,46 @@ int colvarmodule::write_restart_string(std::string &output) int colvarmodule::write_traj_files() { + int error_code = COLVARS_OK; + 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(); - } else { - cv_traj_write_labels = true; - } + std::ostream &cv_traj_os = proxy->output_stream(cv_traj_name, + "colvars trajectory"); + + if (!cv_traj_os) { + return COLVARS_FILE_ERROR; } - // write labels in the traj file every 1000 lines and at first timestep - if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || - cvm::step_relative() == 0 || - cv_traj_write_labels) { - write_traj_label(*cv_traj_os); + // Write labels in the traj file at beginning and then every 1000 lines + if ( (cvm::step_relative() == 0) || cv_traj_write_labels || + ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0) ) { + error_code |= + write_traj_label(cv_traj_os) ? COLVARS_OK : COLVARS_FILE_ERROR; + cv_traj_write_labels = false; + } + + if (cvm::debug()) { + proxy->flush_output_stream(cv_traj_name); } - cv_traj_write_labels = false; if ((cvm::step_absolute() % cv_traj_freq) == 0) { - write_traj(*cv_traj_os); + error_code |= write_traj(cv_traj_os) ? COLVARS_OK : COLVARS_FILE_ERROR; } - if (restart_out_freq && (cv_traj_os != NULL) && - ((cvm::step_absolute() % restart_out_freq) == 0)) { + if (cvm::debug()) { + proxy->flush_output_stream(cv_traj_name); + } + + if (restart_out_freq && ((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); + error_code |= proxy->flush_output_stream(cv_traj_name); } - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return error_code; } @@ -1208,7 +1205,7 @@ int colvarmodule::end_of_step() } -int colvarmodule::setup() +int colvarmodule::update_engine_parameters() { if (this->size() == 0) return cvm::get_error(); for (std::vector::iterator cvi = variables()->begin(); @@ -1221,8 +1218,8 @@ int colvarmodule::setup() colvarmodule::~colvarmodule() { - if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) || - (proxy->smp_thread_id() == 0)) { + if ((proxy->smp_thread_id() < 0) || // not using threads + (proxy->smp_thread_id() == 0)) { // or this is thread 0 reset(); @@ -1250,8 +1247,6 @@ colvarmodule::~colvarmodule() int colvarmodule::reset() { - cvm::log("Resetting the Collective Variables module.\n"); - parse->clear(); // Iterate backwards because we are deleting the elements as we go @@ -1289,30 +1284,34 @@ int colvarmodule::setup_input() // 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 without the suffix - input_is.clear(); + std::istream *input_is = &(proxy->input_stream(restart_in_name, + "restart file/channel", + false)); + if (!*input_is) { + // Try without the suffix ".colvars.state" restart_in_name = proxy->input_prefix(); - input_is.open(restart_in_name.c_str()); + input_is = &(proxy->input_stream(restart_in_name, + "restart file/channel")); + if (!*input_is) { + return COLVARS_FILE_ERROR; + } } - // Now that the file has been opened, clear this for the next round + // Now that the file has been opened, clear this field so that this + // function will not be called twice proxy->input_prefix().clear(); - if (!input_is.good()) { - return cvm::error("Error: in opening input state file \""+ - std::string(restart_in_name)+"\".\n", - COLVARS_FILE_ERROR); - } else { - cvm::log(cvm::line_marker); - cvm::log("Loading state from file \""+restart_in_name+"\".\n"); - read_restart(input_is); - cvm::log(cvm::line_marker); - return cvm::get_error(); - } + cvm::log(cvm::line_marker); + cvm::log("Loading state from file \""+restart_in_name+"\".\n"); + read_restart(*input_is); + cvm::log(cvm::line_marker); + + proxy->close_input_stream(restart_in_name); + + return cvm::get_error(); } + // TODO This could soon be redundant if (proxy->input_buffer() != NULL) { // Read a string buffer char const *buffer = proxy->input_buffer(); @@ -1364,10 +1363,6 @@ int colvarmodule::setup_output() std::string(output_prefix()+".colvars.traj") : std::string("")); - if (cv_traj_freq && cv_traj_name.size()) { - error_code |= open_traj_file(cv_traj_name); - } - for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { @@ -1572,7 +1567,7 @@ int colvarmodule::write_output_files() error_code |= (*bi)->write_state_to_replicas(); } cvm::decrease_depth(); - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return error_code; } @@ -1582,6 +1577,9 @@ int colvarmodule::read_traj(char const *traj_filename, { cvm::log("Opening trajectory file \""+ std::string(traj_filename)+"\".\n"); + // NB: this function is not currently used, but when it will it should + // retain the ability for direct file-based access (in case traj files + // exceed memory) std::ifstream traj_is(traj_filename); while (true) { @@ -1683,42 +1681,6 @@ std::ostream & colvarmodule::write_restart(std::ostream &os) } -int colvarmodule::open_traj_file(std::string const &file_name) -{ - if (cv_traj_os != NULL) { - return COLVARS_OK; - } - - // (re)open trajectory file - if (cv_traj_append) { - cvm::log("Appending to trajectory file \""+file_name+"\".\n"); - cv_traj_os = (cvm::proxy)->output_stream(file_name, std::ios::app); - } else { - cvm::log("Opening trajectory file \""+file_name+"\".\n"); - proxy->backup_file(file_name.c_str()); - cv_traj_os = (cvm::proxy)->output_stream(file_name); - } - - if (cv_traj_os == NULL) { - cvm::error("Error: cannot write to file \""+file_name+"\".\n", - COLVARS_FILE_ERROR); - } - - return cvm::get_error(); -} - - -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; - } - return cvm::get_error(); -} - - std::ostream & colvarmodule::write_traj_label(std::ostream &os) { os.setf(std::ios::scientific, std::ios::floatfield); @@ -1739,10 +1701,6 @@ std::ostream & colvarmodule::write_traj_label(std::ostream &os) } os << "\n"; - if (cvm::debug()) { - proxy->flush_output_stream(&os); - } - cvm::decrease_depth(); return os; } @@ -1768,10 +1726,6 @@ std::ostream & colvarmodule::write_traj(std::ostream &os) } os << "\n"; - if (cvm::debug()) { - proxy->flush_output_stream(&os); - } - cvm::decrease_depth(); return os; } @@ -1780,12 +1734,15 @@ std::ostream & colvarmodule::write_traj(std::ostream &os) void colvarmodule::log(std::string const &message, int min_log_level) { if (cvm::log_level() < min_log_level) return; + + std::string const trailing_newline = (message.size() > 0) ? + (message[message.size()-1] == '\n' ? "" : "\n") : ""; // allow logging when the module is not fully initialized size_t const d = (cvm::main() != NULL) ? depth() : 0; if (d > 0) { - proxy->log((std::string(2*d, ' '))+message); + proxy->log((std::string(2*d, ' ')) + message + trailing_newline); } else { - proxy->log(message); + proxy->log(message + trailing_newline); } } @@ -1854,18 +1811,26 @@ void colvarmodule::clear_error() int colvarmodule::error(std::string const &message, int code) { set_error_bits(code); - proxy->error(message); + + std::string const trailing_newline = (message.size() > 0) ? + (message[message.size()-1] == '\n' ? "" : "\n") : ""; + size_t const d = depth(); + if (d > 0) { + proxy->error((std::string(2*d, ' ')) + message + trailing_newline); + } else { + proxy->error(message + trailing_newline); + } + return get_error(); } int cvm::read_index_file(char const *filename) { - std::ifstream is(filename, std::ios::binary); - if (!is.good()) { - return cvm::error("Error: in opening index file \""+ - std::string(filename)+"\".\n", - COLVARS_FILE_ERROR); + std::istream &is = proxy->input_stream(filename, "index file"); + + if (!is) { + return COLVARS_FILE_ERROR; } else { index_file_names.push_back(std::string(filename)); } @@ -1947,7 +1912,7 @@ int cvm::read_index_file(char const *filename) cvm::to_str((index_groups[i])->size())+" atoms)\n"); } - return COLVARS_OK; + return proxy->close_input_stream(filename); } @@ -2016,9 +1981,10 @@ int cvm::load_coords(char const *file_name, int cvm::load_coords_xyz(char const *filename, std::vector *pos, - cvm::atom_group *atoms) + cvm::atom_group *atoms, + bool keep_open) { - std::ifstream xyz_is(filename); + std::istream &xyz_is = proxy->input_stream(filename, "XYZ file"); unsigned int natoms; char symbol[256]; std::string line; @@ -2028,7 +1994,8 @@ int cvm::load_coords_xyz(char const *filename, std::string(filename)+"\".\n"); if ( ! (xyz_is >> natoms) ) { - return cvm::error(error_msg, COLVARS_INPUT_ERROR); + // Return silent error when reaching the end of multi-frame files + return keep_open ? COLVARS_NO_SUCH_FRAME : cvm::error(error_msg, COLVARS_INPUT_ERROR); } ++xyz_reader_use_count; @@ -2049,6 +2016,12 @@ int cvm::load_coords_xyz(char const *filename, size_t xyz_natoms = 0; if (pos->size() != natoms) { // Use specified indices int next = 0; // indices are zero-based + if (!atoms) { + // In the other branch of this test, reading all positions from the file, + // a valid atom group pointer is not necessary + return cvm::error("Trying to read partial positions with invalid atom group pointer", + COLVARS_BUG_ERROR); + } std::vector::const_iterator index = atoms->sorted_ids().begin(); for ( ; pos_i != pos->end() ; pos_i++, index++) { @@ -2092,7 +2065,11 @@ int cvm::load_coords_xyz(char const *filename, cvm::to_str(pos->size())+".\n", COLVARS_INPUT_ERROR); } - return COLVARS_OK; + if (keep_open) { + return COLVARS_OK; + } else { + return proxy->close_input_stream(filename); + } } @@ -2100,18 +2077,6 @@ int cvm::load_coords_xyz(char const *filename, // Wrappers to proxy functions: these may go in the future -cvm::real cvm::boltzmann() -{ - return proxy->boltzmann(); -} - - -cvm::real cvm::temperature() -{ - return proxy->temperature(); -} - - cvm::real cvm::dt() { return proxy->dt(); diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index ae29e94505..236d432a95 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -40,15 +40,13 @@ You can browse the class hierarchy or the list of source files. #define COLVARS_BUG_ERROR (1<<3) // Inconsistent state indicating bug #define COLVARS_FILE_ERROR (1<<4) #define COLVARS_MEMORY_ERROR (1<<5) -#define COLVARS_NO_SUCH_FRAME (1<<6) // Cannot load the requested frame +#define COLVARS_NO_SUCH_FRAME (1<<6) // Cannot load the requested frame -#include -#include -#include #include #include #include #include +#include class colvarparse; class colvar; @@ -441,9 +439,9 @@ public: // from the proxy that may change during program execution) // No additional parsing is done within these functions - /// (Re)initialize internal data (currently used by LAMMPS) + /// (Re)initialize any internal data affected by changes in the engine /// Also calls setup() member functions of colvars and biases - int setup(); + int update_engine_parameters(); /// (Re)initialize and (re)read the input state file calling read_restart() int setup_input(); @@ -634,12 +632,6 @@ public: // proxy functions - /// \brief Boltmann constant - static real boltzmann(); - - /// \brief Temperature of the simulation (K) - static real temperature(); - /// \brief Time step of MD integrator (fs) static real dt(); @@ -751,7 +743,8 @@ public: /// Load coordinates into an atom group from an XYZ file (assumes Angstroms) int load_coords_xyz(char const *filename, std::vector *pos, - atom_group *atoms); + atom_group *atoms, + bool keep_open = false); /// Frequency for collective variables trajectory output static size_t cv_traj_freq; @@ -766,21 +759,12 @@ public: protected: - /// Configuration file - std::ifstream config_s; - /// Configuration file parser object colvarparse *parse; /// Name of the trajectory file std::string cv_traj_name; - /// Collective variables output trajectory file - std::ostream *cv_traj_os; - - /// Appending to the existing trajectory file? - bool cv_traj_append; - /// Write labels at the next iteration bool cv_traj_write_labels; diff --git a/lib/colvars/colvarmodule_refs.h b/lib/colvars/colvarmodule_refs.h index 0b8791ce34..c3fccf297c 100644 --- a/lib/colvars/colvarmodule_refs.h +++ b/lib/colvars/colvarmodule_refs.h @@ -63,6 +63,22 @@ " url = {https://doi.org/10.1021/ct500874p}\n" "}\n"; + paper_count_[std::string("Ebrahimi2022")] = 0; + paper_url_[std::string("Ebrahimi2022")] = "https://doi.org/10.1021/acs.jctc.1c01235"; + paper_bibtex_[std::string("Ebrahimi2022")] = + "\n" + "@article {Ebrahimi2022,\n" + " author = {Ebrahimi, Mina and H\\'enin, J\\'er\\^ome},\n" + " title = {Symmetry-Adapted Restraints for Binding Free Energy Calculations},\n" + " journal = {Journal of Chemical Theory and Computation},\n" + " volume = {18},\n" + " number = {4},\n" + " pages = {2494-2502},\n" + " year = {2022},\n" + " doi = {10.1021/acs.jctc.1c01235},\n" + " url = {https://doi.org/10.1021/acs.jctc.1c01235}\n" + "}\n"; + paper_count_[std::string("Fiorin2013")] = 0; paper_url_[std::string("Fiorin2013")] = "https://doi.org/10.1080/00268976.2013.813594"; paper_bibtex_[std::string("Fiorin2013")] = @@ -246,20 +262,19 @@ " url = {https://doi.org/10.1063/5.0014475}\n" "}\n"; - paper_count_[std::string("Plimpton1995")] = 0; - paper_url_[std::string("Plimpton1995")] = "https://doi.org/10.1006/jcph.1995.1039"; - paper_bibtex_[std::string("Plimpton1995")] = + paper_count_[std::string("Thompson2022")] = 0; + paper_url_[std::string("Thompson2022")] = "https://doi.org/10.1016/j.cpc.2021.108171"; + paper_bibtex_[std::string("Thompson2022")] = "\n" - "@article{Plimpton1995,\n" - " title = {Fast parallel algorithms for short-range molecular dynamics},\n" - " author = {Plimpton, Steve},\n" - " journal = {J. Comp. Phys.},\n" - " year = {1995},\n" - " volume = {117},\n" - " number = {1},\n" - " pages = {1--19},\n" - " doi = {10.1006/jcph.1995.1039},\n" - " url = {https://doi.org/10.1006/jcph.1995.1039}\n" + "@article{Thompson2022,\n" + " title = {{LAMMPS} - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales},\n" + " author = {Thompson, Aidan P. and Aktulga, H. Metin and Berger, Richard and Bolintineanu, Dan S. and Brown, W. Michael and Crozier, Paul S. and {in't Veld}, Pieter J. and Kohlmeyer, Axel and Moore, Stan G. and Nguyen, Trung Dac and Shan, Ray and Stevens, Mark J. and Tranchida, Julien and Trott, Christian and Plimpton, Steven J.},\n" + " journal = {Comp. Phys. Comm.},\n" + " volume = {271},\n" + " pages = {108171},\n" + " year = {2022},\n" + " doi = {10.1016/j.cpc.2021.108171},\n" + " url = {https://doi.org/10.1016/j.cpc.2021.108171}\n" "}\n"; paper_count_[std::string("Shen2015")] = 0; @@ -329,6 +344,9 @@ feature_count_[std::string("Multiple-walker ABF implementation")] = 0; feature_paper_map_[std::string("Multiple-walker ABF implementation")] = "Comer2014c"; + feature_count_[std::string("Symmetry-adapted RMSD")] = 0; + feature_paper_map_[std::string("Symmetry-adapted RMSD")] = "Ebrahimi2022"; + feature_count_[std::string("Colvars module")] = 0; feature_paper_map_[std::string("Colvars module")] = "Fiorin2013"; @@ -504,7 +522,7 @@ feature_paper_map_[std::string("Scalable center-of-mass computation (NAMD)")] = "Phillips2020"; feature_count_[std::string("LAMMPS engine")] = 0; - feature_paper_map_[std::string("LAMMPS engine")] = "Plimpton1995"; + feature_paper_map_[std::string("LAMMPS engine")] = "Thompson2022"; feature_count_[std::string("distancePairs colvar component")] = 0; feature_paper_map_[std::string("distancePairs colvar component")] = "Shen2015"; diff --git a/lib/colvars/colvarproxy.cpp b/lib/colvars/colvarproxy.cpp index 3d4c9c439b..b09ea6667c 100644 --- a/lib/colvars/colvarproxy.cpp +++ b/lib/colvars/colvarproxy.cpp @@ -7,23 +7,9 @@ // If you wish to distribute your changes, please submit them to the // Colvars repository at GitHub. -// Using access() to check if a file exists (until we can assume C++14/17) -#if !defined(_WIN32) || defined(__CYGWIN__) -#include -#endif -#if defined(_WIN32) -#include -#endif - -#include - -#include -#include -#include - -#if defined(_OPENMP) -#include -#endif +#include +#include +#include #include "colvarmodule.h" #include "colvarproxy.h" @@ -33,205 +19,6 @@ -colvarproxy_system::colvarproxy_system() -{ - angstrom_value = 0.0; - kcal_mol_value = 0.0; - boundaries_type = boundaries_unsupported; - total_force_requested = false; - indirect_lambda_biasing_force = 0.0; - cached_alch_lambda_changed = false; - cached_alch_lambda = -1.0; - reset_pbc_lattice(); -} - - -colvarproxy_system::~colvarproxy_system() {} - - -int colvarproxy_system::set_unit_system(std::string const & /* units */, - bool /* check_only */) -{ - return COLVARS_NOT_IMPLEMENTED; -} - - -cvm::real colvarproxy_system::backend_angstrom_value() -{ - return 1.0; -} - - -cvm::real colvarproxy_system::boltzmann() -{ - return 0.001987191; -} - - -cvm::real colvarproxy_system::temperature() -{ - // TODO define, document and implement a user method to set the value of this - return 300.0; -} - - -cvm::real colvarproxy_system::dt() -{ - // TODO define, document and implement a user method to set the value of this - return 1.0; -} - - -cvm::real colvarproxy_system::rand_gaussian() -{ - // TODO define, document and implement a user method to set the value of this - return 0.0; -} - - -void colvarproxy_system::add_energy(cvm::real /* energy */) {} - - -void colvarproxy_system::request_total_force(bool yesno) -{ - if (yesno == true) - cvm::error("Error: total forces are currently not implemented.\n", - COLVARS_NOT_IMPLEMENTED); -} - - -bool colvarproxy_system::total_forces_enabled() const -{ - return false; -} - - -bool colvarproxy_system::total_forces_same_step() const -{ - return false; -} - - -inline int round_to_integer(cvm::real x) -{ - return int(cvm::floor(x+0.5)); -} - - -void colvarproxy_system::update_pbc_lattice() -{ - // Periodicity is assumed in all directions - - if (boundaries_type == boundaries_unsupported || - boundaries_type == boundaries_non_periodic) { - cvm::error("Error: setting PBC lattice with unsupported boundaries.\n", - COLVARS_BUG_ERROR); - return; - } - - { - cvm::rvector const v = cvm::rvector::outer(unit_cell_y, unit_cell_z); - reciprocal_cell_x = v/(v*unit_cell_x); - } - { - cvm::rvector const v = cvm::rvector::outer(unit_cell_z, unit_cell_x); - reciprocal_cell_y = v/(v*unit_cell_y); - } - { - cvm::rvector const v = cvm::rvector::outer(unit_cell_x, unit_cell_y); - reciprocal_cell_z = v/(v*unit_cell_z); - } -} - - -void colvarproxy_system::reset_pbc_lattice() -{ - unit_cell_x.reset(); - unit_cell_y.reset(); - unit_cell_z.reset(); - reciprocal_cell_x.reset(); - reciprocal_cell_y.reset(); - reciprocal_cell_z.reset(); -} - - -cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) - const -{ - if (boundaries_type == boundaries_unsupported) { - cvm::error("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR); - } - - cvm::rvector diff = (pos2 - pos1); - - if (boundaries_type == boundaries_non_periodic) return diff; - - cvm::real const x_shift = round_to_integer(reciprocal_cell_x*diff); - cvm::real const y_shift = round_to_integer(reciprocal_cell_y*diff); - cvm::real const z_shift = round_to_integer(reciprocal_cell_z*diff); - - diff.x -= x_shift*unit_cell_x.x + y_shift*unit_cell_y.x + - z_shift*unit_cell_z.x; - diff.y -= x_shift*unit_cell_x.y + y_shift*unit_cell_y.y + - z_shift*unit_cell_z.y; - diff.z -= x_shift*unit_cell_x.z + y_shift*unit_cell_y.z + - z_shift*unit_cell_z.z; - - return diff; -} - - -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; -} - - -int colvarproxy_system::get_alch_lambda(cvm::real * /* lambda */) -{ - return cvm::error("Error in get_alch_lambda: alchemical lambda dynamics is not supported by this build.", - COLVARS_NOT_IMPLEMENTED); -} - - -void colvarproxy_system::set_alch_lambda(cvm::real lambda) -{ - cached_alch_lambda = lambda; - cached_alch_lambda_changed = true; -} - - -int colvarproxy_system::send_alch_lambda() -{ - return cvm::error("Error in set_alch_lambda: alchemical lambda dynamics is not supported by this build.", - COLVARS_NOT_IMPLEMENTED); -} - - -int colvarproxy_system::get_dE_dlambda(cvm::real * /* force */) -{ - return cvm::error("Error in get_dE_dlambda: alchemical lambda dynamics is not supported by this build.", - COLVARS_NOT_IMPLEMENTED); -} - - -int colvarproxy_system::apply_force_dE_dlambda(cvm::real* /* force */) -{ - return cvm::error("Error in apply_force_dE_dlambda: function is not implemented by this build.", - COLVARS_NOT_IMPLEMENTED); -} - - -int colvarproxy_system::get_d2E_dlambda2(cvm::real*) -{ - return cvm::error("Error in get_d2E_dlambda2: function is not implemented by this build.", - COLVARS_NOT_IMPLEMENTED); -} - - colvarproxy_atoms::colvarproxy_atoms() { atoms_rms_applied_force_ = atoms_max_applied_force_ = 0.0; @@ -249,7 +36,7 @@ colvarproxy_atoms::~colvarproxy_atoms() int colvarproxy_atoms::reset() { atoms_ids.clear(); - atoms_ncopies.clear(); + atoms_refcount.clear(); atoms_masses.clear(); atoms_charges.clear(); atoms_positions.clear(); @@ -262,7 +49,7 @@ int colvarproxy_atoms::reset() int colvarproxy_atoms::add_atom_slot(int atom_id) { atoms_ids.push_back(atom_id); - atoms_ncopies.push_back(1); + atoms_refcount.push_back(1); atoms_masses.push_back(1.0); atoms_charges.push_back(0.0); atoms_positions.push_back(cvm::rvector(0.0, 0.0, 0.0)); @@ -309,12 +96,22 @@ void colvarproxy_atoms::clear_atom(int index) cvm::error("Error: trying to disable an atom that was not previously requested.\n", COLVARS_INPUT_ERROR); } - if (atoms_ncopies[index] > 0) { - atoms_ncopies[index] -= 1; + if (atoms_refcount[index] > 0) { + atoms_refcount[index] -= 1; } } +size_t colvarproxy_atoms::get_num_active_atoms() const +{ + size_t result = 0; + for (size_t i = 0; i < atoms_refcount.size(); i++) { + if (atoms_refcount[i] > 0) result++; + } + return result; +} + + int colvarproxy_atoms::load_atoms(char const * /* filename */, cvm::atom_group & /* atoms */, std::string const & /* pdb_field */, @@ -382,7 +179,7 @@ colvarproxy_atom_groups::~colvarproxy_atom_groups() int colvarproxy_atom_groups::reset() { atom_groups_ids.clear(); - atom_groups_ncopies.clear(); + atom_groups_refcount.clear(); atom_groups_masses.clear(); atom_groups_charges.clear(); atom_groups_coms.clear(); @@ -395,7 +192,7 @@ int colvarproxy_atom_groups::reset() int colvarproxy_atom_groups::add_atom_group_slot(int atom_group_id) { atom_groups_ids.push_back(atom_group_id); - atom_groups_ncopies.push_back(1); + atom_groups_refcount.push_back(1); atom_groups_masses.push_back(1.0); atom_groups_charges.push_back(0.0); atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0)); @@ -427,12 +224,22 @@ void colvarproxy_atom_groups::clear_atom_group(int index) "that was not previously requested.\n", COLVARS_INPUT_ERROR); } - if (atom_groups_ncopies[index] > 0) { - atom_groups_ncopies[index] -= 1; + if (atom_groups_refcount[index] > 0) { + atom_groups_refcount[index] -= 1; } } +size_t colvarproxy_atom_groups::get_num_active_atom_groups() const +{ + size_t result = 0; + for (size_t i = 0; i < atom_groups_refcount.size(); i++) { + if (atom_groups_refcount[i] > 0) result++; + } + return result; +} + + void colvarproxy_atom_groups::compute_rms_atom_groups_applied_force() { atom_groups_rms_applied_force_ = @@ -454,8 +261,8 @@ colvarproxy_smp::colvarproxy_smp() omp_lock_state = NULL; #if defined(_OPENMP) if (omp_get_thread_num() == 0) { - omp_lock_state = reinterpret_cast(new omp_lock_t); - omp_init_lock(reinterpret_cast(omp_lock_state)); + omp_lock_state = new omp_lock_t; + omp_init_lock(omp_lock_state); } #endif } @@ -466,7 +273,7 @@ colvarproxy_smp::~colvarproxy_smp() #if defined(_OPENMP) if (omp_get_thread_num() == 0) { if (omp_lock_state) { - delete reinterpret_cast(omp_lock_state); + delete omp_lock_state; } } #endif @@ -567,7 +374,7 @@ int colvarproxy_smp::smp_thread_id() #if defined(_OPENMP) return omp_get_thread_num(); #else - return COLVARS_NOT_IMPLEMENTED; + return -1; #endif } @@ -577,7 +384,7 @@ int colvarproxy_smp::smp_num_threads() #if defined(_OPENMP) return omp_get_max_threads(); #else - return COLVARS_NOT_IMPLEMENTED; + return -1; #endif } @@ -585,7 +392,7 @@ int colvarproxy_smp::smp_num_threads() int colvarproxy_smp::smp_lock() { #if defined(_OPENMP) - omp_set_lock(reinterpret_cast(omp_lock_state)); + omp_set_lock(omp_lock_state); #endif return COLVARS_OK; } @@ -594,8 +401,7 @@ int colvarproxy_smp::smp_lock() int colvarproxy_smp::smp_trylock() { #if defined(_OPENMP) - return omp_test_lock(reinterpret_cast(omp_lock_state)) ? - COLVARS_OK : COLVARS_ERROR; + return omp_test_lock(omp_lock_state) ? COLVARS_OK : COLVARS_ERROR; #else return COLVARS_OK; #endif @@ -605,7 +411,7 @@ int colvarproxy_smp::smp_trylock() int colvarproxy_smp::smp_unlock() { #if defined(_OPENMP) - omp_unset_lock(reinterpret_cast(omp_lock_state)); + omp_unset_lock(omp_lock_state); #endif return COLVARS_OK; } @@ -651,151 +457,35 @@ int colvarproxy_script::run_colvar_gradient_callback(std::string const & /* name -colvarproxy_io::colvarproxy_io() -{ - input_buffer_ = NULL; - restart_frequency_engine = 0; -} - - -colvarproxy_io::~colvarproxy_io() {} - - -int colvarproxy_io::get_frame(long int&) -{ - return COLVARS_NOT_IMPLEMENTED; -} - - -int colvarproxy_io::set_frame(long int) -{ - return COLVARS_NOT_IMPLEMENTED; -} - - -int colvarproxy_io::backup_file(char const *filename) -{ - // Simplified version of NAMD_file_exists() - int exit_code; - do { -#if defined(_WIN32) && !defined(__CYGWIN__) - // We could use _access_s here, but it is probably too new - exit_code = _access(filename, 00); -#else - exit_code = access(filename, F_OK); -#endif - } while ((exit_code != 0) && (errno == EINTR)); - if (exit_code != 0) { - if (errno == ENOENT) { - // File does not exist - return COLVARS_OK; - } else { - return cvm::error("Unknown error while checking if file \""+ - std::string(filename)+"\" exists.\n", COLVARS_ERROR); - } - } - - // The file exists, then rename it - if (std::string(filename).rfind(std::string(".colvars.state")) != - std::string::npos) { - return rename_file(filename, (std::string(filename)+".old").c_str()); - } else { - return rename_file(filename, (std::string(filename)+".BAK").c_str()); - } -} - - -int colvarproxy_io::remove_file(char const *filename) -{ - int error_code = COLVARS_OK; -#if defined(_WIN32) && !defined(__CYGWIN__) - // Because the file may be open by other processes, rename it to filename.old - std::string const renamed_file(std::string(filename)+".old"); - // It may still be there from an interrupted run, so remove it to be safe - std::remove(renamed_file.c_str()); - int rename_exit_code = 0; - while ((rename_exit_code = std::rename(filename, - renamed_file.c_str())) != 0) { - if (errno == EINTR) continue; - error_code |= COLVARS_FILE_ERROR; - break; - } - // Ask to remove filename.old, but ignore any errors raised - std::remove(renamed_file.c_str()); -#else - if (std::remove(filename)) { - if (errno != ENOENT) { - error_code |= COLVARS_FILE_ERROR; - } - } -#endif - if (error_code != COLVARS_OK) { - return cvm::error("Error: in removing file \""+std::string(filename)+ - "\".\n.", - error_code); - } - return COLVARS_OK; -} - - -int colvarproxy_io::rename_file(char const *filename, char const *newfilename) -{ - int error_code = COLVARS_OK; -#if defined(_WIN32) && !defined(__CYGWIN__) - // On straight Windows, must remove the destination before renaming it - error_code |= remove_file(newfilename); -#endif - int rename_exit_code = 0; - while ((rename_exit_code = std::rename(filename, newfilename)) != 0) { - if (errno == EINTR) continue; - // Call log() instead of error to allow the next try - cvm::log("Error: in renaming file \""+std::string(filename)+"\" to \""+ - std::string(newfilename)+"\".\n."); - error_code |= COLVARS_FILE_ERROR; - if (errno == EXDEV) continue; - break; - } - return rename_exit_code ? error_code : COLVARS_OK; -} - - - colvarproxy::colvarproxy() { colvars = NULL; + // By default, simulation engines allow to immediately request atoms + engine_ready_ = true; b_simulation_running = true; b_simulation_continuing = false; b_delete_requested = false; version_int = -1; features_hash = 0; + config_queue_ = reinterpret_cast(new std::list >); } colvarproxy::~colvarproxy() { - close_files(); + close_output_streams(); if (colvars != NULL) { delete colvars; colvars = NULL; } + delete reinterpret_cast > *>(config_queue_); } -int colvarproxy::close_files() +bool colvarproxy::io_available() { - if (smp_enabled() == COLVARS_OK && smp_thread_id() > 0) { - // Nothing to do on non-master threads - return COLVARS_OK; - } - std::list::iterator osni = output_stream_names.begin(); - std::list::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; + return (smp_enabled() == COLVARS_OK && smp_thread_id() == 0) || + (smp_enabled() != COLVARS_OK); } @@ -816,12 +506,41 @@ int colvarproxy::request_deletion() } +void colvarproxy::add_config(std::string const &cmd, std::string const &conf) +{ + reinterpret_cast > *>(config_queue_)->push_back(std::make_pair(cmd, conf)); +} + + int colvarproxy::setup() { return COLVARS_OK; } +int colvarproxy::parse_module_config() +{ + int error_code = COLVARS_OK; + // Read any configuration queued up for Colvars + std::list > *config_queue = reinterpret_cast > *>(config_queue_); + while (config_queue->size() > 0) { + std::pair const &p = config_queue->front(); + if (p.first == "config") { + error_code |= colvars->read_config_string(p.second); + } else if (p.first == "configfile") { + error_code |= colvars->read_config_file(p.second.c_str()); + } else { + error_code |= cvm::error(std::string("Error: invalid keyword \"") + + p.first + + std::string("\" in colvarproxy::setup()\n"), + COLVARS_BUG_ERROR); + } + config_queue->pop_front(); + } + return error_code; +} + + int colvarproxy::update_input() { return COLVARS_OK; @@ -874,7 +593,7 @@ void colvarproxy::print_input_atomic_data() cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ "atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ - "atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + "atoms_refcount = "+cvm::to_str(atoms_refcount)+"\n"); cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ "atoms_masses = "+cvm::to_str(atoms_masses)+"\n"); cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ @@ -893,7 +612,7 @@ void colvarproxy::print_input_atomic_data() cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ "atom_groups_ids = "+cvm::to_str(atom_groups_ids)+"\n"); cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ - "atom_groups_ncopies = "+cvm::to_str(atom_groups_ncopies)+"\n"); + "atom_groups_refcount = "+cvm::to_str(atom_groups_refcount)+"\n"); cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ "atom_groups_masses = "+cvm::to_str(atom_groups_masses)+"\n"); cvm::log("Step "+cvm::to_str(cvm::step_absolute())+", "+ @@ -988,103 +707,3 @@ int colvarproxy::get_version_from_string(char const *version_string) } -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", COLVARS_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", - COLVARS_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::iterator osi = output_files.begin(); - std::list::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::iterator osi = output_files.begin(); - std::list::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", COLVARS_BUG_ERROR); -} - - -int colvarproxy::flush_output_streams() -{ - if (smp_enabled() == COLVARS_OK && smp_thread_id() > 0) - return COLVARS_OK; - - std::list::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::iterator osi = output_files.begin(); - std::list::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", COLVARS_BUG_ERROR); -} diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 3ed3768aef..5733b8b66d 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -10,12 +10,11 @@ #ifndef COLVARPROXY_H #define COLVARPROXY_H -#include -#include - #include "colvarmodule.h" #include "colvartypes.h" #include "colvarvalue.h" +#include "colvarproxy_io.h" +#include "colvarproxy_system.h" #include "colvarproxy_tcl.h" #include "colvarproxy_volmaps.h" @@ -38,165 +37,6 @@ class colvarscript; -/// Methods for accessing the simulation system (PBCs, integrator, etc) -class colvarproxy_system { - -public: - - /// Constructor - colvarproxy_system(); - - /// Destructor - virtual ~colvarproxy_system(); - - /// \brief Name of the unit system used internally by Colvars (by default, that of the back-end). - /// Supported depending on the back-end: real (A, kcal/mol), metal (A, eV), electron (Bohr, Hartree), gromacs (nm, kJ/mol) - /// Note: calls to back-end PBC functions assume back-end length unit - /// We use different unit from back-end in VMD bc using PBC functions from colvarproxy base class - /// Colvars internal units are user specified, because the module exchanges info in unknown - /// composite dimensions with user input, while it only exchanges quantities of known - /// dimension with the back-end (length and forces) - std::string units; - - /// \brief Request to set the units used internally by Colvars - virtual int set_unit_system(std::string const &units, bool check_only); - - /// \brief Value of 1 Angstrom in the internal (front-end) Colvars unit for atomic coordinates - /// * defaults to 0. in the base class; derived proxy classes must set it - /// * in VMD proxy, can only be changed when no variables are defined - /// as user-defined values in composite units must be compatible with that system - cvm::real angstrom_value; - - /// \brief Value of 1 Angstrom in the backend's unit for atomic coordinates - virtual cvm::real backend_angstrom_value(); - - /// \brief Value of 1 kcal/mol in the internal Colvars unit for energy - cvm::real kcal_mol_value; - - /// \brief Convert a length from Angstrom to internal - inline cvm::real angstrom_to_internal(cvm::real l) const - { - return l * angstrom_value; - } - - /// \brief Convert a length from internal to Angstrom - inline cvm::real internal_to_angstrom(cvm::real l) const - { - return l / angstrom_value; - } - - // /// \brief Convert a length from back-end unit to internal - // inline cvm::real back_end_to_internal_unit(cvm::real l) { - // if (angstrom_value == 0.) { - // return l / backend_angstrom_value(); - // } - // return l * angstrom_value / backend_angstrom_value(); - // } - - /// \brief Boltzmann constant in internal Colvars units - virtual cvm::real boltzmann(); - - /// \brief Target temperature of the simulation (K units) - virtual cvm::real temperature(); - - /// \brief Time step of the simulation (fs) - virtual cvm::real dt(); - - /// \brief Pseudo-random number with Gaussian distribution - virtual cvm::real rand_gaussian(void); - - /// Pass restraint energy value for current timestep to MD engine - virtual void add_energy(cvm::real energy); - - /// \brief Get the PBC-aware distance vector between two positions - virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) const; - - /// Recompute PBC reciprocal lattice (assumes XYZ periodicity) - void update_pbc_lattice(); - - /// Set the lattice vectors to zero - void reset_pbc_lattice(); - - /// \brief Tell the proxy whether total forces are needed (they may not - /// always be available) - virtual void request_total_force(bool yesno); - - /// Are total forces being used? - virtual bool total_forces_enabled() const; - - /// 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); - - /// Get value of alchemical lambda parameter from back-end (if available) - virtual int get_alch_lambda(cvm::real* lambda); - - /// Set value of alchemical lambda parameter to be sent to back-end at end of timestep - void set_alch_lambda(cvm::real lambda); - - /// Send cached value of alchemical lambda parameter to back-end (if available) - virtual int send_alch_lambda(); - - /// Get energy derivative with respect to lambda (if available) - virtual int get_dE_dlambda(cvm::real* dE_dlambda); - - /// Apply a scalar force on dE_dlambda (back-end distributes it onto atoms) - virtual int apply_force_dE_dlambda(cvm::real* force); - - /// Get energy second derivative with respect to lambda (if available) - virtual int get_d2E_dlambda2(cvm::real* d2E_dlambda2); - - /// Force to be applied onto alch. lambda, propagated from biasing forces on dE_dlambda - cvm::real indirect_lambda_biasing_force; - - /// Get weight factor from accelMD - virtual cvm::real get_accelMD_factor() const { - cvm::error("Error: accessing the reweighting factor of accelerated MD " - "is not yet implemented in the MD engine.\n", - COLVARS_NOT_IMPLEMENTED); - return 1.0; - } - virtual bool accelMD_enabled() const { - return false; - } - -protected: - /// Next value of lambda to be sent to back-end - cvm::real cached_alch_lambda; - - /// Whether lambda has been set and needs to be updated in backend - bool cached_alch_lambda_changed; - - /// 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. - /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS) - /// minimum-image distances are computed by the host engine regardless. - enum Boundaries_type { - boundaries_non_periodic, - boundaries_pbc_ortho, - boundaries_pbc_triclinic, - boundaries_unsupported - }; - - /// Type of boundary conditions - Boundaries_type boundaries_type; - - /// Bravais lattice vectors - cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z; - - /// Reciprocal lattice vectors - cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z; -}; - - /// \brief Container of atomic data for processing by Colvars class colvarproxy_atoms { @@ -229,7 +69,7 @@ public: std::string const &segment_id); /// \brief Used by the atom class destructor: rather than deleting the array slot - /// (costly) set the corresponding atoms_ncopies to zero + /// (costly) set the corresponding atoms_refcount to zero virtual void clear_atom(int index); /// \brief Select atom IDs from a file (usually PDB) \param filename name of @@ -260,37 +100,51 @@ public: /// Clear atomic data int reset(); - /// Get the numeric ID of the given atom (for the program) + /// Get the numeric ID of the given atom + /// \param index Internal index in the Colvars arrays inline int get_atom_id(int index) const { return atoms_ids[index]; } /// Get the mass of the given atom + /// \param index Internal index in the Colvars arrays inline cvm::real get_atom_mass(int index) const { return atoms_masses[index]; } + /// Increase the reference count of the given atom + /// \param index Internal index in the Colvars arrays + inline void increase_refcount(int index) + { + atoms_refcount[index] += 1; + } + /// Get the charge of the given atom + /// \param index Internal index in the Colvars arrays inline cvm::real get_atom_charge(int index) const { return atoms_charges[index]; } /// Read the current position of the given atom + /// \param index Internal index in the Colvars arrays inline cvm::rvector get_atom_position(int index) const { return atoms_positions[index]; } /// Read the current total force of the given atom + /// \param index Internal index in the Colvars arrays inline cvm::rvector get_atom_total_force(int index) const { return atoms_total_forces[index]; } /// Request that this force is applied to the given atom + /// \param index Internal index in the Colvars arrays + /// \param new_force Force to add inline void apply_atom_force(int index, cvm::rvector const &new_force) { atoms_new_colvar_forces[index] += new_force; @@ -310,6 +164,9 @@ public: return &atoms_ids; } + /// Return number of atoms with positive reference count + size_t get_num_active_atoms() const; + inline std::vector const *get_atom_masses() const { return &atoms_masses; @@ -406,7 +263,7 @@ protected: /// within the host program std::vector atoms_ids; /// \brief Keep track of how many times each atom is used by a separate colvar object - std::vector atoms_ncopies; + std::vector atoms_refcount; /// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) std::vector atoms_masses; /// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) @@ -510,6 +367,9 @@ public: return &atom_groups_ids; } + /// Return number of atom groups with positive reference count + size_t get_num_active_atom_groups() const; + inline std::vector *modify_atom_group_masses() { // TODO updated_masses @@ -561,7 +421,7 @@ protected: /// within the host program std::vector atom_groups_ids; /// \brief Keep track of how many times each group is used by a separate cvc - std::vector atom_groups_ncopies; + std::vector atom_groups_refcount; /// \brief Total masses of the atom groups std::vector atom_groups_masses; /// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS) @@ -584,6 +444,12 @@ protected: }; +#if defined(_OPENMP) +#include +#else +struct omp_lock_t; +#endif + /// \brief Methods for SMP parallelization class colvarproxy_smp { @@ -629,7 +495,7 @@ public: protected: /// Lock state for OpenMP - void *omp_lock_state; + omp_lock_t *omp_lock_state; }; @@ -698,108 +564,6 @@ public: }; -/// Methods for data input/output -class colvarproxy_io { - -public: - - /// Constructor - colvarproxy_io(); - - /// Destructor - virtual ~colvarproxy_io(); - - /// \brief Save the current frame number in the argument given - // Returns error code - virtual int get_frame(long int &); - - /// \brief Set the current frame number (as well as colvarmodule::it) - // Returns error code - virtual int set_frame(long int); - - /// \brief Rename the given file, before overwriting it - virtual int backup_file(char const *filename); - - /// \brief Rename the given file, before overwriting it - inline int backup_file(std::string const &filename) - { - return backup_file(filename.c_str()); - } - - /// Remove the given file (on Windows only, rename to filename.old) - virtual int remove_file(char const *filename); - - /// Remove the given file (on Windows only, rename to filename.old) - inline int remove_file(std::string const &filename) - { - return remove_file(filename.c_str()); - } - - /// Rename the given file - virtual int rename_file(char const *filename, char const *newfilename); - - /// Rename the given file - inline int rename_file(std::string const &filename, - std::string const &newfilename) - { - return rename_file(filename.c_str(), newfilename.c_str()); - } - - /// Prefix of the input state file to be read next - inline std::string & input_prefix() - { - return input_prefix_str; - } - - /// 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: - - /// 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 - std::list output_files; - /// \brief Identifiers for output_stream objects: by default, these are the names of the files - std::list output_stream_names; - - /// Buffer from which the input state information may be read - char const *input_buffer_; -}; - - /// \brief Interface between the collective variables module and /// the simulation or analysis program (NAMD, VMD, LAMMPS...). @@ -827,6 +591,8 @@ public: /// Destructor virtual ~colvarproxy(); + virtual bool io_available() /* override */; + /// Request deallocation of the module (currently only implemented by VMD) virtual int request_deletion(); @@ -839,21 +605,30 @@ 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 the module + virtual int parse_module_config(); - /// (Re)initialize required member data after construction + /// (Re)initialize required member data (called after the module) virtual int setup(); - /// \brief Update data required by the colvars module (e.g. cache atom positions) + /// Whether the engine allows to fully initialize Colvars immediately + inline bool engine_ready() const + { + return engine_ready_; + } + + /// Enqueue new configuration text, to be parsed as soon as possible + void add_config(std::string const &cmd, std::string const &conf); + + /// Update data required by Colvars module (e.g. read atom positions) /// /// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these virtual int update_input(); - /// \brief Update data based from the results of a module update (e.g. send forces) + /// Update data based on the results of a Colvars call (e.g. send forces) virtual int update_output(); - /// Carry out operations needed before next step is run + /// Carry out operations needed before next simulation step is run int end_of_step(); /// Print a message to the main log @@ -903,26 +678,11 @@ 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: + /// Whether the engine allows to fully initialize Colvars immediately + bool engine_ready_; + /// Collected error messages std::string error_output; @@ -943,8 +703,10 @@ protected: /// Track which features have been acknowledged during the last run size_t features_hash; - /// Raise when the output stream functions are used on threads other than 0 - void smp_stream_error(); +private: + + /// Queue of config strings or files to be fed to the module + void *config_queue_; }; diff --git a/lib/colvars/colvarproxy_io.cpp b/lib/colvars/colvarproxy_io.cpp new file mode 100644 index 0000000000..225ca40bef --- /dev/null +++ b/lib/colvars/colvarproxy_io.cpp @@ -0,0 +1,310 @@ +// -*- 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. + +// Using access() to check if a file exists (until we can assume C++14/17) +#if !defined(_WIN32) || defined(__CYGWIN__) +#include +#endif +#if defined(_WIN32) +#include +#endif + +#include +#include + +#include +#include +#include +#include + +#include "colvarmodule.h" +#include "colvarproxy_io.h" + + +colvarproxy_io::colvarproxy_io() +{ + input_buffer_ = NULL; + restart_frequency_engine = 0; + input_stream_error_ = new std::istringstream(); + input_stream_error_->setstate(std::ios::badbit); + output_stream_error_ = new std::ostringstream(); + output_stream_error_->setstate(std::ios::badbit); +} + + +colvarproxy_io::~colvarproxy_io() +{ + delete input_stream_error_; + close_input_streams(); + delete output_stream_error_; + close_output_streams(); +} + + +bool colvarproxy_io::io_available() +{ + return false; +} + + +int colvarproxy_io::get_frame(long int&) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_io::set_frame(long int) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_io::backup_file(char const *filename) +{ + // Simplified version of NAMD_file_exists() + int exit_code; + do { +#if defined(_WIN32) && !defined(__CYGWIN__) + // We could use _access_s here, but it is probably too new + exit_code = _access(filename, 00); +#else + exit_code = access(filename, F_OK); +#endif + } while ((exit_code != 0) && (errno == EINTR)); + if (exit_code != 0) { + if (errno == ENOENT) { + // File does not exist + return COLVARS_OK; + } else { + return cvm::error("Unknown error while checking if file \""+ + std::string(filename)+"\" exists.\n", COLVARS_ERROR); + } + } + + // The file exists, then rename it + if (std::string(filename).rfind(std::string(".colvars.state")) != + std::string::npos) { + return rename_file(filename, (std::string(filename)+".old").c_str()); + } else { + return rename_file(filename, (std::string(filename)+".BAK").c_str()); + } +} + + +int colvarproxy_io::remove_file(char const *filename) +{ + int error_code = COLVARS_OK; +#if defined(_WIN32) && !defined(__CYGWIN__) + // Because the file may be open by other processes, rename it to filename.old + std::string const renamed_file(std::string(filename)+".old"); + // It may still be there from an interrupted run, so remove it to be safe + std::remove(renamed_file.c_str()); + int rename_exit_code = 0; + while ((rename_exit_code = std::rename(filename, + renamed_file.c_str())) != 0) { + if (errno == EINTR) continue; + error_code |= COLVARS_FILE_ERROR; + break; + } + // Ask to remove filename.old, but ignore any errors raised + std::remove(renamed_file.c_str()); +#else + if (std::remove(filename)) { + if (errno != ENOENT) { + error_code |= COLVARS_FILE_ERROR; + } + } +#endif + if (error_code != COLVARS_OK) { + return cvm::error("Error: in removing file \""+std::string(filename)+ + "\".\n.", + error_code); + } + return COLVARS_OK; +} + + +int colvarproxy_io::rename_file(char const *filename, char const *newfilename) +{ + int error_code = COLVARS_OK; +#if defined(_WIN32) && !defined(__CYGWIN__) + // On straight Windows, must remove the destination before renaming it + error_code |= remove_file(newfilename); +#endif + int rename_exit_code = 0; + while ((rename_exit_code = std::rename(filename, newfilename)) != 0) { + if (errno == EINTR) continue; + // Call log() instead of error to allow the next try + cvm::log("Error: in renaming file \""+std::string(filename)+"\" to \""+ + std::string(newfilename)+"\".\n."); + error_code |= COLVARS_FILE_ERROR; + if (errno == EXDEV) continue; + break; + } + return rename_exit_code ? error_code : COLVARS_OK; +} + + +std::istream &colvarproxy_io::input_stream(std::string const &input_name, + std::string const description, + bool error_on_fail) +{ + if (!io_available()) { + cvm::error("Error: trying to access an input file/channel " + "from the wrong thread.\n", COLVARS_BUG_ERROR); + return *input_stream_error_; + } + + if (colvarproxy_io::input_stream_exists(input_name)) { + return *(input_streams_[input_name]); + } + + // Using binary to work around differences in line termination conventions + // See https://github.com/Colvars/colvars/commit/8236879f7de4 + input_streams_[input_name] = new std::ifstream(input_name.c_str(), + std::ios::binary); + + if (input_streams_[input_name]->fail() && error_on_fail) { + cvm::error("Error: cannot open "+description+" \""+input_name+"\".\n", + COLVARS_FILE_ERROR); + } + + return *(input_streams_[input_name]); +} + + +bool colvarproxy_io::input_stream_exists(std::string const &input_name) +{ + return (input_streams_.count(input_name) > 0); +} + + +int colvarproxy_io::close_input_stream(std::string const &input_name) +{ + if (colvarproxy_io::input_stream_exists(input_name)) { + delete input_streams_[input_name]; + input_streams_.erase(input_name); + return COLVARS_OK; + } + return cvm::error("Error: input file/channel \""+input_name+ + "\" does not exist.\n", COLVARS_FILE_ERROR); +} + + +int colvarproxy_io::close_input_streams() +{ + for (std::map::iterator ii = input_streams_.begin(); + ii != input_streams_.end(); + ii++) { + delete ii->second; + } + input_streams_.clear(); + return COLVARS_OK; +} + + +std::ostream & colvarproxy_io::output_stream(std::string const &output_name, + std::string const description) +{ + if (cvm::debug()) { + cvm::log("Using colvarproxy_io::output_stream()\n"); + } + + if (!io_available()) { + cvm::error("Error: trying to access an output file/channel " + "from the wrong thread.\n", COLVARS_BUG_ERROR); + return *output_stream_error_; + } + + if (colvarproxy_io::output_stream_exists(output_name)) { + return *(output_streams_[output_name]); + } + + backup_file(output_name.c_str()); + + output_streams_[output_name] = new std::ofstream(output_name.c_str()); + if (!*(output_streams_[output_name])) { + cvm::error("Error: cannot write to "+description+" \""+output_name+"\".\n", + COLVARS_FILE_ERROR); + } + + return *(output_streams_[output_name]); +} + + +bool colvarproxy_io::output_stream_exists(std::string const &output_name) +{ + return (output_streams_.count(output_name) > 0); +} + + +int colvarproxy_io::flush_output_stream(std::string const &output_name) +{ + if (!io_available()) { + // No-op + return COLVARS_OK; + } + + if (colvarproxy_io::output_stream_exists(output_name)) { + (dynamic_cast(output_streams_[output_name]))->flush(); + return COLVARS_OK; + } + + return COLVARS_OK; +} + + +int colvarproxy_io::flush_output_streams() +{ + if (!io_available()) { + return COLVARS_OK; + } + + for (std::map::iterator osi = output_streams_.begin(); + osi != output_streams_.end(); + osi++) { + (dynamic_cast(osi->second))->flush(); + } + + return COLVARS_OK; +} + + +int colvarproxy_io::close_output_stream(std::string const &output_name) +{ + if (!io_available()) { + return cvm::error("Error: trying to access an output file/channel " + "from the wrong thread.\n", COLVARS_BUG_ERROR); + } + + if (colvarproxy_io::output_stream_exists(output_name)) { + (dynamic_cast(output_streams_[output_name]))->close(); + delete output_streams_[output_name]; + output_streams_.erase(output_name); + } + + return COLVARS_OK; +} + + +int colvarproxy_io::close_output_streams() +{ + if (! io_available()) { + return COLVARS_OK; + } + + for (std::map::iterator osi = output_streams_.begin(); + osi != output_streams_.end(); + osi++) { + (dynamic_cast(osi->second))->close(); + } + output_streams_.clear(); + + return COLVARS_OK; +} diff --git a/lib/colvars/colvarproxy_io.h b/lib/colvars/colvarproxy_io.h new file mode 100644 index 0000000000..ee217362d4 --- /dev/null +++ b/lib/colvars/colvarproxy_io.h @@ -0,0 +1,167 @@ +// -*- 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_IO_H +#define COLVARPROXY_IO_H + +#include +#include +#include + + +/// Methods for data input/output +class colvarproxy_io { + +public: + + /// Constructor + colvarproxy_io(); + + /// Destructor + virtual ~colvarproxy_io(); + + /// Ensure that we're on the main thread (derived class will do actual check) + virtual bool io_available(); + + /// \brief Save the current frame number in the argument given + // Returns error code + virtual int get_frame(long int &); + + /// \brief Set the current frame number (as well as colvarmodule::it) + // Returns error code + virtual int set_frame(long int); + + /// \brief Rename the given file, before overwriting it + virtual int backup_file(char const *filename); + + /// \brief Rename the given file, before overwriting it + inline int backup_file(std::string const &filename) + { + return backup_file(filename.c_str()); + } + + /// Remove the given file (on Windows only, rename to filename.old) + virtual int remove_file(char const *filename); + + /// Remove the given file (on Windows only, rename to filename.old) + inline int remove_file(std::string const &filename) + { + return remove_file(filename.c_str()); + } + + /// Rename the given file + virtual int rename_file(char const *filename, char const *newfilename); + + /// Rename the given file + inline int rename_file(std::string const &filename, + std::string const &newfilename) + { + return rename_file(filename.c_str(), newfilename.c_str()); + } + + /// Prefix of the input state file to be read next + inline std::string & input_prefix() + { + return input_prefix_str; + } + + /// 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_; + } + + /// Returns a reference to given input stream, creating it if needed + /// \param input_name File name (later only a handle) + /// \param description Purpose of the file + /// \param error_on_fail Raise error when failing to open (allow testing) + virtual std::istream &input_stream(std::string const &input_name, + std::string const description = "file/channel", + bool error_on_fail = true); + + /// Check if the file/channel is open (without opening it if not) + virtual bool input_stream_exists(std::string const &input_name); + + /// Closes the given input stream + virtual int close_input_stream(std::string const &input_name); + + /// Closes all input streams + virtual int close_input_streams(); + + /// Returns a reference to the named output file/channel (open it if needed) + /// \param output_name File name or identifier + /// \param description Purpose of the file + virtual std::ostream &output_stream(std::string const &output_name, + std::string const description = "file/channel"); + + /// Check if the file/channel is open (without opening it if not) + virtual bool output_stream_exists(std::string const &output_name); + + /// Flushes the given output file/channel + virtual int flush_output_stream(std::string const &output_name); + + /// Flushes all output files/channels + virtual int flush_output_streams(); + + /// Closes the given output file/channel + virtual int close_output_stream(std::string const &output_name); + + /// Close all open files/channels to prevent data loss + virtual int close_output_streams(); + +protected: + + /// 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; + + /// Container of input files/channels indexed by path name + std::map input_streams_; + + /// Object whose reference is returned when read errors occur + std::istream *input_stream_error_; + + /// Currently open output files/channels + std::map output_streams_; + + /// Object whose reference is returned when write errors occur + std::ostream *output_stream_error_; + + /// Buffer from which the input state information may be read + char const *input_buffer_; +}; + + +#endif diff --git a/lib/colvars/colvarproxy_system.cpp b/lib/colvars/colvarproxy_system.cpp new file mode 100644 index 0000000000..0a2769dcc9 --- /dev/null +++ b/lib/colvars/colvarproxy_system.cpp @@ -0,0 +1,203 @@ +// -*- 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 "colvartypes.h" +#include "colvarproxy_system.h" + + + +colvarproxy_system::colvarproxy_system() +{ + angstrom_value_ = 0.0; + kcal_mol_value_ = 0.0; + target_temperature_ = 0.0; + boltzmann_ = 0.001987191; // Default: kcal/mol/K + boundaries_type = boundaries_unsupported; + total_force_requested = false; + indirect_lambda_biasing_force = 0.0; + cached_alch_lambda_changed = false; + cached_alch_lambda = -1.0; + reset_pbc_lattice(); +} + + +colvarproxy_system::~colvarproxy_system() {} + + +int colvarproxy_system::set_unit_system(std::string const & /* units */, + bool /* check_only */) +{ + return COLVARS_NOT_IMPLEMENTED; +} + + +int colvarproxy_system::set_target_temperature(cvm::real T) +{ + target_temperature_ = T; + return COLVARS_OK; +} + + +cvm::real colvarproxy_system::dt() +{ + // TODO define, document and implement a user method to set the value of this + return 1.0; +} + + +cvm::real colvarproxy_system::rand_gaussian() +{ + // TODO define, document and implement a user method to set the value of this + return 0.0; +} + + +void colvarproxy_system::add_energy(cvm::real /* energy */) {} + + +void colvarproxy_system::request_total_force(bool yesno) +{ + if (yesno == true) + cvm::error("Error: total forces are currently not implemented.\n", + COLVARS_NOT_IMPLEMENTED); +} + + +bool colvarproxy_system::total_forces_enabled() const +{ + return false; +} + + +bool colvarproxy_system::total_forces_same_step() const +{ + return false; +} + + +inline int round_to_integer(cvm::real x) +{ + return int(cvm::floor(x+0.5)); +} + + +void colvarproxy_system::update_pbc_lattice() +{ + // Periodicity is assumed in all directions + + if (boundaries_type == boundaries_unsupported || + boundaries_type == boundaries_non_periodic) { + cvm::error("Error: setting PBC lattice with unsupported boundaries.\n", + COLVARS_BUG_ERROR); + return; + } + + { + cvm::rvector const v = cvm::rvector::outer(unit_cell_y, unit_cell_z); + reciprocal_cell_x = v/(v*unit_cell_x); + } + { + cvm::rvector const v = cvm::rvector::outer(unit_cell_z, unit_cell_x); + reciprocal_cell_y = v/(v*unit_cell_y); + } + { + cvm::rvector const v = cvm::rvector::outer(unit_cell_x, unit_cell_y); + reciprocal_cell_z = v/(v*unit_cell_z); + } +} + + +void colvarproxy_system::reset_pbc_lattice() +{ + unit_cell_x.reset(); + unit_cell_y.reset(); + unit_cell_z.reset(); + reciprocal_cell_x.reset(); + reciprocal_cell_y.reset(); + reciprocal_cell_z.reset(); +} + + +cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) + const +{ + if (boundaries_type == boundaries_unsupported) { + cvm::error("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR); + } + + cvm::rvector diff = (pos2 - pos1); + + if (boundaries_type == boundaries_non_periodic) return diff; + + cvm::real const x_shift = round_to_integer(reciprocal_cell_x*diff); + cvm::real const y_shift = round_to_integer(reciprocal_cell_y*diff); + cvm::real const z_shift = round_to_integer(reciprocal_cell_z*diff); + + diff.x -= x_shift*unit_cell_x.x + y_shift*unit_cell_y.x + + z_shift*unit_cell_z.x; + diff.y -= x_shift*unit_cell_x.y + y_shift*unit_cell_y.y + + z_shift*unit_cell_z.y; + diff.z -= x_shift*unit_cell_x.z + y_shift*unit_cell_y.z + + z_shift*unit_cell_z.z; + + return diff; +} + + +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; +} + + +int colvarproxy_system::get_alch_lambda(cvm::real * /* lambda */) +{ + return cvm::error("Error in get_alch_lambda: alchemical lambda dynamics is not supported by this build.", + COLVARS_NOT_IMPLEMENTED); +} + + +void colvarproxy_system::set_alch_lambda(cvm::real lambda) +{ + cached_alch_lambda = lambda; + cached_alch_lambda_changed = true; +} + + +int colvarproxy_system::send_alch_lambda() +{ + return cvm::error("Error in set_alch_lambda: alchemical lambda dynamics is not supported by this build.", + COLVARS_NOT_IMPLEMENTED); +} + + +int colvarproxy_system::get_dE_dlambda(cvm::real * /* force */) +{ + return cvm::error("Error in get_dE_dlambda: alchemical lambda dynamics is not supported by this build.", + COLVARS_NOT_IMPLEMENTED); +} + + +int colvarproxy_system::apply_force_dE_dlambda(cvm::real* /* force */) +{ + return cvm::error("Error in apply_force_dE_dlambda: function is not implemented by this build.", + COLVARS_NOT_IMPLEMENTED); +} + + +int colvarproxy_system::get_d2E_dlambda2(cvm::real*) +{ + return cvm::error("Error in get_d2E_dlambda2: function is not implemented by this build.", + COLVARS_NOT_IMPLEMENTED); +} diff --git a/lib/colvars/colvarproxy_system.h b/lib/colvars/colvarproxy_system.h new file mode 100644 index 0000000000..be3cd346c2 --- /dev/null +++ b/lib/colvars/colvarproxy_system.h @@ -0,0 +1,177 @@ +// -*- 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_SYSTEM_H +#define COLVARPROXY_SYSTEM_H + + +/// Methods for accessing the simulation system (PBCs, integrator, etc) +class colvarproxy_system { + +public: + + /// Constructor + colvarproxy_system(); + + /// Destructor + virtual ~colvarproxy_system(); + + /// \brief Name of the unit system used internally by Colvars (by default, that of the back-end). + /// Supported depending on the back-end: real (A, kcal/mol), metal (A, eV), electron (Bohr, Hartree), gromacs (nm, kJ/mol) + /// Note: calls to back-end PBC functions assume back-end length unit + /// We use different unit from back-end in VMD bc using PBC functions from colvarproxy base class + /// Colvars internal units are user specified, because the module exchanges info in unknown + /// composite dimensions with user input, while it only exchanges quantities of known + /// dimension with the back-end (length and forces) + std::string units; + + /// \brief Request to set the units used internally by Colvars + virtual int set_unit_system(std::string const &units, bool check_only); + + /// \brief Convert a length from Angstrom to internal + inline cvm::real angstrom_to_internal(cvm::real l) const + { + return l * angstrom_value_; + } + + /// \brief Convert a length from internal to Angstrom + inline cvm::real internal_to_angstrom(cvm::real l) const + { + return l / angstrom_value_; + } + + /// Boltzmann constant, with unit the same as energy / K + inline cvm::real boltzmann() const + { + return boltzmann_; + } + + /// Current target temperature of the simulation (K units) + inline cvm::real target_temperature() const + { + return target_temperature_; + } + + /// Set the current target temperature of the simulation (K units) + virtual int set_target_temperature(cvm::real T); + + /// \brief Time step of the simulation (fs) + virtual cvm::real dt(); + + /// \brief Pseudo-random number with Gaussian distribution + virtual cvm::real rand_gaussian(void); + + /// Pass restraint energy value for current timestep to MD engine + virtual void add_energy(cvm::real energy); + + /// \brief Get the PBC-aware distance vector between two positions + virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const; + + /// Recompute PBC reciprocal lattice (assumes XYZ periodicity) + void update_pbc_lattice(); + + /// Set the lattice vectors to zero + void reset_pbc_lattice(); + + /// \brief Tell the proxy whether total forces are needed (they may not + /// always be available) + virtual void request_total_force(bool yesno); + + /// Are total forces being used? + virtual bool total_forces_enabled() const; + + /// 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); + + /// Get value of alchemical lambda parameter from back-end (if available) + virtual int get_alch_lambda(cvm::real* lambda); + + /// Set value of alchemical lambda parameter to be sent to back-end at end of timestep + void set_alch_lambda(cvm::real lambda); + + /// Send cached value of alchemical lambda parameter to back-end (if available) + virtual int send_alch_lambda(); + + /// Get energy derivative with respect to lambda (if available) + virtual int get_dE_dlambda(cvm::real* dE_dlambda); + + /// Apply a scalar force on dE_dlambda (back-end distributes it onto atoms) + virtual int apply_force_dE_dlambda(cvm::real* force); + + /// Get energy second derivative with respect to lambda (if available) + virtual int get_d2E_dlambda2(cvm::real* d2E_dlambda2); + + /// Force to be applied onto alch. lambda, propagated from biasing forces on dE_dlambda + cvm::real indirect_lambda_biasing_force; + + /// Get weight factor from accelMD + virtual cvm::real get_accelMD_factor() const { + cvm::error("Error: accessing the reweighting factor of accelerated MD " + "is not yet implemented in the MD engine.\n", + COLVARS_NOT_IMPLEMENTED); + return 1.0; + } + virtual bool accelMD_enabled() const { + return false; + } + +protected: + + /// Next value of lambda to be sent to back-end + cvm::real cached_alch_lambda; + + /// Whether lambda has been set and needs to be updated in backend + bool cached_alch_lambda_changed; + + /// Boltzmann constant in internal Colvars units + cvm::real boltzmann_; + + /// Most up to date target temperature for the system (in K) + cvm::real target_temperature_; + + /// \brief Value of 1 Angstrom in the internal (front-end) Colvars unit for atomic coordinates + /// * defaults to 0 in the base class; derived proxy classes must set it + /// * in VMD proxy, can only be changed when no variables are defined + /// as user-defined values in composite units must be compatible with that system + cvm::real angstrom_value_; + + /// \brief Value of 1 kcal/mol in the internal Colvars unit for energy + cvm::real kcal_mol_value_; + + /// 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. + /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS) + /// minimum-image distances are computed by the host engine regardless. + enum Boundaries_type { + boundaries_non_periodic, + boundaries_pbc_ortho, + boundaries_pbc_triclinic, + boundaries_unsupported + }; + + /// Type of boundary conditions + Boundaries_type boundaries_type; + + /// Bravais lattice vectors + cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z; + + /// Reciprocal lattice vectors + cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z; +}; + +#endif diff --git a/lib/colvars/colvarproxy_tcl.cpp b/lib/colvars/colvarproxy_tcl.cpp index 700492f0e7..5bf97a0d98 100644 --- a/lib/colvars/colvarproxy_tcl.cpp +++ b/lib/colvars/colvarproxy_tcl.cpp @@ -8,6 +8,7 @@ // Colvars repository at GitHub. #include +#include #include "colvarmodule.h" #include "colvarproxy.h" diff --git a/lib/colvars/colvarproxy_volmaps.cpp b/lib/colvars/colvarproxy_volmaps.cpp index d0ae8b51ff..6c1f11e32e 100644 --- a/lib/colvars/colvarproxy_volmaps.cpp +++ b/lib/colvars/colvarproxy_volmaps.cpp @@ -33,7 +33,7 @@ int colvarproxy_volmaps::reset() clear_volmap(i); } volmaps_ids.clear(); - volmaps_ncopies.clear(); + volmaps_refcount.clear(); volmaps_values.clear(); volmaps_new_colvar_forces.clear(); return COLVARS_OK; @@ -43,7 +43,7 @@ int colvarproxy_volmaps::reset() int colvarproxy_volmaps::add_volmap_slot(int volmap_id) { volmaps_ids.push_back(volmap_id); - volmaps_ncopies.push_back(1); + volmaps_refcount.push_back(1); volmaps_values.push_back(0.0); volmaps_new_colvar_forces.push_back(0.0); return (volmaps_ids.size() - 1); @@ -95,8 +95,8 @@ void colvarproxy_volmaps::clear_volmap(int index) "previously requested.\n", COLVARS_INPUT_ERROR); } - if (volmaps_ncopies[index] > 0) { - volmaps_ncopies[index] -= 1; + if (volmaps_refcount[index] > 0) { + volmaps_refcount[index] -= 1; } } diff --git a/lib/colvars/colvarproxy_volmaps.h b/lib/colvars/colvarproxy_volmaps.h index 6e88ee83f9..f8c9ba8539 100644 --- a/lib/colvars/colvarproxy_volmaps.h +++ b/lib/colvars/colvarproxy_volmaps.h @@ -108,7 +108,7 @@ protected: /// \brief Keep track of how many times each vol map is used by a /// separate colvar object - std::vector volmaps_ncopies; + std::vector volmaps_refcount; /// \brief Current values of the vol maps std::vector volmaps_values; diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index d2a48f8af7..d500c0e5ec 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -1,3 +1,3 @@ #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2022-05-24" +#define COLVARS_VERSION "2023-05-01" #endif diff --git a/lib/colvars/colvarscript_commands.h b/lib/colvars/colvarscript_commands.h index 7f208e5da6..6dd63d82d2 100644 --- a/lib/colvars/colvarscript_commands.h +++ b/lib/colvars/colvarscript_commands.h @@ -134,11 +134,19 @@ CVSCRIPT(cv_config, 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->proxy()->add_config("config", conf); + if (script->proxy()->engine_ready()) { + // Engine allows immediate initialization + if ((script->proxy()->parse_module_config() | + script->proxy()->setup()) == COLVARS_OK) { + return COLVARS_OK; + } else { + script->add_error_msg("Error parsing configuration string"); + return COLVARSCRIPT_ERROR; + } } - script->add_error_msg("Error parsing configuration string"); - return COLVARSCRIPT_ERROR; + // Engine not ready, config will be read during proxy->setup() + return COLVARS_OK; ) CVSCRIPT(cv_configfile, @@ -146,13 +154,20 @@ CVSCRIPT(cv_configfile, 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; + script->obj_to_str(script->get_module_cmd_arg(0, objc, objv)); + script->proxy()->add_config("configfile", std::string(conf_file_name)); + if (script->proxy()->engine_ready()) { + // Engine allows immediate initialization + if ((script->proxy()->parse_module_config() | + script->proxy()->setup()) == COLVARS_OK) { + return COLVARS_OK; + } else { + script->add_error_msg("Error parsing configuration file"); + return COLVARSCRIPT_ERROR; + } } + // Engine not ready, config will be read during proxy->setup() + return COLVARS_OK; ) CVSCRIPT(cv_delete, @@ -309,6 +324,51 @@ CVSCRIPT(cv_getenergy, return COLVARS_OK; ) +CVSCRIPT(cv_getnumactiveatomgroups, + "Get the number of atom groups that currently have positive ref counts\n" + "count : integer - Total number of atom groups", + 0, 0, + "", + script->set_result_int(static_cast(script->proxy()->get_num_active_atom_groups())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getnumactiveatoms, + "Get the number of atoms that currently have positive ref counts\n" + "count : integer - Total number of atoms", + 0, 0, + "", + script->set_result_int(static_cast(script->proxy()->get_num_active_atoms())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getnumatoms, + "Get the number of requested atoms, including those not in use now\n" + "count : integer - Total number of atoms", + 0, 0, + "", + script->set_result_int(static_cast(script->proxy()->get_atom_ids()->size())); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getstepabsolute, + "Get the current step number of the simulation (including restarts)\n" + "step : int - Absolute step number", + 0, 0, + "", + script->set_result_int(cvm::step_absolute()); + return COLVARS_OK; + ) + +CVSCRIPT(cv_getsteprelative, + "Get the current step number from the start of this job\n" + "step : int - Relative step number", + 0, 0, + "", + script->set_result_int(cvm::step_relative()); + return COLVARS_OK; + ) + CVSCRIPT(cv_help, "Get the help string of the Colvars scripting interface\n" "help : string - Help string", @@ -481,6 +541,7 @@ CVSCRIPT(cv_reset, "Delete all internal configuration", 0, 0, "", + cvm::log("Resetting the Collective Variables module."); return script->module()->reset(); ) @@ -516,6 +577,20 @@ CVSCRIPT(cv_savetostring, return script->module()->write_restart_string(script->modify_str_result()); ) +CVSCRIPT(cv_targettemperature, + "Get/set target temperature, overriding what the MD engine provides\n" + "T : float - Current target temperature in K", + 0, 1, + "T : float - New target temperature in K", + char const *Targ = + script->obj_to_str(script->get_module_cmd_arg(0, objc, objv)); + if (Targ == NULL) { + return script->set_result_real(script->proxy()->target_temperature()); + } else { + return script->proxy()->set_target_temperature(strtod(Targ, NULL)); + } + ) + CVSCRIPT(cv_units, "Get or set the current Colvars unit system\n" "units : string - The current unit system", diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index 47117bf2ff..04e007cbc0 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -125,46 +125,16 @@ std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q) std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q) { std::streampos const start_pos = is.tellg(); - - std::string euler(""); - - if ( (is >> euler) && (colvarparse::to_lower_cppstr(euler) == - std::string("euler")) ) { - - // parse the Euler angles - - char sep; - cvm::real phi, theta, psi; - if ( !(is >> sep) || !(sep == '(') || - !(is >> phi) || !(is >> sep) || !(sep == ',') || - !(is >> theta) || !(is >> sep) || !(sep == ',') || - !(is >> psi) || !(is >> sep) || !(sep == ')') ) { - is.clear(); - is.seekg(start_pos, std::ios::beg); - is.setstate(std::ios::failbit); - return is; - } - - q = colvarmodule::quaternion(phi, theta, psi); - - } else { - - // parse the quaternion components - + char sep; + if ( !(is >> sep) || !(sep == '(') || + !(is >> q.q0) || !(is >> sep) || !(sep == ',') || + !(is >> q.q1) || !(is >> sep) || !(sep == ',') || + !(is >> q.q2) || !(is >> sep) || !(sep == ',') || + !(is >> q.q3) || !(is >> sep) || !(sep == ')') ) { + is.clear(); is.seekg(start_pos, std::ios::beg); - char sep; - if ( !(is >> sep) || !(sep == '(') || - !(is >> q.q0) || !(is >> sep) || !(sep == ',') || - !(is >> q.q1) || !(is >> sep) || !(sep == ',') || - !(is >> q.q2) || !(is >> sep) || !(sep == ',') || - !(is >> q.q3) || !(is >> sep) || !(sep == ')') ) { - is.clear(); - is.seekg(start_pos, std::ios::beg); - is.setstate(std::ios::failbit); - return is; - } + is.setstate(std::ios::failbit); } - return is; } @@ -318,15 +288,15 @@ void colvarmodule::rotation::build_correlation_matrix( // build the correlation matrix size_t i; for (i = 0; i < pos1.size(); i++) { - C.xx() += pos1[i].x * pos2[i].x; - C.xy() += pos1[i].x * pos2[i].y; - C.xz() += pos1[i].x * pos2[i].z; - C.yx() += pos1[i].y * pos2[i].x; - C.yy() += pos1[i].y * pos2[i].y; - C.yz() += pos1[i].y * pos2[i].z; - C.zx() += pos1[i].z * pos2[i].x; - C.zy() += pos1[i].z * pos2[i].y; - C.zz() += pos1[i].z * pos2[i].z; + C.xx += pos1[i].x * pos2[i].x; + C.xy += pos1[i].x * pos2[i].y; + C.xz += pos1[i].x * pos2[i].z; + C.yx += pos1[i].y * pos2[i].x; + C.yy += pos1[i].y * pos2[i].y; + C.yz += pos1[i].y * pos2[i].z; + C.zx += pos1[i].z * pos2[i].x; + C.zy += pos1[i].z * pos2[i].y; + C.zz += pos1[i].z * pos2[i].z; } } @@ -335,22 +305,22 @@ void colvarmodule::rotation::compute_overlap_matrix() { // build the "overlap" matrix, whose eigenvectors are stationary // points of the RMSD in the space of rotations - S[0][0] = C.xx() + C.yy() + C.zz(); - S[1][0] = C.yz() - C.zy(); + S[0][0] = C.xx + C.yy + C.zz; + S[1][0] = C.yz - C.zy; S[0][1] = S[1][0]; - S[2][0] = - C.xz() + C.zx() ; + S[2][0] = - C.xz + C.zx ; S[0][2] = S[2][0]; - S[3][0] = C.xy() - C.yx(); + S[3][0] = C.xy - C.yx; S[0][3] = S[3][0]; - S[1][1] = C.xx() - C.yy() - C.zz(); - S[2][1] = C.xy() + C.yx(); + S[1][1] = C.xx - C.yy - C.zz; + S[2][1] = C.xy + C.yx; S[1][2] = S[2][1]; - S[3][1] = C.xz() + C.zx(); + S[3][1] = C.xz + C.zx; S[1][3] = S[3][1]; - S[2][2] = - C.xx() + C.yy() - C.zz(); - S[3][2] = C.yz() + C.zy(); + S[2][2] = - C.xx + C.yy - C.zz; + S[3][2] = C.yz + C.zy; S[2][3] = S[3][2]; - S[3][3] = - C.xx() - C.yy() + C.zz(); + S[3][3] = - C.xx - C.yy + C.zz; } @@ -402,7 +372,6 @@ void colvarmodule::rotation::calc_optimal_rotation( std::vector const &pos1, std::vector const &pos2) { - C.resize(3, 3); C.reset(); build_correlation_matrix(pos1, pos2); diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index 354c216838..13b6a5bdb3 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -868,107 +868,69 @@ public: /// \brief 2-dimensional array of real numbers with three components /// along each dimension (works with colvarmodule::rvector) -class colvarmodule::rmatrix - : public colvarmodule::matrix2d { -private: +class colvarmodule::rmatrix { public: - /// Return the xx element - inline cvm::real & xx() { return (*this)[0][0]; } - /// Return the xy element - inline cvm::real & xy() { return (*this)[0][1]; } - /// Return the xz element - inline cvm::real & xz() { return (*this)[0][2]; } - /// Return the yx element - inline cvm::real & yx() { return (*this)[1][0]; } - /// Return the yy element - inline cvm::real & yy() { return (*this)[1][1]; } - /// Return the yz element - inline cvm::real & yz() { return (*this)[1][2]; } - /// Return the zx element - inline cvm::real & zx() { return (*this)[2][0]; } - /// Return the zy element - inline cvm::real & zy() { return (*this)[2][1]; } - /// Return the zz element - inline cvm::real & zz() { return (*this)[2][2]; } - - /// Return the xx element - inline cvm::real xx() const { return (*this)[0][0]; } - /// Return the xy element - inline cvm::real xy() const { return (*this)[0][1]; } - /// Return the xz element - inline cvm::real xz() const { return (*this)[0][2]; } - /// Return the yx element - inline cvm::real yx() const { return (*this)[1][0]; } - /// Return the yy element - inline cvm::real yy() const { return (*this)[1][1]; } - /// Return the yz element - inline cvm::real yz() const { return (*this)[1][2]; } - /// Return the zx element - inline cvm::real zx() const { return (*this)[2][0]; } - /// Return the zy element - inline cvm::real zy() const { return (*this)[2][1]; } - /// Return the zz element - inline cvm::real zz() const { return (*this)[2][2]; } + cvm::real xx, xy, xz, yx, yy, yz, zx, zy, zz; /// Default constructor inline rmatrix() - : cvm::matrix2d(3, 3) - {} + { + reset(); + } /// Constructor component by component inline rmatrix(cvm::real xxi, cvm::real xyi, cvm::real xzi, cvm::real yxi, cvm::real yyi, cvm::real yzi, cvm::real zxi, cvm::real zyi, cvm::real zzi) - : cvm::matrix2d(3, 3) { - this->xx() = xxi; - this->xy() = xyi; - this->xz() = xzi; - this->yx() = yxi; - this->yy() = yyi; - this->yz() = yzi; - this->zx() = zxi; - this->zy() = zyi; - this->zz() = zzi; + xx = xxi; + xy = xyi; + xz = xzi; + yx = yxi; + yy = yyi; + yz = yzi; + zx = zxi; + zy = zyi; + zz = zzi; } /// Destructor inline ~rmatrix() {} + inline void reset() + { + xx = xy = xz = yx = yy = yz = zx = zy = zz = 0.0; + } + /// Return the determinant inline cvm::real determinant() const { return - ( xx() * (yy()*zz() - zy()*yz())) - - (yx() * (xy()*zz() - zy()*xz())) - + (zx() * (xy()*yz() - yy()*xz())); + ( xx * (yy*zz - zy*yz)) + - (yx * (xy*zz - zy*xz)) + + (zx * (xy*yz - yy*xz)); } inline cvm::rmatrix transpose() const { - return cvm::rmatrix(this->xx(), this->yx(), this->zx(), - this->xy(), this->yy(), this->zy(), - this->xz(), this->yz(), this->zz()); + return cvm::rmatrix(xx, yx, zx, + xy, yy, zy, + xz, yz, zz); } - friend cvm::rvector operator * (cvm::rmatrix const &m, cvm::rvector const &r); - + inline friend cvm::rvector operator * (cvm::rmatrix const &m, + cvm::rvector const &r) + { + return cvm::rvector(m.xx*r.x + m.xy*r.y + m.xz*r.z, + m.yx*r.x + m.yy*r.y + m.yz*r.z, + m.zx*r.x + m.zy*r.y + m.zz*r.z); + } }; -inline cvm::rvector operator * (cvm::rmatrix const &m, - cvm::rvector const &r) -{ - return cvm::rvector(m.xx()*r.x + m.xy()*r.y + m.xz()*r.z, - m.yx()*r.x + m.yy()*r.y + m.yz()*r.z, - m.zx()*r.x + m.zy()*r.y + m.zz()*r.z); -} - - - /// \brief 1-dimensional vector of real numbers with four components and /// a quaternion algebra @@ -1151,11 +1113,6 @@ public: q0-=h.q0; q1-=h.q1; q2-=h.q2; q3-=h.q3; } - /// Promote a 3-vector to a quaternion - static inline cvm::quaternion promote(cvm::rvector const &v) - { - return cvm::quaternion(0.0, v.x, v.y, v.z); - } /// Return the vector component inline cvm::rvector get_vector() const { @@ -1207,7 +1164,8 @@ public: /// reference frame) inline cvm::rvector rotate(cvm::rvector const &v) const { - return ((*this) * promote(v) * ((*this).conjugate())).get_vector(); + return ( (*this) * cvm::quaternion(0.0, v.x, v.y, v.z) * + this->conjugate() ).get_vector(); } /// \brief Rotate Q2 through this quaternion (put it in the rotated @@ -1223,18 +1181,18 @@ public: { cvm::rmatrix R; - R.xx() = q0*q0 + q1*q1 - q2*q2 - q3*q3; - R.yy() = q0*q0 - q1*q1 + q2*q2 - q3*q3; - R.zz() = q0*q0 - q1*q1 - q2*q2 + q3*q3; + R.xx = q0*q0 + q1*q1 - q2*q2 - q3*q3; + R.yy = q0*q0 - q1*q1 + q2*q2 - q3*q3; + R.zz = q0*q0 - q1*q1 - q2*q2 + q3*q3; - R.xy() = 2.0 * (q1*q2 - q0*q3); - R.xz() = 2.0 * (q0*q2 + q1*q3); + R.xy = 2.0 * (q1*q2 - q0*q3); + R.xz = 2.0 * (q0*q2 + q1*q3); - R.yx() = 2.0 * (q0*q3 + q1*q2); - R.yz() = 2.0 * (q2*q3 - q0*q1); + R.yx = 2.0 * (q0*q3 + q1*q2); + R.yz = 2.0 * (q2*q3 - q0*q1); - R.zx() = 2.0 * (q1*q3 - q0*q2); - R.zy() = 2.0 * (q0*q1 + q2*q3); + R.zx = 2.0 * (q1*q3 - q0*q2); + R.zy = 2.0 * (q0*q1 + q2*q3); return R; } diff --git a/src/COLVARS/colvarproxy_lammps.cpp b/src/COLVARS/colvarproxy_lammps.cpp index eb03c14de7..06a2a23ec0 100644 --- a/src/COLVARS/colvarproxy_lammps.cpp +++ b/src/COLVARS/colvarproxy_lammps.cpp @@ -43,9 +43,10 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp, first_timestep=true; previous_step=-1; - t_target=temp; do_exit=false; + engine_ready_ = false; + // set input restart name and strip the extension, if present input_prefix_str = std::string(inp_name ? inp_name : ""); if (input_prefix_str.rfind(".colvars.state") != std::string::npos) @@ -87,7 +88,7 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp, } -void colvarproxy_lammps::init(const char *conf_file) +void colvarproxy_lammps::init() { version_int = get_version_from_string(COLVARPROXY_VERSION); @@ -102,18 +103,13 @@ void colvarproxy_lammps::init(const char *conf_file) my_angstrom = _lmp->force->angstrom; // Front-end unit is the same as back-end - angstrom_value = my_angstrom; + angstrom_value_ = my_angstrom; // my_kcal_mol = _lmp->force->qe2f / 23.060549; // force->qe2f is 1eV expressed in LAMMPS' energy unit (1 if unit is eV, 23 if kcal/mol) - my_boltzmann = _lmp->force->boltz; + boltzmann_ = _lmp->force->boltz; my_timestep = _lmp->update->dt * _lmp->force->femtosecond; - // TODO move one or more of these to setup() if needed - colvars->read_config_file(conf_file); - colvars->setup_input(); - colvars->setup_output(); - if (_lmp->update->ntimestep != 0) { cvm::log("Setting initial step number from LAMMPS: "+ cvm::to_str(_lmp->update->ntimestep)+"\n"); @@ -123,7 +119,7 @@ void colvarproxy_lammps::init(const char *conf_file) if (cvm::debug()) { cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); - cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + cvm::log("atoms_refcount = "+cvm::to_str(atoms_refcount)+"\n"); cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); cvm::log(cvm::line_marker); cvm::log("Info: done initializing the colvars proxy object.\n"); @@ -154,8 +150,12 @@ colvarproxy_lammps::~colvarproxy_lammps() // re-initialize data where needed int colvarproxy_lammps::setup() { + int error_code = colvarproxy::setup(); my_timestep = _lmp->update->dt * _lmp->force->femtosecond; - return colvars->setup(); + error_code |= colvars->update_engine_parameters(); + error_code |= colvars->setup_input(); + error_code |= colvars->setup_output(); + return error_code; } // trigger colvars computation @@ -222,17 +222,19 @@ double colvarproxy_lammps::compute() if (cvm::debug()) { cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); - cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + cvm::log("atoms_refcount = "+cvm::to_str(atoms_refcount)+"\n"); cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); cvm::log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); } - // call the collective variable module - colvars->calc(); + // Call the collective variable module + if (colvars->calc() != COLVARS_OK) { + cvm::error("Error in the collective variables module.\n", COLVARS_ERROR); + } if (cvm::debug()) { cvm::log("atoms_ids = "+cvm::to_str(atoms_ids)+"\n"); - cvm::log("atoms_ncopies = "+cvm::to_str(atoms_ncopies)+"\n"); + cvm::log("atoms_refcount = "+cvm::to_str(atoms_refcount)+"\n"); cvm::log("atoms_positions = "+cvm::to_str(atoms_positions)+"\n"); cvm::log("atoms_new_colvar_forces = "+cvm::to_str(atoms_new_colvar_forces)+"\n"); } @@ -384,7 +386,7 @@ int colvarproxy_lammps::init_atom(int atom_number) for (size_t i = 0; i < atoms_ids.size(); i++) { if (atoms_ids[i] == aid) { // this atom id was already recorded - atoms_ncopies[i] += 1; + atoms_refcount[i] += 1; return i; } } diff --git a/src/COLVARS/colvarproxy_lammps.h b/src/COLVARS/colvarproxy_lammps.h index 338f64fa86..0fc9f1ba12 100644 --- a/src/COLVARS/colvarproxy_lammps.h +++ b/src/COLVARS/colvarproxy_lammps.h @@ -33,7 +33,7 @@ class colvarproxy_lammps : public colvarproxy { LAMMPS_NS::RanPark *_random; // state of LAMMPS properties - double t_target, my_timestep, my_boltzmann, my_angstrom; + double my_timestep, my_angstrom; double bias_energy; int previous_step; @@ -50,7 +50,7 @@ class colvarproxy_lammps : public colvarproxy { colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp, const char *, const char *, const int, const double, MPI_Comm); ~colvarproxy_lammps() override; - void init(const char *); + void init(); int setup() override; // disable default and copy constructor @@ -60,7 +60,7 @@ class colvarproxy_lammps : public colvarproxy { // methods for lammps to move data or trigger actions in the proxy public: - void set_temperature(double t) { t_target = t; }; + bool total_forces_enabled() const override { return total_force_requested; }; bool total_forces_same_step() const override { return true; }; bool want_exit() const { return do_exit; }; @@ -86,10 +86,6 @@ class colvarproxy_lammps : public colvarproxy { // Request to set the units used internally by Colvars int set_unit_system(std::string const &units_in, bool check_only) override; - inline cvm::real backend_angstrom_value() override { return my_angstrom; }; - - inline cvm::real boltzmann() override { return my_boltzmann; }; - inline cvm::real temperature() override { return t_target; }; inline cvm::real dt() override { return my_timestep; @@ -101,8 +97,7 @@ class colvarproxy_lammps : public colvarproxy { void log(std::string const &message) override; void error(std::string const &message) override; - cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) const override; + cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) const override; cvm::real rand_gaussian(void) override { return _random->gaussian(); }; diff --git a/src/COLVARS/colvarproxy_lammps_version.h b/src/COLVARS/colvarproxy_lammps_version.h index b0f8965fb0..4228740554 100644 --- a/src/COLVARS/colvarproxy_lammps_version.h +++ b/src/COLVARS/colvarproxy_lammps_version.h @@ -1,3 +1,3 @@ #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2022-05-09" +#define COLVARPROXY_VERSION "2023-04-12" #endif diff --git a/src/COLVARS/fix_colvars.cpp b/src/COLVARS/fix_colvars.cpp index 0fd7bc9857..baf0209c61 100644 --- a/src/COLVARS/fix_colvars.cpp +++ b/src/COLVARS/fix_colvars.cpp @@ -443,7 +443,9 @@ void FixColvars::one_time_init() } proxy = new colvarproxy_lammps(lmp,inp_name,out_name,rng_seed,t_target,root2root); - proxy->init(conf_file); + proxy->init(); + proxy->add_config("configfile", conf_file); + proxy->parse_module_config(); num_coords = (proxy->modify_atom_positions()->size()); } @@ -677,16 +679,16 @@ void FixColvars::post_force(int /*vflag*/) error->one(FLERR,"Run aborted on request from colvars module.\n"); if (!tstat_fix) { - proxy->set_temperature(0.0); + proxy->set_target_temperature(0.0); } else { int tmp; // get thermostat target temperature from corresponding fix, // if the fix supports extraction. double *tt = (double *) tstat_fix->extract("t_target", tmp); if (tt) - proxy->set_temperature(*tt); + proxy->set_target_temperature(*tt); else - proxy->set_temperature(0.0); + proxy->set_target_temperature(0.0); } }