From 044084847b51f820498af0f8881f7de4bcaeec17 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Tue, 7 Oct 2014 20:30:25 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@12612 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/colvars/Makefile.fermi | 52 +- lib/colvars/Makefile.g++ | 43 +- lib/colvars/Makefile.mingw32-cross | 42 +- lib/colvars/Makefile.mingw64-cross | 42 +- lib/colvars/colvar.cpp | 442 +++++--- lib/colvars/colvar.h | 39 +- lib/colvars/colvaratoms.cpp | 222 ++-- lib/colvars/colvaratoms.h | 6 +- lib/colvars/colvarbias.cpp | 566 +--------- lib/colvars/colvarbias.h | 166 +-- lib/colvars/colvarbias_abf.cpp | 60 +- lib/colvars/colvarbias_abf.h | 2 + lib/colvars/colvarbias_alb.cpp | 42 +- lib/colvars/colvarbias_alb.h | 2 +- lib/colvars/colvarbias_meta.cpp | 77 +- lib/colvars/colvarbias_meta.h | 10 +- lib/colvars/colvarcomp.cpp | 53 +- lib/colvars/colvarcomp.h | 52 +- lib/colvars/colvarcomp_angles.cpp | 18 +- lib/colvars/colvarcomp_coordnums.cpp | 4 +- lib/colvars/colvarcomp_distances.cpp | 98 +- lib/colvars/colvarcomp_protein.cpp | 7 +- lib/colvars/colvarcomp_rotations.cpp | 59 +- lib/colvars/colvargrid.cpp | 2 +- lib/colvars/colvargrid.h | 406 +++---- lib/colvars/colvarmodule.cpp | 1511 ++++++++++---------------- lib/colvars/colvarmodule.h | 159 ++- lib/colvars/colvarparse.cpp | 51 +- lib/colvars/colvarparse.h | 9 +- lib/colvars/colvarproxy.h | 100 +- lib/colvars/colvartypes.h | 29 +- lib/colvars/colvarvalue.cpp | 138 ++- lib/colvars/colvarvalue.h | 370 ++++--- 33 files changed, 2284 insertions(+), 2595 deletions(-) diff --git a/lib/colvars/Makefile.fermi b/lib/colvars/Makefile.fermi index 0309bb108c..1f4b5a5b01 100644 --- a/lib/colvars/Makefile.fermi +++ b/lib/colvars/Makefile.fermi @@ -7,18 +7,20 @@ EXTRAMAKE = Makefile.lammps.empty # ------ SETTINGS ------ CXX = g++ -CXXFLAGS = -O2 -mpc64 -march=native -funroll-loops \ - -fno-rtti -fno-exceptions # -DCOLVARS_DEBUG +CXXFLAGS = -O2 -mpc64 -march=native -funroll-loops -g \ + -fno-rtti -fno-exceptions -Wall -Wno-sign-compare # -DCOLVARS_DEBUG ARCHIVE = ar ARCHFLAG = -rscv SHELL = /bin/sh # ------ DEFINITIONS ------ -SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ - colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ - colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp +SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ + colvarbias_meta.cpp colvarbias_restraint.cpp colvarcomp_angles.cpp \ + colvarcomp_coordnums.cpp colvarcomp.cpp colvarcomp_distances.cpp \ + colvarcomp_protein.cpp colvarcomp_rotations.cpp colvar.cpp colvargrid.cpp \ + colvarmodule.cpp colvarparse.cpp colvarscript.cpp colvartypes.cpp \ + colvarvalue.cpp LIB = libcolvars.a OBJ = $(SRC:.cpp=.o) @@ -26,11 +28,13 @@ EXE = #colvars_standalone # ------ MAKE PROCEDURE ------ -default: $(LIB) $(EXE) +default: $(LIB) $(EXE) Makefile.lammps + +Makefile.lammps: + @cp $(EXTRAMAKE) Makefile.lammps $(LIB): $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) - @cp $(EXTRAMAKE) Makefile.lammps colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) $(CXX) -o $@ $(CXXFLAGS) $^ @@ -49,31 +53,30 @@ colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) # ------ DEPENDENCIES ------ # -colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarproxy_standalone.h colvaratoms.h colvarparse.h colvarvalue.h -colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \ - colvartypes.h colvarproxy.h colvaratoms.h colvarparse.h colvarvalue.h \ - colvarproxy_standalone.h colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarparse.h colvarvalue.h colvaratoms.h colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \ colvarbias.h colvargrid.h +colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarbias_alb.h colvar.h colvarvalue.h colvarparse.h \ + colvarbias_restraint.h colvarbias.h colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarbias.h colvar.h colvarparse.h colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h \ colvarbias_meta.h colvarbias.h colvargrid.h -colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h +colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \ + colvarbias.h colvar.h colvarparse.h colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarcomp.h \ colvaratoms.h -colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarparse.h colvarvalue.h colvaratoms.h \ colvar.h colvarcomp.h +colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h @@ -83,14 +86,23 @@ colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \ colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h +colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ + colvarscript.h colvarbias.h colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ colvargrid.h -colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarparse.h colvarvalue.h colvar.h colvarbias.h colvarbias_meta.h \ - colvargrid.h colvarbias_abf.h +colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarparse.h colvarvalue.h colvar.h colvarbias.h \ + colvarbias_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ + colvarbias_abf.h colvarscript.h colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h +colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias.h +colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarparse.h colvarvalue.h colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h diff --git a/lib/colvars/Makefile.g++ b/lib/colvars/Makefile.g++ index 80422be3a3..6e5d4463ea 100644 --- a/lib/colvars/Makefile.g++ +++ b/lib/colvars/Makefile.g++ @@ -14,10 +14,12 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ -SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ - colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ - colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp colvarbias_alb.cpp +SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ + colvarbias_meta.cpp colvarbias_restraint.cpp colvarcomp_angles.cpp \ + colvarcomp_coordnums.cpp colvarcomp.cpp colvarcomp_distances.cpp \ + colvarcomp_protein.cpp colvarcomp_rotations.cpp colvar.cpp colvargrid.cpp \ + colvarmodule.cpp colvarparse.cpp colvarscript.cpp colvartypes.cpp \ + colvarvalue.cpp LIB = libcolvars.a OBJ = $(SRC:.cpp=.o) @@ -50,34 +52,30 @@ colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) # ------ DEPENDENCIES ------ # -colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarproxy_standalone.h colvaratoms.h colvarparse.h colvarvalue.h -colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \ - colvartypes.h colvarproxy.h colvaratoms.h colvarparse.h colvarvalue.h \ - colvarproxy_standalone.h colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarparse.h colvarvalue.h colvaratoms.h colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \ colvarbias.h colvargrid.h colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \ - colvar.h colvarvalue.h colvarparse.h colvarbias_alb.h \ - colvarbias.h + colvarproxy.h colvarbias_alb.h colvar.h colvarvalue.h colvarparse.h \ + colvarbias_restraint.h colvarbias.h colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarbias.h colvar.h colvarparse.h colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h \ colvarbias_meta.h colvarbias.h colvargrid.h -colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h +colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \ + colvarbias.h colvar.h colvarparse.h colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarcomp.h \ colvaratoms.h -colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarparse.h colvarvalue.h colvaratoms.h \ colvar.h colvarcomp.h +colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h @@ -87,14 +85,23 @@ colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \ colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h +colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ + colvarscript.h colvarbias.h colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ colvargrid.h -colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarparse.h colvarvalue.h colvar.h colvarbias.h colvarbias_meta.h \ - colvargrid.h colvarbias_abf.h colvarbias_alb.h +colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarparse.h colvarvalue.h colvar.h colvarbias.h \ + colvarbias_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ + colvarbias_abf.h colvarscript.h colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h +colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias.h +colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarparse.h colvarvalue.h colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h diff --git a/lib/colvars/Makefile.mingw32-cross b/lib/colvars/Makefile.mingw32-cross index 05bdbf32a5..39b07ca8d7 100644 --- a/lib/colvars/Makefile.mingw32-cross +++ b/lib/colvars/Makefile.mingw32-cross @@ -17,10 +17,12 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ -SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ - colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ - colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp +SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ + colvarbias_meta.cpp colvarbias_restraint.cpp colvarcomp_angles.cpp \ + colvarcomp_coordnums.cpp colvarcomp.cpp colvarcomp_distances.cpp \ + colvarcomp_protein.cpp colvarcomp_rotations.cpp colvar.cpp colvargrid.cpp \ + colvarmodule.cpp colvarparse.cpp colvarscript.cpp colvartypes.cpp \ + colvarvalue.cpp DIR = Obj_mingw32/ LIB = $(DIR)libcolvars.a @@ -58,31 +60,30 @@ $(DIR)%.o: %.cpp # ------ DEPENDENCIES ------ # -$(DIR)colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarproxy_standalone.h colvaratoms.h colvarparse.h colvarvalue.h -$(DIR)colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \ - colvartypes.h colvarproxy.h colvaratoms.h colvarparse.h colvarvalue.h \ - colvarproxy_standalone.h $(DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarparse.h colvarvalue.h colvaratoms.h $(DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \ colvarbias.h colvargrid.h +$(DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarbias_alb.h colvar.h colvarvalue.h colvarparse.h \ + colvarbias_restraint.h colvarbias.h $(DIR)colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarbias.h colvar.h colvarparse.h $(DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h \ colvarbias_meta.h colvarbias.h colvargrid.h -$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h +$(DIR)colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \ + colvarbias.h colvar.h colvarparse.h $(DIR)colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarcomp.h \ colvaratoms.h -$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h $(DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarparse.h colvarvalue.h colvaratoms.h \ colvar.h colvarcomp.h +$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h $(DIR)colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h @@ -92,14 +93,23 @@ $(DIR)colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h $(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h +$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ + colvarscript.h colvarbias.h $(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ colvargrid.h -$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarparse.h colvarvalue.h colvar.h colvarbias.h colvarbias_meta.h \ - colvargrid.h colvarbias_abf.h +$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarparse.h colvarvalue.h colvar.h colvarbias.h \ + colvarbias_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ + colvarbias_abf.h colvarscript.h $(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h +$(DIR)colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias.h +$(DIR)colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarparse.h colvarvalue.h $(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h diff --git a/lib/colvars/Makefile.mingw64-cross b/lib/colvars/Makefile.mingw64-cross index 85d63eae54..5101cf7805 100644 --- a/lib/colvars/Makefile.mingw64-cross +++ b/lib/colvars/Makefile.mingw64-cross @@ -17,10 +17,12 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ -SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ - colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ - colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ - colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp +SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias_alb.cpp colvarbias.cpp \ + colvarbias_meta.cpp colvarbias_restraint.cpp colvarcomp_angles.cpp \ + colvarcomp_coordnums.cpp colvarcomp.cpp colvarcomp_distances.cpp \ + colvarcomp_protein.cpp colvarcomp_rotations.cpp colvar.cpp colvargrid.cpp \ + colvarmodule.cpp colvarparse.cpp colvarscript.cpp colvartypes.cpp \ + colvarvalue.cpp DIR = Obj_mingw64/ LIB = $(DIR)libcolvars.a @@ -57,31 +59,30 @@ $(DIR)%.o: %.cpp # ------ DEPENDENCIES ------ # -$(DIR)colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarproxy_standalone.h colvaratoms.h colvarparse.h colvarvalue.h -$(DIR)colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \ - colvartypes.h colvarproxy.h colvaratoms.h colvarparse.h colvarvalue.h \ - colvarproxy_standalone.h $(DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarparse.h colvarvalue.h colvaratoms.h $(DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \ colvarbias.h colvargrid.h +$(DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarbias_alb.h colvar.h colvarvalue.h colvarparse.h \ + colvarbias_restraint.h colvarbias.h $(DIR)colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarbias.h colvar.h colvarparse.h $(DIR)colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h \ colvarbias_meta.h colvarbias.h colvargrid.h -$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h +$(DIR)colvarbias_restraint.o: colvarbias_restraint.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvarbias_restraint.h \ + colvarbias.h colvar.h colvarparse.h $(DIR)colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarcomp.h \ colvaratoms.h -$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h $(DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarparse.h colvarvalue.h colvaratoms.h \ colvar.h colvarcomp.h +$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h $(DIR)colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h @@ -91,14 +92,23 @@ $(DIR)colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h $(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h +$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ + colvarscript.h colvarbias.h $(DIR)colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h \ colvargrid.h -$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarparse.h colvarvalue.h colvar.h colvarbias.h colvarbias_meta.h \ - colvargrid.h colvarbias_abf.h +$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h \ + colvarproxy.h colvarparse.h colvarvalue.h colvar.h colvarbias.h \ + colvarbias_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ + colvarbias_abf.h colvarscript.h $(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h +$(DIR)colvarscript.o: colvarscript.cpp colvarscript.h colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias.h +$(DIR)colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \ + colvarparse.h colvarvalue.h $(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 8838f55442..96231ec89b 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -1,33 +1,31 @@ -// -*- c++ -*- +/// -*- c++ -*- #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" #include "colvar.h" #include "colvarcomp.h" +#include "colvarscript.h" #include -// XX TODO make the acf code handle forces as well as values and velocities - colvar::colvar (std::string const &conf) { + size_t i, j; cvm::log ("Initializing a new collective variable.\n"); get_keyval (conf, "name", this->name, (std::string ("colvar")+cvm::to_str (cvm::colvars.size()+1))); - for (std::vector::iterator cvi = cvm::colvars.begin(); - cvi < cvm::colvars.end(); - cvi++) { - if ((*cvi)->name == this->name) - cvm::fatal_error ("Error: this colvar cannot have the same name, \""+this->name+ - "\", as another colvar.\n"); + if (cvm::colvar_by_name (this->name) != NULL) { + cvm::error ("Error: this colvar cannot have the same name, \""+this->name+ + "\", as another colvar.\n", + INPUT_ERROR); } // all tasks disabled by default - for (size_t i = 0; i < task_ntot; i++) { + for (i = 0; i < task_ntot; i++) { tasks[i] = false; } @@ -57,19 +55,21 @@ colvar::colvar (std::string const &conf) cvcp->check_keywords (def_conf, def_config_key); \ cvm::decrease_depth(); \ } else { \ - cvm::fatal_error ("Error: in allocating component \"" \ - def_config_key"\".\n"); \ + cvm::error ("Error: in allocating component \"" \ + def_config_key"\".\n", \ + MEMORY_ERROR); \ } \ if ( (cvcp->period != 0.0) || (cvcp->wrap_center != 0.0) ) { \ if ( (cvcp->function_type != std::string ("distance_z")) && \ (cvcp->function_type != std::string ("dihedral")) && \ (cvcp->function_type != std::string ("spin_angle")) ) { \ - cvm::fatal_error ("Error: invalid use of period and/or " \ + cvm::error ("Error: invalid use of period and/or " \ "wrapAround in a \""+ \ std::string (def_config_key)+ \ "\" component.\n"+ \ "Period: "+cvm::to_str(cvcp->period) + \ - " wrapAround: "+cvm::to_str(cvcp->wrap_center));\ + " wrapAround: "+cvm::to_str(cvcp->wrap_center), \ + INPUT_ERROR); \ } \ } \ if ( ! cvcs.back()->name.size()) \ @@ -120,6 +120,8 @@ colvar::colvar (std::string const &conf) initialize_components ("orientation", "orientation", orientation); initialize_components ("orientation " "angle", "orientationAngle",orientation_angle); + initialize_components ("orientation " + "projection", "orientationProj",orientation_proj); initialize_components ("tilt", "tilt", tilt); initialize_components ("spin angle", "spinAngle", spin_angle); @@ -135,12 +137,55 @@ colvar::colvar (std::string const &conf) "inertiaZ", inertia_z); initialize_components ("eigenvector", "eigenvector", eigenvector); - if (!cvcs.size()) - cvm::fatal_error ("Error: no valid components were provided " - "for this collective variable.\n"); + if (!cvcs.size()) { + cvm::error ("Error: no valid components were provided " + "for this collective variable.\n", + INPUT_ERROR); + } cvm::log ("All components initialized.\n"); + // Setup colvar as scripted function of components + if (get_keyval (conf, "scriptedFunction", scripted_function, + "", colvarparse::parse_silent)) { + + enable(task_scripted); + cvm::log("This colvar is a scripted function."); + + std::string type_str; + get_keyval (conf, "scriptedFunctionType", type_str, "scalar"); + + x.type(colvarvalue::type_notset); + for (i = 0; i < colvarvalue::type_all; i++) { + if (type_str == colvarvalue::type_keyword[i]) { + x.type(colvarvalue::Type(i)); + break; + } + } + if (x.type() == colvarvalue::type_notset) { + cvm::error("Could not parse scripted colvar type."); + return; + } + x_reported.type (x.type()); + cvm::log(std::string("Expecting colvar value of type ") + + colvarvalue::type_desc[x.type()]); + + // Build ordered list of component values that will be + // passed to the script + for (i = 1; i <= cvcs.size(); i++) { + for (j = 0; j < cvcs.size(); j++) { + if (cvcs[j]->sup_np == int(i)) { + sorted_cvc_values.push_back(cvcs[j]->p_value()); + break; + } + } + } + if (sorted_cvc_values.size() != cvcs.size()) { + cvm::error("Could not find order numbers for all components" + "in componentExp values."); + return; + } + } // this is set false if any of the components has an exponent // different from 1 in the polynomial @@ -165,7 +210,7 @@ colvar::colvar (std::string const &conf) } // check the available features of each cvc - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { if ((cvcs[i])->b_debug_gradients) enable (task_gradients); @@ -198,22 +243,26 @@ colvar::colvar (std::string const &conf) if (! (cvcs[i])->b_Jacobian_derivative) b_Jacobian_force = false; - for (size_t j = i; j < cvcs.size(); j++) { - if ( (cvcs[i])->type() != (cvcs[j])->type() ) { - cvm::fatal_error ("ERROR: you are definining this collective variable " - "by using components of different types, \""+ - colvarvalue::type_desc[(cvcs[i])->type()]+ - "\" and \""+ - colvarvalue::type_desc[(cvcs[j])->type()]+ - "\". " - "You must use the same type in order to " - " sum them together.\n"); + if (!tasks[task_scripted]) { + // If the combination of components is a scripted function, + // the components may have different types + for (size_t j = i; j < cvcs.size(); j++) { + if ( (cvcs[i])->type() != (cvcs[j])->type() ) { + cvm::log ("ERROR: you are definining this collective variable " + "by using components of different types, \""+ + colvarvalue::type_desc[(cvcs[i])->type()]+ + "\" and \""+ + colvarvalue::type_desc[(cvcs[j])->type()]+ + "\". " + "You must use the same type in order to " + " sum them together.\n"); + cvm::set_error_bits(INPUT_ERROR); + } } } } - - { + if (!tasks[task_scripted]) { colvarvalue::Type const value_type = (cvcs[0])->type(); if (cvm::debug()) cvm::log ("This collective variable is a "+ @@ -225,8 +274,9 @@ colvar::colvar (std::string const &conf) } get_keyval (conf, "width", width, 1.0); - if (width <= 0.0) - cvm::fatal_error ("Error: \"width\" must be positive.\n"); + if (width <= 0.0) { + cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); + } lower_boundary.type (this->type()); lower_wall.type (this->type()); @@ -267,19 +317,21 @@ colvar::colvar (std::string const &conf) // consistency checks for boundaries and walls if (tasks[task_lower_boundary] && tasks[task_upper_boundary]) { if (lower_boundary >= upper_boundary) { - cvm::fatal_error ("Error: the 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"); + cvm::to_str (lower_boundary)+".\n", + INPUT_ERROR); } } if (tasks[task_lower_wall] && tasks[task_upper_wall]) { if (lower_wall >= upper_wall) { - cvm::fatal_error ("Error: the upper wall, "+ + cvm::error ("Error: the upper wall, "+ cvm::to_str (upper_wall)+ ", is not higher than the lower wall, "+ - cvm::to_str (lower_wall)+".\n"); + cvm::to_str (lower_wall)+".\n", + INPUT_ERROR); } if (dist2 (lower_wall, upper_wall) < 1.0E-12) { @@ -292,14 +344,16 @@ colvar::colvar (std::string const &conf) get_keyval (conf, "expandBoundaries", expand_boundaries, false); if (expand_boundaries && periodic_boundaries()) { - cvm::fatal_error ("Error: trying to expand boundaries that already " - "cover a whole period of a periodic colvar.\n"); + cvm::error ("Error: trying to expand boundaries that already " + "cover a whole period of a periodic colvar.\n", + INPUT_ERROR); } if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) { - cvm::fatal_error ("Error: inconsistent configuration " + cvm::error ("Error: inconsistent configuration " "(trying to expand boundaries with both " - "hardLowerBoundary and hardUpperBoundary enabled).\n"); - } + "hardLowerBoundary and hardUpperBoundary enabled).\n", + INPUT_ERROR); +} { bool b_extended_lagrangian; @@ -320,21 +374,24 @@ colvar::colvar (std::string const &conf) const bool found = get_keyval (conf, "extendedTemp", temp, cvm::temperature()); if (temp <= 0.0) { if (found) - cvm::fatal_error ("Error: \"extendedTemp\" must be positive.\n"); + cvm::log ("Error: \"extendedTemp\" must be positive.\n"); else - cvm::fatal_error ("Error: a positive temperature must be provided, either " - "by enabling a thermostat, or through \"extendedTemp\".\n"); + cvm::error ("Error: a positive temperature must be provided, either " + "by enabling a thermostat, or through \"extendedTemp\".\n", + INPUT_ERROR); } - get_keyval (conf, "extendedFluctuation", tolerance, width); - if (tolerance <= 0.0) - cvm::fatal_error ("Error: \"extendedFluctuation\" must be positive.\n"); + get_keyval (conf, "extendedFluctuation", tolerance); + if (tolerance <= 0.0) { + cvm::error("Error: \"extendedFluctuation\" must be positive.\n", INPUT_ERROR); + } ext_force_k = cvm::boltzmann() * temp / (tolerance * tolerance); cvm::log ("Computed extended system force constant: " + cvm::to_str(ext_force_k) + " kcal/mol/U^2"); get_keyval (conf, "extendedTimeConstant", period, 200.0); - if (period <= 0.0) - cvm::fatal_error ("Error: \"extendedTimeConstant\" must be positive.\n"); + if (period <= 0.0) { + cvm::error("Error: \"extendedTimeConstant\" must be positive.\n", INPUT_ERROR); + } ext_mass = (cvm::boltzmann() * temp * period * period) / (4.0 * PI * PI * tolerance * tolerance); cvm::log ("Computed fictitious mass: " + cvm::to_str(ext_mass) + " kcal/mol/(U/fs)^2 (U: colvar unit)"); @@ -348,8 +405,9 @@ colvar::colvar (std::string const &conf) } get_keyval (conf, "extendedLangevinDamping", ext_gamma, 1.0); - if (ext_gamma < 0.0) - cvm::fatal_error ("Error: \"extendedLangevinDamping\" may not be negative.\n"); + if (ext_gamma < 0.0) { + cvm::error("Error: \"extendedLangevinDamping\" may not be negative.\n", INPUT_ERROR); + } if (ext_gamma != 0.0) { enable (task_langevin); ext_gamma *= 1.0e-3; // convert from ps-1 to fs-1 @@ -414,7 +472,15 @@ void colvar::build_atom_list (void) temp_id_list.sort(); temp_id_list.unique(); - atom_ids = std::vector (temp_id_list.begin(), temp_id_list.end()); + + // atom_ids = std::vector (temp_id_list.begin(), temp_id_list.end()); + unsigned int id_i = 0; + std::list::iterator li; + for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) { + atom_ids[id_i] = *li; + id_i++; + } + temp_id_list.clear(); atomic_gradients.resize (atom_ids.size()); @@ -427,7 +493,7 @@ void colvar::build_atom_list (void) } -void colvar::parse_analysis (std::string const &conf) +int colvar::parse_analysis (std::string const &conf) { // if (cvm::debug()) @@ -443,8 +509,9 @@ void colvar::parse_analysis (std::string const &conf) get_keyval (conf, "runAveLength", runave_length, 1000); get_keyval (conf, "runAveStride", runave_stride, 1); - if ((cvm::restart_out_freq % runave_stride) != 0) - cvm::fatal_error ("Error: runAveStride must be commensurate with the restart frequency.\n"); + if ((cvm::restart_out_freq % runave_stride) != 0) { + cvm::error("Error: runAveStride must be commensurate with the restart frequency.\n", INPUT_ERROR); + } std::string runave_outfile; get_keyval (conf, "runAveOutputFile", runave_outfile, @@ -484,30 +551,33 @@ void colvar::parse_analysis (std::string const &conf) acf_type = acf_vel; enable (task_fdiff_velocity); if (acf_colvar_name.size()) - (cvm::colvar_p (acf_colvar_name))->enable (task_fdiff_velocity); + (cvm::colvar_by_name (acf_colvar_name))->enable (task_fdiff_velocity); } else if (acf_type_str == to_lower_cppstr (std::string ("coordinate_p2"))) { acf_type = acf_p2coor; } else { - cvm::fatal_error ("Unknown type of correlation function, \""+ + cvm::log ("Unknown type of correlation function, \""+ acf_type_str+"\".\n"); + cvm::set_error_bits(INPUT_ERROR); } get_keyval (conf, "corrFuncOffset", acf_offset, 0); get_keyval (conf, "corrFuncLength", acf_length, 1000); get_keyval (conf, "corrFuncStride", acf_stride, 1); - if ((cvm::restart_out_freq % acf_stride) != 0) - cvm::fatal_error ("Error: corrFuncStride must be commensurate with the restart frequency.\n"); + if ((cvm::restart_out_freq % acf_stride) != 0) { + cvm::error("Error: corrFuncStride must be commensurate with the restart frequency.\n", INPUT_ERROR); + } get_keyval (conf, "corrFuncNormalize", acf_normalize, true); get_keyval (conf, "corrFuncOutputFile", acf_outfile, std::string (cvm::output_prefix+"."+this->name+ ".corrfunc.dat")); } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvar::enable (colvar::task const &t) +int colvar::enable (colvar::task const &t) { switch (t) { @@ -528,13 +598,17 @@ void colvar::enable (colvar::task const &t) if ( !tasks[task_extended_lagrangian] ) { enable (task_gradients); - if (!b_Jacobian_force) - cvm::fatal_error ("Error: colvar \""+this->name+ - "\" does not have Jacobian forces implemented.\n"); - if (!b_linear) - cvm::fatal_error ("Error: colvar \""+this->name+ + if (!b_Jacobian_force) { + cvm::error ("Error: colvar \""+this->name+ + "\" does not have Jacobian forces implemented.\n", + INPUT_ERROR); + } + if (!b_linear) { + cvm::error ("Error: colvar \""+this->name+ "\" must be defined as a linear combination " - "to calculate the Jacobian force.\n"); + "to calculate the Jacobian force.\n", + INPUT_ERROR); + } if (cvm::debug()) cvm::log ("Enabling calculation of the Jacobian force " "on this colvar.\n"); @@ -545,9 +619,10 @@ void colvar::enable (colvar::task const &t) case task_system_force: if (!tasks[task_extended_lagrangian]) { if (!b_inverse_gradients) { - cvm::fatal_error ("Error: one or more of the components of " + cvm::error ("Error: one or more of the components of " "colvar \""+this->name+ - "\" does not support system force calculation.\n"); + "\" does not support system force calculation.\n", + INPUT_ERROR); } cvm::request_system_force(); } @@ -573,9 +648,11 @@ void colvar::enable (colvar::task const &t) break; case task_grid: - if (this->type() != colvarvalue::type_scalar) - cvm::fatal_error ("Cannot calculate a grid for collective variable, \""+ - this->name+"\", because its value is not a scalar number.\n"); + if (this->type() != colvarvalue::type_scalar) { + cvm::error ("Cannot calculate a grid for collective variable, \""+ + this->name+"\", because its value is not a scalar number.\n", + INPUT_ERROR); + } break; case task_extended_lagrangian: @@ -586,8 +663,9 @@ void colvar::enable (colvar::task const &t) case task_lower_boundary: case task_upper_boundary: if (this->type() != colvarvalue::type_scalar) { - cvm::fatal_error ("Error: this colvar is not a scalar value " - "and cannot produce a grid.\n"); + cvm::error ("Error: this colvar is not a scalar value " + "and cannot produce a grid.\n", + INPUT_ERROR); } break; @@ -597,6 +675,7 @@ void colvar::enable (colvar::task const &t) case task_ntot: case task_langevin: case task_output_energy: + case task_scripted: break; case task_gradients: @@ -605,19 +684,21 @@ void colvar::enable (colvar::task const &t) break; case task_collect_gradients: - if (this->type() != colvarvalue::type_scalar) - cvm::fatal_error ("Collecting atomic gradients for non-scalar collective variable \""+ - this->name+"\" is not yet implemented.\n"); + if (this->type() != colvarvalue::type_scalar) { + cvm::error ("Collecting atomic gradients for non-scalar collective variable \""+ + this->name+"\" is not yet implemented.\n", + INPUT_ERROR); + } + enable (task_gradients); if (atom_ids.size() == 0) { build_atom_list(); } break; - } - tasks[t] = true; + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -665,6 +746,7 @@ void colvar::disable (colvar::task const &t) case task_langevin: case task_output_energy: case task_collect_gradients: + case task_scripted: break; } @@ -689,6 +771,16 @@ colvar::~colvar() for (size_t i = 0; i < cvcs.size(); i++) { delete cvcs[i]; } + + // remove reference to this colvar from the CVM + for (std::vector::iterator cvi = cvm::colvars.begin(); + cvi != cvm::colvars.end(); + ++cvi) { + if ( *cvi == this) { + cvm::colvars.erase (cvi); + break; + } + } } @@ -698,14 +790,15 @@ colvar::~colvar() void colvar::calc() { + size_t i, ig; if (cvm::debug()) cvm::log ("Calculating colvar \""+this->name+"\".\n"); // prepare atom groups for calculation if (cvm::debug()) cvm::log ("Collecting data from atom groups.\n"); - for (size_t i = 0; i < cvcs.size(); i++) { - for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + for (i = 0; i < cvcs.size(); i++) { + for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); atoms.reset_atoms_data(); atoms.read_positions(); @@ -716,15 +809,15 @@ void colvar::calc() } } if (tasks[task_output_velocity]) { - for (size_t i = 0; i < cvcs.size(); i++) { - for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + for (i = 0; i < cvcs.size(); i++) { + for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvcs[i]->atom_groups[ig]->read_velocities(); } } } if (tasks[task_system_force]) { - for (size_t i = 0; i < cvcs.size(); i++) { - for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + for (i = 0; i < cvcs.size(); i++) { + for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvcs[i]->atom_groups[ig]->read_system_forces(); } } @@ -735,35 +828,42 @@ void colvar::calc() if (cvm::debug()) cvm::log ("Calculating colvar components.\n"); x.reset(); - if (x.type() == colvarvalue::type_scalar) { - // polynomial combination allowed - for (size_t i = 0; i < cvcs.size(); i++) { - cvm::increase_depth(); - (cvcs[i])->calc_value(); - cvm::decrease_depth(); - if (cvm::debug()) - cvm::log ("Colvar component no. "+cvm::to_str (i+1)+ - " within colvar \""+this->name+"\" has value "+ - cvm::to_str ((cvcs[i])->value(), - cvm::cv_width, cvm::cv_prec)+".\n"); + // First, update component values + for (i = 0; i < cvcs.size(); i++) { + cvm::increase_depth(); + (cvcs[i])->calc_value(); + cvm::decrease_depth(); + if (cvm::debug()) + cvm::log ("Colvar component no. "+cvm::to_str (i+1)+ + " within colvar \""+this->name+"\" has value "+ + cvm::to_str ((cvcs[i])->value(), + cvm::cv_width, cvm::cv_prec)+".\n"); + } + + // Then combine them appropriately + if (tasks[task_scripted]) { + // cvcs combined by user script + int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x); + if (res == COLVARS_NOT_IMPLEMENTED) { + cvm::error("Scripted colvars are not implemented."); + return; + } + if (res != COLVARS_OK) { + cvm::error("Error running scripted colvar"); + return; + } + } else if (x.type() == colvarvalue::type_scalar) { + // polynomial combination allowed + for (i = 0; i < cvcs.size(); i++) { x += (cvcs[i])->sup_coeff * - ( ((cvcs[i])->sup_np != 1) ? - std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np) : - (cvcs[i])->value().real_value ); + ( ((cvcs[i])->sup_np != 1) ? + std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np) : + (cvcs[i])->value().real_value ); } } else { // only linear combination allowed - - for (size_t i = 0; i < cvcs.size(); i++) { - cvm::increase_depth(); - (cvcs[i])->calc_value(); - cvm::decrease_depth(); - if (cvm::debug()) - cvm::log ("Colvar component no. "+cvm::to_str (i+1)+ - " within colvar \""+this->name+"\" has value "+ - cvm::to_str ((cvcs[i])->value(), - cvm::cv_width, cvm::cv_prec)+".\n"); + for (i = 0; i < cvcs.size(); i++) { x += (cvcs[i])->sup_coeff * (cvcs[i])->value(); } } @@ -777,7 +877,7 @@ void colvar::calc() if (cvm::debug()) cvm::log ("Calculating gradients of colvar \""+this->name+"\".\n"); - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { // calculate the gradients of each component cvm::increase_depth(); @@ -785,7 +885,7 @@ void colvar::calc() // if requested, propagate (via chain rule) the gradients above // to the atoms used to define the roto-translation - for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { if (cvcs[i]->atom_groups[ig]->b_fit_gradients) cvcs[i]->atom_groups[ig]->calc_fit_gradients(); } @@ -797,11 +897,18 @@ void colvar::calc() cvm::log ("Done calculating gradients of colvar \""+this->name+"\".\n"); if (tasks[task_collect_gradients]) { + + if (tasks[task_scripted]) { + cvm::error("Collecting atomic gradients is not implemented for " + "scripted colvars."); + return; + } + // Collect the atomic gradients inside colvar object - for (int a = 0; a < atomic_gradients.size(); a++) { + for (unsigned int a = 0; a < atomic_gradients.size(); a++) { atomic_gradients[a].reset(); } - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real ((cvcs[i])->sup_np) * std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1); @@ -835,14 +942,25 @@ void colvar::calc() if (tasks[task_system_force]) { + if (tasks[task_scripted]) { + // TODO see if this could reasonably be done in a generic way + // from generic inverse gradients + cvm::error("System force is not implemented for " + "scripted colvars."); + return; + } if (cvm::debug()) cvm::log ("Calculating system force of colvar \""+this->name+"\".\n"); ft.reset(); - if(!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) { + // if(!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) { + // Disabled check to allow for explicit system force calculation + // even with extended Lagrangian + + if(cvm::step_relative() > 0) { // get from the cvcs the system forces from the PREVIOUS step - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { (cvcs[i])->calc_force_invgrads(); // linear combination is assumed cvm::increase_depth(); @@ -949,15 +1067,15 @@ cvm::real colvar::update() } if (tasks[task_Jacobian_force]) { - + size_t i; cvm::increase_depth(); - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { (cvcs[i])->calc_Jacobian_derivative(); } cvm::decrease_depth(); fj.reset(); - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { // linear combination is assumed fj += 1.0 / ( cvm::real (cvcs.size()) * cvm::real ((cvcs[i])->sup_coeff) ) * (cvcs[i])->Jacobian_derivative(); @@ -1016,12 +1134,32 @@ cvm::real colvar::update() void colvar::communicate_forces() { + size_t i; if (cvm::debug()) cvm::log ("Communicating forces from colvar \""+this->name+"\".\n"); - if (x.type() == colvarvalue::type_scalar) { + if (tasks[task_scripted]) { + std::vector func_grads(cvcs.size()); + int res = cvm::proxy->run_colvar_gradient_callback(scripted_function, sorted_cvc_values, func_grads); - for (size_t i = 0; i < cvcs.size(); i++) { + if (res == COLVARS_NOT_IMPLEMENTED) { + cvm::error("Colvar gradient scripts are not implemented."); + return; + } + if (res != COLVARS_OK) { + cvm::error("Error running colvar gradient script"); + return; + } + + for (i = 0; i < cvcs.size(); i++) { + cvm::increase_depth(); + // Note: we need a dot product here + (cvcs[i])->apply_force (f * func_grads[i]); + cvm::decrease_depth(); + } + } else if (x.type() == colvarvalue::type_scalar) { + + for (i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff * cvm::real ((cvcs[i])->sup_np) * @@ -1032,7 +1170,7 @@ void colvar::communicate_forces() } else { - for (size_t i = 0; i < cvcs.size(); i++) { + for (i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff); cvm::decrease_depth(); @@ -1052,10 +1190,11 @@ void colvar::communicate_forces() bool colvar::periodic_boundaries (colvarvalue const &lb, colvarvalue const &ub) const { if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) { - cvm::fatal_error ("Error: requesting to histogram the " + cvm::log ("Error: requesting to histogram the " "collective variable \""+this->name+"\", but a " "pair of lower and upper boundaries must be " "defined.\n"); + cvm::set_error_bits(INPUT_ERROR); } if (period > 0.0) { @@ -1071,7 +1210,7 @@ bool colvar::periodic_boundaries (colvarvalue const &lb, colvarvalue const &ub) bool colvar::periodic_boundaries() const { if ( (!tasks[task_lower_boundary]) || (!tasks[task_upper_boundary]) ) { - cvm::fatal_error ("Error: requesting to histogram the " + cvm::error ("Error: requesting to histogram the " "collective variable \""+this->name+"\", but a " "pair of lower and upper boundaries must be " "defined.\n"); @@ -1132,11 +1271,11 @@ std::istream & colvar::read_restart (std::istream &is) if ( (get_keyval (conf, "name", check_name, std::string (""), colvarparse::parse_silent)) && (check_name != name) ) { - cvm::fatal_error ("Error: the state file does not match the " + cvm::error ("Error: the state file does not match the " "configuration file, at colvar \""+name+"\".\n"); } if (check_name.size() == 0) { - cvm::fatal_error ("Error: Collective variable in the " + cvm::error ("Error: Collective variable in the " "restart file without any identifier.\n"); } } @@ -1383,7 +1522,7 @@ std::ostream & colvar::write_traj (std::ostream &os) return os; } -void colvar::write_output_files() +int colvar::write_output_files() { if (cvm::b_analysis) { @@ -1391,8 +1530,9 @@ void colvar::write_output_files() cvm::log ("Writing acf to file \""+acf_outfile+"\".\n"); std::ofstream acf_os (acf_outfile.c_str()); - if (! acf_os.good()) - cvm::fatal_error ("Cannot open file \""+acf_outfile+"\".\n"); + if (! acf_os.good()) { + cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR); + } write_acf (acf_os); acf_os.close(); } @@ -1401,6 +1541,7 @@ void colvar::write_output_files() runave_os.close(); } } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -1436,7 +1577,7 @@ inline void history_incr (std::list< std::list > &history } -void colvar::calc_acf() +int colvar::calc_acf() { // using here an acf_stride-long list of vectors for either // coordinates (acf_x_history) or velocities (acf_v_history); each vector can @@ -1444,17 +1585,19 @@ void colvar::calc_acf() // representation but separated by acf_stride in the time series; // the pointer to each vector is changed at every step - if (! (acf_x_history.size() || acf_v_history.size()) ) { + if (acf_x_history.empty() && acf_v_history.empty()) { // first-step operations colvar *cfcv = (acf_colvar_name.size() ? - cvm::colvar_p (acf_colvar_name) : + cvm::colvar_by_name (acf_colvar_name) : this); - if (cfcv->type() != this->type()) - cvm::fatal_error ("Error: correlation function between \""+cfcv->name+ + if (cfcv->type() != this->type()) { + cvm::error ("Error: correlation function between \""+cfcv->name+ "\" and \""+this->name+"\" cannot be calculated, " - "because their value types are different.\n"); + "because their value types are different.\n", + INPUT_ERROR); + } acf_nframes = 0; cvm::log ("Colvar \""+this->name+"\": initializing ACF calculation.\n"); @@ -1462,11 +1605,12 @@ void colvar::calc_acf() if (acf.size() < acf_length+1) acf.resize (acf_length+1, 0.0); + size_t i; switch (acf_type) { case acf_vel: // allocate space for the velocities history - for (size_t i = 0; i < acf_stride; i++) { + for (i = 0; i < acf_stride; i++) { acf_v_history.push_back (std::list()); } acf_v_history_p = acf_v_history.begin(); @@ -1475,7 +1619,7 @@ void colvar::calc_acf() case acf_coor: case acf_p2coor: // allocate space for the coordinates history - for (size_t i = 0; i < acf_stride; i++) { + for (i = 0; i < acf_stride; i++) { acf_x_history.push_back (std::list()); } acf_x_history_p = acf_x_history.begin(); @@ -1488,7 +1632,7 @@ void colvar::calc_acf() } else { colvar *cfcv = (acf_colvar_name.size() ? - cvm::colvar_p (acf_colvar_name) : + cvm::colvar_by_name (acf_colvar_name) : this); switch (acf_type) { @@ -1532,10 +1676,11 @@ void colvar::calc_acf() // set it for the next step x_old = x; } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvar::calc_vel_acf (std::list &v_list, +int colvar::calc_vel_acf (std::list &v_list, colvarvalue const &v) { // loop over stored velocities and add to the ACF, but only the @@ -1545,10 +1690,11 @@ void colvar::calc_vel_acf (std::list &v_list, std::vector::iterator acf_i = acf.begin(); for (size_t i = 0; i < acf_offset; i++) - vs_i++; + ++vs_i; // current vel with itself - *(acf_i++) += v.norm2(); + *(acf_i) += v.norm2(); + ++acf_i; // inner products of previous velocities with current (acf_i and // vs_i are updated) @@ -1556,6 +1702,7 @@ void colvar::calc_vel_acf (std::list &v_list, acf_nframes++; } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -1568,7 +1715,7 @@ void colvar::calc_coor_acf (std::list &x_list, std::vector::iterator acf_i = acf.begin(); for (size_t i = 0; i < acf_offset; i++) - xs_i++; + ++xs_i; *(acf_i++) += x.norm2(); @@ -1589,7 +1736,7 @@ void colvar::calc_p2coor_acf (std::list &x_list, std::vector::iterator acf_i = acf.begin(); for (size_t i = 0; i < acf_offset; i++) - xs_i++; + ++xs_i; // value of P2(0) = 1 *(acf_i++) += 1.0; @@ -1617,7 +1764,7 @@ void colvar::write_acf (std::ostream &os) cvm::real const acf_norm = acf.front() / cvm::real (acf_nframes); std::vector::iterator acf_i; size_t it = acf_offset; - for (acf_i = acf.begin(); acf_i != acf.end(); acf_i++) { + for (acf_i = acf.begin(); acf_i != acf.end(); ++acf_i) { os << std::setw (cvm::it_width) << acf_stride * (it++) << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) @@ -1630,7 +1777,7 @@ void colvar::write_acf (std::ostream &os) void colvar::calc_runave() { - if (!x_history.size()) { + if (x_history.empty()) { runave.type (x.type()); runave.reset(); @@ -1653,8 +1800,9 @@ void colvar::calc_runave() if ((*x_history_p).size() >= runave_length-1) { runave = x; - for (std::list::iterator xs_i = (*x_history_p).begin(); - xs_i != (*x_history_p).end(); xs_i++) { + std::list::iterator xs_i; + for (xs_i = (*x_history_p).begin(); + xs_i != (*x_history_p).end(); ++xs_i) { runave += (*xs_i); } runave *= 1.0 / cvm::real (runave_length); @@ -1662,8 +1810,8 @@ void colvar::calc_runave() runave_variance = 0.0; runave_variance += this->dist2 (x, runave); - for (std::list::iterator xs_i = (*x_history_p).begin(); - xs_i != (*x_history_p).end(); xs_i++) { + for (xs_i = (*x_history_p).begin(); + xs_i != (*x_history_p).end(); ++xs_i) { runave_variance += this->dist2 (x, (*xs_i)); } runave_variance *= 1.0 / cvm::real (runave_length-1); diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 6f6692e904..896613d7da 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -1,4 +1,4 @@ -// -*- c++ -*- +/// -*- c++ -*- #ifndef COLVAR_H #define COLVAR_H @@ -148,6 +148,8 @@ public: task_runave, /// \brief Compute time correlation function task_corrfunc, + /// \brief Value and gradients computed by user script + task_scripted, /// \brief Number of possible tasks task_ntot }; @@ -155,6 +157,8 @@ public: /// Tasks performed by this colvar bool tasks[task_ntot]; + /// List of biases that depend on this colvar + std::vector biases; protected: @@ -283,7 +287,7 @@ public: colvar (std::string const &conf); /// Enable the specified task - void enable (colvar::task const &t); + int enable (colvar::task const &t); /// Disable the specified task void disable (colvar::task const &t); @@ -355,7 +359,7 @@ public: /// Read the analysis tasks - void parse_analysis (std::string const &conf); + int parse_analysis (std::string const &conf); /// Perform analysis tasks void analyse(); @@ -374,7 +378,7 @@ public: std::ostream & write_restart (std::ostream &os); /// Write output files (if defined, e.g. in analysis mode) - void write_output_files(); + int write_output_files(); protected: @@ -430,7 +434,7 @@ protected: acf_type_e acf_type; /// \brief Velocity ACF, scalar product between v(0) and v(t) - void calc_vel_acf (std::list &v_history, + int calc_vel_acf (std::list &v_history, colvarvalue const &v); /// \brief Coordinate ACF, scalar product between x(0) and x(t) @@ -444,7 +448,7 @@ protected: colvarvalue const &x); /// Calculate the auto-correlation function (ACF) - void calc_acf(); + int calc_acf(); /// Save the ACF to a file void write_acf (std::ostream &os); @@ -485,6 +489,7 @@ public: class h_bond; class rmsd; class orientation_angle; + class orientation_proj; class tilt; class spin_angle; class gyration; @@ -509,6 +514,14 @@ protected: /// in all cvcs (called when enabling task_collect_gradients) void build_atom_list (void); +private: + /// Name of scripted function to be used + std::string scripted_function; + + /// Current cvc values in the order requested by script + /// when using scriptedFunction + std::vector sorted_cvc_values; + public: /// \brief Sorted array of (zero-based) IDs for all atoms involved std::vector atom_ids; @@ -523,20 +536,6 @@ public: } }; - -inline colvar * cvm::colvar_p (std::string const &name) -{ - for (std::vector::iterator cvi = cvm::colvars.begin(); - cvi != cvm::colvars.end(); - cvi++) { - if ((*cvi)->name == name) { - return (*cvi); - } - } - return NULL; -} - - inline colvarvalue::Type colvar::type() const { return x.type(); diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 11b5b8e14c..f56871d614 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include "colvarmodule.h" #include "colvarparse.h" #include "colvaratoms.h" @@ -44,8 +46,8 @@ cvm::atom_group::atom_group (std::vector const &atoms) cvm::atom_group::atom_group() : b_dummy (false), b_center (false), b_rotate (false), - b_fit_gradients (false), ref_pos_group (NULL), - noforce (false) + b_user_defined_fit (false), b_fit_gradients (false), + ref_pos_group (NULL), noforce (false) { total_mass = 0.0; } @@ -63,7 +65,7 @@ cvm::atom_group::~atom_group() void cvm::atom_group::add_atom (cvm::atom const &a) { if (b_dummy) { - cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n"); + cvm::error ("Error: cannot add atoms to a dummy group.\n", INPUT_ERROR); } else { this->push_back (a); total_mass += a.mass; @@ -83,7 +85,7 @@ void cvm::atom_group::reset_mass(std::string &name, int i, int j) " atoms: total mass = "+cvm::to_str (this->total_mass)+".\n"); } -void cvm::atom_group::parse (std::string const &conf, +int cvm::atom_group::parse (std::string const &conf, char const *key) { std::string group_conf; @@ -98,9 +100,10 @@ void cvm::atom_group::parse (std::string const &conf, save_delimiters = true; if (group_conf.size() == 0) { - cvm::fatal_error ("Error: atom group \""+ - std::string (key)+"\" is set, but " - "has no definition.\n"); + cvm::error ("Error: atom group \""+std::string (key)+ + "\" is set, but has no definition.\n", + INPUT_ERROR); + return COLVARS_ERROR; } cvm::increase_depth(); @@ -135,9 +138,11 @@ void cvm::atom_group::parse (std::string const &conf, for (size_t i = 0; i < atom_indexes.size(); i++) { this->push_back (cvm::atom (atom_indexes[i])); } - } else - cvm::fatal_error ("Error: no numbers provided for \"" - "atomNumbers\".\n"); + } else { + cvm::error ("Error: no numbers provided for \"" + "atomNumbers\".\n", INPUT_ERROR); + return COLVARS_ERROR; + } atom_indexes.clear(); } @@ -147,13 +152,15 @@ void cvm::atom_group::parse (std::string const &conf, // use an index group from the index file read globally std::list::iterator names_i = cvm::index_group_names.begin(); std::list >::iterator index_groups_i = cvm::index_groups.begin(); - for ( ; names_i != cvm::index_group_names.end() ; names_i++, index_groups_i++) { + for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) { if (*names_i == index_group_name) break; } if (names_i == cvm::index_group_names.end()) { - cvm::fatal_error ("Error: could not find index group "+ - index_group_name+" among those provided by the index file.\n"); + cvm::error ("Error: could not find index group "+ + index_group_name+" among those provided by the index file.\n", + INPUT_ERROR); + return COLVARS_ERROR; } this->reserve (index_groups_i->size()); for (size_t i = 0; i < index_groups_i->size(); i++) { @@ -183,20 +190,19 @@ void cvm::atom_group::parse (std::string const &conf, } } - cvm::fatal_error ("Error: no valid definition for \"" - "atomNumbersRange\", \""+ - range_conf+"\".\n"); + cvm::error ("Error: no valid definition for \"atomNumbersRange\", \""+ + range_conf+"\".\n", INPUT_ERROR); } } std::vector psf_segids; get_keyval (group_conf, "psfSegID", psf_segids, std::vector (), mode); for (std::vector::iterator psii = psf_segids.begin(); - psii < psf_segids.end(); psii++) { + psii < psf_segids.end(); ++psii) { if ( (psii->size() == 0) || (psii->size() > 4) ) { - cvm::fatal_error ("Error: invalid segmend identifier provided, \""+ - (*psii)+"\".\n"); + cvm::error ("Error: invalid segmend identifier provided, \""+ + (*psii)+"\".\n", INPUT_ERROR); } } @@ -210,8 +216,8 @@ void cvm::atom_group::parse (std::string const &conf, range_count++; if (range_count > psf_segids.size()) { - cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than " - "values of \"psfSegID\".\n"); + cvm::error ("Error: more instances of \"atomNameResidueRange\" than " + "values of \"psfSegID\".\n", INPUT_ERROR); } std::string const &psf_segid = psf_segids.size() ? *psii : std::string (""); @@ -231,17 +237,17 @@ void cvm::atom_group::parse (std::string const &conf, } range_conf = ""; } else { - cvm::fatal_error ("Error: cannot parse definition for \"" - "atomNameResidueRange\", \""+ - range_conf+"\".\n"); + cvm::error ("Error: cannot parse definition for \"" + "atomNameResidueRange\", \""+ + range_conf+"\".\n"); } } else { - cvm::fatal_error ("Error: atomNameResidueRange with empty definition.\n"); + cvm::error ("Error: atomNameResidueRange with empty definition.\n"); } if (psf_segid.size()) - psii++; + ++psii; } } @@ -252,24 +258,25 @@ void cvm::atom_group::parse (std::string const &conf, std::string atoms_col; if (!get_keyval (group_conf, "atomsCol", atoms_col, std::string (""), mode)) { - cvm::fatal_error ("Error: parameter atomsCol is required if atomsFile is set.\n"); + cvm::error ("Error: parameter atomsCol is required if atomsFile is set.\n", + INPUT_ERROR); } double atoms_col_value; bool const atoms_col_value_defined = get_keyval (group_conf, "atomsColValue", atoms_col_value, 0.0, mode); - if (atoms_col_value_defined && (!atoms_col_value)) - cvm::fatal_error ("Error: atomsColValue, " - "if provided, must be non-zero.\n"); + if (atoms_col_value_defined && (!atoms_col_value)) { + cvm::error ("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR); + } cvm::load_atoms (atoms_file_name.c_str(), *this, atoms_col, atoms_col_value); } } for (std::vector::iterator a1 = this->begin(); - a1 != this->end(); a1++) { + a1 != this->end(); ++a1) { std::vector::iterator a2 = a1; ++a2; - for ( ; a2 != this->end(); a2++) { + for ( ; a2 != this->end(); ++a2) { if (a1->id == a2->id) { if (cvm::debug()) cvm::log ("Discarding doubly counted atom with number "+ @@ -287,14 +294,15 @@ void cvm::atom_group::parse (std::string const &conf, } else b_dummy = false; - if (b_dummy && (this->size())) - cvm::fatal_error ("Error: cannot set up group \""+ - std::string (key)+"\" as a dummy atom " - "and provide it with atom definitions.\n"); + if (b_dummy && (this->size())) { + cvm::error ("Error: cannot set up group \""+ + std::string (key)+"\" as a dummy atom " + "and provide it with atom definitions.\n", INPUT_ERROR); + } #if (! defined (COLVARS_STANDALONE)) if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) { - cvm::fatal_error ("Error: no atoms defined for atom group \""+ + cvm::error ("Error: no atoms defined for atom group \""+ std::string (key)+"\".\n"); } #endif @@ -329,17 +337,17 @@ void cvm::atom_group::parse (std::string const &conf, if (b_center || b_rotate) { if (b_dummy) - cvm::fatal_error ("Error: centerReference or rotateReference " - "cannot be defined for a dummy atom.\n"); + cvm::error ("Error: centerReference or rotateReference " + "cannot be defined for a dummy atom.\n"); if (key_lookup (group_conf, "refPositionsGroup")) { // instead of this group, define another group to compute the fit if (ref_pos_group) { - cvm::fatal_error ("Error: the atom group \""+ - std::string (key)+"\" has already a reference group " - "for the rototranslational fit, which was communicated by the " - "colvar component. You should not use refPositionsGroup " - "in this case.\n"); + cvm::error ("Error: the atom group \""+ + std::string (key)+"\" has already a reference group " + "for the rototranslational fit, which was communicated by the " + "colvar component. You should not use refPositionsGroup " + "in this case.\n"); } cvm::log ("Within atom group \""+std::string (key)+"\":\n"); ref_pos_group = new atom_group (group_conf, "refPositionsGroup"); @@ -357,8 +365,8 @@ void cvm::atom_group::parse (std::string const &conf, if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) { if (ref_pos.size()) { - cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and " - "\"refPositions\" at the same time.\n"); + cvm::error ("Error: cannot specify \"refPositionsFile\" and " + "\"refPositions\" at the same time.\n"); } std::string ref_pos_col; @@ -368,8 +376,8 @@ void cvm::atom_group::parse (std::string const &conf, // if provided, use PDB column to select coordinates bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode); if (found && !ref_pos_col_value) - cvm::fatal_error ("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); + cvm::error ("Error: refPositionsColValue, " + "if provided, must be non-zero.\n"); } else { // if not, rely on existing atom indices for the group group_for_fit->create_sorted_ids(); @@ -384,13 +392,13 @@ void cvm::atom_group::parse (std::string const &conf, if (b_rotate) { if (ref_pos.size() != group_for_fit->size()) - cvm::fatal_error ("Error: the number of reference positions provided ("+ - cvm::to_str (ref_pos.size())+ - ") does not match the number of atoms within \""+ - std::string (key)+ - "\" ("+cvm::to_str (group_for_fit->size())+ - "): to perform a rotational fit, "+ - "these numbers should be equal.\n"); + cvm::error ("Error: the number of reference positions provided ("+ + cvm::to_str (ref_pos.size())+ + ") does not match the number of atoms within \""+ + std::string (key)+ + "\" ("+cvm::to_str (group_for_fit->size())+ + "): to perform a rotational fit, "+ + "these numbers should be equal.\n", INPUT_ERROR); } // save the center of geometry of ref_pos and subtract it @@ -398,7 +406,7 @@ void cvm::atom_group::parse (std::string const &conf, } else { #if (! defined (COLVARS_STANDALONE)) - cvm::fatal_error ("Error: no reference positions provided.\n"); + cvm::error ("Error: no reference positions provided.\n"); #endif } @@ -431,29 +439,39 @@ void cvm::atom_group::parse (std::string const &conf, cvm::to_str (this->total_mass)+".\n"); cvm::decrease_depth(); + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void cvm::atom_group::create_sorted_ids (void) +int cvm::atom_group::create_sorted_ids (void) { // Only do the work if the vector is not yet populated if (sorted_ids.size()) - return; + return COLVARS_OK; std::list temp_id_list; - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + cvm::atom_iter ai; + for (ai = this->begin(); ai != this->end(); ai++) { temp_id_list.push_back (ai->id); } temp_id_list.sort(); temp_id_list.unique(); if (temp_id_list.size() != this->size()) { - cvm::fatal_error ("Error: duplicate atom IDs in atom group? (found " + - cvm::to_str(temp_id_list.size()) + - " unique atom IDs instead of" + - cvm::to_str(this->size()) + ").\n"); + cvm::error ("Error: duplicate atom IDs in atom group? (found " + + cvm::to_str(temp_id_list.size()) + + " unique atom IDs instead of" + + cvm::to_str(this->size()) + ").\n"); + return COLVARS_ERROR; } - sorted_ids = std::vector (temp_id_list.begin(), temp_id_list.end()); - return; + sorted_ids = std::vector (temp_id_list.size()); + unsigned int id_i = 0; + std::list::iterator li; + for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) { + sorted_ids[id_i] = *li; + id_i++; + } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } void cvm::atom_group::center_ref_pos() @@ -464,13 +482,12 @@ void cvm::atom_group::center_ref_pos() // This is called either by atom_group::parse or by CVCs that set // reference positions (eg. RMSD, eigenvector) ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0); - std::vector::iterator pi = ref_pos.begin(); - for ( ; pi != ref_pos.end(); pi++) { + std::vector::iterator pi; + for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) { ref_pos_cog += *pi; } ref_pos_cog /= (cvm::real) ref_pos.size(); - for (std::vector::iterator pi = ref_pos.begin(); - pi != ref_pos.end(); pi++) { + for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) { (*pi) -= ref_pos_cog; } } @@ -679,14 +696,15 @@ void cvm::atom_group::calc_fit_gradients() std::vector cvm::atom_group::positions() const { - if (b_dummy) - cvm::fatal_error ("Error: positions are not available " + if (b_dummy) { + cvm::error ("Error: positions are not available " "from a dummy atom group.\n"); + } std::vector x (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator xi = x.begin(); - for ( ; ai != this->end(); xi++, ai++) { + for ( ; ai != this->end(); ++xi, ++ai) { *xi = ai->pos; } return x; @@ -694,14 +712,15 @@ std::vector cvm::atom_group::positions() const std::vector cvm::atom_group::positions_shifted (cvm::rvector const &shift) const { - if (b_dummy) - cvm::fatal_error ("Error: positions are not available " - "from a dummy atom group.\n"); + if (b_dummy) { + cvm::error ("Error: positions are not available " + "from a dummy atom group.\n"); + } std::vector x (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator xi = x.begin(); - for ( ; ai != this->end(); xi++, ai++) { + for ( ; ai != this->end(); ++xi, ++ai) { *xi = (ai->pos + shift); } return x; @@ -709,9 +728,10 @@ std::vector cvm::atom_group::positions_shifted (cvm::rvector cons std::vector cvm::atom_group::velocities() const { - if (b_dummy) - cvm::fatal_error ("Error: velocities are not available " - "from a dummy atom group.\n"); + if (b_dummy) { + cvm::error ("Error: velocities are not available " + "from a dummy atom group.\n"); + } std::vector v (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); @@ -724,14 +744,15 @@ std::vector cvm::atom_group::velocities() const std::vector cvm::atom_group::system_forces() const { - if (b_dummy) - cvm::fatal_error ("Error: system forces are not available " - "from a dummy atom group.\n"); + if (b_dummy) { + cvm::error ("Error: system forces are not available " + "from a dummy atom group.\n"); + } std::vector f (this->size(), 0.0); cvm::atom_const_iter ai = this->begin(); std::vector::iterator fi = f.begin(); - for ( ; ai != this->end(); fi++, ai++) { + for ( ; ai != this->end(); ++fi, ++ai) { *fi = ai->system_force; } return f; @@ -739,9 +760,10 @@ std::vector cvm::atom_group::system_forces() const cvm::rvector cvm::atom_group::system_force() const { - if (b_dummy) - cvm::fatal_error ("Error: system forces are not available " - "from a dummy atom group.\n"); + if (b_dummy) { + cvm::error ("Error: system forces are not available " + "from a dummy atom group.\n"); + } cvm::rvector f (0.0); for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { @@ -756,9 +778,11 @@ void cvm::atom_group::apply_colvar_force (cvm::real const &force) if (b_dummy) return; - if (noforce) - cvm::fatal_error ("Error: sending a force to a group that has " - "\"enableForces\" set to off.\n"); + if (noforce) { + cvm::error ("Error: sending a force to a group that has " + "\"enableForces\" set to off.\n"); + return; + } if (b_rotate) { @@ -803,9 +827,11 @@ void cvm::atom_group::apply_force (cvm::rvector const &force) if (b_dummy) return; - if (noforce) - cvm::fatal_error ("Error: sending a force to a group that has " - "\"disableForces\" defined.\n"); + if (noforce) { + cvm::error ("Error: sending a force to a group that has " + "\"disableForces\" defined.\n"); + return; + } if (b_rotate) { @@ -831,12 +857,12 @@ void cvm::atom_group::apply_forces (std::vector const &forces) return; if (noforce) - cvm::fatal_error ("Error: sending a force to a group that has " - "\"disableForces\" defined.\n"); + cvm::error ("Error: sending a force to a group that has " + "\"disableForces\" defined.\n"); if (forces.size() != this->size()) { - cvm::fatal_error ("Error: trying to apply an array of forces to an atom " - "group which does not have the same length.\n"); + cvm::error ("Error: trying to apply an array of forces to an atom " + "group which does not have the same length.\n"); } if (b_rotate) { @@ -844,7 +870,7 @@ void cvm::atom_group::apply_forces (std::vector const &forces) cvm::rotation const rot_inv = rot.inverse(); cvm::atom_iter ai = this->begin(); std::vector::const_iterator fi = forces.begin(); - for ( ; ai != this->end(); fi++, ai++) { + for ( ; ai != this->end(); ++fi, ++ai) { ai->apply_force (rot_inv.rotate (*fi)); } @@ -852,7 +878,7 @@ void cvm::atom_group::apply_forces (std::vector const &forces) cvm::atom_iter ai = this->begin(); std::vector::const_iterator fi = forces.begin(); - for ( ; ai != this->end(); fi++, ai++) { + for ( ; ai != this->end(); ++fi, ++ai) { ai->apply_force (*fi); } } diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index aac8c857db..663871ffa6 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -1,4 +1,4 @@ -// -*- c++ -*- +/// -*- c++ -*- #ifndef COLVARATOMS_H #define COLVARATOMS_H @@ -137,7 +137,7 @@ public: std::vector sorted_ids; /// Allocates and populates the sorted list of atom ids - void create_sorted_ids (void); + int create_sorted_ids (void); /// \brief When updating atomic coordinates, translate them to align with the @@ -193,7 +193,7 @@ public: /// \brief Initialize the group by looking up its configuration /// string in conf and parsing it - void parse (std::string const &conf, + int parse (std::string const &conf, char const *key); /// \brief Initialize the group after a temporary vector of atoms diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 328bc6a662..d88611b9b9 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarbias.h" @@ -27,12 +29,10 @@ colvarbias::colvarbias (std::string const &conf, char const *key) get_keyval (conf, "name", name, key_str+cvm::to_str (rank)); - for (std::vector::iterator bi = cvm::biases.begin(); - bi != cvm::biases.end(); - bi++) { - if ((*bi)->name == this->name) - cvm::fatal_error ("Error: this bias cannot have the same name, \""+this->name+ - "\", as another bias.\n"); + if (cvm::bias_by_name (this->name) != NULL) { + cvm::error ("Error: this bias cannot have the same name, \""+this->name+ + "\", as another bias.\n", INPUT_ERROR); + return; } // lookup the associated colvars @@ -43,7 +43,8 @@ colvarbias::colvarbias (std::string const &conf, char const *key) } } if (!colvars.size()) { - cvm::fatal_error ("Error: no collective variables specified.\n"); + cvm::error ("Error: no collective variables specified.\n"); + return; } get_keyval (conf, "outputEnergy", b_output_energy, false); @@ -54,19 +55,45 @@ colvarbias::colvarbias() : colvarparse(), has_data (false) {} +colvarbias::~colvarbias() +{ + // Remove references to this bias from colvars + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); + ++cvi) { + for (std::vector::iterator bi = (*cvi)->biases.begin(); + bi != (*cvi)->biases.end(); + ++bi) { + if ( *bi == this) { + (*cvi)->biases.erase (bi); + break; + } + } + } + // ...and from the colvars module + for (std::vector::iterator bi = cvm::biases.begin(); + bi != cvm::biases.end(); + ++bi) { + if ( *bi == this) { + cvm::biases.erase (bi); + break; + } + } +} void colvarbias::add_colvar (std::string const &cv_name) { - if (colvar *cvp = cvm::colvar_p (cv_name)) { + if (colvar *cvp = cvm::colvar_by_name (cv_name)) { cvp->enable (colvar::task_gradients); if (cvm::debug()) cvm::log ("Applying this bias to collective variable \""+ cvp->name+"\".\n"); colvars.push_back (cvp); colvar_forces.push_back (colvarvalue (cvp->type())); + cvp->biases.push_back (this); // add back-reference to this bias to colvar } else { - cvm::fatal_error ("Error: cannot find a colvar named \""+ - cv_name+"\".\n"); + cvm::error ("Error: cannot find a colvar named \""+ + cv_name+"\".\n"); } } @@ -86,13 +113,13 @@ void colvarbias::communicate_forces() void colvarbias::change_configuration(std::string const &conf) { - cvm::fatal_error ("Error: change_configuration() not implemented.\n"); + cvm::error ("Error: change_configuration() not implemented.\n"); } cvm::real colvarbias::energy_difference(std::string const &conf) { - cvm::fatal_error ("Error: energy_difference() not implemented.\n"); + cvm::error ("Error: energy_difference() not implemented.\n"); return 0.; } @@ -115,518 +142,3 @@ std::ostream & colvarbias::write_traj (std::ostream &os) << bias_energy; return os; } - - - - -colvarbias_restraint::colvarbias_restraint (std::string const &conf, - char const *key) - : colvarbias (conf, key), - target_nstages (0), - target_nsteps (0) -{ - get_keyval (conf, "forceConstant", force_k, 1.0); - - // get the initial restraint centers - colvar_centers.resize (colvars.size()); - colvar_centers_raw.resize (colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].type (colvars[i]->type()); - colvar_centers_raw[i].type (colvars[i]->type()); - } - if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) { - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].apply_constraints(); - colvar_centers_raw[i] = colvar_centers[i]; - } - } else { - colvar_centers.clear(); - cvm::fatal_error ("Error: must define the initial centers of the restraints.\n"); - } - - if (colvar_centers.size() != colvars.size()) - cvm::fatal_error ("Error: number of centers does not match " - "that of collective variables.\n"); - - if (get_keyval (conf, "targetCenters", target_centers, colvar_centers)) { - b_chg_centers = true; - for (size_t i = 0; i < target_centers.size(); i++) { - target_centers[i].apply_constraints(); - } - } else { - b_chg_centers = false; - target_centers.clear(); - } - - if (get_keyval (conf, "targetForceConstant", target_force_k, 0.0)) { - if (b_chg_centers) - cvm::fatal_error ("Error: cannot specify both targetCenters and targetForceConstant.\n"); - - starting_force_k = force_k; - b_chg_force_k = true; - - get_keyval (conf, "targetEquilSteps", target_equil_steps, 0); - - get_keyval (conf, "lambdaSchedule", lambda_schedule, lambda_schedule); - if (lambda_schedule.size()) { - // There is one more lambda-point than stages - target_nstages = lambda_schedule.size() - 1; - } - } else { - b_chg_force_k = false; - } - - if (b_chg_centers || b_chg_force_k) { - get_keyval (conf, "targetNumSteps", target_nsteps, 0); - if (!target_nsteps) - cvm::fatal_error ("Error: targetNumSteps must be non-zero.\n"); - - if (get_keyval (conf, "targetNumStages", target_nstages, target_nstages) && - lambda_schedule.size()) { - cvm::fatal_error ("Error: targetNumStages and lambdaSchedule are incompatible.\n"); - } - - if (target_nstages) { - // This means that either numStages of lambdaSchedule has been provided - stage = 0; - restraint_FE = 0.0; - } - - if (get_keyval (conf, "targetForceExponent", force_k_exp, 1.0)) { - if (! b_chg_force_k) - cvm::log ("Warning: not changing force constant: targetForceExponent will be ignored\n"); - if (force_k_exp < 1.0) - cvm::log ("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n"); - } - } - - get_keyval (conf, "outputCenters", b_output_centers, false); - get_keyval (conf, "outputAccumulatedWork", b_output_acc_work, false); - acc_work = 0.0; - - if (cvm::debug()) - cvm::log ("Done initializing a new restraint bias.\n"); -} - -colvarbias_restraint::~colvarbias_restraint () -{ - if (cvm::n_rest_biases > 0) - cvm::n_rest_biases -= 1; -} - - -void colvarbias_restraint::change_configuration (std::string const &conf) -{ - get_keyval (conf, "forceConstant", force_k, force_k); - if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) { - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].apply_constraints(); - colvar_centers_raw[i] = colvar_centers[i]; - } - } -} - - -cvm::real colvarbias_restraint::energy_difference (std::string const &conf) -{ - std::vector alt_colvar_centers; - cvm::real alt_force_k; - cvm::real alt_bias_energy = 0.0; - - get_keyval (conf, "forceConstant", alt_force_k, force_k); - - alt_colvar_centers.resize (colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - alt_colvar_centers[i].type (colvars[i]->type()); - } - if (get_keyval (conf, "centers", alt_colvar_centers, colvar_centers)) { - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].apply_constraints(); - } - } - - for (size_t i = 0; i < colvars.size(); i++) { - alt_bias_energy += restraint_potential(restraint_convert_k(alt_force_k, colvars[i]->width), - colvars[i], - alt_colvar_centers[i]); - } - - return alt_bias_energy - bias_energy; -} - - -cvm::real colvarbias_restraint::update() -{ - bias_energy = 0.0; - - if (cvm::debug()) - cvm::log ("Updating the restraint bias \""+this->name+"\".\n"); - - // Setup first stage of staged variable force constant calculation - if (b_chg_force_k && target_nstages && cvm::step_absolute() == 0) { - cvm::real lambda; - if (lambda_schedule.size()) { - lambda = lambda_schedule[0]; - } else { - lambda = 0.0; - } - force_k = starting_force_k + (target_force_k - starting_force_k) - * std::pow (lambda, force_k_exp); - cvm::log ("Restraint " + this->name + ", stage " + - cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); - cvm::log ("Setting force constant to " + cvm::to_str (force_k)); - } - - if (b_chg_centers) { - - if (!centers_incr.size()) { - // if this is the first calculation, calculate the advancement - // at each simulation step (or stage, if applicable) - // (take current stage into account: it can be non-zero - // if we are restarting a staged calculation) - centers_incr.resize (colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - centers_incr[i].type (colvars[i]->type()); - centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) / - cvm::real ( target_nstages ? (target_nstages - stage) : - (target_nsteps - cvm::step_absolute())); - } - if (cvm::debug()) - cvm::log ("Center increment for the restraint bias \""+ - this->name+"\": "+cvm::to_str (centers_incr)+" at stage "+cvm::to_str (stage)+ ".\n"); - - } - - if (target_nstages) { - if ((cvm::step_relative() > 0) - && (cvm::step_absolute() % target_nsteps) == 0 - && stage < target_nstages) { - - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers_raw[i] += centers_incr[i]; - colvar_centers[i] = colvar_centers_raw[i]; - colvars[i]->wrap(colvar_centers[i]); - colvar_centers[i].apply_constraints(); - } - stage++; - cvm::log ("Moving restraint stage " + cvm::to_str(stage) + - " : setting centers to " + cvm::to_str (colvar_centers) + - " at step " + cvm::to_str (cvm::step_absolute())); - } - } else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) { - // move the restraint centers in the direction of the targets - // (slow growth) - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers_raw[i] += centers_incr[i]; - colvar_centers[i] = colvar_centers_raw[i]; - colvars[i]->wrap(colvar_centers[i]); - colvar_centers[i].apply_constraints(); - } - } - - if (cvm::debug()) - cvm::log ("Current centers for the restraint bias \""+ - this->name+"\": "+cvm::to_str (colvar_centers)+".\n"); - } - - if (b_chg_force_k) { - // Coupling parameter, between 0 and 1 - cvm::real lambda; - - if (target_nstages) { - // TI calculation: estimate free energy derivative - // need current lambda - if (lambda_schedule.size()) { - lambda = lambda_schedule[stage]; - } else { - lambda = cvm::real(stage) / cvm::real(target_nstages); - } - - if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) { - // Start averaging after equilibration period, if requested - - // Square distance normalized by square colvar width - cvm::real dist_sq = 0.0; - for (size_t i = 0; i < colvars.size(); i++) { - dist_sq += colvars[i]->dist2 (colvars[i]->value(), colvar_centers[i]) - / (colvars[i]->width * colvars[i]->width); - } - - restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0) - * (target_force_k - starting_force_k) * dist_sq; - } - - // Finish current stage... - if (cvm::step_absolute() % target_nsteps == 0 && - cvm::step_absolute() > 0) { - - cvm::log ("Lambda= " + cvm::to_str (lambda) + " dA/dLambda= " - + cvm::to_str (restraint_FE / cvm::real(target_nsteps - target_equil_steps))); - - // ...and move on to the next one - if (stage < target_nstages) { - - restraint_FE = 0.0; - stage++; - if (lambda_schedule.size()) { - lambda = lambda_schedule[stage]; - } else { - lambda = cvm::real(stage) / cvm::real(target_nstages); - } - force_k = starting_force_k + (target_force_k - starting_force_k) - * std::pow (lambda, force_k_exp); - cvm::log ("Restraint " + this->name + ", stage " + - cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); - cvm::log ("Setting force constant to " + cvm::to_str (force_k)); - } - } - } else if (cvm::step_absolute() <= target_nsteps) { - // update force constant (slow growth) - lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps); - force_k = starting_force_k + (target_force_k - starting_force_k) - * std::pow (lambda, force_k_exp); - } - } - - if (cvm::debug()) - cvm::log ("Done updating the restraint bias \""+this->name+"\".\n"); - - // Force and energy calculation - for (size_t i = 0; i < colvars.size(); i++) { - colvar_forces[i] = -restraint_force(restraint_convert_k(force_k, colvars[i]->width), - colvars[i], - colvar_centers[i]); - bias_energy += restraint_potential(restraint_convert_k(force_k, colvars[i]->width), - colvars[i], - colvar_centers[i]); - } - - if (b_output_acc_work) { - if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) { - for (size_t i = 0; i < colvars.size(); i++) { - // project forces on the calculated increments at this step - acc_work += colvar_forces[i] * centers_incr[i]; - } - } - } - - if (cvm::debug()) - cvm::log ("Current forces for the restraint bias \""+ - this->name+"\": "+cvm::to_str (colvar_forces)+".\n"); - - return bias_energy; -} - - -std::istream & colvarbias_restraint::read_restart (std::istream &is) -{ - size_t const start_pos = is.tellg(); - - cvm::log ("Restarting restraint bias \""+ - this->name+"\".\n"); - - std::string key, brace, conf; - if ( !(is >> key) || !(key == "restraint" || key == "harmonic") || - !(is >> brace) || !(brace == "{") || - !(is >> colvarparse::read_block ("configuration", conf)) ) { - - cvm::log ("Error: in reading restart configuration for restraint bias \""+ - this->name+"\" at position "+ - cvm::to_str (is.tellg())+" in stream.\n"); - is.clear(); - is.seekg (start_pos, std::ios::beg); - is.setstate (std::ios::failbit); - return is; - } - -// int id = -1; - std::string name = ""; -// if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) && -// (id != this->id) ) || - if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) && - (name != this->name) ) - cvm::fatal_error ("Error: in the restart file, the " - "\"restraint\" block has a wrong name\n"); -// if ( (id == -1) && (name == "") ) { - if (name.size() == 0) { - cvm::fatal_error ("Error: \"restraint\" block in the restart file " - "has no identifiers.\n"); - } - - if (b_chg_centers) { -// cvm::log ("Reading the updated restraint centers from the restart.\n"); - if (!get_keyval (conf, "centers", colvar_centers)) - cvm::fatal_error ("Error: restraint centers are missing from the restart.\n"); - if (!get_keyval (conf, "centers_raw", colvar_centers_raw)) - cvm::fatal_error ("Error: \"raw\" restraint centers are missing from the restart.\n"); - } - - if (b_chg_force_k) { -// cvm::log ("Reading the updated force constant from the restart.\n"); - if (!get_keyval (conf, "forceConstant", force_k)) - cvm::fatal_error ("Error: force constant is missing from the restart.\n"); - } - - if (target_nstages) { -// cvm::log ("Reading current stage from the restart.\n"); - if (!get_keyval (conf, "stage", stage)) - cvm::fatal_error ("Error: current stage is missing from the restart.\n"); - } - - if (b_output_acc_work) { - if (!get_keyval (conf, "accumulatedWork", acc_work)) - cvm::fatal_error ("Error: accumulatedWork is missing from the restart.\n"); - } - - is >> brace; - if (brace != "}") { - cvm::fatal_error ("Error: corrupt restart information for restraint bias \""+ - this->name+"\": no matching brace at position "+ - cvm::to_str (is.tellg())+" in the restart file.\n"); - is.setstate (std::ios::failbit); - } - return is; -} - - -std::ostream & colvarbias_restraint::write_restart (std::ostream &os) -{ - os << "restraint {\n" - << " configuration {\n" - // << " id " << this->id << "\n" - << " name " << this->name << "\n"; - - if (b_chg_centers) { - os << " centers "; - for (size_t i = 0; i < colvars.size(); i++) { - os << " " << colvar_centers[i]; - } - os << "\n"; - os << " centers_raw "; - for (size_t i = 0; i < colvars.size(); i++) { - os << " " << colvar_centers_raw[i]; - } - os << "\n"; - } - - if (b_chg_force_k) { - os << " forceConstant " - << std::setprecision (cvm::en_prec) - << std::setw (cvm::en_width) << force_k << "\n"; - } - - if (target_nstages) { - os << " stage " << std::setw (cvm::it_width) - << stage << "\n"; - } - - if (b_output_acc_work) { - os << " accumulatedWork " << acc_work << "\n"; - } - - os << " }\n" - << "}\n\n"; - - return os; -} - - -std::ostream & colvarbias_restraint::write_traj_label (std::ostream &os) -{ - os << " "; - - if (b_output_energy) - os << " E_" - << cvm::wrap_string (this->name, cvm::en_width-2); - - if (b_output_centers) - for (size_t i = 0; i < colvars.size(); i++) { - size_t const this_cv_width = (colvars[i]->value()).output_width (cvm::cv_width); - os << " x0_" - << cvm::wrap_string (colvars[i]->name, this_cv_width-3); - } - - if (b_output_acc_work) - os << " W_" - << cvm::wrap_string (this->name, cvm::en_width-2); - - return os; -} - - -std::ostream & colvarbias_restraint::write_traj (std::ostream &os) -{ - os << " "; - - if (b_output_energy) - os << " " - << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) - << bias_energy; - - if (b_output_centers) - for (size_t i = 0; i < colvars.size(); i++) { - os << " " - << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) - << colvar_centers[i]; - } - - if (b_output_acc_work) - os << " " - << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) - << acc_work; - - return os; -} - -colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(std::string const &conf, char const *key) : - colvarbias_restraint(conf, key) { - for (size_t i = 0; i < colvars.size(); i++) { - if (colvars[i]->width != 1.0) - cvm::log ("The force constant for colvar \""+colvars[i]->name+ - "\" will be rescaled to "+ - cvm::to_str (restraint_convert_k(force_k, colvars[i]->width))+ - " according to the specified width.\n"); - } -} - -cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k, colvar* x, const colvarvalue &xcenter) const -{ - return 0.5 * k * x->dist2(x->value(), xcenter); -} - -colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k, colvar* x, const colvarvalue &xcenter) const -{ - return 0.5 * k * x->dist2_lgrad(x->value(), xcenter); -} - -cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k, cvm::real dist_measure) const -{ - return k / (dist_measure * dist_measure); -} - - -colvarbias_restraint_linear::colvarbias_restraint_linear(std::string const &conf, char const *key) : - colvarbias_restraint(conf, key) { - for (size_t i = 0; i < colvars.size(); i++) { - if (colvars[i]->width != 1.0) - cvm::log ("The force constant for colvar \""+colvars[i]->name+ - "\" will be rescaled to "+ - cvm::to_str (restraint_convert_k(force_k, colvars[i]->width))+ - " according to the specified width.\n"); - } -} - -cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k, colvar* x, const colvarvalue &xcenter) const -{ - return k * (x->value() - xcenter); -} - -colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k, colvar* x, const colvarvalue &xcenter) const -{ - return k; -} - -cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k, cvm::real dist_measure) const -{ - return k / dist_measure; -} diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index c4e6493151..c3624784d1 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #ifndef COLVARBIAS_H #define COLVARBIAS_H @@ -9,9 +11,6 @@ class colvarbias : public colvarparse { public: - /// Numeric id of this bias - int id; - /// Name of this bias std::string name; @@ -22,6 +21,8 @@ public: /// Return bias energy virtual cvm::real update() = 0; + // TODO: move update_bias here (share with metadynamics) + /// Load new configuration - force constant and/or centers only virtual void change_configuration(std::string const &conf); @@ -44,7 +45,7 @@ public: colvarbias(); /// Destructor - virtual inline ~colvarbias() {} + virtual ~colvarbias(); /// Read the bias configuration from a restart file virtual std::istream & read_restart (std::istream &is) = 0; @@ -58,7 +59,9 @@ public: /// Output quantities such as the bias energy to the trajectory file virtual std::ostream & write_traj (std::ostream &os); - + inline cvm::real get_energy () { + return bias_energy; + } protected: /// \brief Pointers to collective variables to which the bias is @@ -81,157 +84,4 @@ protected: }; - -/// \brief Bias restraint, optionally moving towards a target -/// (implementation of \link colvarbias \endlink) -class colvarbias_restraint : public colvarbias { - -public: - - /// Retrieve colvar values and calculate their biasing forces - virtual cvm::real update(); - - /// Load new configuration - force constant and/or centers only - virtual void change_configuration(std::string const &conf); - - /// Calculate change in energy from using alternate configuration - virtual cvm::real energy_difference(std::string const &conf); - - /// Read the bias configuration from a restart file - virtual std::istream & read_restart (std::istream &is); - - /// Write the bias configuration to a restart file - virtual std::ostream & write_restart (std::ostream &os); - - /// Write a label to the trajectory file (comment line) - virtual std::ostream & write_traj_label (std::ostream &os); - - /// Output quantities such as the bias energy to the trajectory file - virtual std::ostream & write_traj (std::ostream &os); - - /// \brief Constructor - colvarbias_restraint (std::string const &conf, char const *key); - - /// Destructor - virtual ~colvarbias_restraint(); - - -protected: - - /// \brief Potential function - virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0; - - /// \brief Force function - virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0; - - ///\brief Unit scaling - virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0; - - /// \brief Restraint centers - std::vector colvar_centers; - - /// \brief Restraint centers without wrapping or constraints applied - std::vector colvar_centers_raw; - - /// \brief Moving target? - bool b_chg_centers; - - /// \brief New restraint centers - std::vector target_centers; - - /// \brief Amplitude of the restraint centers' increment at each step - /// (or stage) towards the new values (calculated from target_nsteps) - std::vector centers_incr; - - /// Whether to write the current restraint centers to the trajectory file - bool b_output_centers; - - /// Whether to write the current accumulated work to the trajectory file - bool b_output_acc_work; - - /// \brief Accumulated work - cvm::real acc_work; - - /// \brief Restraint force constant - cvm::real force_k; - - /// \brief Changing force constant? - bool b_chg_force_k; - - /// \brief Restraint force constant (target value) - cvm::real target_force_k; - - /// \brief Restraint force constant (starting value) - cvm::real starting_force_k; - - /// \brief Lambda-schedule for custom varying force constant - std::vector lambda_schedule; - - /// \brief Exponent for varying the force constant - cvm::real force_k_exp; - - /// \brief Intermediate quantity to compute the restraint free energy - /// (in TI, would be the accumulating FE derivative) - cvm::real restraint_FE; - - - /// \brief Equilibration steps for restraint FE calculation through TI - cvm::real target_equil_steps; - - /// \brief Number of stages over which to perform the change - /// If zero, perform a continuous change - int target_nstages; - - /// \brief Number of current stage of the perturbation - int stage; - - /// \brief Number of steps required to reach the target force constant - /// or restraint centers - size_t target_nsteps; -}; - -/// \brief Harmonic bias restraint -/// (implementation of \link colvarbias_restraint \endlink) -class colvarbias_restraint_harmonic : public colvarbias_restraint { - -public: - colvarbias_restraint_harmonic(std::string const &conf, char const *key); - -protected: /// \brief Potential function - virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const; - - /// \brief Force function - virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const; - - ///\brief Unit scaling - virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const; - -}; - -/// \brief Linear bias restraint -/// (implementation of \link colvarbias_restraint \endlink) -class colvarbias_restraint_linear : public colvarbias_restraint { - -public: - colvarbias_restraint_linear(std::string const &conf, char const *key); - -protected: /// \brief Potential function - virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const; - - /// \brief Force function - virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const; - - ///\brief Unit scaling - virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const; - -}; - - #endif - - - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index dcbfba6acd..cf5825650c 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + /******************************************************************************** * Implementation of the ABF and histogram biases * ********************************************************************************/ @@ -13,6 +15,7 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key) gradients (NULL), samples (NULL) { + // TODO relax this in case of VMD plugin if (cvm::temperature() == 0.0) cvm::log ("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n"); @@ -44,13 +47,13 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key) // ************* checking the associated colvars ******************* if (colvars.size() == 0) { - cvm::fatal_error ("Error: no collective variables specified for the ABF bias.\n"); + cvm::error ("Error: no collective variables specified for the ABF bias.\n"); } for (size_t i = 0; i < colvars.size(); i++) { if (colvars[i]->type() != colvarvalue::type_scalar) { - cvm::fatal_error ("Error: ABF bias can only use scalar-type variables.\n"); + cvm::error ("Error: ABF bias can only use scalar-type variables.\n"); } colvars[i]->enable (colvar::task_gradients); @@ -84,11 +87,11 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key) if (get_keyval (conf, "maxForce", max_force)) { if (max_force.size() != colvars.size()) { - cvm::fatal_error ("Error: Number of parameters to maxForce does not match number of colvars."); + cvm::error ("Error: Number of parameters to maxForce does not match number of colvars."); } for (size_t i=0; iwrite_multicol (samples_os); samples_os.close (); if (!append) cvm::backup_file (gradients_out_name.c_str()); gradients_os.open (gradients_out_name.c_str(), mode); - if (!gradients_os.good()) cvm::fatal_error ("Error opening ABF gradient file " + gradients_out_name + " for writing"); + if (!gradients_os.good()) cvm::error ("Error opening ABF gradient file " + gradients_out_name + " for writing"); gradients->write_multicol (gradients_os); gradients_os.close (); @@ -257,7 +260,7 @@ void colvarbias_abf::write_gradients_samples (const std::string &prefix, bool ap std::ofstream pmf_os; // Do numerical integration and output a PMF pmf_os.open (pmf_out_name.c_str(), mode); - if (!pmf_os.good()) cvm::fatal_error ("Error opening pmf file " + pmf_out_name + " for writing"); + if (!pmf_os.good()) cvm::error ("Error opening pmf file " + pmf_out_name + " for writing"); gradients->write_1D_integral (pmf_os); pmf_os << std::endl; pmf_os.close (); @@ -278,13 +281,13 @@ void colvarbias_abf::read_gradients_samples () cvm::log ("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name); is.open (samples_in_name.c_str()); - if (!is.good()) cvm::fatal_error ("Error opening ABF samples file " + samples_in_name + " for reading"); + if (!is.good()) 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.good()) cvm::fatal_error ("Error opening ABF gradient file " + gradients_in_name + " for reading"); + if (!is.good()) cvm::error ("Error opening ABF gradient file " + gradients_in_name + " for reading"); gradients->read_multicol (is, true); is.close (); } @@ -319,7 +322,7 @@ std::ostream & colvarbias_abf::write_restart (std::ostream& os) std::istream & colvarbias_abf::read_restart (std::istream& is) { if ( input_prefix.size() > 0 ) { - cvm::fatal_error ("ERROR: cannot provide both inputPrefix and restart information (colvarsInput)"); + cvm::error ("ERROR: cannot provide both inputPrefix and restart information (colvarsInput)"); } size_t const start_pos = is.tellg(); @@ -343,10 +346,10 @@ std::istream & colvarbias_abf::read_restart (std::istream& is) std::string name = ""; if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) && (name != this->name) ) - cvm::fatal_error ("Error: in the restart file, the " + cvm::error ("Error: in the restart file, the " "\"abf\" block has wrong name (" + name + ")\n"); if ( name == "" ) { - cvm::fatal_error ("Error: \"abf\" block in the restart file has no name.\n"); + cvm::error ("Error: \"abf\" block in the restart file has no name.\n"); } if ( !(is >> key) || !(key == "samples")) { @@ -359,7 +362,10 @@ std::istream & colvarbias_abf::read_restart (std::istream& is) return is; } if (! samples->read_raw (is)) { - samples->read_raw_error(); + is.clear(); + is.seekg (start_pos, std::ios::beg); + is.setstate (std::ios::failbit); + return is; } if ( !(is >> key) || !(key == "gradient")) { @@ -372,12 +378,15 @@ std::istream & colvarbias_abf::read_restart (std::istream& is) return is; } if (! gradients->read_raw (is)) { - gradients->read_raw_error(); + is.clear(); + is.seekg (start_pos, std::ios::beg); + is.setstate (std::ios::failbit); + return is; } is >> brace; if (brace != "}") { - cvm::fatal_error ("Error: corrupt restart information for ABF bias \""+ + cvm::error ("Error: corrupt restart information for ABF bias \""+ this->name+"\": no matching brace at position "+ cvm::to_str (is.tellg())+" in the restart file.\n"); is.setstate (std::ios::failbit); @@ -397,7 +406,7 @@ colvarbias_histogram::colvarbias_histogram (std::string const &conf, char const get_keyval (conf, "outputfreq", output_freq, cvm::restart_out_freq); if ( output_freq == 0 ) { - cvm::fatal_error ("User required histogram with zero output frequency"); + cvm::error ("User required histogram with zero output frequency"); } grid = new colvar_grid_count (colvars); @@ -440,7 +449,7 @@ cvm::real colvarbias_histogram::update() if (cvm::debug()) cvm::log ("Histogram bias trying to write grid to disk"); grid_os.open (out_name.c_str()); - if (!grid_os.good()) cvm::fatal_error ("Error opening histogram file " + out_name + " for writing"); + if (!grid_os.good()) cvm::error ("Error opening histogram file " + out_name + " for writing"); grid->write_multicol (grid_os); grid_os.close (); } @@ -472,15 +481,15 @@ std::istream & colvarbias_histogram::read_restart (std::istream& is) std::string name = ""; if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) && (name != this->name) ) - cvm::fatal_error ("Error: in the restart file, the " + cvm::error ("Error: in the restart file, the " "\"histogram\" block has a wrong name: different system?\n"); if ( (id == -1) && (name == "") ) { - cvm::fatal_error ("Error: \"histogram\" block in the restart file " + cvm::error ("Error: \"histogram\" block in the restart file " "has no name.\n"); } if ( !(is >> key) || !(key == "grid")) { - cvm::log ("Error: in reading restart configuration for histogram \""+ + cvm::error ("Error: in reading restart configuration for histogram \""+ this->name+"\" at position "+ cvm::to_str (is.tellg())+" in stream.\n"); is.clear(); @@ -489,14 +498,17 @@ std::istream & colvarbias_histogram::read_restart (std::istream& is) return is; } if (! grid->read_raw (is)) { - grid->read_raw_error(); + is.clear(); + is.seekg (start_pos, std::ios::beg); + is.setstate (std::ios::failbit); + return is; } is >> brace; if (brace != "}") { - cvm::fatal_error ("Error: corrupt restart information for ABF bias \""+ - this->name+"\": no matching brace at position "+ - cvm::to_str (is.tellg())+" in the restart file.\n"); + cvm::error ("Error: corrupt restart information for ABF bias \""+ + this->name+"\": no matching brace at position "+ + cvm::to_str (is.tellg())+" in the restart file.\n"); is.setstate (std::ios::failbit); } return is; diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index e82d2f08ac..6873a8f9b1 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + /************************************************************************ * Headers for the ABF and histogram biases * ************************************************************************/ diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 41e01d9695..867203a8ab 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -13,6 +13,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : colvarbias(conf, key), update_calls(0), b_equilibration(true) { + size_t i; // get the initial restraint centers colvar_centers.resize (colvars.size()); @@ -29,7 +30,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : coupling_rate.resize(colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { + for (i = 0; i < colvars.size(); i++) { colvar_centers[i].type (colvars[i]->type()); //zero moments means[i] = ssd[i] = 0; @@ -39,7 +40,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : } if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) { - for (size_t i = 0; i < colvars.size(); i++) { + for (i = 0; i < colvars.size(); i++) { colvar_centers[i].apply_constraints(); } } else { @@ -52,7 +53,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : "that of collective variables.\n"); if(!get_keyval (conf, "UpdateFrequency", update_freq, 0)) - cvm::fatal_error("Error: must set updateFrequency for apadtive linear bias.\n"); + cvm::fatal_error("Error: must set updateFrequency for adaptive linear bias.\n"); //we split the time between updating and equilibrating update_freq /= 2; @@ -67,21 +68,21 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : //initial guess if(!get_keyval (conf, "forceConstant", set_coupling, set_coupling)) - for(size_t i =0 ; i < colvars.size(); i++) + for(i =0 ; i < colvars.size(); i++) set_coupling[i] = 0.; //how we're going to increase to that point - for(size_t i = 0; i < colvars.size(); i++) + for(i = 0; i < colvars.size(); i++) coupling_rate[i] = (set_coupling[i] - current_coupling[i]) / update_freq; if(!get_keyval (conf, "forceRange", max_coupling_range, max_coupling_range)) { //set to default - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { if(cvm::temperature() > 0) - max_coupling_range[i] = 3 * cvm::temperature() * cvm::boltzmann(); + max_coupling_range[i] = 3 * cvm::temperature() * cvm::boltzmann(); else - max_coupling_range[i] = 3 * cvm::boltzmann(); + max_coupling_range[i] = 3 * cvm::boltzmann(); } } @@ -89,7 +90,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : if(!get_keyval (conf, "rateMax", max_coupling_rate, max_coupling_rate)) { //set to default - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { max_coupling_rate[i] = max_coupling_range[i] / (10 * update_freq); } } @@ -123,9 +124,9 @@ cvm::real colvarbias_alb::update() { bool finished_equil_flag = 1; cvm::real delta; for (size_t i = 0; i < colvars.size(); i++) { - colvar_forces[i] = -restraint_force(restraint_convert_k(current_coupling[i], colvars[i]->width), - colvars[i], - colvar_centers[i]); + colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(current_coupling[i], colvars[i]->width), + colvars[i], + colvar_centers[i]); bias_energy += restraint_potential(restraint_convert_k(current_coupling[i], colvars[i]->width), colvars[i], colvar_centers[i]); @@ -193,7 +194,7 @@ cvm::real colvarbias_alb::update() { ssd[i] = 0; //stochastic if we do that update or not - if(colvars.size() == 1 || rand() < RAND_MAX / colvars.size()) { + if(colvars.size() == 1 || rand() < RAND_MAX / ((int) colvars.size())) { coupling_accum[i] += step_size * step_size; current_coupling[i] = set_coupling[i]; set_coupling[i] += max_coupling_range[i] / sqrt(coupling_accum[i]) * step_size; @@ -294,37 +295,38 @@ std::ostream & colvarbias_alb::write_restart (std::ostream &os) << " configuration {\n" << " name " << this->name << "\n"; os << " setCoupling "; - for(size_t i = 0; i < colvars.size(); i++) { + size_t i; + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << set_coupling[i] << "\n"; } os << " currentCoupling "; - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << current_coupling[i] << "\n"; } os << " maxCouplingRange "; - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << max_coupling_range[i] << "\n"; } os << " couplingRate "; - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << coupling_rate[i] << "\n"; } os << " couplingAccum "; - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << coupling_accum[i] << "\n"; } os << " mean "; - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << means[i] << "\n"; } os << " ssd "; - for(size_t i = 0; i < colvars.size(); i++) { + for(i = 0; i < colvars.size(); i++) { os << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) << ssd[i] << "\n"; } diff --git a/lib/colvars/colvarbias_alb.h b/lib/colvars/colvarbias_alb.h index 20ad20cf75..a79a9861e1 100644 --- a/lib/colvars/colvarbias_alb.h +++ b/lib/colvars/colvarbias_alb.h @@ -2,7 +2,7 @@ #define COLVARBIAS_ALB_H #include "colvar.h" -#include "colvarbias.h" +#include "colvarbias_restraint.h" class colvarbias_alb : public colvarbias { diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index cbdd2c20f6..b49aee87c0 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include #include #include @@ -64,7 +66,8 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) get_keyval (conf, "rebinGrids", rebin_grids, false); expand_grids = false; - for (size_t i = 0; i < colvars.size(); i++) { + size_t i; + for (i = 0; i < colvars.size(); i++) { if (colvars[i]->expand_boundaries) { expand_grids = true; cvm::log ("Metadynamics bias \""+this->name+"\""+ @@ -79,7 +82,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent); get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false); - for (size_t i = 0; i < colvars.size(); i++) { + for (i = 0; i < colvars.size(); i++) { colvars[i]->enable (colvar::task_grid); } @@ -91,6 +94,9 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) dump_fes = false; dump_fes_save = false; dump_replica_fes = false; + + hills_energy = NULL; + hills_energy_gradients = NULL; } if (comm != single_replica) { @@ -144,7 +150,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) replica_state_file = (std::string (pwd)+std::string (PATHSEP)+ cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state"); - delete pwd; + delete[] pwd; } // now register this replica @@ -316,7 +322,7 @@ colvarbias_meta::delete_hill (hill_iter &h) "")+".\n"); } - if (use_grids && hills_off_grid.size()) { + if (use_grids && !hills_off_grid.empty()) { for (hill_iter hoff = hills_off_grid.begin(); hoff != hills_off_grid.end(); hoff++) { if (*h == *hoff) { @@ -458,9 +464,14 @@ cvm::real colvarbias_meta::update() case single_replica: if (well_tempered) { - std::vector curr_bin = hills_energy->get_colvars_index(); - cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin); - cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); + cvm::real hills_energy_sum_here = 0.0; + if (use_grids) { + std::vector curr_bin = hills_energy->get_colvars_index(); + hills_energy_sum_here = hills_energy->value(curr_bin); + } else { + calc_hills (new_hills_begin, hills.end(), hills_energy_sum_here); + } + cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); create_hill (hill ((hill_weight*exp_weight), colvars, hill_width)); } else { create_hill (hill (hill_weight, colvars, hill_width)); @@ -469,9 +480,14 @@ cvm::real colvarbias_meta::update() case multiple_replicas: if (well_tempered) { - std::vector curr_bin = hills_energy->get_colvars_index(); - cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin); - cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); + cvm::real hills_energy_sum_here = 0.0; + if (use_grids) { + std::vector curr_bin = hills_energy->get_colvars_index(); + hills_energy_sum_here = hills_energy->value(curr_bin); + } else { + calc_hills (new_hills_begin, hills.end(), hills_energy_sum_here); + } + cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); create_hill (hill ((hill_weight*exp_weight), colvars, hill_width, replica_id)); } else { create_hill (hill (hill_weight, colvars, hill_width, replica_id)); @@ -680,10 +696,11 @@ void colvarbias_meta::calc_hills_force (size_t const &i, // were already saved with their types matching those in the // colvars) + hill_iter h; switch (colvars[i]->type()) { case colvarvalue::type_scalar: - for (hill_iter h = h_first; h != h_last; h++) { + for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; @@ -695,7 +712,8 @@ void colvarbias_meta::calc_hills_force (size_t const &i, case colvarvalue::type_vector: case colvarvalue::type_unitvector: - for (hill_iter h = h_first; h != h_last; h++) { + case colvarvalue::type_unitvectorderiv: + for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; @@ -705,8 +723,8 @@ void colvarbias_meta::calc_hills_force (size_t const &i, } break; - case colvarvalue::type_quaternion: - for (hill_iter h = h_first; h != h_last; h++) { + case colvarvalue::type_quaternionderiv: + for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; @@ -717,6 +735,7 @@ void colvarbias_meta::calc_hills_force (size_t const &i, break; case colvarvalue::type_notset: + case colvarvalue::type_all: break; } } @@ -756,8 +775,8 @@ void colvarbias_meta::project_hills (colvarbias_meta::hill_iter h_first, for ( ; (he->index_ok (he_ix)) && (hg->index_ok (hg_ix)); count++) { - - for (size_t i = 0; i < colvars.size(); i++) { + size_t i; + for (i = 0; i < colvars.size(); i++) { colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i); } @@ -766,7 +785,7 @@ void colvarbias_meta::project_hills (colvarbias_meta::hill_iter h_first, calc_hills (h_first, h_last, hills_energy_here, colvar_values); he->acc_value (he_ix, hills_energy_here); - for (size_t i = 0; i < colvars.size(); i++) { + for (i = 0; i < colvars.size(); i++) { hills_forces_here[i].reset(); calc_hills_force (i, h_first, h_last, hills_forces_here, colvar_values); colvar_forces_scalar[i] = hills_forces_here[i].real_value; @@ -935,10 +954,6 @@ void colvarbias_meta::update_replicas_registry() } } } - - // continue the cycle - new_replica_file = ""; - new_replica = ""; } else { cvm::fatal_error ("Error: cannot read the replicas registry file \""+ replicas_registry+"\".\n"); @@ -1279,7 +1294,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) } } - bool const existing_hills = (hills.size() > 0); + bool const existing_hills = !hills.empty(); size_t const old_hills_size = hills.size(); hill_iter old_hills_end = hills.end(); hill_iter old_hills_off_grid_end = hills_off_grid.end(); @@ -1299,7 +1314,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) cvm::log ("Read "+cvm::to_str (hills.size())+ " hills in addition to the grids.\n"); } else { - if (hills.size()) + if (!hills.empty()) cvm::log ("Read "+cvm::to_str (hills.size())+" hills.\n"); } @@ -1314,7 +1329,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) colvar_grid_gradient *new_hills_energy_gradients = new colvar_grid_gradient (colvars); - if (!grids_from_restart_file || (keep_hills && hills.size())) { + if (!grids_from_restart_file || (keep_hills && !hills.empty())) { // if there are hills, recompute the new grids from them cvm::log ("Rebinning the energy and forces grids from "+ cvm::to_str (hills.size())+" hills (this may take a while)...\n"); @@ -1337,12 +1352,12 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) // assuming that some boundaries have expanded, eliminate those // off-grid hills that aren't necessary any more - if (hills.size()) + if (!hills.empty()) recount_hills_off_grid (hills.begin(), hills.end(), hills_energy); } if (use_grids) { - if (hills_off_grid.size()) { + if (!hills_off_grid.empty()) { cvm::log (cvm::to_str (hills_off_grid.size())+" hills are near the " "grid boundaries: they will be computed analytically " "and saved to the state files.\n"); @@ -1656,15 +1671,16 @@ std::string colvarbias_meta::hill::output_traj() os.setf (std::ios::scientific, std::ios::floatfield); + size_t i; os << " "; - for (size_t i = 0; i < centers.size(); i++) { + for (i = 0; i < centers.size(); i++) { os << " "; os << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << centers[i]; } os << " "; - for (size_t i = 0; i < widths.size(); i++) { + for (i = 0; i < widths.size(); i++) { os << " "; os << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) << widths[i]; @@ -1692,8 +1708,9 @@ std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h) if (h.replica.size()) os << " replicaID " << h.replica << "\n"; + size_t i; os << " centers "; - for (size_t i = 0; i < (h.centers).size(); i++) { + for (i = 0; i < (h.centers).size(); i++) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) @@ -1702,7 +1719,7 @@ std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h) os << "\n"; os << " widths "; - for (size_t i = 0; i < (h.widths).size(); i++) { + for (i = 0; i < (h.widths).size(); i++) { os << " " << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 8c3eb91003..42faac9d5c 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #ifndef COLVARBIAS_META_H #define COLVARBIAS_META_H @@ -224,7 +226,7 @@ protected: std::ofstream replica_hills_os; /// Position within replica_hills_file (when reading it) - size_t replica_hills_file_pos; + int replica_hills_file_pos; }; @@ -417,9 +419,3 @@ public: #endif - - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index f509d8fcc3..f333f5db3e 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include "colvarmodule.h" #include "colvarvalue.h" #include "colvar.h" @@ -78,18 +80,14 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group) // this function should work for any scalar variable: // the only difference will be the name of the atom group (here, "group") - // collect into a vector for convenience - std::vector gradients (group.size()); - for (size_t i = 0; i < group.size(); i++) { - gradients[i] = group[i].grad; - } + if (group.b_dummy) return; cvm::rotation const rot_0 = group.rot; cvm::rotation const rot_inv = group.rot.inverse(); cvm::real const x_0 = x.real_value; - cvm::log ("gradients = "+cvm::to_str (gradients)+"\n"); + // cvm::log ("gradients = "+cvm::to_str (gradients)+"\n"); // it only makes sense to debug the fit gradients // when the fitting group is the same as this group @@ -131,7 +129,8 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group) cvm::log ("Atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n"); cvm::log ("dx(actual) = "+cvm::to_str (x_1 - x_0, 21, 14)+"\n"); - cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ? + //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ? + cvm::real const dx_pred = (group.fit_gradients.size()) ? (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) : (cvm::debug_gradients_step_size * atom_grad[id]); cvm::log ("dx(interp) = "+cvm::to_str (dx_pred, @@ -141,7 +140,45 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group) std::fabs (x_1 - x_0), 12, 5)+ ".\n"); - } } + +/* + * The code below is WIP + */ +// if (group.ref_pos_group != NULL) { +// cvm::atom_group &ref = *group.ref_pos_group; +// group.calc_fit_gradients(); +// +// for (size_t ia = 0; ia < ref.size(); ia++) { +// +// for (size_t id = 0; id < 3; id++) { +// // (re)read original positions +// group.read_positions(); +// ref.read_positions(); +// // change one coordinate +// ref[ia].pos[id] += cvm::debug_gradients_step_size; +// group.calc_apply_roto_translation(); +// calc_value(); +// cvm::real const x_1 = x.real_value; +// cvm::log ("refPosGroup atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n"); +// cvm::log ("dx(actual) = "+cvm::to_str (x_1 - x_0, +// 21, 14)+"\n"); +// //cvm::real const dx_pred = (group.fit_gradients.size() && (group.ref_pos_group == NULL)) ? +// // cvm::real const dx_pred = (group.fit_gradients.size()) ? +// // (cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) : +// // (cvm::debug_gradients_step_size * atom_grad[id]); +// cvm::real const dx_pred = cvm::debug_gradients_step_size * ref.fit_gradients[ia][id]; +// cvm::log ("dx(interp) = "+cvm::to_str (dx_pred, +// 21, 14)+"\n"); +// cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+ +// cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) / +// std::fabs (x_1 - x_0), +// 12, 5)+ +// ".\n"); +// } +// } +// } + + return; } diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 6a4315d9d1..253c1d8de7 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #ifndef COLVARCOMP_H #define COLVARCOMP_H @@ -153,6 +155,9 @@ public: /// \brief Return the previously calculated value virtual colvarvalue value() const; + /// \brief Return const pointer to the previously calculated value + virtual const colvarvalue *p_value() const; + /// \brief Return the previously calculated system force virtual colvarvalue system_force() const; @@ -261,6 +266,11 @@ inline colvarvalue colvar::cvc::value() const return x; } +inline const colvarvalue * colvar::cvc::p_value() const +{ + return &x; +} + inline colvarvalue colvar::cvc::system_force() const { return ft; @@ -296,8 +306,8 @@ inline cvm::real colvar::cvc::compare (colvarvalue const &x1, if (this->type() == colvarvalue::type_scalar) { return cvm::real (x1 - x2); } else { - cvm::fatal_error ("Error: you requested an operation which requires " - "comparison between two non-scalar values.\n"); + cvm::error ("Error: you requested an operation which requires " + "comparison between two non-scalar values.\n"); return 0.0; } } @@ -1050,6 +1060,31 @@ public: }; +/// \brief Colvar component: cosine of the angle of rotation with respect to a set +/// of reference coordinates (colvarvalue::type_scalar type, range +/// [-1:1]) +class colvar::orientation_proj + : public colvar::orientation +{ +public: + + orientation_proj (std::string const &conf); + orientation_proj(); + virtual inline ~orientation_proj() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force (colvarvalue const &force); + virtual cvm::real dist2 (colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_lgrad (colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_rgrad (colvarvalue const &x1, + colvarvalue const &x2) const; + virtual cvm::real compare (colvarvalue const &x1, + colvarvalue const &x2) const; +}; + + /// \brief Colvar component: projection of the orientation vector onto /// a predefined axis (colvarvalue::type_scalar type, range [-1:1]) class colvar::tilt @@ -1199,6 +1234,7 @@ public: simple_scalar_dist_functions (inertia_z) simple_scalar_dist_functions (rmsd) simple_scalar_dist_functions (orientation_angle) + simple_scalar_dist_functions (orientation_proj) simple_scalar_dist_functions (tilt) simple_scalar_dist_functions (eigenvector) // simple_scalar_dist_functions (alpha_dihedrals) @@ -1379,7 +1415,7 @@ inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1, inline cvm::real colvar::distance_vec::compare (colvarvalue const &x1, colvarvalue const &x2) const { - cvm::fatal_error ("Error: cannot compare() two distance vectors.\n"); + cvm::error ("Error: cannot compare() two distance vectors.\n"); return 0.0; } @@ -1404,7 +1440,7 @@ inline colvarvalue colvar::distance_dir::dist2_rgrad (colvarvalue const &x1, inline cvm::real colvar::distance_dir::compare (colvarvalue const &x1, colvarvalue const &x2) const { - cvm::fatal_error ("Error: cannot compare() two distance directions.\n"); + cvm::error ("Error: cannot compare() two distance directions.\n"); return 0.0; } @@ -1431,15 +1467,9 @@ inline colvarvalue colvar::orientation::dist2_rgrad (colvarvalue const &x1, inline cvm::real colvar::orientation::compare (colvarvalue const &x1, colvarvalue const &x2) const { - cvm::fatal_error ("Error: cannot compare() two quaternions.\n"); + cvm::error ("Error: cannot compare() two quaternions.\n"); return 0.0; } #endif - - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index bc2aee0598..7733242405 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -1,3 +1,4 @@ +/// -*- c++ -*- #include "colvarmodule.h" #include "colvar.h" @@ -6,7 +7,6 @@ #include - colvar::angle::angle (std::string const &conf) : cvc (conf) { @@ -75,6 +75,7 @@ void colvar::angle::calc_value() void colvar::angle::calc_gradients() { + size_t i; cvm::real const cos_theta = (r21*r23)/(r21l*r23l); cvm::real const dxdcos = -1.0 / std::sqrt (1.0 - cos_theta*cos_theta); @@ -84,17 +85,17 @@ void colvar::angle::calc_gradients() dxdr3 = (180.0/PI) * dxdcos * (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l ); - for (size_t i = 0; i < group1.size(); i++) { + for (i = 0; i < group1.size(); i++) { group1[i].grad = (group1[i].mass/group1.total_mass) * (dxdr1); } - for (size_t i = 0; i < group2.size(); i++) { + for (i = 0; i < group2.size(); i++) { group2[i].grad = (group2[i].mass/group2.total_mass) * (dxdr1 + dxdr3) * (-1.0); } - for (size_t i = 0; i < group3.size(); i++) { + for (i = 0; i < group3.size(); i++) { group3[i].grad = (group3[i].mass/group3.total_mass) * (dxdr3); } @@ -307,16 +308,17 @@ void colvar::dihedral::calc_gradients() +dsindB.y*r34.x - dsindB.x*r34.y); } - for (size_t i = 0; i < group1.size(); i++) + size_t i; + for (i = 0; i < group1.size(); i++) group1[i].grad = (group1[i].mass/group1.total_mass) * (-f1); - for (size_t i = 0; i < group2.size(); i++) + for (i = 0; i < group2.size(); i++) group2[i].grad = (group2[i].mass/group2.total_mass) * (-f2 + f1); - for (size_t i = 0; i < group3.size(); i++) + for (i = 0; i < group3.size(); i++) group3[i].grad = (group3[i].mass/group3.total_mass) * (-f3 + f2); - for (size_t i = 0; i < group4.size(); i++) + for (i = 0; i < group4.size(); i++) group4[i].grad = (group4[i].mass/group4.total_mass) * (f3); } diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index d08a5e505c..8b6448419f 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include #include "colvarmodule.h" @@ -276,7 +278,7 @@ colvar::h_bond::h_bond() colvar::h_bond::~h_bond() { - for (int i=0; i #include "colvarmodule.h" @@ -40,6 +42,7 @@ colvar::distance::distance() b_inverse_gradients = true; b_Jacobian_derivative = true; b_1site_force = false; + b_no_PBC = false; x.type (colvarvalue::type_scalar); } @@ -143,8 +146,9 @@ colvar::distance_z::distance_z (std::string const &conf) b_periodic = true; if ((wrap_center != 0.0) && (period == 0.0)) { - cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component," - " but its period has not been set.\n"); + cvm::error ("Error: wrapAround was defined in a distanceZ component," + " but its period has not been set.\n"); + return; } parse_group (conf, "main", main); @@ -162,8 +166,10 @@ colvar::distance_z::distance_z (std::string const &conf) cvm::log ("Warning: explicit axis definition will be ignored!"); } else { if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) { - if (axis.norm2() == 0.0) - cvm::fatal_error ("Axis vector is zero!"); + if (axis.norm2() == 0.0) { + cvm::error ("Axis vector is zero!"); + return; + } if (axis.norm2() != 1.0) { axis = axis.unit(); cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n"); @@ -234,6 +240,15 @@ void colvar::distance_z::calc_gradients() cvm::position_distance (main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis )); } } + + if (b_debug_gradients) { + cvm::log ("Debugging gradients for group main:\n"); + debug_gradients (main); + cvm::log ("Debugging gradients for group ref1:\n"); + debug_gradients (ref1); + cvm::log ("Debugging gradients for group ref2:\n"); + debug_gradients (ref2); + } } void colvar::distance_z::calc_force_invgrads() @@ -422,16 +437,19 @@ colvar::distance_inv::distance_inv (std::string const &conf) function_type = "distance_inv"; get_keyval (conf, "exponent", exponent, 6); if (exponent%2) { - cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n"); + cvm::error ("Error: odd exponent provided, can only use even ones.\n"); + return; } if (exponent <= 0) { - cvm::fatal_error ("Error: negative or zero exponent provided.\n"); + cvm::error ("Error: negative or zero exponent provided.\n"); + return; } for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { if (ai1->id == ai2->id) - cvm::fatal_error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n"); + cvm::error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n"); + return; } } @@ -640,8 +658,10 @@ colvar::inertia_z::inertia_z (std::string const &conf) { function_type = "inertia_z"; if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) { - if (axis.norm2() == 0.0) - cvm::fatal_error ("Axis vector is zero!"); + if (axis.norm2() == 0.0) { + cvm::error ("Axis vector is zero!"); + return; + } if (axis.norm2() != 1.0) { axis = axis.unit(); cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n"); @@ -700,8 +720,10 @@ colvar::rmsd::rmsd (std::string const &conf) parse_group (conf, "atoms", atoms); atom_groups.push_back (&atoms); - if (atoms.b_dummy) - cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom."); + if (atoms.b_dummy) { + cvm::error ("Error: \"atoms\" cannot be a dummy atom."); + return; + } if (atoms.ref_pos_group != NULL) { cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: " @@ -715,10 +737,11 @@ colvar::rmsd::rmsd (std::string const &conf) if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) { cvm::log ("Using reference positions from configuration file to calculate the variable.\n"); if (ref_pos.size() != atoms.size()) { - cvm::fatal_error ("Error: the number of reference positions provided ("+ - cvm::to_str (ref_pos.size())+ - ") does not match the number of atoms of group \"atoms\" ("+ - cvm::to_str (atoms.size())+").\n"); + cvm::error ("Error: the number of reference positions provided ("+ + cvm::to_str (ref_pos.size())+ + ") does not match the number of atoms of group \"atoms\" ("+ + cvm::to_str (atoms.size())+").\n"); + return; } } { @@ -726,8 +749,9 @@ colvar::rmsd::rmsd (std::string const &conf) if (get_keyval (conf, "refPositionsFile", ref_pos_file, std::string (""))) { if (ref_pos.size()) { - cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and " + cvm::error ("Error: cannot specify \"refPositionsFile\" and " "\"refPositions\" at the same time.\n"); + return; } std::string ref_pos_col; @@ -737,8 +761,9 @@ colvar::rmsd::rmsd (std::string const &conf) // if provided, use PDB column to select coordinates bool found = get_keyval (conf, "refPositionsColValue", ref_pos_col_value, 0.0); if (found && !ref_pos_col_value) - cvm::fatal_error ("Error: refPositionsColValue, " + cvm::error ("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); + return; } else { // if not, rely on existing atom indices for the group atoms.create_sorted_ids(); @@ -897,25 +922,30 @@ colvar::eigenvector::eigenvector (std::string const &conf) if (b_inline) { cvm::log ("Using reference positions from input file.\n"); if (ref_pos.size() != atoms.size()) { - cvm::fatal_error ("Error: reference positions do not " + cvm::error ("Error: reference positions do not " "match the number of requested atoms.\n"); + return; } } std::string file_name; if (get_keyval (conf, "refPositionsFile", file_name)) { - if (b_inline) - cvm::fatal_error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n"); + if (b_inline) { + cvm::error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n"); + return; + } std::string file_col; double file_col_value; if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) { // use PDB flags if column is provided bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0); - if (found && !file_col_value) - cvm::fatal_error ("Error: refPositionsColValue, " + if (found && !file_col_value) { + cvm::error ("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); + return; + } } else { // if not, use atom indices atoms.create_sorted_ids(); @@ -965,17 +995,20 @@ colvar::eigenvector::eigenvector (std::string const &conf) std::string file_name; if (get_keyval (conf, "vectorFile", file_name)) { - if (b_inline) - cvm::fatal_error ("Error: vector and vectorFile cannot be specified at the same time.\n"); + if (b_inline) { + cvm::error ("Error: vector and vectorFile cannot be specified at the same time.\n"); + return; + } std::string file_col; double file_col_value; if (get_keyval (conf, "vectorCol", file_col, std::string (""))) { // use PDB flags if column is provided bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0); - if (found && !file_col_value) - cvm::fatal_error ("Error: vectorColValue, " - "if provided, must be non-zero.\n"); + if (found && !file_col_value) { + cvm::error ("Error: vectorColValue, if provided, must be non-zero.\n"); + return; + } } else { // if not, use atom indices atoms.create_sorted_ids(); @@ -987,13 +1020,14 @@ colvar::eigenvector::eigenvector (std::string const &conf) } if (!ref_pos.size() || !eigenvec.size()) { - cvm::fatal_error ("Error: both reference coordinates" + cvm::error ("Error: both reference coordinates" "and eigenvector must be defined.\n"); + return; } cvm::atom_pos eig_center (0.0, 0.0, 0.0); - for (size_t i = 0; i < atoms.size(); i++) { - eig_center += eigenvec[i]; + for (size_t eil = 0; eil < atoms.size(); eil++) { + eig_center += eigenvec[eil]; } eig_center *= 1.0 / atoms.size(); cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n"); @@ -1038,8 +1072,8 @@ colvar::eigenvector::eigenvector (std::string const &conf) // for inverse gradients eigenvec_invnorm2 = 0.0; - for (size_t i = 0; i < atoms.size(); i++) { - eigenvec_invnorm2 += eigenvec[i].norm2(); + for (size_t ein = 0; ein < atoms.size(); ein++) { + eigenvec_invnorm2 += eigenvec[ein].norm2(); } eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2; diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index 1750782b7a..5021f60598 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include #include "colvarmodule.h" @@ -162,10 +164,11 @@ void colvar::alpha_angles::calc_value() void colvar::alpha_angles::calc_gradients() { - for (size_t i = 0; i < theta.size(); i++) + size_t i; + for (i = 0; i < theta.size(); i++) (theta[i])->calc_gradients(); - for (size_t i = 0; i < hb.size(); i++) + for (i = 0; i < hb.size(); i++) (hb[i])->calc_gradients(); } diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index b77c563878..0c52b989bf 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include #include "colvarmodule.h" @@ -57,11 +59,12 @@ colvar::orientation::orientation (std::string const &conf) "assumed that each atom is the closest " "periodic image to the center of geometry.\n"); cvm::rvector cog (0.0, 0.0, 0.0); - for (size_t i = 0; i < ref_pos.size(); i++) { + size_t i; + for (i = 0; i < ref_pos.size(); i++) { cog += ref_pos[i]; } cog /= cvm::real (ref_pos.size()); - for (size_t i = 0; i < ref_pos.size(); i++) { + for (i = 0; i < ref_pos.size(); i++) { ref_pos[i] -= cog; } @@ -167,10 +170,61 @@ void colvar::orientation_angle::calc_gradients() for (size_t ia = 0; ia < atoms.size(); ia++) { atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); } + if (b_debug_gradients) { + cvm::log ("Debugging orientationAngle component gradients:\n"); + debug_gradients (atoms); + } } void colvar::orientation_angle::apply_force (colvarvalue const &force) +{ + cvm::real const &fw = force.real_value; + if (!atoms.noforce) { + atoms.apply_colvar_force (fw); + } +} + + + +colvar::orientation_proj::orientation_proj (std::string const &conf) + : orientation (conf) +{ + function_type = "orientation_proj"; + x.type (colvarvalue::type_scalar); +} + + +colvar::orientation_proj::orientation_proj() + : orientation() +{ + function_type = "orientation_proj"; + x.type (colvarvalue::type_scalar); +} + + +void colvar::orientation_proj::calc_value() +{ + atoms_cog = atoms.center_of_geometry(); + rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog)); + x.real_value = 2.0 * (rot.q).q0 * (rot.q).q0 - 1.0; +} + + +void colvar::orientation_proj::calc_gradients() +{ + cvm::real const dxdq0 = 2.0 * 2.0 * (rot.q).q0; + for (size_t ia = 0; ia < atoms.size(); ia++) { + atoms[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]); + } + if (b_debug_gradients) { + cvm::log ("Debugging orientationProj component gradients:\n"); + debug_gradients (atoms); + } +} + + +void colvar::orientation_proj::apply_force (colvarvalue const &force) { cvm::real const &fw = force.real_value; @@ -180,6 +234,7 @@ void colvar::orientation_angle::apply_force (colvarvalue const &force) } + colvar::tilt::tilt (std::string const &conf) : orientation (conf) { diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index a300b91340..1035cc5ac5 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -1,4 +1,4 @@ -// -*- c++ -*- +/// -*- c++ -*- #include "colvarmodule.h" #include "colvarvalue.h" diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 1095e0b845..c3a61dcaa0 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -1,4 +1,4 @@ -// -*- c++ -*- +/// -*- c++ -*- #ifndef COLVARGRID_H #define COLVARGRID_H @@ -55,8 +55,10 @@ protected: for (size_t i = 0; i < nd; i++) { addr += ix[i]*nxc[i]; if (cvm::debug()) { - if (ix[i] >= nx[i]) - cvm::fatal_error ("Error: exceeding bounds in colvar_grid.\n"); + if (ix[i] >= nx[i]) { + cvm::error ("Error: exceeding bounds in colvar_grid.\n"); + return 0; + } } } return addr; @@ -135,9 +137,11 @@ public: nt = mult; for (int i = nd-1; i >= 0; i--) { - if (nx_i[i] <= 0) - cvm::fatal_error ("Error: providing an invalid number of points, "+ + if (nx_i[i] <= 0) { + cvm::error ("Error: providing an invalid number of points, "+ cvm::to_str (nx_i[i])+".\n"); + return; + } nxc[i] = nt; nt *= nx[i]; } @@ -220,20 +224,23 @@ public: for (size_t i = 0; i < cv.size(); i++) { if (cv[i]->type() != colvarvalue::type_scalar) { - cvm::fatal_error ("Colvar grids can only be automatically " - "constructed for scalar variables. " - "ABF and histogram can not be used; metadynamics " - "can be used with useGrids disabled.\n"); + cvm::error ("Colvar grids can only be automatically " + "constructed for scalar variables. " + "ABF and histogram can not be used; metadynamics " + "can be used with useGrids disabled.\n"); + return; } if (cv[i]->width <= 0.0) { - cvm::fatal_error ("Tried to initialize a grid on a " + cvm::error ("Tried to initialize a grid on a " "variable with negative or zero width.\n"); + return; } if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) { - cvm::fatal_error ("Tried to initialize a grid on a " - "variable with undefined boundaries.\n"); + cvm::error ("Tried to initialize a grid on a " + "variable with undefined boundaries.\n"); + return; } widths.push_back (cv[i]->width); @@ -303,8 +310,9 @@ public: ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result } else { if (ix[i] < 0 || ix[i] >= nx[i]) - cvm::fatal_error ("Trying to wrap illegal index vector (non-PBC): " - + cvm::to_str (ix)); + cvm::error ("Trying to wrap illegal index vector (non-PBC): " + + cvm::to_str (ix)); + return; } } } @@ -437,9 +445,11 @@ public: /// of whether it fits or not. void map_grid (colvar_grid const &other_grid) { - if (other_grid.multiplicity() != this->multiplicity()) - cvm::fatal_error ("Error: trying to merge two grids with values of " + if (other_grid.multiplicity() != this->multiplicity()) { + cvm::error ("Error: trying to merge two grids with values of " "different multiplicity.\n"); + return; + } std::vector const &gb = this->lower_boundaries; std::vector const &gw = this->widths; @@ -479,9 +489,11 @@ public: void add_grid (colvar_grid const &other_grid, cvm::real scale_factor = 1.0) { - if (other_grid.multiplicity() != this->multiplicity()) - cvm::fatal_error ("Error: trying to sum togetehr two grids with values of " + if (other_grid.multiplicity() != this->multiplicity()) { + cvm::error ("Error: trying to sum togetehr two grids with values of " "different multiplicity.\n"); + return; + } if (scale_factor != 1.0) for (size_t i = 0; i < data.size(); i++) { data[i] += scale_factor * other_grid.data[i]; @@ -570,25 +582,26 @@ public: /// \brief Write the grid parameters (number of colvars, boundaries, width and number of points) std::ostream & write_params (std::ostream &os) { + size_t i; os << "grid_parameters {\n n_colvars " << nd << "\n"; os << " lower_boundaries "; - for (size_t i = 0; i < nd; i++) + for (i = 0; i < nd; i++) os << " " << lower_boundaries[i]; os << "\n"; os << " upper_boundaries "; - for (size_t i = 0; i < nd; i++) + for (i = 0; i < nd; i++) os << " " << upper_boundaries[i]; os << "\n"; os << " widths "; - for (size_t i = 0; i < nd; i++) + for (i = 0; i < nd; i++) os << " " << widths[i]; os << "\n"; os << " sizes "; - for (size_t i = 0; i < nd; i++) + for (i = 0; i < nd; i++) os << " " << nx[i]; os << "\n"; @@ -605,12 +618,14 @@ public: { size_t nd_in = 0; colvarparse::get_keyval (conf, "n_colvars", nd_in, nd, colvarparse::parse_silent); - if (nd_in != nd) - cvm::fatal_error ("Error: trying to read data for a grid " - "that contains a different number of colvars ("+ - cvm::to_str (nd_in)+") than the grid defined " - "in the configuration file ("+cvm::to_str (nd)+ - ").\n"); + if (nd_in != nd) { + cvm::error ("Error: trying to read data for a grid " + "that contains a different number of colvars ("+ + cvm::to_str (nd_in)+") than the grid defined " + "in the configuration file ("+cvm::to_str (nd)+ + ").\n"); + return false; + } } colvarparse::get_keyval (conf, "lower_boundaries", @@ -653,8 +668,9 @@ public: upper_boundaries[i])) > 1.0E-10) || (std::sqrt (cv[i]->dist2 (cv[i]->width, widths[i])) > 1.0E-10) ) { - cvm::fatal_error ("Error: restart information for a grid is " - "inconsistent with that of its colvars.\n"); + cvm::error ("Error: restart information for a grid is " + "inconsistent with that of its colvars.\n"); + return; } } } @@ -675,9 +691,10 @@ public: (std::fabs (other_grid.widths[i] - widths[i]) > 1.0E-10) || (data.size() != other_grid.data.size()) ) { - cvm::fatal_error ("Error: inconsistency between " - "two grids that are supposed to be equal, " - "aside from the data stored.\n"); + cvm::error ("Error: inconsistency between " + "two grids that are supposed to be equal, " + "aside from the data stored.\n"); + return; } } } @@ -710,169 +727,167 @@ public: return os; } -/// \brief Read data written by colvar_grid::write_raw() -std::istream & read_raw (std::istream &is) -{ - size_t const start_pos = is.tellg(); - - for (std::vector ix = new_index(); index_ok (ix); incr (ix)) { - for (size_t imult = 0; imult < mult; imult++) { - T new_value; - if (is >> new_value) { - value_input (ix, new_value, imult); - } else { - is.clear(); - is.seekg (start_pos, std::ios::beg); - is.setstate (std::ios::failbit); - return is; - } - } - } - - has_data = true; - return is; -} - -/// \brief To be called after colvar_grid::read_raw() returns an error -void read_raw_error() -{ - cvm::fatal_error ("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n"); -} - -/// \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) -{ - std::streamsize const w = os.width(); - std::streamsize const p = os.precision(); - - // Data in the header: nColvars, then for each - // xiMin, dXi, nPoints, periodic - - 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"; - } - - 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 (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 - - std::string hash; - cvm::real lower, width, x; - size_t n, periodic; - bool remap; - std::vector new_value; - std::vector nx_read; - std::vector bin; - - if ( cv.size() != nd ) { - cvm::fatal_error ("Cannot read grid file: missing reference to colvars."); - } - - if ( !(is >> hash) || (hash != "#") ) { - cvm::fatal_error ("Error reading grid at position "+ - cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n"); - } - - is >> n; - if ( n != nd ) { - cvm::fatal_error ("Error reading grid: wrong number of collective variables.\n"); - } - - 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::fatal_error ("Error reading grid at position "+ - cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n"); - } - - is >> lower >> width >> nx_read[i] >> periodic; - - - if ( (std::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) || - (std::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; + /// \brief Read data written by colvar_grid::write_raw() + std::istream & read_raw (std::istream &is) + { + size_t const start_pos = is.tellg(); + for (std::vector ix = new_index(); index_ok (ix); incr (ix)) { 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); + T new_value; + if (is >> new_value) { + value_input (ix, new_value, imult); + } else { + is.clear(); + is.seekg (start_pos, std::ios::beg); + is.setstate (std::ios::failbit); + cvm::error ("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n"); + return is; } } } - } else { - // do not re-grid the data but assume the same grid is used + + has_data = true; + 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) + { + std::streamsize const w = os.width(); + std::streamsize const p = os.precision(); + + // Data in the header: nColvars, then for each + // xiMin, dXi, nPoints, periodic + + 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"; + } + for (std::vector ix = new_index(); index_ok (ix); incr (ix) ) { - for (size_t i = 0; i < nd; i++ ) { - is >> x; + + 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++) { - is >> new_value[imult]; - value_input (ix, new_value[imult], imult, add); + os << " " + << std::setw (w) << std::setprecision (p) + << value_output (ix, imult); } + os << "\n"; } } - has_data = true; - return is; -} + /// \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 + + std::string hash; + cvm::real lower, width, x; + size_t n, periodic; + bool remap; + std::vector new_value; + std::vector nx_read; + std::vector bin; + + if ( cv.size() != nd ) { + cvm::error ("Cannot read grid file: missing reference to colvars."); + return is; + } + + if ( !(is >> hash) || (hash != "#") ) { + cvm::error ("Error reading grid at position "+ + cvm::to_str (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 (is.tellg())+" in stream (read \"" + hash + "\")\n"); + return is; + } + + is >> lower >> width >> nx_read[i] >> periodic; + + + if ( (std::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) || + (std::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; + } }; @@ -981,8 +996,11 @@ public: { cvm::real A0, A1; std::vector ix; - if (nd != 2) cvm::fatal_error ("Finite differences available in dimension 2 only."); - for (int n = 0; n < nd; n++) { + if (nd != 2) { + cvm::error ("Finite differences available in dimension 2 only."); + return grad; + } + for (unsigned int n = 0; n < nd; n++) { ix = ix0; A0 = data[address (ix)]; ix[n]++; wrap (ix); @@ -1001,9 +1019,11 @@ public: virtual cvm::real value_output (std::vector const &ix, size_t const &imult = 0) { - if (imult > 0) - cvm::fatal_error ("Error: trying to access a component " - "larger than 1 in a scalar data grid.\n"); + if (imult > 0) { + cvm::error ("Error: trying to access a component " + "larger than 1 in a scalar data grid.\n"); + return 0.; + } if (samples) return (samples->value (ix) > 0) ? (data[address (ix)] / cvm::real (samples->value (ix))) : @@ -1020,9 +1040,11 @@ public: size_t const &imult = 0, bool add = false) { - if (imult > 0) - cvm::fatal_error ("Error: trying to access a component " - "larger than 1 in a scalar data grid.\n"); + if (imult > 0) { + cvm::error ("Error: trying to access a component " + "larger than 1 in a scalar data grid.\n"); + return; + } if (add) { if (samples) data[address (ix)] += new_value * samples->new_count (ix); diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 3420618605..3b85273a60 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include #include "colvarmodule.h" #include "colvarparse.h" @@ -7,183 +9,164 @@ #include "colvarbias_alb.h" #include "colvarbias_meta.h" #include "colvarbias_abf.h" +#include "colvarbias_restraint.h" +#include "colvarscript.h" - -colvarmodule::colvarmodule (char const *config_filename, - colvarproxy *proxy_in) +colvarmodule::colvarmodule (colvarproxy *proxy_in) { // pointer to the proxy object if (proxy == NULL) { proxy = proxy_in; parse = new colvarparse(); } else { - cvm::fatal_error ("Error: trying to allocate the collective " + // TODO relax this error to handle multiple molecules in VMD + // once the module is not static anymore + cvm::error ("Error: trying to allocate the collective " "variable module twice.\n"); + return; } - cvm::log (cvm::line_marker); cvm::log ("Initializing the collective variables module, version "+ - cvm::to_str(COLVARS_VERSION)+".\n"); + cvm::to_str (COLVARS_VERSION)+".\n"); + + // set initial default values // "it_restart" will be set by the input state file, if any; // "it" should be updated by the proxy - it = it_restart = 0; - it_restart_from_state_file = true; + colvarmodule::it = colvarmodule::it_restart = 0; + colvarmodule::it_restart_from_state_file = true; + + colvarmodule::use_scripted_forces = false; + + colvarmodule::b_analysis = false; + colvarmodule::debug_gradients_step_size = 1.0e-07; + colvarmodule::rotation::crossing_threshold = 1.0e-02; + + colvarmodule::cv_traj_freq = 100; + colvarmodule::restart_out_freq = proxy->restart_frequency(); + + // by default overwrite the existing trajectory file + colvarmodule::cv_traj_append = false; +} + + +int colvarmodule::config_file (char const *config_filename) +{ + cvm::log (cvm::line_marker); + cvm::log ("Reading new configuration from file \""+ + std::string (config_filename)+"\":\n"); // open the configfile config_s.open (config_filename); - if (!config_s) - cvm::fatal_error ("Error: in opening configuration file \""+ - std::string (config_filename)+"\".\n"); + if (!config_s) { + cvm::error ("Error: in opening configuration file \""+ + std::string (config_filename)+"\".\n", + FILE_ERROR); + return COLVARS_ERROR; + } // read the config file into a string std::string conf = ""; - { - std::string line; - while (colvarparse::getline_nocomments (config_s, line)) - conf.append (line+"\n"); - // don't need the stream any more - config_s.close(); + std::string line; + while (colvarparse::getline_nocomments (config_s, line)) { + conf.append (line+"\n"); + } + config_s.close(); + + return config (conf); +} + + +int colvarmodule::config_string (std::string const &config_str) +{ + cvm::log (cvm::line_marker); + cvm::log ("Reading new configuration:\n"); + std::istringstream config_s (config_str); + + // strip the comments away + std::string conf = ""; + std::string line; + while (colvarparse::getline_nocomments (config_s, line)) { + conf.append (line+"\n"); + } + return config (conf); +} + +int colvarmodule::config (std::string &conf) +{ + int error_code = 0; + + // parse global options + error_code |= parse_global_params (conf); + + if (error_code != COLVARS_OK) { + set_error_bits(INPUT_ERROR); + return COLVARS_ERROR; } + // parse the options for collective variables + error_code |= parse_colvars (conf); + + // parse the options for biases + error_code |= parse_biases (conf); + + // done parsing known keywords, check that all keywords found were valid ones + error_code |= parse->check_keywords (conf, "colvarmodule"); + + if (error_code != COLVARS_OK) { + set_error_bits(INPUT_ERROR); + return COLVARS_ERROR; + } + + cvm::log (cvm::line_marker); + cvm::log ("Collective variables module (re)initialized.\n"); + cvm::log (cvm::line_marker); + + // configuration might have changed, better redo the labels + write_traj_label(cv_traj_os); + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + + +int colvarmodule::parse_global_params (std::string const &conf) +{ std::string index_file_name; if (parse->get_keyval (conf, "indexFile", index_file_name)) { read_index_file (index_file_name.c_str()); } - parse->get_keyval (conf, "analysis", b_analysis, false); + parse->get_keyval (conf, "analysis", b_analysis, b_analysis); - parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07, + parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, + debug_gradients_step_size, colvarparse::parse_silent); parse->get_keyval (conf, "eigenvalueCrossingThreshold", - colvarmodule::rotation::crossing_threshold, 1.0e-02, + colvarmodule::rotation::crossing_threshold, + colvarmodule::rotation::crossing_threshold, colvarparse::parse_silent); - parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, 100); - parse->get_keyval (conf, "colvarsRestartFrequency", restart_out_freq, - proxy->restart_frequency()); + parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq); + parse->get_keyval (conf, "colvarsRestartFrequency", + restart_out_freq, restart_out_freq); - // by default overwrite the existing trajectory file - parse->get_keyval (conf, "colvarsTrajAppend", cv_traj_append, false); + // if this is true when initializing, it means + // we are continuing after a reset(): default to true + parse->get_keyval (conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append); - // input restart file - restart_in_name = proxy->input_prefix().size() ? - std::string (proxy->input_prefix()+".colvars.state") : - std::string ("") ; + parse->get_keyval (conf, "scriptedColvarForces", use_scripted_forces, false, + colvarparse::parse_silent); - // output restart file - restart_out_name = proxy->restart_output_prefix().size() ? - std::string (proxy->restart_output_prefix()+".colvars.state") : - std::string (""); - - if (restart_out_name.size()) - cvm::log ("The restart output state file will be \""+restart_out_name+"\".\n"); - - output_prefix = proxy->output_prefix(); - - cvm::log ("The final output state file will be \""+ - (output_prefix.size() ? - std::string (output_prefix+".colvars.state") : - std::string ("colvars.state"))+"\".\n"); - - cv_traj_name = - (output_prefix.size() ? - std::string (output_prefix+".colvars.traj") : - std::string ("colvars.traj")); - - if (cv_traj_freq) { - // open trajectory file - if (cv_traj_append) { - cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+ - "\".\n"); - cv_traj_os.open (cv_traj_name.c_str(), std::ios::app); - } else { - cvm::log ("Writing to colvar trajectory file \""+cv_traj_name+ - "\".\n"); - proxy->backup_file (cv_traj_name.c_str()); - cv_traj_os.open (cv_traj_name.c_str(), std::ios::out); - } - cv_traj_os.setf (std::ios::scientific, std::ios::floatfield); + if (use_scripted_forces && !proxy->force_script_defined) { + cvm::fatal_error("User script for scripted colvars forces not found."); } - // parse the options for collective variables - init_colvars (conf); - - // parse the options for biases - init_biases (conf); - - // done with the parsing, check that all keywords are valid - parse->check_keywords (conf, "colvarmodule"); - cvm::log (cvm::line_marker); - - // read the restart configuration, if available - if (restart_in_name.size()) { - // read the restart file - std::ifstream input_is (restart_in_name.c_str()); - if (!input_is.good()) - fatal_error ("Error: in opening restart file \""+ - std::string (restart_in_name)+"\".\n"); - else { - cvm::log ("Restarting from file \""+restart_in_name+"\".\n"); - read_restart (input_is); - cvm::log (cvm::line_marker); - } - } - - // check if it is possible to save output configuration - if ((!output_prefix.size()) && (!restart_out_name.size())) { - cvm::fatal_error ("Error: neither the final output state file or " - "the output restart file could be defined, exiting.\n"); - } - - cvm::log ("Collective variables module initialized.\n"); - cvm::log (cvm::line_marker); + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -std::istream & colvarmodule::read_restart (std::istream &is) -{ - { - // read global restart information - std::string restart_conf; - if (is >> colvarparse::read_block ("configuration", restart_conf)) { - if (it_restart_from_state_file) { - parse->get_keyval (restart_conf, "step", - it_restart, (size_t) 0, - colvarparse::parse_silent); - it = it_restart; - } - } - is.clear(); - } - - // colvars restart - cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { - if ( !((*cvi)->read_restart (is)) ) - cvm::fatal_error ("Error: in reading restart configuration for collective variable \""+ - (*cvi)->name+"\".\n"); - } - - // biases restart - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { - if (!((*bi)->read_restart (is))) - fatal_error ("Error: in reading restart configuration for bias \""+ - (*bi)->name+"\".\n"); - } - cvm::decrease_depth(); - - return is; -} - - - -void colvarmodule::init_colvars (std::string const &conf) +int colvarmodule::parse_colvars (std::string const &conf) { if (cvm::debug()) cvm::log ("Initializing the collective variables.\n"); @@ -196,27 +179,37 @@ void colvarmodule::init_colvars (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); colvars.push_back (new colvar (colvar_conf)); - (colvars.back())->check_keywords (colvar_conf, "colvar"); + if (cvm::get_error()) { + delete colvars.back(); + colvars.pop_back(); + return COLVARS_ERROR; + } + if ((colvars.back())->check_keywords (colvar_conf, "colvar") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); } else { - cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n"); + cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } colvar_conf = ""; } - - if (!colvars.size()) - cvm::fatal_error ("Error: no collective variables defined.\n"); + if (!colvars.size()) { + cvm::log ("Warning: no collective variables defined.\n"); + } if (colvars.size()) cvm::log (cvm::line_marker); cvm::log ("Collective variables initialized, "+ cvm::to_str (colvars.size())+ " in total.\n"); + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvarmodule::init_biases (std::string const &conf) +int colvarmodule::parse_biases (std::string const &conf) { if (cvm::debug()) cvm::log ("Initializing the collective variables biases.\n"); @@ -230,11 +223,19 @@ void colvarmodule::init_biases (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_abf (abf_conf, "abf")); - (biases.back())->check_keywords (abf_conf, "abf"); + if (cvm::get_error()) { + delete biases.back(); + biases.pop_back(); + return COLVARS_ERROR; + } + if ((biases.back())->check_keywords (abf_conf, "abf") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); n_abf_biases++; } else { - cvm::log ("Warning: \"abf\" keyword found without configuration.\n"); + cvm::error("Error: \"abf\" keyword found without configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } abf_conf = ""; } @@ -249,11 +250,19 @@ void colvarmodule::init_biases (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_restraint_harmonic (harm_conf, "harmonic")); - (biases.back())->check_keywords (harm_conf, "harmonic"); + if (cvm::get_error()) { + delete biases.back(); + biases.pop_back(); + return COLVARS_ERROR; + } + if ((biases.back())->check_keywords (harm_conf, "harmonic") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); n_rest_biases++; } else { - cvm::log ("Warning: \"harmonic\" keyword found without configuration.\n"); + cvm::error("Error: \"harmonic\" keyword found without configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } harm_conf = ""; } @@ -268,11 +277,19 @@ void colvarmodule::init_biases (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_restraint_linear (lin_conf, "linear")); - (biases.back())->check_keywords (lin_conf, "linear"); + if (cvm::get_error()) { + delete biases.back(); + biases.pop_back(); + return COLVARS_ERROR; + } + if ((biases.back())->check_keywords (lin_conf, "linear") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); n_rest_biases++; } else { - cvm::log ("Warning: \"linear\" keyword found without configuration.\n"); + cvm::error("Error: \"linear\" keyword found without configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } lin_conf = ""; } @@ -287,11 +304,19 @@ void colvarmodule::init_biases (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_alb (alb_conf, "ALB")); - (biases.back())->check_keywords (alb_conf, "ALB"); + if (cvm::get_error()) { + delete biases.back(); + biases.pop_back(); + return COLVARS_ERROR; + } + if ((biases.back())->check_keywords (alb_conf, "ALB") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); n_rest_biases++; } else { - cvm::log ("Warning: \"ALB\" keyword found without configuration.\n"); + cvm::error("Error: \"ALB\" keyword found without configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } alb_conf = ""; } @@ -307,11 +332,19 @@ void colvarmodule::init_biases (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_histogram (histo_conf, "histogram")); - (biases.back())->check_keywords (histo_conf, "histogram"); + if (cvm::get_error()) { + delete biases.back(); + biases.pop_back(); + return COLVARS_ERROR; + } + if ((biases.back())->check_keywords (histo_conf, "histogram") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); n_histo_biases++; } else { - cvm::log ("Warning: \"histogram\" keyword found without configuration.\n"); + cvm::error("Error: \"histogram\" keyword found without configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } histo_conf = ""; } @@ -326,57 +359,87 @@ void colvarmodule::init_biases (std::string const &conf) cvm::log (cvm::line_marker); cvm::increase_depth(); biases.push_back (new colvarbias_meta (meta_conf, "metadynamics")); - (biases.back())->check_keywords (meta_conf, "metadynamics"); + if (cvm::get_error()) { + delete biases.back(); + biases.pop_back(); + return COLVARS_ERROR; + } + if ((biases.back())->check_keywords (meta_conf, "metadynamics") != COLVARS_OK) { + return COLVARS_ERROR; + } cvm::decrease_depth(); n_meta_biases++; } else { - cvm::log ("Warning: \"metadynamics\" keyword found without configuration.\n"); + cvm::error("Error: \"metadynamics\" keyword found without configuration.\n", INPUT_ERROR); + return COLVARS_ERROR; } meta_conf = ""; } } - if (biases.size()) + if (use_scripted_forces) { + cvm::log (cvm::line_marker); + cvm::increase_depth(); + cvm::log("User forces script will be run at each bias update."); + cvm::decrease_depth(); + } + + if (biases.size() || use_scripted_forces) cvm::log (cvm::line_marker); cvm::log ("Collective variables biases initialized, "+ cvm::to_str (biases.size())+" in total.\n"); + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvarmodule::change_configuration (std::string const &bias_name, - std::string const &conf) -{ - cvm::increase_depth(); - int found = 0; +colvarbias * colvarmodule::bias_by_name(std::string const &name) { for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { - if ( (*bi)->name == bias_name ) { - ++found; - (*bi)->change_configuration (conf); + if ((*bi)->name == name) { + return (*bi); } } - if (found < 1) cvm::fatal_error ("Error: bias not found.\n"); - if (found > 1) cvm::fatal_error ("Error: duplicate bias name.\n"); + return NULL; +} + + +colvar *colvarmodule::colvar_by_name(std::string const &name) { + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); + cvi++) { + if ((*cvi)->name == name) { + return (*cvi); + } + } + return NULL; +} + + +int colvarmodule::change_configuration (std::string const &bias_name, + std::string const &conf) +{ + // This is deprecated; supported strategy is to delete the bias + // and parse the new config + cvm::increase_depth(); + colvarbias *b; + b = bias_by_name (bias_name); + if (b == NULL) { cvm::error ("Error: bias not found: " + bias_name); } + b->change_configuration (conf); cvm::decrease_depth(); + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::string colvarmodule::read_colvar(std::string const &name) { cvm::increase_depth(); - int found = 0; + colvar *c; std::stringstream ss; - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { - if ( (*cvi)->name == name ) { - ++found; - ss << (*cvi)->value(); - } - } - if ( found < 1 ) cvm::fatal_error ("Error: colvar not found"); - if ( found > 1 ) cvm::fatal_error ("Error: duplicate colvar"); + c = colvar_by_name (name); + if (c == NULL) { cvm::fatal_error ("Error: colvar not found: " + name); } + ss << c->value(); cvm::decrease_depth(); return ss.str(); } @@ -386,27 +449,24 @@ cvm::real colvarmodule::energy_difference (std::string const &bias_name, std::string const &conf) { cvm::increase_depth(); + colvarbias *b; cvm::real energy_diff = 0.; - int found = 0; - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { - if ( (*bi)->name == bias_name ) { - ++found; - energy_diff = (*bi)->energy_difference (conf); - } - } - if (found < 1) cvm::fatal_error ("Error: bias not found.\n"); - if (found > 1) cvm::fatal_error ("Error: duplicate bias name.\n"); + b = bias_by_name (bias_name); + if (b == NULL) { cvm::fatal_error ("Error: bias not found: " + bias_name); } + energy_diff = b->energy_difference (conf); cvm::decrease_depth(); return energy_diff; } -void colvarmodule::calc() { + +int colvarmodule::calc() { cvm::real total_bias_energy = 0.0; cvm::real total_colvar_energy = 0.0; + std::vector::iterator cvi; + std::vector::iterator bi; + if (cvm::debug()) { cvm::log (cvm::line_marker); cvm::log ("Collective variables module, step no. "+ @@ -417,10 +477,11 @@ void colvarmodule::calc() { if (cvm::debug()) cvm::log ("Calculating collective variables.\n"); cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->calc(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } } cvm::decrease_depth(); @@ -429,10 +490,11 @@ void colvarmodule::calc() { if (cvm::debug() && biases.size()) cvm::log ("Updating collective variable biases.\n"); cvm::increase_depth(); - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { + for (bi = biases.begin(); bi != biases.end(); bi++) { total_bias_energy += (*bi)->update(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } } cvm::decrease_depth(); @@ -440,16 +502,31 @@ void colvarmodule::calc() { if (cvm::debug() && biases.size()) cvm::log ("Collecting forces from all biases.\n"); cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->reset_bias_force(); } - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { + for (bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->communicate_forces(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } } + + // Run user force script, if provided, + // potentially adding scripted forces to the colvars + if (use_scripted_forces) { + int res; + res = proxy->run_force_callback(); + if (res == COLVARS_NOT_IMPLEMENTED) { + cvm::error("Colvar forces scripts are not implemented."); + return COLVARS_ERROR; + } + if (res != COLVARS_OK) { + cvm::error("Error running user colvar forces script"); + return COLVARS_ERROR; + } + } + cvm::decrease_depth(); if (cvm::b_analysis) { @@ -457,29 +534,33 @@ void colvarmodule::calc() { if (cvm::debug() && biases.size()) cvm::log ("Perform runtime analyses.\n"); cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->analyse(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } } - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { + for (bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->analyse(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } } cvm::decrease_depth(); } - // sum up the forces for each colvar and integrate any internal - // equation of motion + // sum up the forces for each colvar, including wall forces + // and integrate any internal + // equation of motion (extended system) if (cvm::debug()) cvm::log ("Updating the internal degrees of freedom " "of colvars (if they have any).\n"); cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { total_colvar_energy += (*cvi)->update(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } } cvm::decrease_depth(); proxy->add_energy (total_bias_energy + total_colvar_energy); @@ -489,93 +570,45 @@ void colvarmodule::calc() { if (cvm::debug()) cvm::log ("Communicating forces from the colvars to the atoms.\n"); cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { - if ((*cvi)->tasks[colvar::task_gradients]) + for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { + if ((*cvi)->tasks[colvar::task_gradients]) { (*cvi)->communicate_forces(); + if (cvm::get_error()) { + return COLVARS_ERROR; + } + } } cvm::decrease_depth(); // write restart file, if needed - if (restart_out_freq && restart_out_name.size() && !cvm::b_analysis) { + if (restart_out_freq && restart_out_name.size()) { if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log ("Writing the state file \""+ restart_out_name+"\".\n"); proxy->backup_file (restart_out_name.c_str()); restart_out_os.open (restart_out_name.c_str()); - restart_out_os.setf (std::ios::scientific, std::ios::floatfield); if (!write_restart (restart_out_os)) - cvm::fatal_error ("Error: in writing restart file.\n"); + cvm::error ("Error: in writing restart file.\n"); restart_out_os.close(); } } // write trajectory file, if needed - if (cv_traj_freq) { + if (cv_traj_freq && cv_traj_name.size()) { - if (cvm::debug()) - cvm::log ("Writing trajectory file.\n"); - - // (re)open trajectory file if (!cv_traj_os.good()) { - if (cv_traj_append) { - cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+ - "\".\n"); - cv_traj_os.open (cv_traj_name.c_str(), std::ios::app); - } else { - cvm::log ("Overwriting colvar trajectory file \""+cv_traj_name+ - "\".\n"); - proxy->backup_file (cv_traj_name.c_str()); - cv_traj_os.open (cv_traj_name.c_str(), std::ios::out); - } - cv_traj_os.setf (std::ios::scientific, std::ios::floatfield); + open_traj_file (cv_traj_name); } - // write labels in the traj file every 1000 lines and at first ts - cvm::increase_depth(); + // 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_os << "# " << cvm::wrap_string ("step", cvm::it_width-2) - << " "; - if (cvm::debug()) - cv_traj_os.flush(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { - (*cvi)->write_traj_label (cv_traj_os); - } - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { - (*bi)->write_traj_label (cv_traj_os); - } - cv_traj_os << "\n"; - if (cvm::debug()) - cv_traj_os.flush(); + write_traj_label (cv_traj_os); } - cvm::decrease_depth(); - // write collective variable values to the traj file - cvm::increase_depth(); if ((cvm::step_absolute() % cv_traj_freq) == 0) { - cv_traj_os << std::setw (cvm::it_width) << it - << " "; - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { - (*cvi)->write_traj (cv_traj_os); - } - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { - (*bi)->write_traj (cv_traj_os); - } - cv_traj_os << "\n"; - if (cvm::debug()) - cv_traj_os.flush(); + write_traj (cv_traj_os); } - cvm::decrease_depth(); if (restart_out_freq) { // flush the trajectory file if we are at the restart frequency @@ -587,10 +620,12 @@ void colvarmodule::calc() { } } } // end if (cv_traj_freq) + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvarmodule::analyze() +int colvarmodule::analyze() { if (cvm::debug()) { cvm::log ("colvarmodule::analyze(), step = "+cvm::to_str (it)+".\n"); @@ -617,23 +652,36 @@ void colvarmodule::analyze() cvm::decrease_depth(); } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvarmodule::setup() +int colvarmodule::setup() { // loop over all components of all colvars to reset masses of all groups for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->setup(); } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + } colvarmodule::~colvarmodule() { + reset(); + delete parse; + proxy = NULL; +} + +int colvarmodule::reset() +{ + if (cvm::debug()) + cvm::log ("colvars::reset() called.\n"); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { delete *cvi; + cvi--; } colvars.clear(); @@ -641,19 +689,144 @@ colvarmodule::~colvarmodule() bi != biases.end(); bi++) { delete *bi; + bi--; } biases.clear(); - if (cv_traj_os.good()) { - cv_traj_os.close(); - } + index_groups.clear(); + index_group_names.clear(); - delete parse; - proxy = NULL; + // Do not close file here, as we might not be done with it yet. + cv_traj_os.flush(); + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void colvarmodule::write_output_files() +int colvarmodule::setup_input() +{ + // name of input state file + restart_in_name = proxy->input_prefix().size() ? + std::string (proxy->input_prefix()+".colvars.state") : + std::string ("") ; + + // read the restart configuration, if available + if (restart_in_name.size()) { + // read the restart file + std::ifstream input_is (restart_in_name.c_str()); + if (!input_is.good()) { + cvm::error ("Error: in opening restart file \""+ + std::string (restart_in_name)+"\".\n", + FILE_ERROR); + return COLVARS_ERROR; + } else { + cvm::log ("Restarting from file \""+restart_in_name+"\".\n"); + read_restart (input_is); + if (cvm::get_error() != COLVARS_OK) { + return COLVARS_ERROR; + } + cvm::log (cvm::line_marker); + } + } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + +} + + +int colvarmodule::setup_output() +{ + // output state file (restart) + restart_out_name = proxy->restart_output_prefix().size() ? + std::string (proxy->restart_output_prefix()+".colvars.state") : + std::string (""); + + if (restart_out_name.size()) { + cvm::log ("The restart output state file will be \""+restart_out_name+"\".\n"); + } + + output_prefix = proxy->output_prefix(); + if (output_prefix.size()) { + cvm::log ("The final output state file will be \""+ + (output_prefix.size() ? + std::string (output_prefix+".colvars.state") : + std::string ("colvars.state"))+"\".\n"); + // cvm::log (cvm::line_marker); + } + + cv_traj_name = + (output_prefix.size() ? + std::string (output_prefix+".colvars.traj") : + std::string ("")); + + if (cv_traj_freq && cv_traj_name.size()) { + // open trajectory file + if (cv_traj_append) { + cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+ + "\".\n"); + cv_traj_os.open (cv_traj_name.c_str(), std::ios::app); + } else { + cvm::log ("Writing to colvar trajectory file \""+cv_traj_name+ + "\".\n"); + proxy->backup_file (cv_traj_name.c_str()); + cv_traj_os.open (cv_traj_name.c_str(), std::ios::out); + } + cv_traj_os.setf (std::ios::scientific, std::ios::floatfield); + } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + +} + + +std::istream & colvarmodule::read_restart (std::istream &is) +{ + { + // read global restart information + std::string restart_conf; + if (is >> colvarparse::read_block ("configuration", restart_conf)) { + if (it_restart_from_state_file) { + parse->get_keyval (restart_conf, "step", + it_restart, (size_t) 0, + colvarparse::parse_silent); + it = it_restart; + } + } + is.clear(); + } + + // colvars restart + cvm::increase_depth(); + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); + cvi++) { + if ( !((*cvi)->read_restart (is)) ) { + cvm::error ("Error: in reading restart configuration for collective variable \""+ + (*cvi)->name+"\".\n", + INPUT_ERROR); + } + } + + // biases restart + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + if (!((*bi)->read_restart (is))) + cvm::error ("Error: in reading restart configuration for bias \""+ + (*bi)->name+"\".\n", + INPUT_ERROR); + } + cvm::decrease_depth(); + + return is; +} + + +int colvarmodule::backup_file (char const *filename) +{ + return proxy->backup_file (filename); +} + + +int colvarmodule::write_output_files() { // if this is a simulation run (i.e. not a postprocessing), output data // must be written to be able to restart the simulation @@ -678,11 +851,12 @@ void colvarmodule::write_output_files() // do not close to avoid problems with multiple NAMD runs cv_traj_os.flush(); + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -bool colvarmodule::read_traj (char const *traj_filename, +int colvarmodule::read_traj (char const *traj_filename, size_t traj_read_begin, size_t traj_read_end) { @@ -725,19 +899,21 @@ bool colvarmodule::read_traj (char const *traj_filename, if ( (traj_read_end > traj_read_begin) && (it > traj_read_end) ) { std::cerr << "\n"; - cvm::log ("Reached the end of the trajectory, " - "read_end = "+cvm::to_str (traj_read_end)+"\n"); - return false; + cvm::error ("Reached the end of the trajectory, " + "read_end = "+cvm::to_str (traj_read_end)+"\n", + FILE_ERROR); + return COLVARS_ERROR; } for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->read_traj (is)) { - cvm::log ("Error: in reading colvar \""+(*cvi)->name+ + cvm::error ("Error: in reading colvar \""+(*cvi)->name+ "\" from trajectory file \""+ - std::string (traj_filename)+"\".\n"); - return false; + std::string (traj_filename)+"\".\n", + FILE_ERROR); + return COLVARS_ERROR; } } @@ -745,13 +921,13 @@ bool colvarmodule::read_traj (char const *traj_filename, } } } - - return true; + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::ostream & colvarmodule::write_restart (std::ostream &os) { + os.setf (std::ios::scientific, std::ios::floatfield); os << "configuration {\n" << " step " << std::setw (it_width) << it << "\n" @@ -775,6 +951,87 @@ std::ostream & colvarmodule::write_restart (std::ostream &os) return os; } +int colvarmodule::open_traj_file (std::string const &file_name) +{ + // (re)open trajectory file + if (cv_traj_append) { + cvm::log ("Appending to colvar trajectory file \""+file_name+ + "\".\n"); + cv_traj_os.open (file_name.c_str(), std::ios::app); + } else { + cvm::log ("Overwriting colvar trajectory file \""+file_name+ + "\".\n"); + proxy->backup_file (file_name.c_str()); + cv_traj_os.open (file_name.c_str(), std::ios::out); + } + if (!cv_traj_os.good()) { + cvm::error ("Error: cannot write to file \""+file_name+"\".\n", + FILE_ERROR); + } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + +int colvarmodule::close_traj_file() +{ + if (cv_traj_os.good()) { + cv_traj_os.close(); + } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + +std::ostream & colvarmodule::write_traj_label (std::ostream &os) +{ + if (!os.good()) { + cvm::error("Cannot write to trajectory file."); + return os; + } + os.setf (std::ios::scientific, std::ios::floatfield); + + os << "# " << cvm::wrap_string ("step", cvm::it_width-2) + << " "; + + cvm::increase_depth(); + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); + cvi++) { + (*cvi)->write_traj_label (os); + } + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + (*bi)->write_traj_label (os); + } + os << "\n"; + if (cvm::debug()) + os.flush(); + cvm::decrease_depth(); + return os; +} + +std::ostream & colvarmodule::write_traj (std::ostream &os) +{ + os.setf (std::ios::scientific, std::ios::floatfield); + + os << std::setw (cvm::it_width) << it + << " "; + + cvm::increase_depth(); + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); + cvi++) { + (*cvi)->write_traj (os); + } + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + (*bi)->write_traj (os); + } + os << "\n"; + if (cvm::debug()) + os.flush(); + cvm::decrease_depth(); + return os; +} void cvm::log (std::string const &message) @@ -795,8 +1052,15 @@ void cvm::decrease_depth() if (depth) depth--; } +void cvm::error (std::string const &message, int code) +{ + set_error_bits(code); + proxy->error (message); +} + void cvm::fatal_error (std::string const &message) { + set_error_bits(FATAL_ERROR); proxy->fatal_error (message); } @@ -806,25 +1070,35 @@ void cvm::exit (std::string const &message) } -void cvm::read_index_file (char const *filename) +int cvm::read_index_file (char const *filename) { std::ifstream is (filename, std::ios::binary); if (!is.good()) - fatal_error ("Error: in opening index file \""+ - std::string (filename)+"\".\n"); - // std::list::iterator names_i = cvm::index_group_names.begin(); - // std::list >::iterator lists_i = cvm::index_groups.begin(); + cvm::error ("Error: in opening index file \""+ + std::string (filename)+"\".\n", + FILE_ERROR); + while (is.good()) { char open, close; std::string group_name; if ( (is >> open) && (open == '[') && (is >> group_name) && (is >> close) && (close == ']') ) { + for (std::list::iterator names_i = index_group_names.begin(); + names_i != index_group_names.end(); + names_i++) { + if (*names_i == group_name) { + cvm::error ("Error: the group name \""+group_name+ + "\" appears in multiple index files.\n", + FILE_ERROR); + } + } cvm::index_group_names.push_back (group_name); cvm::index_groups.push_back (std::vector ()); } else { - cvm::fatal_error ("Error: in parsing index file \""+ - std::string (filename)+"\".\n"); + cvm::error ("Error: in parsing index file \""+ + std::string (filename)+"\".\n", + INPUT_ERROR); } int atom_number = 1; @@ -847,23 +1121,23 @@ void cvm::read_index_file (char const *filename) cvm::log ("The following index groups were read from the index file \""+ std::string (filename)+"\":\n"); - std::list::iterator names_i = cvm::index_group_names.begin(); - std::list >::iterator lists_i = cvm::index_groups.begin(); - for ( ; names_i != cvm::index_group_names.end() ; names_i++, lists_i++) { + std::list::iterator names_i = index_group_names.begin(); + std::list >::iterator lists_i = index_groups.begin(); + for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) { cvm::log (" "+(*names_i)+" ("+cvm::to_str (lists_i->size())+" atoms).\n"); } - + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } -void cvm::load_atoms (char const *file_name, +int cvm::load_atoms (char const *file_name, std::vector &atoms, std::string const &pdb_field, double const pdb_field_value) { - proxy->load_atoms (file_name, atoms, pdb_field, pdb_field_value); + return proxy->load_atoms (file_name, atoms, pdb_field, pdb_field_value); } -void cvm::load_coords (char const *file_name, +int cvm::load_coords (char const *file_name, std::vector &pos, const std::vector &indices, std::string const &pdb_field, @@ -876,31 +1150,28 @@ void cvm::load_coords (char const *file_name, std::string const ext (strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name); if (colvarparse::to_lower_cppstr (ext) == std::string (".xyz")) { if ( pdb_field.size() > 0 ) { - cvm::fatal_error ("Error: PDB column may not be specified for XYZ coordinate file.\n"); + cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR); + return COLVARS_ERROR; } - cvm::load_coords_xyz (file_name, pos, indices); + return cvm::load_coords_xyz (file_name, pos, indices); } else { - proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value); + return proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value); } } -void cvm::backup_file (char const *filename) -{ - proxy->backup_file (filename); -} -void cvm::load_coords_xyz (char const *filename, +int cvm::load_coords_xyz (char const *filename, std::vector &pos, const std::vector &indices) { std::ifstream xyz_is (filename); - int natoms; + unsigned int natoms; char symbol[256]; std::string line; if ( ! (xyz_is >> natoms) ) { - cvm::fatal_error ("Error: cannot parse XYZ file " - + std::string (filename) + ".\n"); + cvm::error ("Error: cannot parse XYZ file " + + std::string (filename) + ".\n", INPUT_ERROR); } // skip comment line std::getline (xyz_is, line); @@ -926,6 +1197,7 @@ void cvm::load_coords_xyz (char const *filename, xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; } } + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } // static pointers @@ -940,6 +1212,7 @@ colvarproxy *colvarmodule::proxy = NULL; // static runtime data cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03; +int colvarmodule::errorCode = 0; size_t colvarmodule::it = 0; size_t colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; @@ -953,7 +1226,6 @@ std::list > colvarmodule::index_groups; // file name prefixes std::string colvarmodule::output_prefix = ""; -std::string colvarmodule::input_prefix = ""; std::string colvarmodule::restart_in_name = ""; @@ -965,612 +1237,3 @@ size_t const colvarmodule::en_prec = 14; size_t const colvarmodule::en_width = 21; std::string const colvarmodule::line_marker = "----------------------------------------------------------------------\n"; - - - - - - -std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v) -{ - std::streamsize const w = os.width(); - std::streamsize const p = os.precision(); - - os.width (2); - os << "( "; - os.width (w); os.precision (p); - os << v.x << " , "; - os.width (w); os.precision (p); - os << v.y << " , "; - os.width (w); os.precision (p); - os << v.z << " )"; - return os; -} - - -std::istream & operator >> (std::istream &is, colvarmodule::rvector &v) -{ - size_t const start_pos = is.tellg(); - char sep; - if ( !(is >> sep) || !(sep == '(') || - !(is >> v.x) || !(is >> sep) || !(sep == ',') || - !(is >> v.y) || !(is >> sep) || !(sep == ',') || - !(is >> v.z) || !(is >> sep) || !(sep == ')') ) { - is.clear(); - is.seekg (start_pos, std::ios::beg); - is.setstate (std::ios::failbit); - return is; - } - return is; -} - - - -std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q) -{ - std::streamsize const w = os.width(); - std::streamsize const p = os.precision(); - - os.width (2); - os << "( "; - os.width (w); os.precision (p); - os << q.q0 << " , "; - os.width (w); os.precision (p); - os << q.q1 << " , "; - os.width (w); os.precision (p); - os << q.q2 << " , "; - os.width (w); os.precision (p); - os << q.q3 << " )"; - return os; -} - - -std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q) -{ - size_t 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 - - 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; - } - } - - return is; -} - - -cvm::quaternion -cvm::quaternion::position_derivative_inner (cvm::rvector const &pos, - cvm::rvector const &vec) const -{ - cvm::quaternion result (0.0, 0.0, 0.0, 0.0); - - - result.q0 = 2.0 * pos.x * q0 * vec.x - +2.0 * pos.y * q0 * vec.y - +2.0 * pos.z * q0 * vec.z - - -2.0 * pos.y * q3 * vec.x - +2.0 * pos.z * q2 * vec.x - - +2.0 * pos.x * q3 * vec.y - -2.0 * pos.z * q1 * vec.y - - -2.0 * pos.x * q2 * vec.z - +2.0 * pos.y * q1 * vec.z; - - - result.q1 = +2.0 * pos.x * q1 * vec.x - -2.0 * pos.y * q1 * vec.y - -2.0 * pos.z * q1 * vec.z - - +2.0 * pos.y * q2 * vec.x - +2.0 * pos.z * q3 * vec.x - - +2.0 * pos.x * q2 * vec.y - -2.0 * pos.z * q0 * vec.y - - +2.0 * pos.x * q3 * vec.z - +2.0 * pos.y * q0 * vec.z; - - - result.q2 = -2.0 * pos.x * q2 * vec.x - +2.0 * pos.y * q2 * vec.y - -2.0 * pos.z * q2 * vec.z - - +2.0 * pos.y * q1 * vec.x - +2.0 * pos.z * q0 * vec.x - - +2.0 * pos.x * q1 * vec.y - +2.0 * pos.z * q3 * vec.y - - -2.0 * pos.x * q0 * vec.z - +2.0 * pos.y * q3 * vec.z; - - - result.q3 = -2.0 * pos.x * q3 * vec.x - -2.0 * pos.y * q3 * vec.y - +2.0 * pos.z * q3 * vec.z - - -2.0 * pos.y * q0 * vec.x - +2.0 * pos.z * q1 * vec.x - - +2.0 * pos.x * q0 * vec.y - +2.0 * pos.z * q2 * vec.y - - +2.0 * pos.x * q1 * vec.z - +2.0 * pos.y * q2 * vec.z; - - return result; -} - - - - - - -// Calculate the optimal rotation between two groups, and implement it -// as a quaternion. Uses the method documented in: Coutsias EA, -// Seok C, Dill KA. Using quaternions to calculate RMSD. J Comput -// Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254 - -void colvarmodule::rotation::build_matrix (std::vector const &pos1, - std::vector const &pos2, - matrix2d &S) -{ - cvm::rmatrix C; - - // build the correlation matrix - C.reset(); - for (size_t 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; - } - - // 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][1] = S[1][0]; - S[2][0] = - C.xz() + C.zx() ; - S[0][2] = S[2][0]; - 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][2] = S[2][1]; - 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][3] = S[3][2]; - S[3][3] = - C.xx() - C.yy() + C.zz(); - - // if (cvm::debug()) { - // for (size_t i = 0; i < 4; i++) { - // std::string line (""); - // for (size_t j = 0; j < 4; j++) { - // line += std::string (" S["+cvm::to_str (i)+ - // "]["+cvm::to_str (j)+"] ="+cvm::to_str (S[i][j])); - // } - // cvm::log (line+"\n"); - // } - // } -} - - -void colvarmodule::rotation::diagonalize_matrix (matrix2d &S, - cvm::real S_eigval[4], - matrix2d &S_eigvec) -{ - // diagonalize - int jac_nrot = 0; - jacobi (S, 4, S_eigval, S_eigvec, &jac_nrot); - eigsrt (S_eigval, S_eigvec, 4); - // jacobi saves eigenvectors by columns - transpose (S_eigvec, 4); - - // normalize eigenvectors - for (size_t ie = 0; ie < 4; ie++) { - cvm::real norm2 = 0.0; - for (size_t i = 0; i < 4; i++) norm2 += std::pow (S_eigvec[ie][i], int (2)); - cvm::real const norm = std::sqrt (norm2); - for (size_t i = 0; i < 4; i++) S_eigvec[ie][i] /= norm; - } -} - - -// Calculate the rotation, plus its derivatives - -void colvarmodule::rotation::calc_optimal_rotation -(std::vector const &pos1, - std::vector const &pos2) -{ - matrix2d S; - matrix2d S_backup; - cvm::real S_eigval[4]; - matrix2d S_eigvec; - -// if (cvm::debug()) { -// cvm::atom_pos cog1 (0.0, 0.0, 0.0); -// for (size_t i = 0; i < pos1.size(); i++) { -// cog1 += pos1[i]; -// } -// cog1 /= cvm::real (pos1.size()); -// cvm::atom_pos cog2 (0.0, 0.0, 0.0); -// for (size_t i = 0; i < pos2.size(); i++) { -// cog2 += pos2[i]; -// } -// cog2 /= cvm::real (pos1.size()); -// cvm::log ("calc_optimal_rotation: centers of geometry are: "+ -// cvm::to_str (cog1, cvm::cv_width, cvm::cv_prec)+ -// " and "+cvm::to_str (cog2, cvm::cv_width, cvm::cv_prec)+".\n"); -// } - - build_matrix (pos1, pos2, S); - S_backup = S; - - if (cvm::debug()) { - if (b_debug_gradients) { - cvm::log ("S = "+cvm::to_str (cvm::to_str (S_backup), cvm::cv_width, cvm::cv_prec)+"\n"); - } - } - - diagonalize_matrix (S, S_eigval, S_eigvec); - - // eigenvalues and eigenvectors - cvm::real const &L0 = S_eigval[0]; - cvm::real const &L1 = S_eigval[1]; - cvm::real const &L2 = S_eigval[2]; - cvm::real const &L3 = S_eigval[3]; - cvm::real const *Q0 = S_eigvec[0]; - cvm::real const *Q1 = S_eigvec[1]; - cvm::real const *Q2 = S_eigvec[2]; - cvm::real const *Q3 = S_eigvec[3]; - - lambda = L0; - q = cvm::quaternion (Q0); - - if (q_old.norm2() > 0.0) { - q.match (q_old); - if (q_old.inner (q) < (1.0 - crossing_threshold)) { - cvm::log ("Warning: one molecular orientation has changed by more than "+ - cvm::to_str (crossing_threshold)+": discontinuous rotation ?\n"); - } - } - q_old = q; - - if (cvm::debug()) { - if (b_debug_gradients) { - cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+ - ", Q0 = "+cvm::to_str (cvm::quaternion (Q0), cvm::cv_width, cvm::cv_prec)+ - ", Q0*Q0 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q0)), cvm::cv_width, cvm::cv_prec)+ - "\n"); - cvm::log ("L1 = "+cvm::to_str (L1, cvm::cv_width, cvm::cv_prec)+ - ", Q1 = "+cvm::to_str (cvm::quaternion (Q1), cvm::cv_width, cvm::cv_prec)+ - ", Q0*Q1 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q1)), cvm::cv_width, cvm::cv_prec)+ - "\n"); - cvm::log ("L2 = "+cvm::to_str (L2, cvm::cv_width, cvm::cv_prec)+ - ", Q2 = "+cvm::to_str (cvm::quaternion (Q2), cvm::cv_width, cvm::cv_prec)+ - ", Q0*Q2 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q2)), cvm::cv_width, cvm::cv_prec)+ - "\n"); - cvm::log ("L3 = "+cvm::to_str (L3, cvm::cv_width, cvm::cv_prec)+ - ", Q3 = "+cvm::to_str (cvm::quaternion (Q3), cvm::cv_width, cvm::cv_prec)+ - ", Q0*Q3 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q3)), cvm::cv_width, cvm::cv_prec)+ - "\n"); - } - } - - // calculate derivatives of L0 and Q0 with respect to each atom in - // either group; note: if dS_1 is a null vector, nothing will be - // calculated - for (size_t ia = 0; ia < dS_1.size(); ia++) { - - cvm::real const &a2x = pos2[ia].x; - cvm::real const &a2y = pos2[ia].y; - cvm::real const &a2z = pos2[ia].z; - - matrix2d &ds_1 = dS_1[ia]; - - // derivative of the S matrix - ds_1.reset(); - ds_1[0][0] = cvm::rvector ( a2x, a2y, a2z); - ds_1[1][0] = cvm::rvector ( 0.0, a2z, -a2y); - ds_1[0][1] = ds_1[1][0]; - ds_1[2][0] = cvm::rvector (-a2z, 0.0, a2x); - ds_1[0][2] = ds_1[2][0]; - ds_1[3][0] = cvm::rvector ( a2y, -a2x, 0.0); - ds_1[0][3] = ds_1[3][0]; - ds_1[1][1] = cvm::rvector ( a2x, -a2y, -a2z); - ds_1[2][1] = cvm::rvector ( a2y, a2x, 0.0); - ds_1[1][2] = ds_1[2][1]; - ds_1[3][1] = cvm::rvector ( a2z, 0.0, a2x); - ds_1[1][3] = ds_1[3][1]; - ds_1[2][2] = cvm::rvector (-a2x, a2y, -a2z); - ds_1[3][2] = cvm::rvector ( 0.0, a2z, a2y); - ds_1[2][3] = ds_1[3][2]; - ds_1[3][3] = cvm::rvector (-a2x, -a2y, a2z); - - cvm::rvector &dl0_1 = dL0_1[ia]; - vector1d &dq0_1 = dQ0_1[ia]; - - // matrix multiplications; derivatives of L_0 and Q_0 are - // calculated using Hellmann-Feynman theorem (i.e. exploiting the - // fact that the eigenvectors Q_i form an orthonormal basis) - - dl0_1.reset(); - for (size_t i = 0; i < 4; i++) { - for (size_t j = 0; j < 4; j++) { - dl0_1 += Q0[i] * ds_1[i][j] * Q0[j]; - } - } - - dq0_1.reset(); - for (size_t p = 0; p < 4; p++) { - for (size_t i = 0; i < 4; i++) { - for (size_t j = 0; j < 4; j++) { - dq0_1[p] += - (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] + - (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] + - (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p]; - } - } - } - } - - // do the same for the second group - for (size_t ia = 0; ia < dS_2.size(); ia++) { - - cvm::real const &a1x = pos1[ia].x; - cvm::real const &a1y = pos1[ia].y; - cvm::real const &a1z = pos1[ia].z; - - matrix2d &ds_2 = dS_2[ia]; - - ds_2.reset(); - ds_2[0][0] = cvm::rvector ( a1x, a1y, a1z); - ds_2[1][0] = cvm::rvector ( 0.0, -a1z, a1y); - ds_2[0][1] = ds_2[1][0]; - ds_2[2][0] = cvm::rvector ( a1z, 0.0, -a1x); - ds_2[0][2] = ds_2[2][0]; - ds_2[3][0] = cvm::rvector (-a1y, a1x, 0.0); - ds_2[0][3] = ds_2[3][0]; - ds_2[1][1] = cvm::rvector ( a1x, -a1y, -a1z); - ds_2[2][1] = cvm::rvector ( a1y, a1x, 0.0); - ds_2[1][2] = ds_2[2][1]; - ds_2[3][1] = cvm::rvector ( a1z, 0.0, a1x); - ds_2[1][3] = ds_2[3][1]; - ds_2[2][2] = cvm::rvector (-a1x, a1y, -a1z); - ds_2[3][2] = cvm::rvector ( 0.0, a1z, a1y); - ds_2[2][3] = ds_2[3][2]; - ds_2[3][3] = cvm::rvector (-a1x, -a1y, a1z); - - cvm::rvector &dl0_2 = dL0_2[ia]; - vector1d &dq0_2 = dQ0_2[ia]; - - dl0_2.reset(); - for (size_t i = 0; i < 4; i++) { - for (size_t j = 0; j < 4; j++) { - dl0_2 += Q0[i] * ds_2[i][j] * Q0[j]; - } - } - - dq0_2.reset(); - for (size_t p = 0; p < 4; p++) { - for (size_t i = 0; i < 4; i++) { - for (size_t j = 0; j < 4; j++) { - dq0_2[p] += - (Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] + - (Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] + - (Q3[i] * ds_2[i][j] * Q0[j]) / (L0-L3) * Q3[p]; - } - } - } - - if (cvm::debug()) { - - if (b_debug_gradients) { - - matrix2d S_new; - cvm::real S_new_eigval[4]; - matrix2d S_new_eigvec; - - // make an infitesimal move along each cartesian coordinate of - // this atom, and solve again the eigenvector problem - for (size_t comp = 0; comp < 3; comp++) { - - S_new = S_backup; - // diagonalize the new overlap matrix - for (size_t i = 0; i < 4; i++) { - for (size_t j = 0; j < 4; j++) { - S_new[i][j] += - colvarmodule::debug_gradients_step_size * ds_2[i][j][comp]; - } - } - -// cvm::log ("S_new = "+cvm::to_str (cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n"); - - diagonalize_matrix (S_new, S_new_eigval, S_new_eigvec); - - cvm::real const &L0_new = S_new_eigval[0]; - cvm::real const *Q0_new = S_new_eigvec[0]; - - cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size; - cvm::quaternion const q0 (Q0); - cvm::quaternion const DQ0 (dq0_2[0][comp] * colvarmodule::debug_gradients_step_size, - dq0_2[1][comp] * colvarmodule::debug_gradients_step_size, - dq0_2[2][comp] * colvarmodule::debug_gradients_step_size, - dq0_2[3][comp] * colvarmodule::debug_gradients_step_size); - - cvm::log ( "|(l_0+dl_0) - l_0^new|/l_0 = "+ - cvm::to_str (std::fabs (L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+ - ", |(q_0+dq_0) - q_0^new| = "+ - cvm::to_str ((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+ - "\n"); - } - } - } - } -} - - - -// Numerical Recipes routine for diagonalization - -#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau); \ - a[k][l]=h+s*(g-h*tau); -void jacobi(cvm::real **a, int n, cvm::real d[], cvm::real **v, int *nrot) -{ - int j,iq,ip,i; - cvm::real tresh,theta,tau,t,sm,s,h,g,c; - - std::vector b (n, 0.0); - std::vector z (n, 0.0); - - for (ip=0;ip 4 && (cvm::real)(std::fabs(d[ip])+g) == (cvm::real)std::fabs(d[ip]) - && (cvm::real)(std::fabs(d[iq])+g) == (cvm::real)std::fabs(d[iq])) - a[ip][iq]=0.0; - else if (std::fabs(a[ip][iq]) > tresh) { - h=d[iq]-d[ip]; - if ((cvm::real)(std::fabs(h)+g) == (cvm::real)std::fabs(h)) - t=(a[ip][iq])/h; - else { - theta=0.5*h/(a[ip][iq]); - t=1.0/(std::fabs(theta)+std::sqrt(1.0+theta*theta)); - if (theta < 0.0) t = -t; - } - c=1.0/std::sqrt(1+t*t); - s=t*c; - tau=s/(1.0+c); - h=t*a[ip][iq]; - z[ip] -= h; - z[iq] += h; - d[ip] -= h; - d[iq] += h; - a[ip][iq]=0.0; - for (j=0;j<=ip-1;j++) { - ROTATE(a,j,ip,j,iq) - } - for (j=ip+1;j<=iq-1;j++) { - ROTATE(a,ip,j,j,iq) - } - for (j=iq+1;j= p) p=d[k=j]; - if (k != i) { - d[k]=d[i]; - d[i]=p; - for (j=0;j #include @@ -34,6 +49,7 @@ class colvarparse; class colvar; class colvarbias; class colvarproxy; +class colvarscript; /// \brief Collective variables module (main class) @@ -55,6 +71,8 @@ private: public: friend class colvarproxy; + // TODO colvarscript should be unaware of colvarmodule's internals + friend class colvarscript; /// Defining an abstract real number allows to switch precision typedef double real; @@ -85,6 +103,22 @@ public: typedef std::vector::iterator atom_iter; typedef std::vector::const_iterator atom_const_iter; + /// Module-wide error state + /// see constants at the top of this file + static int errorCode; + static inline void set_error_bits(int code) + { + errorCode |= code; + } + static inline int get_error() + { + return errorCode; + } + static inline void clear_error() + { + errorCode = 0; + } + /// Current step number static size_t it; /// Starting step number for this run @@ -115,9 +149,6 @@ public: /// Prefix for all output files for this run static std::string output_prefix; - /// Prefix for files from a previous run (including restart/output) - static std::string input_prefix; - /// input restart file name (determined from input_prefix) static std::string restart_in_name; @@ -158,26 +189,79 @@ public: /// \brief Constructor \param config_name Configuration file name /// \param restart_name (optional) Restart file name - colvarmodule (char const *config_name, - colvarproxy *proxy_in); + colvarmodule (colvarproxy *proxy); /// Destructor ~colvarmodule(); - /// Initialize collective variables - void init_colvars (std::string const &conf); + /// Actual function called by the destructor + int reset(); - /// Initialize collective variable biases - void init_biases (std::string const &conf); + /// Open a config file, load its contents, and pass it to config_string() + int config_file (char const *config_file_name); - /// Re-initialize data at the beginning of a run. For use with - /// MD codes that can change system parameters like atom masses - /// between run commands. - void setup(); + /// \brief Parse a config string assuming it is a complete configuration + /// (i.e. calling all parse functions) + int config_string (std::string const &conf); + + /// \brief Parse a "clean" config string (no comments) + int config (std::string &conf); + + + // Parse functions (setup internal data based on a string) + + /// Parse the few module's global parameters + int parse_global_params (std::string const &conf); + + /// Parse and initialize collective variables + int parse_colvars (std::string const &conf); + + /// Parse and initialize collective variable biases + int parse_biases (std::string const &conf); + + + // "Setup" functions (change internal data based on related data + // 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) + /// Also calls setup() member functions of colvars and biases + int setup(); + + /// (Re)initialize and (re)read the input state file calling read_restart() + int setup_input(); + + /// (Re)initialize the output trajectory and state file (does not write it yet) + int setup_output(); + + /// Read the input restart file + std::istream & read_restart (std::istream &is); + /// Write the output restart file + std::ostream & write_restart (std::ostream &os); + + /// Open a trajectory file if requested (and leave it open) + int open_traj_file (std::string const &file_name); + /// Close it + int close_traj_file(); + /// Write in the trajectory file + std::ostream & write_traj (std::ostream &os); + /// Write explanatory labels in the trajectory file + std::ostream & write_traj_label (std::ostream &os); + + /// Write all FINAL output files + int write_output_files(); + /// Backup a file before writing it + static int backup_file (char const *filename); + + /// Look up a bias by name; returns NULL if not found + static colvarbias * bias_by_name(std::string const &name); + + /// Look up a colvar by name; returns NULL if not found + static colvar * colvar_by_name(std::string const &name); /// Load new configuration for the given bias - /// currently works for harmonic (force constant and/or centers) - void change_configuration (std::string const &bias_name, std::string const &conf); + int change_configuration (std::string const &bias_name, std::string const &conf); /// Read a colvar value std::string read_colvar(std::string const &name); @@ -187,27 +271,16 @@ public: real energy_difference (std::string const &bias_name, std::string const &conf); /// Calculate collective variables and biases - void calc(); - /// Read the input restart file - std::istream & read_restart (std::istream &is); - /// Write the output restart file - std::ostream & write_restart (std::ostream &os); - /// Write all output files (called by the proxy) - void write_output_files(); - /// \brief Call colvarproxy::backup_file() - static void backup_file (char const *filename); + int calc(); /// Perform analysis - void analyze(); + int analyze(); /// \brief Read a collective variable trajectory (post-processing /// only, not called at runtime) - bool read_traj (char const *traj_filename, + int read_traj (char const *traj_filename, size_t traj_read_begin, size_t traj_read_end); - /// Get the pointer of a colvar from its name (returns NULL if not found) - static colvar * colvar_p (std::string const &name); - /// Quick conversion of an object to a string template static std::string to_str (T const &x, size_t const &width = 0, @@ -267,6 +340,9 @@ public: /// Print a message to the main log and exit with error code static void fatal_error (std::string const &message); + /// Print a message to the main log and set global error code + static void error (std::string const &message, int code = GENERAL_ERROR); + /// Print a message to the main log and exit normally static void exit (std::string const &message); @@ -307,21 +383,21 @@ public: static std::list > index_groups; /// \brief Read a Gromacs .ndx file - static void read_index_file (char const *filename); + static int read_index_file (char const *filename); /// \brief Create atoms from a file \param filename name of the file /// (usually a PDB) \param atoms array of the atoms to be allocated /// \param pdb_field (optiona) if "filename" is a PDB file, use this /// field to determine which are the atoms to be set - static void load_atoms (char const *filename, + static int load_atoms (char const *filename, std::vector &atoms, std::string const &pdb_field, double const pdb_field_value = 0.0); /// \brief Load the coordinates for a group of atoms from a file /// (PDB or XYZ) - static void load_coords (char const *filename, + static int load_coords (char const *filename, std::vector &pos, const std::vector &indices, std::string const &pdb_field, @@ -329,7 +405,7 @@ public: /// \brief Load the coordinates for a group of atoms from an /// XYZ file - static void load_coords_xyz (char const *filename, + static int load_coords_xyz (char const *filename, std::vector &pos, const std::vector &indices); @@ -366,15 +442,18 @@ protected: /// Output restart file std::ofstream restart_out_os; + /// \brief Counter for the current depth in the object hierarchy (useg e.g. in output + static size_t depth; + + /// Use scripted colvars forces? + bool use_scripted_forces; + +public: /// \brief Pointer to the proxy object, used to retrieve atomic data /// from the hosting program; it is static in order to be accessible /// from static functions in the colvarmodule class static colvarproxy *proxy; - /// \brief Counter for the current depth in the object hierarchy (useg e.g. in outpu - static size_t depth; - -public: /// Increase the depth (number of indentations in the output) static void increase_depth(); @@ -489,9 +568,3 @@ inline cvm::real cvm::rand_gaussian (void) } #endif - - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 5ccaf98294..486a8f01a4 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -1,14 +1,13 @@ +/// -*- c++ -*- + #include #include - - #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" - // space & tab std::string const colvarparse::white_space = " \t"; @@ -51,8 +50,8 @@ size_t colvarparse::dummy_pos = 0; if (is >> x) \ value = x; \ else \ - cvm::fatal_error ("Error: in parsing \""+ \ - std::string (key)+"\".\n"); \ + cvm::error ("Error: in parsing \""+ \ + std::string (key)+"\".\n", INPUT_ERROR); \ if (parse_mode != parse_silent) { \ cvm::log ("# "+std::string (key)+" = "+ \ cvm::to_str (value)+"\n"); \ @@ -60,8 +59,8 @@ size_t colvarparse::dummy_pos = 0; } else { \ \ if (b_found_any) \ - cvm::fatal_error ("Error: improper or missing value " \ - "for \""+std::string (key)+"\".\n"); \ + cvm::error ("Error: improper or missing value " \ + "for \""+std::string (key)+"\".\n", INPUT_ERROR); \ value = def_value; \ if (parse_mode != parse_silent) { \ cvm::log ("# "+std::string (key)+" = \""+ \ @@ -112,9 +111,9 @@ size_t colvarparse::dummy_pos = 0; cvm::fatal_error ("Error: in parsing \""+ \ std::string (key)+"\".\n"); \ if (data_count > 1) \ - cvm::fatal_error ("Error: multiple values " \ - "are not allowed for keyword \""+ \ - std::string (key)+"\".\n"); \ + cvm::error ("Error: multiple values " \ + "are not allowed for keyword \""+ \ + std::string (key)+"\".\n", INPUT_ERROR); \ if (parse_mode != parse_silent) { \ cvm::log ("# "+std::string (key)+" = "+ \ cvm::to_str (value)+"\n"); \ @@ -122,8 +121,8 @@ size_t colvarparse::dummy_pos = 0; } else { \ \ if (b_found_any) \ - cvm::fatal_error ("Error: improper or missing value " \ - "for \""+std::string (key)+"\".\n"); \ + cvm::error ("Error: improper or missing value " \ + "for \""+std::string (key)+"\".\n", INPUT_ERROR); \ value = def_value; \ if (parse_mode != parse_silent) { \ cvm::log ("# "+std::string (key)+" = "+ \ @@ -189,8 +188,8 @@ size_t colvarparse::dummy_pos = 0; if (is >> x) \ values[i] = x; \ else \ - cvm::fatal_error ("Error: in parsing \""+ \ - std::string (key)+"\".\n"); \ + cvm::error ("Error: in parsing \""+ \ + std::string (key)+"\".\n", INPUT_ERROR); \ } \ } \ \ @@ -202,8 +201,8 @@ size_t colvarparse::dummy_pos = 0; } else { \ \ if (b_found_any) \ - cvm::fatal_error ("Error: improper or missing values for \""+ \ - std::string (key)+"\".\n"); \ + cvm::error ("Error: improper or missing values for \""+ \ + std::string (key)+"\".\n", INPUT_ERROR); \ \ for (size_t i = 0; i < values.size(); i++) \ values[i] = def_values[ (i > def_values.size()) ? 0 : i ]; \ @@ -267,7 +266,6 @@ bool colvarparse::get_keyval (std::string const &conf, std::string (key)+"\".\n"); if (data.size()) { - std::istringstream is (data); if ( (data == std::string ("on")) || (data == std::string ("yes")) || (data == std::string ("true")) ) { @@ -349,7 +347,7 @@ void colvarparse::strip_values (std::string &conf) } -void colvarparse::check_keywords (std::string &conf, char const *key) +int colvarparse::check_keywords (std::string &conf, char const *key) { if (cvm::debug()) cvm::log ("Configuration string for \""+std::string (key)+ @@ -383,10 +381,17 @@ void colvarparse::check_keywords (std::string &conf, char const *key) break; } } - if (!found_keyword) - cvm::fatal_error ("Error: keyword \""+uk+"\" is not supported, " + if (!found_keyword) { + cvm::log ("Error: keyword \""+uk+"\" is not supported, " "or not recognized in this context.\n"); + cvm::set_error_bits(INPUT_ERROR); + return COLVARS_ERROR; + } } + allowed_keywords.clear(); + data_begin_pos.clear(); + data_end_pos.clear(); + return COLVARS_OK; } @@ -635,9 +640,3 @@ bool colvarparse::brace_check (std::string const &conf, else return true; } - - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index d3b3968a79..8ea4ca853c 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #ifndef COLVARPARSE_H #define COLVARPARSE_H @@ -125,7 +127,7 @@ public: /// \brief Check that all the keywords within "conf" are in the list /// of allowed keywords; this will invoke strip_values() first and /// then loop over all words - void check_keywords (std::string &conf, char const *key); + int check_keywords (std::string &conf, char const *key); /// \brief Return a lowercased copy of the string @@ -194,8 +196,3 @@ public: #endif - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index c1416f13b7..89c229dde3 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -1,24 +1,46 @@ +/// -*- c++ -*- + #ifndef COLVARPROXY_H #define COLVARPROXY_H +#ifndef COLVARPROXY_VERSION +#define COLVARPROXY_VERSION "2014-09-19" +#endif + #include "colvarmodule.h" +#include "colvarvalue.h" -/// \brief Interface class between the collective variables module and -/// the simulation program +// return values for the frame() routine +#define COLVARS_NO_SUCH_FRAME -1 +#define COLVARS_NOT_IMPLEMENTED -2 + +// forward declarations +class colvarscript; + +/// \brief Interface between the collective variables module and +/// the simulation or analysis program. +/// This is the base class: each program is supported by a derived class. +/// Only pure virtual functions ("= 0") must be reimplemented in a new interface. class colvarproxy { public: - /// Pointer to the instance of colvarmodule + /// Pointer to the main object colvarmodule *colvars; + /// Default constructor + inline colvarproxy() : script (NULL) {} + /// Default destructor virtual inline ~colvarproxy() {} + /// (Re)initialize member data after construction + virtual void setup() {} + // **************** SYSTEM-WIDE PHYSICAL QUANTITIES **************** @@ -38,19 +60,38 @@ public: /// \brief Pseudo-random number with Gaussian distribution virtual cvm::real rand_gaussian (void) = 0; + /// \brief Get the current frame number + virtual int frame() { return COLVARS_NOT_IMPLEMENTED; } + + /// \brief Set the current frame number + // return 0 on success, -1 on failure + virtual int frame (int) { return COLVARS_NOT_IMPLEMENTED; } + // **************** SIMULATION PARAMETERS **************** + /// \brief Prefix to be used for input files (restarts, not /// configuration) - virtual std::string input_prefix() = 0; + std::string input_prefix_str, output_prefix_str, restart_output_prefix_str; + + inline std::string input_prefix() + { + return input_prefix_str; + } /// \brief Prefix to be used for output restart files - virtual std::string restart_output_prefix() = 0; + inline std::string restart_output_prefix() + { + return restart_output_prefix_str; + } /// \brief Prefix to be used for output files (final system /// configuration) - virtual std::string output_prefix() = 0; + inline std::string output_prefix() + { + return output_prefix_str; + } /// \brief Restarts will be fritten each time this number of steps has passed virtual size_t restart_frequency() = 0; @@ -91,13 +132,39 @@ public: void select_closest_images (std::vector &pos, cvm::atom_pos const &ref_pos); + // **************** SCRIPTING INTERFACE **************** + /// Pointer to the scripting interface object + /// (does not need to be allocated in a new interface) + colvarscript *script; + + /// is a user force script defined? + bool force_script_defined; + + /// Do we have a scripting interface? + bool have_scripts; + + /// Run a user-defined colvar forces script + virtual int run_force_callback() { return COLVARS_NOT_IMPLEMENTED; } + + virtual int run_colvar_callback(std::string const &name, + std::vector const &cvcs, + colvarvalue &value) + { return COLVARS_NOT_IMPLEMENTED; } + + virtual int run_colvar_gradient_callback(std::string const &name, + std::vector const &cvcs, + std::vector &gradient) + { return COLVARS_NOT_IMPLEMENTED; } // **************** INPUT/OUTPUT **************** /// Print a message to the main log virtual void log (std::string const &message) = 0; + /// Print a message to the main log and let the rest of the program handle the error + virtual void error (std::string const &message) = 0; + /// Print a message to the main log and exit with error code virtual void fatal_error (std::string const &message) = 0; @@ -109,23 +176,23 @@ public: /// from "filename" will be appended \param pdb_field (optiona) if /// "filename" is a PDB file, use this field to determine which are /// the atoms to be set - virtual void load_atoms (char const *filename, + virtual int load_atoms (char const *filename, std::vector &atoms, - std::string const pdb_field, - double const pdb_field_value = 0.0) {} + std::string const &pdb_field, + double const pdb_field_value = 0.0) = 0; /// \brief Load the coordinates for a group of atoms from a file /// (usually a PDB); if "pos" is already allocated, the number of its /// elements must match the number of atoms in "filename" - virtual void load_coords (char const *filename, + virtual int load_coords (char const *filename, std::vector &pos, const std::vector &indices, - std::string const pdb_field, + std::string const &pdb_field, double const pdb_field_value = 0.0) = 0; /// \brief Rename the given file, before overwriting it - virtual void backup_file (char const *filename) {} - + virtual int backup_file (char const *filename) + { return COLVARS_NOT_IMPLEMENTED; } }; @@ -133,7 +200,7 @@ inline void colvarproxy::select_closest_images (std::vector &pos, cvm::atom_pos const &ref_pos) { for (std::vector::iterator pi = pos.begin(); - pi != pos.end(); pi++) { + pi != pos.end(); ++pi) { select_closest_image (*pi, ref_pos); } } @@ -146,8 +213,3 @@ inline cvm::real colvarproxy::position_dist2 (cvm::atom_pos const &pos1, #endif - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index 011f22cbdf..0d97740622 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #ifndef COLVARTYPES_H #define COLVARTYPES_H @@ -35,7 +37,7 @@ public: {} /// \brief Set all components to a scalar value - inline void set (cvm::real const value = 0.0) { + inline void set (cvm::real const &value = 0.0) { x = y = z = value; } @@ -539,13 +541,13 @@ inline cvm::rvector operator * (cvm::rmatrix const &m, /// Numerical recipes diagonalization -void jacobi (cvm::real **a, int n, cvm::real d[], cvm::real **v, int *nrot); +void jacobi (cvm::real **a, cvm::real d[], cvm::real **v, int *nrot); /// Eigenvector sort -void eigsrt (cvm::real d[], cvm::real **v, int n); +void eigsrt (cvm::real d[], cvm::real **v); /// Transpose the matrix -void transpose (cvm::real **v, int n); +void transpose (cvm::real **v); @@ -603,7 +605,7 @@ public: } /// \brief Set all components to a scalar - inline void set (cvm::real const value = 0.0) + inline void set (cvm::real const &value = 0.0) { q0 = q1 = q2 = q3 = value; } @@ -645,7 +647,7 @@ public: case 3: return this->q3; default: - cvm::fatal_error ("Error: incorrect quaternion component.\n"); + cvm::error ("Error: incorrect quaternion component.\n"); return q0; } } @@ -662,9 +664,9 @@ public: case 3: return this->q3; default: - cvm::fatal_error ("Error: trying to access a quaternion " - "component which is not between 0 and 3.\n"); - return this->q0; + cvm::error ("Error: trying to access a quaternion " + "component which is not between 0 and 3.\n"); + return 0.0; } } @@ -943,9 +945,9 @@ public: /// Constructor after a quaternion inline rotation (cvm::quaternion const &qi) - : b_debug_gradients (false) + : q (qi), + b_debug_gradients (false) { - q = qi; } /// Constructor after an axis of rotation and an angle (in radians) @@ -1091,8 +1093,3 @@ protected: #endif - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index 18a3bd00c5..fa209ee19d 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #include #include "colvarmodule.h" @@ -5,48 +7,60 @@ -std::string const colvarvalue::type_desc[colvarvalue::type_quaternion+1] = +std::string const colvarvalue::type_desc[colvarvalue::type_all+1] = { "not_set", "scalar number", "3-dimensional vector", "3-dimensional unit vector", - "4-dimensional unit vector" }; + "3-dimensional tangent vector", + "4-dimensional unit quaternion", + "4-dimensional tangent vector", + }; -size_t const colvarvalue::dof_num[ colvarvalue::type_quaternion+1] = - { 0, 1, 3, 2, 3 }; +std::string const colvarvalue::type_keyword[colvarvalue::type_all+1] = + { "not_set", + "scalar", + "vector", + "unit_vector", + "", + "unit_quaternion", + "", + }; + +size_t const colvarvalue::dof_num[ colvarvalue::type_all+1] = + { 0, 1, 3, 2, 2, 3, 3 }; void colvarvalue::undef_op() const { - cvm::fatal_error ("Error: Undefined operation on a colvar of type \""+ - colvarvalue::type_desc[this->value_type]+"\".\n"); + cvm::error ("Error: Undefined operation on a colvar of type \""+ + colvarvalue::type_desc[this->value_type]+"\".\n"); } void colvarvalue::error_rside (colvarvalue::Type const &vt) const { - cvm::fatal_error ("Trying to assign a colvar value with type \""+ - type_desc[this->value_type]+"\" to one with type \""+ - type_desc[vt]+"\".\n"); + cvm::error("Trying to assign a colvar value with type \""+ + type_desc[this->value_type]+"\" to one with type \""+ + type_desc[vt]+"\".\n"); } -void colvarvalue::error_lside -(colvarvalue::Type const &vt) const +void colvarvalue::error_lside(colvarvalue::Type const &vt) const { - cvm::fatal_error ("Trying to use a colvar value with type \""+ - type_desc[vt]+"\" as one of type \""+ - type_desc[this->value_type]+"\".\n"); + cvm::error("Trying to use a colvar value with type \""+ + type_desc[vt]+"\" as one of type \""+ + type_desc[this->value_type]+"\".\n"); } -void colvarvalue::inner_opt (colvarvalue const &x, - std::vector::iterator &xv, - std::vector::iterator const &xv_end, - std::vector::iterator &inner) +void colvarvalue::inner_opt(colvarvalue const &x, + std::vector::iterator &xv, + std::vector::iterator const &xv_end, + std::vector::iterator &inner) { // doing type check only once, here - colvarvalue::check_types (x, *xv); + colvarvalue::check_types(x, *xv); std::vector::iterator &xvi = xv; std::vector::iterator &ii = inner; @@ -59,11 +73,13 @@ void colvarvalue::inner_opt (colvarvalue const &x, break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: while (xvi != xv_end) { *(ii++) += (xvi++)->rvector_value * x.rvector_value; } break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: while (xvi != xv_end) { *(ii++) += ((xvi++)->quaternion_value).cosine (x.quaternion_value); } @@ -73,13 +89,13 @@ void colvarvalue::inner_opt (colvarvalue const &x, }; } -void colvarvalue::inner_opt (colvarvalue const &x, - std::list::iterator &xv, - std::list::iterator const &xv_end, - std::vector::iterator &inner) +void colvarvalue::inner_opt(colvarvalue const &x, + std::list::iterator &xv, + std::list::iterator const &xv_end, + std::vector::iterator &inner) { // doing type check only once, here - colvarvalue::check_types (x, *xv); + colvarvalue::check_types(x, *xv); std::list::iterator &xvi = xv; std::vector::iterator &ii = inner; @@ -92,11 +108,13 @@ void colvarvalue::inner_opt (colvarvalue const &x, break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: while (xvi != xv_end) { *(ii++) += (xvi++)->rvector_value * x.rvector_value; } break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: while (xvi != xv_end) { *(ii++) += ((xvi++)->quaternion_value).cosine (x.quaternion_value); } @@ -107,21 +125,22 @@ void colvarvalue::inner_opt (colvarvalue const &x, } -void colvarvalue::p2leg_opt (colvarvalue const &x, - std::vector::iterator &xv, - std::vector::iterator const &xv_end, - std::vector::iterator &inner) +void colvarvalue::p2leg_opt(colvarvalue const &x, + std::vector::iterator &xv, + std::vector::iterator const &xv_end, + std::vector::iterator &inner) { // doing type check only once, here - colvarvalue::check_types (x, *xv); + colvarvalue::check_types(x, *xv); std::vector::iterator &xvi = xv; std::vector::iterator &ii = inner; switch (x.value_type) { case colvarvalue::type_scalar: - cvm::fatal_error ("Error: cannot calculate Legendre polynomials " - "for scalar variables.\n"); + cvm::error("Error: cannot calculate Legendre polynomials " + "for scalar variables.\n"); + return; break; case colvarvalue::type_vector: while (xvi != xv_end) { @@ -133,12 +152,14 @@ void colvarvalue::p2leg_opt (colvarvalue const &x, } break; case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: while (xvi != xv_end) { cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value; *(ii++) += 1.5*cosine*cosine - 0.5; } break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: while (xvi != xv_end) { cvm::real const cosine = (xvi++)->quaternion_value.cosine (x.quaternion_value); *(ii++) += 1.5*cosine*cosine - 0.5; @@ -149,21 +170,21 @@ void colvarvalue::p2leg_opt (colvarvalue const &x, }; } -void colvarvalue::p2leg_opt (colvarvalue const &x, - std::list::iterator &xv, - std::list::iterator const &xv_end, - std::vector::iterator &inner) +void colvarvalue::p2leg_opt(colvarvalue const &x, + std::list::iterator &xv, + std::list::iterator const &xv_end, + std::vector::iterator &inner) { // doing type check only once, here - colvarvalue::check_types (x, *xv); + colvarvalue::check_types(x, *xv); std::list::iterator &xvi = xv; std::vector::iterator &ii = inner; switch (x.value_type) { case colvarvalue::type_scalar: - cvm::fatal_error ("Error: cannot calculate Legendre polynomials " - "for scalar variables.\n"); + cvm::error("Error: cannot calculate Legendre polynomials " + "for scalar variables.\n"); break; case colvarvalue::type_vector: while (xvi != xv_end) { @@ -175,12 +196,14 @@ void colvarvalue::p2leg_opt (colvarvalue const &x, } break; case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: while (xvi != xv_end) { cvm::real const cosine = (xvi++)->rvector_value * x.rvector_value; *(ii++) += 1.5*cosine*cosine - 0.5; } break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: while (xvi != xv_end) { cvm::real const cosine = (xvi++)->quaternion_value.cosine (x.quaternion_value); *(ii++) += 1.5*cosine*cosine - 0.5; @@ -196,12 +219,17 @@ std::ostream & operator << (std::ostream &os, colvarvalue const &x) { switch (x.type()) { case colvarvalue::type_scalar: - os << x.real_value; break; + os << x.real_value; + break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: - os << x.rvector_value; break; + case colvarvalue::type_unitvectorderiv: + os << x.rvector_value; + break; case colvarvalue::type_quaternion: - os << x.quaternion_value; break; + case colvarvalue::type_quaternionderiv: + os << x.quaternion_value; + break; case colvarvalue::type_notset: os << "not set"; break; } @@ -211,7 +239,9 @@ std::ostream & operator << (std::ostream &os, colvarvalue const &x) std::ostream & operator << (std::ostream &os, std::vector const &v) { - for (size_t i = 0; i < v.size(); i++) os << v[i]; + for (size_t i = 0; i < v.size(); i++) { + os << v[i]; + } return os; } @@ -219,8 +249,9 @@ std::ostream & operator << (std::ostream &os, std::vector const &v) std::istream & operator >> (std::istream &is, colvarvalue &x) { if (x.type() == colvarvalue::type_notset) { - cvm::fatal_error ("Trying to read from a stream a colvarvalue, " - "which has not yet been assigned a data type.\n"); + cvm::error("Trying to read from a stream a colvarvalue, " + "which has not yet been assigned a data type.\n"); + return is; } switch (x.type()) { @@ -228,30 +259,39 @@ std::istream & operator >> (std::istream &is, colvarvalue &x) is >> x.real_value; break; case colvarvalue::type_vector: - case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: is >> x.rvector_value; break; + case colvarvalue::type_unitvector: + is >> x.rvector_value; + x.apply_constraints(); + break; case colvarvalue::type_quaternion: + is >> x.quaternion_value; + x.apply_constraints(); + break; + case colvarvalue::type_quaternionderiv: is >> x.quaternion_value; break; default: x.undef_op(); } - x.apply_constraints(); return is; } -size_t colvarvalue::output_width (size_t const &real_width) const +size_t colvarvalue::output_width(size_t const &real_width) const { switch (this->value_type) { case colvarvalue::type_scalar: return real_width; case colvarvalue::type_vector: case colvarvalue::type_unitvector: - return cvm::rvector::output_width (real_width); + case colvarvalue::type_unitvectorderiv: + return cvm::rvector::output_width(real_width); case colvarvalue::type_quaternion: - return cvm::quaternion::output_width (real_width); + case colvarvalue::type_quaternionderiv: + return cvm::quaternion::output_width(real_width); case colvarvalue::type_notset: default: return 0; diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index e262e7febe..20c960259c 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -1,3 +1,5 @@ +/// -*- c++ -*- + #ifndef COLVARVALUE_H #define COLVARVALUE_H @@ -72,15 +74,22 @@ public: type_vector, /// 3-dimensional unit vector, implemented as \link colvarmodule::rvector \endlink type_unitvector, + /// 3-dimensional vector that is a derivative of a unitvector + type_unitvectorderiv, /// 4-dimensional unit vector representing a rotation, implemented as \link colvarmodule::quaternion \endlink - type_quaternion + type_quaternion, + /// 4-dimensional vector that is a derivative of a quaternion + type_quaternionderiv, + /// bogus type to hold the size of the enum + type_all }; /// Runtime description of value types - std::string static const type_desc[colvarvalue::type_quaternion+1]; - + std::string static const type_desc[colvarvalue::type_all+1]; + /// User keywords for specifying value types + std::string static const type_keyword[colvarvalue::type_all+1]; /// Number of degrees of freedom for each type - size_t static const dof_num[ colvarvalue::type_quaternion+1]; + size_t static const dof_num[ colvarvalue::type_all+1]; /// \brief Real data member cvm::real real_value; @@ -102,43 +111,43 @@ public: /// \brief Default constructor: this class defaults to a scalar /// number and always behaves like it unless you change its type inline colvarvalue() - : real_value (0.0), value_type (type_scalar) + : real_value(0.0), value_type(type_scalar) {} /// Constructor from a type specification - inline colvarvalue (Type const &vti) - : value_type (vti) + inline colvarvalue(Type const &vti) + : value_type(vti) { reset(); } /// Copy constructor from real base type - inline colvarvalue (cvm::real const &x) - : real_value (x), value_type (type_scalar) + inline colvarvalue(cvm::real const &x) + : real_value(x), value_type(type_scalar) {} /// \brief Copy constructor from rvector base type (Note: this sets /// automatically a type \link type_vector \endlink , if you want a /// \link type_unitvector \endlink you must set it explicitly) - inline colvarvalue (cvm::rvector const &v) - : rvector_value (v), value_type (type_vector) + inline colvarvalue(cvm::rvector const &v) + : rvector_value(v), value_type(type_vector) {} /// \brief Copy constructor from rvector base type (additional /// argument to make possible to choose a \link type_unitvector /// \endlink - inline colvarvalue (cvm::rvector const &v, Type const &vti) - : rvector_value (v), value_type (vti) + inline colvarvalue(cvm::rvector const &v, Type const &vti) + : rvector_value(v), value_type(vti) {} /// \brief Copy constructor from quaternion base type - inline colvarvalue (cvm::quaternion const &q) - : quaternion_value (q), value_type (type_quaternion) + inline colvarvalue(cvm::quaternion const &q) + : quaternion_value(q), value_type(type_quaternion) {} /// Copy constructor from another \link colvarvalue \endlink - inline colvarvalue (colvarvalue const &x) - : value_type (x.value_type) + inline colvarvalue(colvarvalue const &x) + : value_type(x.value_type) { reset(); @@ -160,18 +169,15 @@ public: } - /// Set to the null value for the data type currently defined void reset(); - /// \brief If the variable has constraints (e.g. unitvector or - /// quaternion), transform it to satisfy them; use it when the \link - /// colvarvalue \endlink is not calculated from \link cvc - /// \endlink objects, but manipulated by you + /// quaternion), transform it to satisfy them; this function needs + /// to be called only when the \link colvarvalue \endlink + /// is calculated outside of \link cvc \endlink objects void apply_constraints(); - /// Get the current type inline Type type() const { @@ -179,7 +185,7 @@ public: } /// Set the type explicitly - inline void type (Type const &vti) + inline void type(Type const &vti) { reset(); value_type = vti; @@ -187,20 +193,24 @@ public: } /// Set the type after another \link colvarvalue \endlink - inline void type (colvarvalue const &x) + inline void type(colvarvalue const &x) { reset(); value_type = x.value_type; reset(); } + /// Make the type a derivative of the original type + /// (constraints do not apply on time derivatives of vector values) + inline void is_derivative(); + /// Square norm of this colvarvalue cvm::real norm2() const; /// Norm of this colvarvalue inline cvm::real norm() const { - return std::sqrt (this->norm2()); + return std::sqrt(this->norm2()); } /// \brief Return the value whose scalar product with this value is @@ -208,10 +218,10 @@ public: inline colvarvalue inverse() const; /// Square distance between this \link colvarvalue \endlink and another - cvm::real dist2 (colvarvalue const &x2) const; + cvm::real dist2(colvarvalue const &x2) const; /// Derivative with respect to this \link colvarvalue \endlink of the square distance - colvarvalue dist2_grad (colvarvalue const &x2) const; + colvarvalue dist2_grad(colvarvalue const &x2) const; /// Assignment operator (type of x is checked) colvarvalue & operator = (colvarvalue const &x); @@ -225,50 +235,57 @@ public: // Casting operator inline operator cvm::real() const { - if (value_type != type_scalar) error_rside (type_scalar); + if (value_type != type_scalar) { + error_rside(type_scalar); + } return real_value; } // Casting operator inline operator cvm::rvector() const { - if ( (value_type != type_vector) && (value_type != type_unitvector)) - error_rside (type_vector); + if ((value_type != type_vector) && + (value_type != type_unitvector) && + (value_type != type_unitvectorderiv)) { + error_rside(type_vector); + } return rvector_value; } // Casting operator inline operator cvm::quaternion() const { - if (value_type != type_quaternion) error_rside (type_quaternion); + if ((value_type != type_quaternion) && + (value_type != type_quaternionderiv)) { + error_rside(type_quaternion); + } return quaternion_value; } /// Special case when the variable is a real number, and all operations are defined - inline bool is_scalar() + inline bool is_scalar() const { return (value_type == type_scalar); } - /// Ensure that the two types are the same within a binary operator - void static check_types (colvarvalue const &x1, colvarvalue const &x2); + void static check_types(colvarvalue const &x1, colvarvalue const &x2); /// Undefined operation void undef_op() const; /// Trying to assign this \link colvarvalue \endlink object to /// another object set with a different type - void error_lside (Type const &vt) const; + void error_lside(Type const &vt) const; /// Trying to assign another \link colvarvalue \endlink object set /// with a different type to this object - void error_rside (Type const &vt) const; + void error_rside(Type const &vt) const; /// Give the number of characters required to output this /// colvarvalue, given the current type assigned and the number of /// characters for a real number - size_t output_width (size_t const &real_width) const; + size_t output_width(size_t const &real_width) const; // optimized routines for operations with an array; xv and inner as @@ -277,32 +294,32 @@ public: /// \brief Optimized routine for the inner product of one collective /// variable with an array - void static inner_opt (colvarvalue const &x, - std::vector::iterator &xv, - std::vector::iterator const &xv_end, - std::vector::iterator &inner); + void static inner_opt(colvarvalue const &x, + std::vector::iterator &xv, + std::vector::iterator const &xv_end, + std::vector::iterator &inner); /// \brief Optimized routine for the inner product of one collective /// variable with an array - void static inner_opt (colvarvalue const &x, - std::list::iterator &xv, - std::list::iterator const &xv_end, - std::vector::iterator &inner); + void static inner_opt(colvarvalue const &x, + std::list::iterator &xv, + std::list::iterator const &xv_end, + std::vector::iterator &inner); /// \brief Optimized routine for the second order Legendre /// polynomial, (3cos^2(w)-1)/2, of one collective variable with an /// array - void static p2leg_opt (colvarvalue const &x, - std::vector::iterator &xv, - std::vector::iterator const &xv_end, - std::vector::iterator &inner); + void static p2leg_opt(colvarvalue const &x, + std::vector::iterator &xv, + std::vector::iterator const &xv_end, + std::vector::iterator &inner); /// \brief Optimized routine for the second order Legendre /// polynomial of one collective variable with an array - void static p2leg_opt (colvarvalue const &x, - std::list::iterator &xv, - std::list::iterator const &xv_end, - std::vector::iterator &inner); + void static p2leg_opt(colvarvalue const &x, + std::list::iterator &xv, + std::list::iterator const &xv_end, + std::vector::iterator &inner); }; @@ -312,14 +329,16 @@ inline void colvarvalue::reset() { switch (value_type) { case colvarvalue::type_scalar: - real_value = cvm::real (0.0); + real_value = cvm::real(0.0); break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: - rvector_value = cvm::rvector (0.0, 0.0, 0.0); + case colvarvalue::type_unitvectorderiv: + rvector_value = cvm::rvector(0.0, 0.0, 0.0); break; case colvarvalue::type_quaternion: - quaternion_value = cvm::quaternion (0.0, 0.0, 0.0, 0.0); + case colvarvalue::type_quaternionderiv: + quaternion_value = cvm::quaternion(0.0, 0.0, 0.0, 0.0); break; case colvarvalue::type_notset: default: @@ -332,14 +351,36 @@ inline void colvarvalue::apply_constraints() { switch (value_type) { case colvarvalue::type_scalar: - break; case colvarvalue::type_vector: - break; + case colvarvalue::type_unitvectorderiv: + case colvarvalue::type_quaternionderiv: + break; case colvarvalue::type_unitvector: - rvector_value /= std::sqrt (rvector_value.norm2()); + rvector_value /= std::sqrt(rvector_value.norm2()); break; case colvarvalue::type_quaternion: - quaternion_value /= std::sqrt (quaternion_value.norm2()); + quaternion_value /= std::sqrt(quaternion_value.norm2()); + break; + case colvarvalue::type_notset: + default: + break; + } +} + + +inline void colvarvalue::is_derivative() +{ + switch (value_type) { + case colvarvalue::type_scalar: + case colvarvalue::type_vector: + case colvarvalue::type_unitvectorderiv: + case colvarvalue::type_quaternionderiv: + break; + case colvarvalue::type_unitvector: + type(colvarvalue::type_unitvectorderiv); + break; + case colvarvalue::type_quaternion: + type(colvarvalue::type_quaternionderiv); break; case colvarvalue::type_notset: default: @@ -356,8 +397,10 @@ inline cvm::real colvarvalue::norm2() const return (this->real_value)*(this->real_value); case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: return (this->rvector_value).norm2(); case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: return (this->quaternion_value).norm2(); case colvarvalue::type_notset: default: @@ -370,16 +413,19 @@ inline colvarvalue colvarvalue::inverse() const { switch (value_type) { case colvarvalue::type_scalar: - return colvarvalue (1.0/real_value); + return colvarvalue(1.0/real_value); case colvarvalue::type_vector: - return colvarvalue (cvm::rvector (1.0/rvector_value.x, - 1.0/rvector_value.y, - 1.0/rvector_value.z)); + case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: + return colvarvalue(cvm::rvector(1.0/rvector_value.x, + 1.0/rvector_value.y, + 1.0/rvector_value.z)); case colvarvalue::type_quaternion: - return colvarvalue (cvm::quaternion (1.0/quaternion_value.q0, - 1.0/quaternion_value.q1, - 1.0/quaternion_value.q2, - 1.0/quaternion_value.q3)); + case colvarvalue::type_quaternionderiv: + return colvarvalue(cvm::quaternion(1.0/quaternion_value.q0, + 1.0/quaternion_value.q1, + 1.0/quaternion_value.q2, + 1.0/quaternion_value.q3)); case colvarvalue::type_notset: default: undef_op(); @@ -387,12 +433,11 @@ inline colvarvalue colvarvalue::inverse() const return colvarvalue(); } - inline colvarvalue & colvarvalue::operator = (colvarvalue const &x) { if (this->value_type != type_notset) if (this->value_type != x.value_type) - error_lside (x.value_type); + error_lside(x.value_type); this->value_type = x.value_type; @@ -402,9 +447,11 @@ inline colvarvalue & colvarvalue::operator = (colvarvalue const &x) break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: this->rvector_value = x.rvector_value; break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: this->quaternion_value = x.quaternion_value; break; case colvarvalue::type_notset: @@ -418,7 +465,7 @@ inline colvarvalue & colvarvalue::operator = (colvarvalue const &x) inline void colvarvalue::operator += (colvarvalue const &x) { if (colvarvalue::type_checking()) - colvarvalue::check_types (*this, x); + colvarvalue::check_types(*this, x); switch (this->value_type) { case colvarvalue::type_scalar: @@ -426,9 +473,11 @@ inline void colvarvalue::operator += (colvarvalue const &x) break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: this->rvector_value += x.rvector_value; break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: this->quaternion_value += x.quaternion_value; break; case colvarvalue::type_notset: @@ -440,16 +489,21 @@ inline void colvarvalue::operator += (colvarvalue const &x) inline void colvarvalue::operator -= (colvarvalue const &x) { if (colvarvalue::type_checking()) - colvarvalue::check_types (*this, x); + colvarvalue::check_types(*this, x); switch (value_type) { case colvarvalue::type_scalar: - real_value -= x.real_value; break; + real_value -= x.real_value; + break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: - rvector_value -= x.rvector_value; break; + case colvarvalue::type_unitvectorderiv: + rvector_value -= x.rvector_value; + break; case colvarvalue::type_quaternion: - quaternion_value -= x.quaternion_value; break; + case colvarvalue::type_quaternionderiv: + quaternion_value -= x.quaternion_value; + break; case colvarvalue::type_notset: default: undef_op(); @@ -463,10 +517,11 @@ inline void colvarvalue::operator *= (cvm::real const &a) real_value *= a; break; case colvarvalue::type_vector: - case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: rvector_value *= a; break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: quaternion_value *= a; break; case colvarvalue::type_notset: @@ -482,8 +537,10 @@ inline void colvarvalue::operator /= (cvm::real const &a) real_value /= a; break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: rvector_value /= a; break; case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: quaternion_value /= a; break; case colvarvalue::type_notset: default: @@ -498,22 +555,24 @@ inline colvarvalue operator + (colvarvalue const &x1, colvarvalue const &x2) { if (colvarvalue::type_checking()) - colvarvalue::check_types (x1, x2); + colvarvalue::check_types(x1, x2); switch (x1.value_type) { case colvarvalue::type_scalar: - return colvarvalue (x1.real_value + x2.real_value); + return colvarvalue(x1.real_value + x2.real_value); case colvarvalue::type_vector: - return colvarvalue (x1.rvector_value + x2.rvector_value); + return colvarvalue(x1.rvector_value + x2.rvector_value); case colvarvalue::type_unitvector: - return colvarvalue (x1.rvector_value + x2.rvector_value, - colvarvalue::type_unitvector); + case colvarvalue::type_unitvectorderiv: + return colvarvalue(x1.rvector_value + x2.rvector_value, + colvarvalue::type_unitvector); case colvarvalue::type_quaternion: - return colvarvalue (x1.quaternion_value + x2.quaternion_value); + case colvarvalue::type_quaternionderiv: + return colvarvalue(x1.quaternion_value + x2.quaternion_value); case colvarvalue::type_notset: default: x1.undef_op(); - return colvarvalue (colvarvalue::type_notset); + return colvarvalue(colvarvalue::type_notset); }; } @@ -521,21 +580,23 @@ inline colvarvalue operator - (colvarvalue const &x1, colvarvalue const &x2) { if (colvarvalue::type_checking()) - colvarvalue::check_types (x1, x2); + colvarvalue::check_types(x1, x2); switch (x1.value_type) { case colvarvalue::type_scalar: - return colvarvalue (x1.real_value - x2.real_value); + return colvarvalue(x1.real_value - x2.real_value); case colvarvalue::type_vector: - return colvarvalue (x1.rvector_value - x2.rvector_value); + return colvarvalue(x1.rvector_value - x2.rvector_value); case colvarvalue::type_unitvector: - return colvarvalue (x1.rvector_value - x2.rvector_value, - colvarvalue::type_unitvector); + case colvarvalue::type_unitvectorderiv: + return colvarvalue(x1.rvector_value - x2.rvector_value, + colvarvalue::type_unitvector); case colvarvalue::type_quaternion: - return colvarvalue (x1.quaternion_value - x2.quaternion_value); + case colvarvalue::type_quaternionderiv: + return colvarvalue(x1.quaternion_value - x2.quaternion_value); default: x1.undef_op(); - return colvarvalue (colvarvalue::type_notset); + return colvarvalue(colvarvalue::type_notset); }; } @@ -547,18 +608,20 @@ inline colvarvalue operator * (cvm::real const &a, { switch (x.value_type) { case colvarvalue::type_scalar: - return colvarvalue (a * x.real_value); + return colvarvalue(a * x.real_value); case colvarvalue::type_vector: - return colvarvalue (a * x.rvector_value); + return colvarvalue(a * x.rvector_value); case colvarvalue::type_unitvector: - return colvarvalue (a * x.rvector_value, - colvarvalue::type_unitvector); + case colvarvalue::type_unitvectorderiv: + return colvarvalue(a * x.rvector_value, + colvarvalue::type_unitvector); case colvarvalue::type_quaternion: - return colvarvalue (a * x.quaternion_value); + case colvarvalue::type_quaternionderiv: + return colvarvalue(a * x.quaternion_value); case colvarvalue::type_notset: default: x.undef_op(); - return colvarvalue (colvarvalue::type_notset); + return colvarvalue(colvarvalue::type_notset); } } @@ -567,18 +630,20 @@ inline colvarvalue operator * (colvarvalue const &x, { switch (x.value_type) { case colvarvalue::type_scalar: - return colvarvalue (x.real_value * a); + return colvarvalue(x.real_value * a); case colvarvalue::type_vector: - return colvarvalue (x.rvector_value * a); + return colvarvalue(x.rvector_value * a); case colvarvalue::type_unitvector: - return colvarvalue (x.rvector_value * a, - colvarvalue::type_unitvector); + case colvarvalue::type_unitvectorderiv: + return colvarvalue(x.rvector_value * a, + colvarvalue::type_unitvector); case colvarvalue::type_quaternion: - return colvarvalue (x.quaternion_value * a); + case colvarvalue::type_quaternionderiv: + return colvarvalue(x.quaternion_value * a); case colvarvalue::type_notset: default: x.undef_op(); - return colvarvalue (colvarvalue::type_notset); + return colvarvalue(colvarvalue::type_notset); } } @@ -587,18 +652,20 @@ inline colvarvalue operator / (colvarvalue const &x, { switch (x.value_type) { case colvarvalue::type_scalar: - return colvarvalue (x.real_value / a); + return colvarvalue(x.real_value / a); case colvarvalue::type_vector: - return colvarvalue (x.rvector_value / a); + return colvarvalue(x.rvector_value / a); case colvarvalue::type_unitvector: - return colvarvalue (x.rvector_value / a, - colvarvalue::type_unitvector); + case colvarvalue::type_unitvectorderiv: + return colvarvalue(x.rvector_value / a, + colvarvalue::type_unitvector); case colvarvalue::type_quaternion: - return colvarvalue (x.quaternion_value / a); + case colvarvalue::type_quaternionderiv: + return colvarvalue(x.quaternion_value / a); case colvarvalue::type_notset: default: x.undef_op(); - return colvarvalue (colvarvalue::type_notset); + return colvarvalue(colvarvalue::type_notset); } } @@ -609,18 +676,20 @@ inline cvm::real operator * (colvarvalue const &x1, colvarvalue const &x2) { if (colvarvalue::type_checking()) - colvarvalue::check_types (x1, x2); + colvarvalue::check_types(x1, x2); switch (x1.value_type) { case colvarvalue::type_scalar: return (x1.real_value * x2.real_value); case colvarvalue::type_vector: case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: return (x1.rvector_value * x2.rvector_value); case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: // the "*" product is the quaternion product, here the inner // member function is used instead - return (x1.quaternion_value.inner (x2.quaternion_value)); + return (x1.quaternion_value.inner(x2.quaternion_value)); case colvarvalue::type_notset: default: x1.undef_op(); @@ -630,10 +699,10 @@ inline cvm::real operator * (colvarvalue const &x1, -inline cvm::real colvarvalue::dist2 (colvarvalue const &x2) const +inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const { if (colvarvalue::type_checking()) - colvarvalue::check_types (*this, x2); + colvarvalue::check_types(*this, x2); switch (this->value_type) { case colvarvalue::type_scalar: @@ -641,12 +710,14 @@ inline cvm::real colvarvalue::dist2 (colvarvalue const &x2) const case colvarvalue::type_vector: return (this->rvector_value - x2.rvector_value).norm2(); case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: // angle between (*this) and x2 is the distance - return std::acos (this->rvector_value * x2.rvector_value) * std::acos (this->rvector_value * x2.rvector_value); + return std::acos(this->rvector_value * x2.rvector_value) * std::acos(this->rvector_value * x2.rvector_value); case colvarvalue::type_quaternion: + case colvarvalue::type_quaternionderiv: // angle between (*this) and x2 is the distance, the quaternion // object has it implemented internally - return this->quaternion_value.dist2 (x2.quaternion_value); + return this->quaternion_value.dist2(x2.quaternion_value); case colvarvalue::type_notset: default: this->undef_op(); @@ -655,10 +726,10 @@ inline cvm::real colvarvalue::dist2 (colvarvalue const &x2) const } -inline colvarvalue colvarvalue::dist2_grad (colvarvalue const &x2) const +inline colvarvalue colvarvalue::dist2_grad(colvarvalue const &x2) const { if (colvarvalue::type_checking()) - colvarvalue::check_types (*this, x2); + colvarvalue::check_types(*this, x2); switch (this->value_type) { case colvarvalue::type_scalar: @@ -666,62 +737,61 @@ inline colvarvalue colvarvalue::dist2_grad (colvarvalue const &x2) const case colvarvalue::type_vector: return 2.0 * (this->rvector_value - x2.rvector_value); case colvarvalue::type_unitvector: + case colvarvalue::type_unitvectorderiv: { cvm::rvector const &v1 = this->rvector_value; cvm::rvector const &v2 = x2.rvector_value; cvm::real const cos_t = v1 * v2; - cvm::real const sin_t = std::sqrt (1.0 - cos_t*cos_t); - return colvarvalue ( 2.0 * sin_t * - cvm::rvector ( (-1.0) * sin_t * v2.x + - cos_t/sin_t * (v1.x - cos_t*v2.x), - (-1.0) * sin_t * v2.y + - cos_t/sin_t * (v1.y - cos_t*v2.y), - (-1.0) * sin_t * v2.z + - cos_t/sin_t * (v1.z - cos_t*v2.z) - ), - colvarvalue::type_unitvector ); + cvm::real const sin_t = std::sqrt(1.0 - cos_t*cos_t); + return colvarvalue( 2.0 * sin_t * + cvm::rvector((-1.0) * sin_t * v2.x + + cos_t/sin_t * (v1.x - cos_t*v2.x), + (-1.0) * sin_t * v2.y + + cos_t/sin_t * (v1.y - cos_t*v2.y), + (-1.0) * sin_t * v2.z + + cos_t/sin_t * (v1.z - cos_t*v2.z) + ), + colvarvalue::type_unitvector ); } case colvarvalue::type_quaternion: - return this->quaternion_value.dist2_grad (x2.quaternion_value); + case colvarvalue::type_quaternionderiv: + return this->quaternion_value.dist2_grad(x2.quaternion_value); case colvarvalue::type_notset: default: this->undef_op(); - return colvarvalue (colvarvalue::type_notset); + return colvarvalue(colvarvalue::type_notset); }; } -inline void colvarvalue::check_types (colvarvalue const &x1, - colvarvalue const &x2) +inline void colvarvalue::check_types(colvarvalue const &x1, + colvarvalue const &x2) { if (x1.value_type != x2.value_type) { - cvm::log ("x1 type = "+cvm::to_str (x1.value_type)+ - ", x2 type = "+cvm::to_str (x2.value_type)+"\n"); - cvm::fatal_error ("Performing an operation between two colvar " - "values with different types, \""+ - colvarvalue::type_desc[x1.value_type]+ - "\" and \""+ - colvarvalue::type_desc[x2.value_type]+ - "\".\n"); + if (((x1.value_type == type_unitvector) && + (x2.value_type == type_unitvectorderiv)) || + ((x2.value_type == type_unitvector) && + (x1.value_type == type_unitvectorderiv)) || + ((x1.value_type == type_quaternion) && + (x2.value_type == type_quaternionderiv)) || + ((x2.value_type == type_quaternion) && + (x1.value_type == type_quaternionderiv))) { + return; + } + cvm::error("Performing an operation between two colvar " + "values with different types, \""+ + colvarvalue::type_desc[x1.value_type]+ + "\" and \""+ + colvarvalue::type_desc[x2.value_type]+ + "\".\n"); } } - - std::ostream & operator << (std::ostream &os, colvarvalue const &x); std::ostream & operator << (std::ostream &os, std::vector const &v); std::istream & operator >> (std::istream &is, colvarvalue &x); - - - #endif - - -// Emacs -// Local Variables: -// mode: C++ -// End: