From f539e43b2253f14930e4c3fd79403f8ee55dee82 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 30 Apr 2015 14:09:42 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13444 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/colvars/Makefile.fermi | 21 +- lib/colvars/Makefile.g++ | 35 +- lib/colvars/Makefile.mingw32-cross | 20 +- lib/colvars/Makefile.mingw64-cross | 20 +- lib/colvars/colvar.cpp | 5 +- lib/colvars/colvar.h | 2 +- lib/colvars/colvaratoms.cpp | 6 +- lib/colvars/colvaratoms.h | 4 - lib/colvars/colvarbias.h | 14 +- lib/colvars/colvarbias_abf.cpp | 205 +--------- lib/colvars/colvarbias_abf.h | 36 +- lib/colvars/colvarbias_meta.cpp | 573 ++++++++++++++------------- lib/colvars/colvarbias_meta.h | 2 + lib/colvars/colvarbias_restraint.h | 2 +- lib/colvars/colvarcomp_distances.cpp | 37 +- lib/colvars/colvarcomp_rotations.cpp | 4 +- lib/colvars/colvarmodule.cpp | 45 ++- lib/colvars/colvarmodule.h | 14 +- lib/colvars/colvarparse.cpp | 2 + lib/colvars/colvarparse.h | 2 + lib/colvars/colvarscript.cpp | 13 +- lib/colvars/colvartypes.cpp | 1 - 22 files changed, 438 insertions(+), 625 deletions(-) diff --git a/lib/colvars/Makefile.fermi b/lib/colvars/Makefile.fermi index 1f4b5a5b01..2923b564e1 100644 --- a/lib/colvars/Makefile.fermi +++ b/lib/colvars/Makefile.fermi @@ -1,4 +1,4 @@ -# library build makefile for colvars module +# library build -*- makefile -*- for colvars module # which file will be copied to Makefile.lammps @@ -16,11 +16,11 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ 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 + colvarbias_histogram.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) @@ -63,6 +63,9 @@ colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.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_histogram.o: colvarbias_histogram.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias_histogram.h colvarbias.h colvargrid.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 @@ -93,9 +96,9 @@ 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_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ - colvarbias_abf.h colvarscript.h + colvarproxy.h colvarvalue.h colvarparse.h colvar.h colvarbias.h \ + colvarbias_abf.h colvargrid.h colvarbias_alb.h colvarbias_restraint.h \ + colvarbias_histogram.h colvarbias_meta.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 \ diff --git a/lib/colvars/Makefile.g++ b/lib/colvars/Makefile.g++ index 6e5d4463ea..8f66b0c743 100644 --- a/lib/colvars/Makefile.g++ +++ b/lib/colvars/Makefile.g++ @@ -15,11 +15,11 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ 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 + colvarbias_histogram.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) @@ -53,15 +53,18 @@ colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) # ------ DEPENDENCIES ------ # colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarparse.h colvarvalue.h colvaratoms.h + colvarvalue.h colvarparse.h colvaratoms.h colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ - colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \ + colvarproxy.h colvarvalue.h colvar.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 \ + colvarproxy.h colvarvalue.h colvarbias_alb.h colvar.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_histogram.o: colvarbias_histogram.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias_histogram.h colvarbias.h colvargrid.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 @@ -69,10 +72,10 @@ 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 \ + 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 \ + colvartypes.h colvarproxy.h colvarvalue.h colvarparse.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 @@ -92,16 +95,16 @@ 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_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ - colvarbias_abf.h colvarscript.h + colvarproxy.h colvarvalue.h colvarparse.h colvar.h colvarbias.h \ + colvarbias_abf.h colvargrid.h colvarbias_alb.h colvarbias_restraint.h \ + colvarbias_histogram.h colvarbias_meta.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.h colvarproxy.h colvarvalue.h colvarbias.h colvar.h \ + colvarparse.h colvartypes.o: colvartypes.cpp colvarmodule.h colvartypes.h colvarproxy.h \ - colvarparse.h colvarvalue.h + colvarvalue.h colvarparse.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 39b07ca8d7..381f4455d5 100644 --- a/lib/colvars/Makefile.mingw32-cross +++ b/lib/colvars/Makefile.mingw32-cross @@ -18,11 +18,11 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ 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 + colvarbias_histogram.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 @@ -70,6 +70,9 @@ $(DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.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_histogram.o: colvarbias_histogram.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias_histogram.h colvarbias.h colvargrid.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 @@ -100,9 +103,9 @@ $(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_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ - colvarbias_abf.h colvarscript.h + colvarproxy.h colvarvalue.h colvarparse.h colvar.h colvarbias.h \ + colvarbias_abf.h colvargrid.h colvarbias_alb.h colvarbias_restraint.h \ + colvarbias_histogram.h colvarbias_meta.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 \ @@ -118,4 +121,3 @@ $(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h clean: -rm $(DIR)*.o *~ $(LIB) -rmdir $(DIR) - diff --git a/lib/colvars/Makefile.mingw64-cross b/lib/colvars/Makefile.mingw64-cross index 5101cf7805..bcb23c6475 100644 --- a/lib/colvars/Makefile.mingw64-cross +++ b/lib/colvars/Makefile.mingw64-cross @@ -18,11 +18,11 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ 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 + colvarbias_histogram.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 @@ -41,6 +41,7 @@ Makefile.lammps: $(LIB): $(DIR) $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) + @cp $(EXTRAMAKE) Makefile.lammps $(DIR)colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) $(CXX) -o $@ $(CXXFLAGS) $^ @@ -69,6 +70,9 @@ $(DIR)colvarbias_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.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_histogram.o: colvarbias_histogram.cpp colvarmodule.h \ + colvartypes.h colvarproxy.h colvarvalue.h colvar.h colvarparse.h \ + colvarbias_histogram.h colvarbias.h colvargrid.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 @@ -99,9 +103,9 @@ $(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_alb.h colvarbias_restraint.h colvarbias_meta.h colvargrid.h \ - colvarbias_abf.h colvarscript.h + colvarproxy.h colvarvalue.h colvarparse.h colvar.h colvarbias.h \ + colvarbias_abf.h colvargrid.h colvarbias_alb.h colvarbias_restraint.h \ + colvarbias_histogram.h colvarbias_meta.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 \ diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 80c04f9f8b..b60c685395 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -20,9 +20,10 @@ colvar::colvar(std::string const &conf) (std::string("colvar")+cvm::to_str(cvm::colvars.size()+1))); if (cvm::colvar_by_name(this->name) != NULL) { - cvm::error("Error: this colvar cannot have the same name, \""+this->name+ + cvm::error("Error: this colvar cannot have the same name, \""+this->name+ "\", as another colvar.\n", INPUT_ERROR); + return; } // all tasks disabled by default @@ -1616,7 +1617,7 @@ int colvar::write_output_files() // ******************** ANALYSIS FUNCTIONS ******************** -void colvar::analyse() +void colvar::analyze() { if (tasks[task_runave]) { calc_runave(); diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 2a9abb9d47..beadc52f51 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -361,7 +361,7 @@ public: /// Read the analysis tasks int parse_analysis(std::string const &conf); /// Perform analysis tasks - void analyse(); + void analyze(); /// Read the value from a collective variable trajectory file diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 99f63c5af6..1f0c97516a 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -331,8 +331,6 @@ int cvm::atom_group::parse(std::string const &conf, } else { get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent); } - - get_keyval(group_conf, "weights", weights, weights, colvarparse::parse_silent); } // FITTING OPTIONS @@ -380,12 +378,12 @@ int cvm::atom_group::parse(std::string const &conf, } std::string ref_pos_col; - double ref_pos_col_value; + double ref_pos_col_value=0.0; if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""), mode)) { // 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) + if (found && ref_pos_col_value == 0.0) cvm::error("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); } else { diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index d950ffb7ac..ef3bbcb39a 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -139,10 +139,6 @@ public: /// Allocates and populates the sorted list of atom ids int create_sorted_ids(void); - /// List of user-defined weights to be used by certain CVCs - std::vector weights; - - /// \brief When updating atomic coordinates, translate them to align with the /// center of mass of the reference coordinates bool b_center; diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index 03674cc1fd..e4cc89f463 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -39,15 +39,12 @@ public: virtual int replica_share(); /// Perform analysis tasks - virtual inline void analyse() {} + virtual void analyze() {} /// Send forces to the collective variables void communicate_forces(); /// \brief Constructor - /// - /// The constructor of the colvarbias base class is protected, so - /// that it can only be called from inherited classes colvarbias(std::string const &conf, char const *key); /// Default constructor @@ -65,9 +62,18 @@ public: /// Write a label to the trajectory file (comment line) virtual std::ostream & write_traj_label(std::ostream &os); + /// (Re)initialize the output files (does not write them yet) + virtual int setup_output() { return COLVARS_OK; } + /// Output quantities such as the bias energy to the trajectory file virtual std::ostream & write_traj(std::ostream &os); + /// Write output files (if defined, e.g. in analysis mode) + virtual int write_output_files() + { + return COLVARS_OK; + } + inline cvm::real get_energy() { return bias_energy; } diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 10490152d7..2d9202f4b3 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -1,9 +1,5 @@ /// -*- c++ -*- -/******************************************************************************** - * Implementation of the ABF and histogram biases * - ********************************************************************************/ - #include "colvarmodule.h" #include "colvar.h" #include "colvarbias_abf.h" @@ -47,7 +43,7 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) // shared ABF get_keyval(conf, "shared", shared_on, false); if (shared_on) { - if (!cvm::replica_enabled || cvm::replica_num() <= 1) + if (!cvm::replica_enabled() || cvm::replica_num() <= 1) cvm::error("Error: shared ABF requires more than one replica."); else cvm::log("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n"); @@ -552,202 +548,3 @@ std::istream & colvarbias_abf::read_restart(std::istream& is) } return is; } - - - - -/// Histogram "bias" constructor - -colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *key) - : colvarbias(conf, key), - grid(NULL), out_name("") -{ - get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq); - /// with VMD, this may not be an error - // if ( output_freq == 0 ) { - // cvm::error("User required histogram with zero output frequency"); - // } - - { - colvar_array_size = 1; - bool colvar_array = false; - if (get_keyval(conf, "sumVectorColvars", colvar_array, colvar_array)) { - size_t i; - for (i = 0; i < colvars.size(); i++) { - if (colvars[i]->value().type() == colvarvalue::type_vector) { - if (colvar_array_size == 1) { - colvar_array_size = colvars[i]->value().size(); - } else { - if (colvar_array_size != colvars[i]->value().size()) { - cvm::error("Error: trying to combine vector colvars of different lengths.\n", INPUT_ERROR); - } - } - } - } - } - } - - grid = new colvar_grid_count(); - { - std::string grid_conf; - if (key_lookup(conf, "grid", grid_conf)) { - grid->parse_params(grid_conf); - } else { - grid->init_from_colvars(colvars); - } - } - - bin.assign(colvars.size(), 0); - - cvm::log("Finished histogram setup.\n"); -} - -/// Destructor -colvarbias_histogram::~colvarbias_histogram() -{ - if (grid_os.is_open()) - grid_os.close(); - - if (grid) { - delete grid; - grid = NULL; - } - - if (cvm::n_histo_biases > 0) - cvm::n_histo_biases -= 1; -} - -/// Update the grid -cvm::real colvarbias_histogram::update() -{ - - if (cvm::debug()) { - cvm::log("Updating histogram bias " + this->name); - } - - // At the first timestep, we need to assign out_name since - // output_prefix is unset during the constructor - - if (cvm::step_relative() == 0) { - out_name = cvm::output_prefix + "." + this->name + ".dat"; - cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\""); - } - - - bin.assign(colvars.size(), 0); - - { - // update indices for all scalar values - size_t i; - for (i = 0; i < colvars.size(); i++) { - if (colvars[i]->value().type() == colvarvalue::type_scalar) { - bin[i] = grid->value_to_bin_scalar(colvars[i]->value(), i); - } - } - } - - if (colvar_array_size > 1) { - // update indices for all vector/array values - size_t iv, i; - for (iv = 0; iv < colvar_array_size; iv++) { - for (i = 0; i < colvars.size(); i++) { - if (colvars[i]->value().type() == colvarvalue::type_vector) { - bin[i] = grid->value_to_bin_scalar(colvars[i]->value().vector1d_value[iv], i); - } - } - if (grid->index_ok(bin)) { - // Only within bounds of the grid... - grid->incr_count(bin); - } - } - } else { - if (grid->index_ok(bin)) { - // Only within bounds of the grid... - grid->incr_count(bin); - } - } - - if (output_freq && (cvm::step_absolute() % output_freq) == 0) { - if (cvm::debug()) cvm::log("Histogram bias trying to write grid to disk"); - - grid_os.open(out_name.c_str()); - if (!grid_os.is_open()) cvm::error("Error opening histogram file " + out_name + " for writing"); - grid->write_multicol(grid_os); - grid_os.close(); - } - - return 0.0; // no bias energy for histogram -} - - -std::istream & colvarbias_histogram::read_restart(std::istream& is) -{ - size_t const start_pos = is.tellg(); - - cvm::log("Restarting collective variable histogram \""+ - this->name+"\".\n"); - std::string key, brace, conf; - - if ( !(is >> key) || !(key == "histogram") || - !(is >> brace) || !(brace == "{") || - !(is >> colvarparse::read_block("configuration", conf)) ) { - cvm::log("Error: in reading restart configuration for histogram \""+ - 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, "name", name, std::string(""), colvarparse::parse_silent)) && - (name != this->name) ) - cvm::error("Error: in the restart file, the " - "\"histogram\" block has a wrong name: different system?\n"); - if ( (id == -1) && (name == "") ) { - cvm::error("Error: \"histogram\" block in the restart file " - "has no name.\n"); - } - - if ( !(is >> key) || !(key == "grid")) { - cvm::error("Error: in reading restart configuration for histogram \""+ - 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; - } - if (! grid->read_raw(is)) { - is.clear(); - is.seekg(start_pos, std::ios::beg); - is.setstate(std::ios::failbit); - return is; - } - - is >> brace; - if (brace != "}") { - 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; -} - -std::ostream & colvarbias_histogram::write_restart(std::ostream& os) -{ - os << "histogram {\n" - << " configuration {\n" - << " name " << this->name << "\n"; - os << " }\n"; - - os << "grid\n"; - grid->write_raw(os, 8); - - os << "}\n\n"; - - return os; -} diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index c1e3a1fb31..eb02d1f55e 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -1,7 +1,4 @@ // -*- c++ -*- -/************************************************************************ - * Headers for the ABF and histogram biases * - ************************************************************************/ #ifndef COLVARBIAS_ABF_H #define COLVARBIAS_ABF_H @@ -40,7 +37,7 @@ private: bool hide_Jacobian; size_t full_samples; size_t min_samples; - /// frequency for updating output files (default: same as restartFreq?) + /// frequency for updating output files int output_freq; /// Write combined files with a history of all output data? bool b_history_files; @@ -90,34 +87,5 @@ private: std::ostream& write_restart(std::ostream&); }; - -/// Histogram "bias" (does as the name says) -class colvarbias_histogram : public colvarbias { - -public: - - colvarbias_histogram(std::string const &conf, char const *key); - ~colvarbias_histogram(); - - cvm::real update(); - -protected: - - /// n-dim histogram - colvar_grid_count *grid; - std::vector bin; - std::string out_name; - - int output_freq; - - /// If one or more of the variables are \link type_vector \endlink, treat them as arrays of this length - size_t colvar_array_size; - - void write_grid(); - cvm::ofstream grid_os; /// Stream for writing grid to disk - - std::istream& read_restart(std::istream&); - std::ostream& write_restart(std::ostream&); -}; - #endif + diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 0276359f0d..1fa804d303 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -44,7 +44,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) get_keyval(conf, "hillWeight", hill_weight, 0.01); if (hill_weight == 0.0) cvm::log("Warning: hillWeight has been set to zero, " - "this bias will have no effect.\n"); + "this bias will have no effect.\n"); get_keyval(conf, "newHillFrequency", new_hill_freq, 1000); @@ -71,9 +71,9 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) if (colvars[i]->expand_boundaries) { expand_grids = true; cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": Will expand grids when the colvar \""+ - colvars[i]->name+"\" approaches its boundaries.\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": Will expand grids when the colvar \""+ + colvars[i]->name+"\" approaches its boundaries.\n"); } } @@ -103,9 +103,9 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) if (expand_grids) cvm::fatal_error("Error: expandBoundaries is not supported when " - "using more than one replicas; please allocate " - "wide enough boundaries for each colvar" - "ahead of time.\n"); + "using more than one replicas; please allocate " + "wide enough boundaries for each colvar" + "ahead of time.\n"); if (get_keyval(conf, "dumpPartialFreeEnergyFile", dump_replica_fes, false)) { if (dump_replica_fes && (! dump_fes)) { @@ -115,118 +115,26 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) get_keyval(conf, "replicaID", replica_id, std::string("")); if (!replica_id.size()) - cvm::fatal_error("Error: replicaID must be defined " - "when using more than one replica.\n"); + cvm::error("Error: replicaID must be defined " + "when using more than one replica.\n", INPUT_ERROR); get_keyval(conf, "replicasRegistry", - replicas_registry_file, - (this->name+".replicas.registry.txt")); + replicas_registry_file, + (this->name+".replicas.registry.txt")); get_keyval(conf, "replicaUpdateFrequency", - replica_update_freq, new_hill_freq); + replica_update_freq, new_hill_freq); if (keep_hills) cvm::log("Warning: in metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": keepHills with more than one replica can lead to a very " - "large amount input/output and slow down your calculations. " - "Please consider disabling it.\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": keepHills with more than one replica can lead to a very " + "large amount of input/output and slow down your calculations. " + "Please consider disabling it.\n"); - - { - // TODO: one may want to specify the path manually for intricated filesystems? - char *pwd = new char[3001]; - if (GETCWD(pwd, 3000) == NULL) - cvm::fatal_error("Error: cannot get the path of the current working directory.\n"); - replica_list_file = - (std::string(pwd)+std::string(PATHSEP)+ - this->name+"."+replica_id+".files.txt"); - // replica_hills_file and replica_state_file are those written - // by the current replica; within the mirror biases, they are - // those by another replica - replica_hills_file = - (std::string(pwd)+std::string(PATHSEP)+ - cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills"); - replica_state_file = - (std::string(pwd)+std::string(PATHSEP)+ - cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state"); - delete[] pwd; - } - - // now register this replica - - // first check that it isn't already there - bool registered_replica = false; - std::ifstream reg_is(replicas_registry_file.c_str()); - if (reg_is.good()) { // the file may not be there yet - std::string existing_replica(""); - std::string existing_replica_file(""); - while ((reg_is >> existing_replica) && existing_replica.size() && - (reg_is >> existing_replica_file) && existing_replica_file.size()) { - if (existing_replica == replica_id) { - // this replica was already registered - replica_list_file = existing_replica_file; - reg_is.close(); - registered_replica = true; - break; - } - } - reg_is.close(); - } - - // if this replica was not included yet, we should generate a - // new record for it: but first, we write this replica's files, - // for the others to read - - // open the "hills" buffer file - replica_hills_os.open(replica_hills_file.c_str()); - if (!replica_hills_os.good()) - cvm::fatal_error("Error: in opening file \""+ - replica_hills_file+"\" for writing.\n"); - replica_hills_os.setf(std::ios::scientific, std::ios::floatfield); - - // write the state file (so that there is always one available) - write_replica_state_file(); - // schedule to read the state files of the other replicas - for (size_t ir = 0; ir < replicas.size(); ir++) { - (replicas[ir])->replica_state_file_in_sync = false; - } - - // if we're running without grids, use a growing list of "hills" files - // otherwise, just one state file and one "hills" file as buffer - std::ofstream list_os(replica_list_file.c_str(), - (use_grids ? std::ios::trunc : std::ios::app)); - if (! list_os.good()) - cvm::fatal_error("Error: in opening file \""+ - replica_list_file+"\" for writing.\n"); - list_os << "stateFile " << replica_state_file << "\n"; - list_os << "hillsFile " << replica_hills_file << "\n"; - list_os.close(); - - // finally, if add a new record for this replica to the registry - if (! registered_replica) { - std::ofstream reg_os(replicas_registry_file.c_str(), std::ios::app); - if (! reg_os.good()) - cvm::fatal_error("Error: in opening file \""+ - replicas_registry_file+"\" for writing.\n"); - reg_os << replica_id << " " << replica_list_file << "\n"; - reg_os.close(); - } } get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false); - if (b_hills_traj) { - std::string const traj_file_name(cvm::output_prefix+ - ".colvars."+this->name+ - ( (comm != single_replica) ? - ("."+replica_id) : - ("") )+ - ".hills.traj"); - hills_traj_os.open(traj_file_name.c_str()); - if (!hills_traj_os.good()) - cvm::fatal_error("Error: in opening hills output file \"" + - traj_file_name + "\".\n"); - } // for well-tempered metadynamics get_keyval(conf, "wellTempered", well_tempered, false); @@ -241,7 +149,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) if (cvm::debug()) cvm::log("Done initializing the metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); save_delimiters = false; } @@ -259,10 +167,10 @@ colvarbias_meta::~colvarbias_meta() hills_energy_gradients = NULL; } - if (replica_hills_os.good()) + if (replica_hills_os.is_open()) replica_hills_os.close(); - if (hills_traj_os.good()) + if (hills_traj_os.is_open()) hills_traj_os.close(); if (cvm::n_meta_biases > 0) @@ -298,7 +206,7 @@ colvarbias_meta::create_hill(colvarbias_meta::hill const &h) } // output to trajectory (if specified) - if (hills_traj_os.good()) { + if (hills_traj_os.is_open()) { hills_traj_os << (hills.back()).output_traj(); if (cvm::debug()) { hills_traj_os.flush(); @@ -315,11 +223,11 @@ colvarbias_meta::delete_hill(hill_iter &h) { if (cvm::debug()) { cvm::log("Deleting hill from the metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ", with step number "+ - cvm::to_str(h->it)+(h->replica.size() ? - ", replica id \""+h->replica : - "")+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ", with step number "+ + cvm::to_str(h->it)+(h->replica.size() ? + ", replica id \""+h->replica : + "")+".\n"); } if (use_grids && !hills_off_grid.empty()) { @@ -332,7 +240,7 @@ colvarbias_meta::delete_hill(hill_iter &h) } } - if (hills_traj_os.good()) { + if (hills_traj_os.is_open()) { // output to the trajectory hills_traj_os << "# DELETED this hill: " << (hills.back()).output_traj() @@ -349,7 +257,7 @@ cvm::real colvarbias_meta::update() { if (cvm::debug()) cvm::log("Updating the metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); if (use_grids) { @@ -360,9 +268,9 @@ cvm::real colvarbias_meta::update() // first of all, expand the grids, if specified if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": current coordinates on the grid: "+ - cvm::to_str(curr_bin)+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": current coordinates on the grid: "+ + cvm::to_str(curr_bin)+".\n"); bool changed_grids = false; int const min_buffer = @@ -393,10 +301,10 @@ cvm::real colvarbias_meta::update() changed_lb = true; cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": new lower boundary for colvar \""+ - colvars[i]->name+"\", at "+ - cvm::to_str(new_lower_boundaries[i])+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": new lower boundary for colvar \""+ + colvars[i]->name+"\", at "+ + cvm::to_str(new_lower_boundaries[i])+".\n"); } if (!colvars[i]->hard_upper_boundary) @@ -407,10 +315,10 @@ cvm::real colvarbias_meta::update() changed_ub = true; cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": new upper boundary for colvar \""+ - colvars[i]->name+"\", at "+ - cvm::to_str(new_upper_boundaries[i])+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": new upper boundary for colvar \""+ + colvars[i]->name+"\", at "+ + cvm::to_str(new_upper_boundaries[i])+".\n"); } if (changed_lb || changed_ub) @@ -447,7 +355,7 @@ cvm::real colvarbias_meta::update() curr_bin = hills_energy->get_colvars_index(); if (cvm::debug()) cvm::log("Coordinates on the new grid: "+ - cvm::to_str(curr_bin)+".\n"); + cvm::to_str(curr_bin)+".\n"); } } } @@ -457,8 +365,8 @@ cvm::real colvarbias_meta::update() if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": adding a new hill at step "+cvm::to_str(cvm::step_absolute())+".\n"); switch (comm) { @@ -492,12 +400,12 @@ cvm::real colvarbias_meta::update() } else { create_hill(hill(hill_weight, colvars, hill_width, replica_id)); } - if (replica_hills_os.good()) { + if (replica_hills_os.is_open()) { replica_hills_os << hills.back(); } else { cvm::fatal_error("Error: in metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - " while writing hills for the other replicas.\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + " while writing hills for the other replicas.\n"); } break; } @@ -510,7 +418,8 @@ cvm::real colvarbias_meta::update() if ((cvm::step_absolute() % replica_update_freq) == 0) { update_replicas_registry(); // empty the output buffer - replica_hills_os.flush(); + if (replica_hills_os.is_open()) + replica_hills_os.flush(); read_replica_files(); } @@ -536,7 +445,7 @@ cvm::real colvarbias_meta::update() if ((cvm::step_absolute() % grids_freq) == 0) { // map the most recent gaussians to the grids project_hills(new_hills_begin, hills.end(), - hills_energy, hills_energy_gradients); + hills_energy, hills_energy_gradients); new_hills_begin = hills.end(); // TODO: we may want to condense all into one replicas array, @@ -544,9 +453,9 @@ cvm::real colvarbias_meta::update() if (comm == multiple_replicas) { for (size_t ir = 0; ir < replicas.size(); ir++) { replicas[ir]->project_hills(replicas[ir]->new_hills_begin, - replicas[ir]->hills.end(), - replicas[ir]->hills_energy, - replicas[ir]->hills_energy_gradients); + replicas[ir]->hills.end(), + replicas[ir]->hills_energy, + replicas[ir]->hills_energy_gradients); replicas[ir]->new_hills_begin = replicas[ir]->hills.end(); } } @@ -555,9 +464,9 @@ cvm::real colvarbias_meta::update() std::vector curr_bin = hills_energy->get_colvars_index(); if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": current coordinates on the grid: "+ - cvm::to_str(curr_bin)+".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": current coordinates on the grid: "+ + cvm::to_str(curr_bin)+".\n"); if (hills_energy->index_ok(curr_bin)) { @@ -590,13 +499,13 @@ cvm::real colvarbias_meta::update() if (comm == multiple_replicas) for (size_t ir = 0; ir < replicas.size(); ir++) { calc_hills(replicas[ir]->hills_off_grid.begin(), - replicas[ir]->hills_off_grid.end(), - bias_energy); + replicas[ir]->hills_off_grid.end(), + bias_energy); for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills_force(ic, - replicas[ir]->hills_off_grid.begin(), - replicas[ir]->hills_off_grid.end(), - colvar_forces); + replicas[ir]->hills_off_grid.begin(), + replicas[ir]->hills_off_grid.end(), + colvar_forces); } } } @@ -612,27 +521,27 @@ cvm::real colvarbias_meta::update() if (cvm::debug()) cvm::log("Hills energy = "+cvm::to_str(bias_energy)+ - ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); + ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": adding the forces from the other replicas.\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": adding the forces from the other replicas.\n"); if (comm == multiple_replicas) for (size_t ir = 0; ir < replicas.size(); ir++) { calc_hills(replicas[ir]->new_hills_begin, - replicas[ir]->hills.end(), - bias_energy); + replicas[ir]->hills.end(), + bias_energy); for (size_t ic = 0; ic < colvars.size(); ic++) { calc_hills_force(ic, - replicas[ir]->new_hills_begin, - replicas[ir]->hills.end(), - colvar_forces); + replicas[ir]->new_hills_begin, + replicas[ir]->hills.end(), + colvar_forces); } if (cvm::debug()) cvm::log("Hills energy = "+cvm::to_str(bias_energy)+ - ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); + ", hills forces = "+cvm::to_str(colvar_forces)+".\n"); } return bias_energy; @@ -640,9 +549,9 @@ cvm::real colvarbias_meta::update() void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, - colvarbias_meta::hill_iter h_last, - cvm::real &energy, - std::vector const &colvar_values) + colvarbias_meta::hill_iter h_last, + cvm::real &energy, + std::vector const &colvar_values) { std::vector curr_values(colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { @@ -683,10 +592,10 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, void colvarbias_meta::calc_hills_force(size_t const &i, - colvarbias_meta::hill_iter h_first, - colvarbias_meta::hill_iter h_last, - std::vector &forces, - std::vector const &values) + colvarbias_meta::hill_iter h_first, + colvarbias_meta::hill_iter h_last, + std::vector &forces, + std::vector const &values) { // Retrieve the value of the colvar colvarvalue const x(values.size() ? values[i] : colvars[i]->value()); @@ -758,15 +667,15 @@ void colvarbias_meta::calc_hills_force(size_t const &i, // ********************************************************************** void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, - colvarbias_meta::hill_iter h_last, - colvar_grid_scalar *he, - colvar_grid_gradient *hg, - bool print_progress) + colvarbias_meta::hill_iter h_last, + colvar_grid_scalar *he, + colvar_grid_gradient *hg, + bool print_progress) { if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": projecting hills.\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": projecting hills.\n"); // TODO: improve it by looping over a small subgrid instead of the whole grid @@ -862,8 +771,8 @@ void colvarbias_meta::project_hills(colvarbias_meta::hill_iter h_first, void colvarbias_meta::recount_hills_off_grid(colvarbias_meta::hill_iter h_first, - colvarbias_meta::hill_iter h_last, - colvar_grid_scalar *he) + colvarbias_meta::hill_iter h_last, + colvar_grid_scalar *he) { hills_off_grid.clear(); @@ -886,20 +795,20 @@ void colvarbias_meta::update_replicas_registry() { if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ": updating the list of replicas, currently containing "+ - cvm::to_str(replicas.size())+" elements.\n"); + ": updating the list of replicas, currently containing "+ + cvm::to_str(replicas.size())+" elements.\n"); { // copy the whole file into a string for convenience std::string line(""); std::ifstream reg_file(replicas_registry_file.c_str()); - if (reg_file.good()) { + if (reg_file.is_open()) { replicas_registry.clear(); while (colvarparse::getline_nocomments(reg_file, line)) replicas_registry.append(line+"\n"); } else { - cvm::fatal_error("Error: failed to open file \""+replicas_registry_file+ - "\" for reading.\n"); + cvm::error("Error: failed to open file \""+replicas_registry_file+ + "\" for reading.\n", FILE_ERROR); } } @@ -916,9 +825,9 @@ void colvarbias_meta::update_replicas_registry() // this is the record for this same replica, skip it if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": skipping this replica's own record: \""+ - new_replica+"\", \""+new_replica_file+"\"\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": skipping this replica's own record: \""+ + new_replica+"\", \""+new_replica_file+"\"\n"); new_replica_file.clear(); new_replica.clear(); continue; @@ -930,9 +839,9 @@ void colvarbias_meta::update_replicas_registry() // this replica was already added if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": skipping a replica already loaded, \""+ - (replicas[ir])->replica_id+"\".\n"); + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": skipping a replica already loaded, \""+ + (replicas[ir])->replica_id+"\".\n"); already_loaded = true; break; } @@ -941,7 +850,7 @@ void colvarbias_meta::update_replicas_registry() if (!already_loaded) { // add this replica to the registry cvm::log("Metadynamics bias \""+this->name+"\""+ - ": accessing replica \""+new_replica+"\".\n"); + ": accessing replica \""+new_replica+"\".\n"); replicas.push_back(new colvarbias_meta()); (replicas.back())->replica_id = new_replica; (replicas.back())->replica_list_file = new_replica_file; @@ -968,15 +877,15 @@ void colvarbias_meta::update_replicas_registry() } } else { cvm::fatal_error("Error: cannot read the replicas registry file \""+ - replicas_registry+"\".\n"); + replicas_registry+"\".\n"); } // now (re)read the list file of each replica for (size_t ir = 0; ir < replicas.size(); ir++) { if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ": reading the list file for replica \""+ - (replicas[ir])->replica_id+"\".\n"); + ": reading the list file for replica \""+ + (replicas[ir])->replica_id+"\".\n"); std::ifstream list_is((replicas[ir])->replica_list_file.c_str()); std::string key; @@ -988,17 +897,17 @@ void colvarbias_meta::update_replicas_registry() !(key == std::string("hillsFile")) || !(list_is >> new_hills_file)) { cvm::log("Metadynamics bias \""+this->name+"\""+ - ": failed to read the file \""+ - (replicas[ir])->replica_list_file+"\": will try again after "+ - cvm::to_str(replica_update_freq)+" steps.\n"); + ": failed to read the file \""+ + (replicas[ir])->replica_list_file+"\": will try again after "+ + cvm::to_str(replica_update_freq)+" steps.\n"); (replicas[ir])->update_status++; } else { (replicas[ir])->update_status = 0; if (new_state_file != (replicas[ir])->replica_state_file) { cvm::log("Metadynamics bias \""+this->name+"\""+ - ": replica \""+(replicas[ir])->replica_id+ - "\" has supplied a new state file, \""+new_state_file+ - "\".\n"); + ": replica \""+(replicas[ir])->replica_id+ + "\" has supplied a new state file, \""+new_state_file+ + "\".\n"); (replicas[ir])->replica_state_file_in_sync = false; (replicas[ir])->replica_state_file = new_state_file; (replicas[ir])->replica_hills_file = new_hills_file; @@ -1009,7 +918,7 @@ void colvarbias_meta::update_replicas_registry() if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\": the list of replicas contains "+ - cvm::to_str(replicas.size())+" elements.\n"); + cvm::to_str(replicas.size())+" elements.\n"); } @@ -1027,15 +936,15 @@ void colvarbias_meta::read_replica_files() (! (replicas[ir])->replica_state_file_in_sync) ) { cvm::log("Metadynamics bias \""+this->name+"\""+ - ": reading the state of replica \""+ - (replicas[ir])->replica_id+"\" from file \""+ - (replicas[ir])->replica_state_file+"\".\n"); + ": reading the state of replica \""+ + (replicas[ir])->replica_id+"\" from file \""+ + (replicas[ir])->replica_state_file+"\".\n"); std::ifstream is((replicas[ir])->replica_state_file.c_str()); if (! (replicas[ir])->read_restart(is)) { cvm::log("Reading from file \""+(replicas[ir])->replica_state_file+ - "\" failed or incomplete: will try again in "+ - cvm::to_str(replica_update_freq)+" steps.\n"); + "\" failed or incomplete: will try again in "+ + cvm::to_str(replica_update_freq)+" steps.\n"); } else { // state file has been read successfully (replicas[ir])->replica_state_file_in_sync = true; @@ -1049,19 +958,19 @@ void colvarbias_meta::read_replica_files() if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ": checking for new hills from replica \""+ - (replicas[ir])->replica_id+"\" in the file \""+ - (replicas[ir])->replica_hills_file+"\".\n"); + ": checking for new hills from replica \""+ + (replicas[ir])->replica_id+"\" in the file \""+ + (replicas[ir])->replica_hills_file+"\".\n"); // read hills from the other replicas' files; for each file, resume // the position recorded previously std::ifstream is((replicas[ir])->replica_hills_file.c_str()); - if (is.good()) { + if (is.is_open()) { // try to resume the previous position is.seekg((replicas[ir])->replica_hills_file_pos, std::ios::beg); - if (!is.good()){ + if (!is.is_open()){ // if fail (the file may have been overwritten), reset this // position is.clear(); @@ -1073,26 +982,26 @@ void colvarbias_meta::read_replica_files() // and record the failure (replicas[ir])->update_status++; cvm::log("Failed to read the file \""+(replicas[ir])->replica_hills_file+ - "\" at the previous position: will try again in "+ - cvm::to_str(replica_update_freq)+" steps.\n"); + "\" at the previous position: will try again in "+ + cvm::to_str(replica_update_freq)+" steps.\n"); } else { while ((replicas[ir])->read_hill(is)) { // if (cvm::debug()) - cvm::log("Metadynamics bias \""+this->name+"\""+ - ": received a hill from replica \""+ - (replicas[ir])->replica_id+ - "\" at step "+ - cvm::to_str(((replicas[ir])->hills.back()).it)+".\n"); + cvm::log("Metadynamics bias \""+this->name+"\""+ + ": received a hill from replica \""+ + (replicas[ir])->replica_id+ + "\" at step "+ + cvm::to_str(((replicas[ir])->hills.back()).it)+".\n"); } is.clear(); // store the position for the next read (replicas[ir])->replica_hills_file_pos = is.tellg(); if (cvm::debug()) cvm::log("Metadynamics bias \""+this->name+"\""+ - ": stopped reading file \""+(replicas[ir])->replica_hills_file+ - "\" at position "+ - cvm::to_str((replicas[ir])->replica_hills_file_pos)+".\n"); + ": stopped reading file \""+(replicas[ir])->replica_hills_file+ + "\" at position "+ + cvm::to_str((replicas[ir])->replica_hills_file_pos)+".\n"); // test whether this is the end of the file is.seekg(0, std::ios::end); @@ -1105,8 +1014,8 @@ void colvarbias_meta::read_replica_files() } else { cvm::log("Failed to read the file \""+(replicas[ir])->replica_hills_file+ - "\": will try again in "+ - cvm::to_str(replica_update_freq)+" steps.\n"); + "\": will try again in "+ + cvm::to_str(replica_update_freq)+" steps.\n"); (replicas[ir])->update_status++; // cvm::fatal_error ("Error: cannot read from file \""+ // (replicas[ir])->replica_hills_file+"\".\n"); @@ -1118,11 +1027,11 @@ void colvarbias_meta::read_replica_files() if ((replicas[ir])->update_status > 3*n_flush) { // TODO: suspend the calculation? cvm::log("WARNING: in metadynamics bias \""+this->name+"\""+ - " failed to read completely the output of replica \""+ - (replicas[ir])->replica_id+ - "\" after more than "+ - cvm::to_str((replicas[ir])->update_status * replica_update_freq)+ - " steps. Ensure that it is still running.\n"); + " failed to read completely the output of replica \""+ + (replicas[ir])->replica_id+ + "\" after more than "+ + cvm::to_str((replicas[ir])->update_status * replica_update_freq)+ + " steps. Ensure that it is still running.\n"); } } } @@ -1141,7 +1050,7 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) // if using a multiple replicas scheme, output messages // are printed before and after calling this function cvm::log("Restarting metadynamics bias \""+this->name+"\""+ - ".\n"); + ".\n"); } std::string key, brace, conf; if ( !(is >> key) || !(key == "metadynamics") || @@ -1150,11 +1059,11 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) if (comm == single_replica) cvm::log("Error: in reading restart configuration for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - (replica_state_file_in_sync ? ("at position "+ - cvm::to_str(start_pos)+ - " in the state file") : "")+".\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + (replica_state_file_in_sync ? ("at position "+ + cvm::to_str(start_pos)+ + " in the state file") : "")+".\n"); is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); @@ -1163,26 +1072,26 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) std::string name = ""; if ( colvarparse::get_keyval(conf, "name", name, - std::string(""), colvarparse::parse_silent) && + std::string(""), colvarparse::parse_silent) && (name != this->name) ) cvm::fatal_error("Error: in the restart file, the " - "\"metadynamics\" block has a different name: different system?\n"); + "\"metadynamics\" block has a different name: different system?\n"); if (name.size() == 0) { cvm::fatal_error("Error: \"metadynamics\" block within the restart file " - "has no identifiers.\n"); + "has no identifiers.\n"); } if (comm != single_replica) { std::string replica = ""; if (colvarparse::get_keyval(conf, "replicaID", replica, - std::string(""), colvarparse::parse_silent) && + std::string(""), colvarparse::parse_silent) && (replica != this->replica_id)) cvm::fatal_error("Error: in the restart file, the " - "\"metadynamics\" block has a different replicaID: different system?\n"); + "\"metadynamics\" block has a different replicaID: different system?\n"); colvarparse::get_keyval(conf, "step", state_file_step, - cvm::step_absolute(), colvarparse::parse_silent); + cvm::step_absolute(), colvarparse::parse_silent); } bool grids_from_restart_file = use_grids; @@ -1206,8 +1115,8 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) if (has_data) { if (cvm::debug()) cvm::log("Backupping grids for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); hills_energy_backup = hills_energy; hills_energy_gradients_backup = hills_energy_gradients; hills_energy = new colvar_grid_scalar(colvars); @@ -1234,14 +1143,14 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) if (!rebin_grids) { if (hills_energy_backup == NULL) cvm::fatal_error("Error: couldn't read the free energy grid for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - "; if useGrids was off when the state file was written, " - "enable rebinGrids now to regenerate the grids.\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + "; if useGrids was off when the state file was written, " + "enable rebinGrids now to regenerate the grids.\n"); else { if (comm == single_replica) cvm::log("Error: couldn't read the free energy grid for metadynamics bias \""+ - this->name+"\".\n"); + this->name+"\".\n"); delete hills_energy; delete hills_energy_gradients; hills_energy = hills_energy_backup; @@ -1272,14 +1181,14 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) if (!rebin_grids) { if (hills_energy_backup == NULL) cvm::fatal_error("Error: couldn't read the free energy gradients grid for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - "; if useGrids was off when the state file was written, " - "enable rebinGrids now to regenerate the grids.\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + "; if useGrids was off when the state file was written, " + "enable rebinGrids now to regenerate the grids.\n"); else { if (comm == single_replica) cvm::log("Error: couldn't read the free energy gradients grid for metadynamics bias \""+ - this->name+"\".\n"); + this->name+"\".\n"); delete hills_energy; delete hills_energy_gradients; hills_energy = hills_energy_backup; @@ -1292,8 +1201,8 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) if (cvm::debug()) cvm::log("Successfully read new grids for bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); if (hills_energy_backup != NULL) { // now that we have successfully updated the grids, delete the @@ -1315,16 +1224,16 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) while (read_hill(is)) { if (cvm::debug()) cvm::log("Read a previously saved hill under the " - "metadynamics bias \""+ - this->name+"\", created at step "+ - cvm::to_str((hills.back()).it)+".\n"); + "metadynamics bias \""+ + this->name+"\", created at step "+ + cvm::to_str((hills.back()).it)+".\n"); } is.clear(); new_hills_begin = hills.end(); if (grids_from_restart_file) { if (hills.size() > old_hills_size) cvm::log("Read "+cvm::to_str(hills.size())+ - " hills in addition to the grids.\n"); + " hills in addition to the grids.\n"); } else { if (!hills.empty()) cvm::log("Read "+cvm::to_str(hills.size())+" hills.\n"); @@ -1344,15 +1253,15 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) 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"); + cvm::to_str(hills.size())+" hills (this may take a while)...\n"); project_hills(hills.begin(), hills.end(), - new_hills_energy, new_hills_energy_gradients, true); + new_hills_energy, new_hills_energy_gradients, true); cvm::log("rebinning done.\n"); } else { // otherwise, use the grids in the restart file cvm::log("Rebinning the energy and forces grids " - "from the grids in the restart file.\n"); + "from the grids in the restart file.\n"); new_hills_energy->map_grid(*hills_energy); new_hills_energy_gradients->map_grid(*hills_energy_gradients); } @@ -1371,18 +1280,18 @@ std::istream & colvarbias_meta::read_restart(std::istream& is) if (use_grids) { 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"); + "grid boundaries: they will be computed analytically " + "and saved to the state files.\n"); } } is >> brace; if (brace != "}") { cvm::log("Incomplete restart information for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ - ": no closing brace at position "+ - cvm::to_str(is.tellg())+" in the file.\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ + ": no closing brace at position "+ + cvm::to_str(is.tellg())+" in the file.\n"); is.setstate(std::ios::failbit); return is; } @@ -1424,8 +1333,8 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) if (h_it <= state_file_step) { if (cvm::debug()) cvm::log("Skipping a hill older than the state file for metadynamics bias \""+ - this->name+"\""+ - ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); + this->name+"\""+ + ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n"); return is; } @@ -1449,16 +1358,16 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) std::vector h_widths(colvars.size()); get_keyval(data, "widths", h_widths, - std::vector (colvars.size(), (std::sqrt(2.0 * PI) / 2.0)), - parse_silent); + std::vector (colvars.size(), (std::sqrt(2.0 * PI) / 2.0)), + parse_silent); std::string h_replica = ""; if (comm != single_replica) { get_keyval(data, "replicaID", h_replica, replica_id, parse_silent); if (h_replica != replica_id) cvm::fatal_error("Error: trying to read a hill created by replica \""+h_replica+ - "\" for replica \""+replica_id+ - "\"; did you swap output files?\n"); + "\" for replica \""+replica_id+ + "\"; did you swap output files?\n"); } hill_iter const hills_end = hills.end(); @@ -1490,6 +1399,111 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) // output functions // ********************************************************************** + +int colvarbias_meta::setup_output() +{ + + if (comm == multiple_replicas) { + + // TODO: one may want to specify the path manually for intricated filesystems? + char *pwd = new char[3001]; + if (GETCWD(pwd, 3000) == NULL) + cvm::fatal_error("Error: cannot get the path of the current working directory.\n"); + replica_list_file = + (std::string(pwd)+std::string(PATHSEP)+ + this->name+"."+replica_id+".files.txt"); + // replica_hills_file and replica_state_file are those written + // by the current replica; within the mirror biases, they are + // those by another replica + replica_hills_file = + (std::string(pwd)+std::string(PATHSEP)+ + cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills"); + replica_state_file = + (std::string(pwd)+std::string(PATHSEP)+ + cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state"); + delete[] pwd; + + // now register this replica + + // first check that it isn't already there + bool registered_replica = false; + std::ifstream reg_is(replicas_registry_file.c_str()); + if (reg_is.is_open()) { // the file may not be there yet + std::string existing_replica(""); + std::string existing_replica_file(""); + while ((reg_is >> existing_replica) && existing_replica.size() && + (reg_is >> existing_replica_file) && existing_replica_file.size()) { + if (existing_replica == replica_id) { + // this replica was already registered + replica_list_file = existing_replica_file; + reg_is.close(); + registered_replica = true; + break; + } + } + reg_is.close(); + } + + // if this replica was not included yet, we should generate a + // new record for it: but first, we write this replica's files, + // for the others to read + + // open the "hills" buffer file + if (!replica_hills_os.is_open()) { + cvm::backup_file(replica_hills_file.c_str()); + replica_hills_os.open(replica_hills_file.c_str()); + if (!replica_hills_os.is_open()) + cvm::error("Error: in opening file \""+ + replica_hills_file+"\" for writing.\n", FILE_ERROR); + replica_hills_os.setf(std::ios::scientific, std::ios::floatfield); + } + + // write the state file (so that there is always one available) + write_replica_state_file(); + // schedule to read the state files of the other replicas + for (size_t ir = 0; ir < replicas.size(); ir++) { + (replicas[ir])->replica_state_file_in_sync = false; + } + + // if we're running without grids, use a growing list of "hills" files + // otherwise, just one state file and one "hills" file as buffer + std::ofstream list_os(replica_list_file.c_str(), + (use_grids ? std::ios::trunc : std::ios::app)); + if (! list_os.is_open()) + cvm::fatal_error("Error: in opening file \""+ + replica_list_file+"\" for writing.\n"); + list_os << "stateFile " << replica_state_file << "\n"; + list_os << "hillsFile " << replica_hills_file << "\n"; + list_os.close(); + + // finally, if add a new record for this replica to the registry + if (! registered_replica) { + std::ofstream reg_os(replicas_registry_file.c_str(), std::ios::app); + if (! reg_os.is_open()) + cvm::error("Error: in opening file \""+ + replicas_registry_file+"\" for writing.\n", FILE_ERROR); + reg_os << replica_id << " " << replica_list_file << "\n"; + reg_os.close(); + } + } + + if (b_hills_traj) { + std::string const traj_file_name(cvm::output_prefix+ + ".colvars."+this->name+ + ( (comm != single_replica) ? + ("."+replica_id) : + ("") )+ + ".hills.traj"); + hills_traj_os.open(traj_file_name.c_str()); + if (!hills_traj_os.is_open()) + cvm::error("Error: in opening hills output file \"" + + traj_file_name+"\".\n", FILE_ERROR); + } + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + + std::ostream & colvarbias_meta::write_restart(std::ostream& os) { os << "metadynamics {\n" @@ -1505,7 +1519,7 @@ std::ostream & colvarbias_meta::write_restart(std::ostream& os) // this is a very good time to project hills, if you haven't done // it already! project_hills(new_hills_begin, hills.end(), - hills_energy, hills_energy_gradients); + hills_energy, hills_energy_gradients); new_hills_begin = hills.end(); // write down the grids to the restart file @@ -1579,10 +1593,10 @@ void colvarbias_meta::write_pmf() } { std::string const fes_file_name(fes_file_name_prefix + - ((comm != single_replica) ? ".partial" : "") + - (dump_fes_save ? - "."+cvm::to_str(cvm::step_absolute()) : "") + - ".pmf"); + ((comm != single_replica) ? ".partial" : "") + + (dump_fes_save ? + "."+cvm::to_str(cvm::step_absolute()) : "") + + ".pmf"); cvm::backup_file(fes_file_name.c_str()); cvm::ofstream fes_os(fes_file_name.c_str()); pmf->write_multicol(fes_os); @@ -1604,9 +1618,9 @@ void colvarbias_meta::write_pmf() pmf->multiply_constant(well_temper_scale); } std::string const fes_file_name(fes_file_name_prefix + - (dump_fes_save ? - "."+cvm::to_str(cvm::step_absolute()) : "") + - ".pmf"); + (dump_fes_save ? + "."+cvm::to_str(cvm::step_absolute()) : "") + + ".pmf"); cvm::backup_file(fes_file_name.c_str()); cvm::ofstream fes_os(fes_file_name.c_str()); pmf->write_multicol(fes_os); @@ -1624,9 +1638,9 @@ void colvarbias_meta::write_replica_state_file() // is duplicated code, that could be cleaned up later cvm::backup_file(replica_state_file.c_str()); cvm::ofstream rep_state_os(replica_state_file.c_str()); - if (!rep_state_os.good()) + if (!rep_state_os.is_open()) cvm::fatal_error("Error: in opening file \""+ - replica_state_file+"\" for writing.\n"); + replica_state_file+"\" for writing.\n"); rep_state_os.setf(std::ios::scientific, std::ios::floatfield); rep_state_os << "\n" @@ -1668,10 +1682,11 @@ void colvarbias_meta::write_replica_state_file() // reopen the hills file replica_hills_os.close(); + cvm::backup_file(replica_hills_file.c_str()); replica_hills_os.open(replica_hills_file.c_str()); - if (!replica_hills_os.good()) + if (!replica_hills_os.is_open()) cvm::fatal_error("Error: in opening file \""+ - replica_hills_file+"\" for writing.\n"); + replica_hills_file+"\" for writing.\n"); replica_hills_os.setf(std::ios::scientific, std::ios::floatfield); } diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index e2f175ad92..d4c521ae62 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -42,6 +42,8 @@ public: virtual std::ostream & write_restart(std::ostream &os); + virtual int setup_output(); + virtual void write_pmf(); class hill; diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index 91f37181d6..3148fe0138 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -110,7 +110,7 @@ protected: /// \brief Number of steps required to reach the target force constant /// or restraint centers - size_t target_nsteps; + long target_nsteps; }; /// \brief Harmonic bias restraint diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index f253aeffc6..ad8c9c605c 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -849,12 +849,12 @@ colvar::rmsd::rmsd(std::string const &conf) } std::string ref_pos_col; - double ref_pos_col_value; + double ref_pos_col_value=0.0; if (get_keyval(conf, "refPositionsCol", ref_pos_col, std::string(""))) { // 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) { + if (found && ref_pos_col_value==0.0) { cvm::error("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); return; @@ -1043,11 +1043,11 @@ colvar::eigenvector::eigenvector(std::string const &conf) } std::string file_col; - double file_col_value; + double file_col_value=0.0; 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) { + if (found && file_col_value==0.0) { cvm::error("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); return; @@ -1107,11 +1107,11 @@ colvar::eigenvector::eigenvector(std::string const &conf) } std::string file_col; - double file_col_value; + double file_col_value=0.0; 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) { + if (found && file_col_value==0.0) { cvm::error("Error: vectorColValue, if provided, must be non-zero.\n"); return; } @@ -1324,14 +1324,6 @@ void colvar::cartesian::calc_value() x.vector1d_value[dim*ia + j] = atoms[ia].pos[axes[j]]; } } - - if (atoms.weights.size()) { - for (ia = 0; ia < atoms.size(); ia++) { - for (j = 0; j < dim; j++) { - x.vector1d_value[dim*ia + j] *= atoms.weights[ia]; - } - } - } } @@ -1349,20 +1341,11 @@ void colvar::cartesian::apply_force(colvarvalue const &force) size_t ia, j; if (!atoms.noforce) { cvm::rvector f; - if (atoms.weights.size()) { - for (ia = 0; ia < atoms.size(); ia++) { - for (j = 0; j < dim; j++) { - f[axes[j]] = force.vector1d_value[dim*ia + j] / atoms.weights[ia]; - } - atoms[ia].apply_force(f); - } - } else { - for (ia = 0; ia < atoms.size(); ia++) { - for (j = 0; j < dim; j++) { - f[axes[j]] = force.vector1d_value[dim*ia + j]; - } - atoms[ia].apply_force(f); + for (ia = 0; ia < atoms.size(); ia++) { + for (j = 0; j < dim; j++) { + f[axes[j]] = force.vector1d_value[dim*ia + j]; } + atoms[ia].apply_force(f); } } } diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index a4c5f2187d..4bc0dcff0e 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -33,11 +33,11 @@ colvar::orientation::orientation(std::string const &conf) if (get_keyval(conf, "refPositionsFile", file_name)) { std::string file_col; - double file_col_value; + double file_col_value=0.0; 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) + if (found && file_col_value==0.0) cvm::fatal_error("Error: refPositionsColValue, " "if provided, must be non-zero.\n"); } else { diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index e91ed1e6cc..99a32ca8ff 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -8,9 +8,10 @@ #include "colvarproxy.h" #include "colvar.h" #include "colvarbias.h" -#include "colvarbias_alb.h" -#include "colvarbias_meta.h" #include "colvarbias_abf.h" +#include "colvarbias_alb.h" +#include "colvarbias_histogram.h" +#include "colvarbias_meta.h" #include "colvarbias_restraint.h" #include "colvarscript.h" @@ -523,13 +524,13 @@ int colvarmodule::calc() { cvm::log("Perform runtime analyses.\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->analyse(); + (*cvi)->analyze(); if (cvm::get_error()) { return COLVARS_ERROR; } } for (bi = biases.begin(); bi != biases.end(); bi++) { - (*bi)->analyse(); + (*bi)->analyze(); if (cvm::get_error()) { return COLVARS_ERROR; } @@ -627,7 +628,7 @@ int colvarmodule::analyze() cvi != colvars.end(); cvi++) { cvm::increase_depth(); - (*cvi)->analyse(); + (*cvi)->analyze(); cvm::decrease_depth(); } @@ -636,7 +637,7 @@ int colvarmodule::analyze() bi != biases.end(); bi++) { cvm::increase_depth(); - (*bi)->analyse(); + (*bi)->analyze(); cvm::decrease_depth(); } @@ -717,13 +718,15 @@ int colvarmodule::setup_input() cvm::log(cvm::line_marker); } } - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup_output() { + int error_code = 0; + // output state file (restart) restart_out_name = proxy->restart_output_prefix().size() ? std::string(proxy->restart_output_prefix()+".colvars.state") : @@ -748,7 +751,17 @@ int colvarmodule::setup_output() std::string("")); if (cv_traj_freq && cv_traj_name.size()) { - open_traj_file(cv_traj_name); + error_code |= open_traj_file(cv_traj_name); + } + + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + error_code |= (*bi)->setup_output(); + } + + if (error_code != COLVARS_OK || cvm::get_error()) { + set_error_bits(FILE_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -828,6 +841,14 @@ int colvarmodule::write_output_files() } cvm::decrease_depth(); + cvm::increase_depth(); + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + (*bi)->write_output_files(); + } + cvm::decrease_depth(); + if (cv_traj_os.is_open()) { // do not close to avoid problems with multiple NAMD runs cv_traj_os.flush(); @@ -839,8 +860,8 @@ int colvarmodule::write_output_files() int colvarmodule::read_traj(char const *traj_filename, - size_t traj_read_begin, - size_t traj_read_end) + long traj_read_begin, + long traj_read_end) { cvm::log("Opening trajectory file \""+ std::string(traj_filename)+"\".\n"); @@ -1204,8 +1225,8 @@ 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; +long colvarmodule::it = 0; +long colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; size_t colvarmodule::depth = 0; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 07134c2db5..25d2ce817b 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -4,7 +4,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2015-03-15" +#define COLVARS_VERSION "2015-04-22" #endif #ifndef COLVARS_DEBUG @@ -122,19 +122,19 @@ public: } /// Current step number - static size_t it; + static long it; /// Starting step number for this run - static size_t it_restart; + static long it_restart; /// Return the current step number from the beginning of this run - static inline size_t step_relative() + static inline long step_relative() { return it - it_restart; } /// Return the current step number from the beginning of the whole /// calculation - static inline size_t step_absolute() + static inline long step_absolute() { return it; } @@ -302,8 +302,8 @@ public: /// \brief Read a collective variable trajectory (post-processing /// only, not called at runtime) int read_traj(char const *traj_filename, - size_t traj_read_begin, - size_t traj_read_end); + long traj_read_begin, + long traj_read_end); /// Quick conversion of an object to a string template static std::string to_str(T const &x, diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 84777149c0..b564ec331b 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -227,6 +227,7 @@ size_t colvarparse::dummy_pos = 0; _get_keyval_scalar_(int); _get_keyval_scalar_(size_t); +_get_keyval_scalar_(long); _get_keyval_scalar_string_(std::string); _get_keyval_scalar_(cvm::real); _get_keyval_scalar_(cvm::rvector); @@ -238,6 +239,7 @@ _get_keyval_scalar_(colvarvalue); _get_keyval_vector_(int); _get_keyval_vector_(size_t); +_get_keyval_vector_(long); _get_keyval_vector_(std::string); _get_keyval_vector_(cvm::real); _get_keyval_vector_(cvm::rvector); diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index 5f7ba26ba3..5a753780a8 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -113,6 +113,7 @@ public: _get_keyval_scalar_proto_(int, (int)0); _get_keyval_scalar_proto_(size_t, (size_t)0); + _get_keyval_scalar_proto_(long, 0); _get_keyval_scalar_proto_(std::string, std::string("")); _get_keyval_scalar_proto_(cvm::real, (cvm::real)0.0); _get_keyval_scalar_proto_(cvm::rvector, cvm::rvector()); @@ -130,6 +131,7 @@ public: _get_keyval_vector_proto_(int, 0); _get_keyval_vector_proto_(size_t, 0); + _get_keyval_vector_proto_(long, 0); _get_keyval_vector_proto_(std::string, std::string("")); _get_keyval_vector_proto_(cvm::real, 0.0); _get_keyval_vector_proto_(cvm::rvector, cvm::rvector()); diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index b985ea5d2d..380f9d6fda 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -128,7 +128,18 @@ int colvarscript::run(int argc, char const *argv[]) { } } - /// TODO Write an output state file? (Useful for testing) + /// Save to an output state file + if (cmd == "save") { + if (argc < 3) { + result = "Missing arguments"; + return COLVARSCRIPT_ERROR; + } + proxy->output_prefix_str = argv[2]; + int error = 0; + error |= colvars->setup_output(); + error |= colvars->write_output_files(); + return error ? COLVARSCRIPT_ERROR : COLVARSCRIPT_OK; + } /// Print the values that would go on colvars.traj if (cmd == "printframelabels") { diff --git a/lib/colvars/colvartypes.cpp b/lib/colvars/colvartypes.cpp index 8e0aa57dfc..ad297a6ac9 100644 --- a/lib/colvars/colvartypes.cpp +++ b/lib/colvars/colvartypes.cpp @@ -501,7 +501,6 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector co cvm::quaternion 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,