diff --git a/doc/PDF/colvars-refman-lammps.pdf b/doc/PDF/colvars-refman-lammps.pdf index c828ddbefb..706617bcd5 100644 Binary files a/doc/PDF/colvars-refman-lammps.pdf and b/doc/PDF/colvars-refman-lammps.pdf differ 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..dd59f60eff 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 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..5f71f312e3 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" @@ -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_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp new file mode 100644 index 0000000000..29f3c9bd9f --- /dev/null +++ b/lib/colvars/colvarbias_histogram.cpp @@ -0,0 +1,226 @@ +/// -*- c++ -*- + +#include "colvarmodule.h" +#include "colvar.h" +#include "colvarbias_histogram.h" + +/// Histogram "bias" constructor + +colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *key) + : colvarbias(conf, key), + grid(NULL), out_name("") +{ + get_keyval(conf, "outputFile", out_name, std::string("")); + 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 = 0; + { + size_t i; + bool colvar_array = false; + get_keyval(conf, "gatherVectorColvars", colvar_array, colvar_array); + + if (colvar_array) { + for (i = 0; i < colvars.size(); i++) { // should be all vector + if (colvars[i]->value().type() != colvarvalue::type_vector) { + cvm::error("Error: used gatherVectorColvars with non-vector colvar.\n", INPUT_ERROR); + return; + } + if (i == 0) { + colvar_array_size = colvars[i]->value().size(); + if (colvar_array_size < 1) { + cvm::error("Error: vector variable has dimension less than one.\n", INPUT_ERROR); + return; + } + } else { + if (colvar_array_size != colvars[i]->value().size()) { + cvm::error("Error: trying to combine vector colvars of different lengths.\n", INPUT_ERROR); + return; + } + } + } + } else { + for (i = 0; i < colvars.size(); i++) { // should be all scalar + if (colvars[i]->value().type() != colvarvalue::type_scalar) { + cvm::error("Error: only scalar colvars are supported when gatherVectorColvars is off.\n", INPUT_ERROR); + return; + } + } + } + } + + if (colvar_array_size > 0) { + weights.assign(colvar_array_size, 1.0); + get_keyval(conf, "weights", weights, weights, colvarparse::parse_silent); + } + + grid = new colvar_grid_scalar(); + + { + std::string grid_conf; + if (key_lookup(conf, "grid", grid_conf)) { + grid->parse_params(grid_conf); + } else { + grid->init_from_colvars(colvars); + } + } + + 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); + } + + // assign a valid bin size + bin.assign(colvars.size(), 0); + + if (out_name.size() == 0) { + // 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 + "\""); + } + } + + if (colvar_array_size == 0) { + // update indices for scalar values + size_t i; + for (i = 0; i < colvars.size(); i++) { + bin[i] = grid->value_to_bin_scalar(colvars[i]->value(), i); + } + + if (grid->index_ok(bin)) { + grid->acc_value(bin, 1.0); + } + } else { + // update indices for vector/array values + size_t iv, i; + for (iv = 0; iv < colvar_array_size; iv++) { + for (i = 0; i < colvars.size(); i++) { + bin[i] = grid->value_to_bin_scalar(colvars[i]->value().vector1d_value[iv], i); + } + + if (grid->index_ok(bin)) { + grid->acc_value(bin, weights[iv]); + } + } + } + + if (output_freq && (cvm::step_absolute() % output_freq) == 0) { + write_output_files(); + } + return 0.0; // no bias energy for histogram +} + + +int colvarbias_histogram::write_output_files() +{ + if (out_name.size()) { + cvm::log("Writing the histogram file \""+out_name+"\".\n"); + + grid_os.open(out_name.c_str()); + if (!grid_os.is_open()) { + cvm::error("Error opening histogram file " + out_name + " for writing"); + } + // TODO add return code here + grid->write_multicol(grid_os); + grid_os.close(); + } + return COLVARS_OK; +} + + +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_histogram.h b/lib/colvars/colvarbias_histogram.h new file mode 100644 index 0000000000..b9d07151d3 --- /dev/null +++ b/lib/colvars/colvarbias_histogram.h @@ -0,0 +1,47 @@ +// -*- c++ -*- + +#ifndef COLVARBIAS_HISTOGRAM_H +#define COLVARBIAS_HISTOGRAM_H + +#include +#include +#include +#include + +#include "colvarbias.h" +#include "colvargrid.h" + +/// 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(); + + int write_output_files(); + +protected: + + /// n-dim histogram + colvar_grid_scalar *grid; + std::vector bin; + std::string out_name; + + size_t 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; + /// If colvar_array_size is larger than 1, weigh each one by this number before accumulating the histogram + std::vector weights; + + 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..a1426f71ee 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -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/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..b6f060b495 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-03" #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/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index 94c3bacdd4..74f7193784 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -15,7 +15,7 @@ #include #ifndef COLVARPROXY_VERSION -#define COLVARPROXY_VERSION "2015-04-01" +#define COLVARPROXY_VERSION "2015-04-02" #endif /* struct for packed data communication of coordinates and forces. */