From 067833cf074c9783f0dfafb1d8d40571bf390109 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 19 Feb 2015 18:20:52 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13126 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/colvars/colvar.cpp | 8 +- lib/colvars/colvar.h | 2 +- lib/colvars/colvaratoms.cpp | 9 +- lib/colvars/colvaratoms.h | 3 + lib/colvars/colvarbias_abf.cpp | 115 +++++++++--- lib/colvars/colvarbias_abf.h | 9 +- lib/colvars/colvarbias_meta.cpp | 18 +- lib/colvars/colvarbias_meta.h | 4 +- lib/colvars/colvarcomp.h | 4 +- lib/colvars/colvarcomp_angles.cpp | 9 - lib/colvars/colvarcomp_coordnums.cpp | 79 ++++---- lib/colvars/colvarcomp_distances.cpp | 63 +++++-- lib/colvars/colvarcomp_protein.cpp | 12 +- lib/colvars/colvargrid.cpp | 8 +- lib/colvars/colvargrid.h | 261 +++++++++++++++++---------- lib/colvars/colvarmodule.cpp | 152 ++++++++-------- lib/colvars/colvarmodule.h | 23 ++- lib/colvars/colvarparse.cpp | 36 ++-- lib/colvars/colvarparse.h | 52 +++--- lib/colvars/colvarproxy.h | 89 +++++++-- lib/colvars/colvartypes.h | 1 - 21 files changed, 597 insertions(+), 360 deletions(-) diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 852cc7984e..8f0850c0ff 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -177,7 +177,7 @@ colvar::colvar(std::string const &conf) cvm::error("Could not parse scripted colvar type."); return; } - x_reported.type (x.type()); + x_reported.type(x.type()); cvm::log(std::string("Expecting colvar value of type ") + colvarvalue::type_desc(x.type())); @@ -987,7 +987,9 @@ void colvar::calc() } } - if (tasks[task_system_force]) { + if (tasks[task_system_force] && !tasks[task_extended_lagrangian]) { + // If extended Lagrangian is enabled, system force calculation is trivial + // and done together with integration of the extended coordinate. if (tasks[task_scripted]) { // TODO see if this could reasonably be done in a generic way @@ -1594,7 +1596,7 @@ int colvar::write_output_files() if (acf.size()) { cvm::log("Writing acf to file \""+acf_outfile+"\".\n"); - std::ofstream acf_os(acf_outfile.c_str()); + cvm::ofstream acf_os(acf_outfile.c_str()); if (! acf_os.good()) { cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR); } diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index bb7d70d09c..2a9abb9d47 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -457,7 +457,7 @@ protected: /// Timesteps to skip between two values in the running average series size_t runave_stride; /// Name of the file to write the running average - std::ofstream runave_os; + cvm::ofstream runave_os; /// Current value of the running average colvarvalue runave; /// Current value of the square deviation from the running average diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index bd8c50ee34..99f63c5af6 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -315,21 +315,24 @@ int cvm::atom_group::parse(std::string const &conf, #endif if (!b_dummy) { + + // calculate total mass (TODO: this is the step that most needs deferred re-initialization) this->total_mass = 0.0; for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { this->total_mass += ai->mass; } - } - if (!b_dummy) { + // whether these atoms will ever receive forces or not bool enable_forces = true; // disableForces is deprecated if (get_keyval(group_conf, "enableForces", enable_forces, true, mode)) { noforce = !enable_forces; } else { - get_keyval(group_conf, "disableForces", noforce, false, mode); + 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 5a44d106cb..d950ffb7ac 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -139,6 +139,9 @@ 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 diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 543be6a95f..d9a56794fd 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -76,6 +76,8 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) if (!colvars[i]->tasks[colvar::task_extended_lagrangian]) { // request computation of Jacobian force + // ultimately, will be done regardless of extended Lagrangian + // and colvar should then just report zero Jacobian force colvars[i]->enable(colvar::task_Jacobian_force); // request Jacobian force as part as system force @@ -116,7 +118,11 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) force = new cvm::real [colvars.size()]; // Construct empty grids based on the colvars - samples = new colvar_grid_count (colvars); + if (cvm::debug()) { + cvm::log("Allocating count and free energy gradient grids.\n"); + } + + samples = new colvar_grid_count(colvars); gradients = new colvar_grid_gradient(colvars); gradients->samples = samples; samples->has_parent_data = true; @@ -124,7 +130,7 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) // For shared ABF, we store a second set of grids. // This used to be only if "shared" was defined, // but now we allow calling share externally (e.g. from Tcl). - last_samples = new colvar_grid_count (colvars); + last_samples = new colvar_grid_count(colvars); last_gradients = new colvar_grid_gradient(colvars); last_gradients->samples = last_samples; last_samples->has_parent_data = true; @@ -260,9 +266,9 @@ cvm::real colvarbias_abf::update() } if (b_history_files && (cvm::step_absolute() % history_freq) == 0) { cvm::log("ABFHISTORYFILE "+cvm::to_str(cvm::step_absolute())); - // append to existing file only if cvm::step_absolute() > 0 + // file already exists iff cvm::step_relative() > 0 // otherwise, backup and replace - write_gradients_samples(output_prefix + ".hist", (cvm::step_absolute() > 0)); + write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0)); } if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) { @@ -365,28 +371,32 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app std::string gradients_out_name = prefix + ".grad"; std::ios::openmode mode = (append ? std::ios::app : std::ios::out); - std::ofstream samples_os; - std::ofstream gradients_os; + cvm::ofstream samples_os; + cvm::ofstream gradients_os; if (!append) cvm::backup_file(samples_out_name.c_str()); samples_os.open(samples_out_name.c_str(), mode); - if (!samples_os.good()) cvm::error("Error opening ABF samples file " + samples_out_name + " for writing"); + if (!samples_os.is_open()) { + cvm::error("Error opening ABF samples file " + samples_out_name + " for writing"); + } samples->write_multicol(samples_os); samples_os.close(); if (!append) cvm::backup_file(gradients_out_name.c_str()); gradients_os.open(gradients_out_name.c_str(), mode); - if (!gradients_os.good()) cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing"); + if (!gradients_os.is_open()) { + cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing"); + } gradients->write_multicol(gradients_os); gradients_os.close(); if (colvars.size() == 1) { std::string pmf_out_name = prefix + ".pmf"; if (!append) cvm::backup_file(pmf_out_name.c_str()); - std::ofstream pmf_os; + cvm::ofstream pmf_os; // Do numerical integration and output a PMF pmf_os.open(pmf_out_name.c_str(), mode); - if (!pmf_os.good()) cvm::error("Error opening pmf file " + pmf_out_name + " for writing"); + if (!pmf_os.is_open()) cvm::error("Error opening pmf file " + pmf_out_name + " for writing"); gradients->write_1D_integral(pmf_os); pmf_os << std::endl; pmf_os.close(); @@ -428,13 +438,13 @@ void colvarbias_abf::read_gradients_samples() cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name); is.open(samples_in_name.c_str()); - if (!is.good()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading"); + if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading"); samples->read_multicol(is, true); is.close(); is.clear(); is.open(gradients_in_name.c_str()); - if (!is.good()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading"); + if (!is.is_open()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading"); gradients->read_multicol(is, true); is.close(); } @@ -550,13 +560,41 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const * : colvarbias(conf, key), grid(NULL), out_name("") { - get_keyval(conf, "outputfreq", output_freq, cvm::restart_out_freq); + 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"); + // } - 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); + } } - grid = new colvar_grid_count (colvars); bin.assign(colvars.size(), 0); cvm::log("Finished histogram setup.\n"); @@ -565,7 +603,8 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const * /// Destructor colvarbias_histogram::~colvarbias_histogram() { - if (grid_os.good()) grid_os.close(); + if (grid_os.is_open()) + grid_os.close(); if (grid) { delete grid; @@ -579,7 +618,10 @@ colvarbias_histogram::~colvarbias_histogram() /// Update the grid cvm::real colvarbias_histogram::update() { - if (cvm::debug()) cvm::log("Updating Grid bias " + this->name); + + 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 @@ -589,22 +631,49 @@ cvm::real colvarbias_histogram::update() cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\""); } - for (size_t i=0; icurrent_bin_scalar(i); + + 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 ( grid->index_ok(bin) ) { // Only within bounds of the grid... - grid->incr_count(bin); + 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.good()) cvm::error("Error opening histogram file " + out_name + " for writing"); + 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 } diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index f9f6cb702b..c1e3a1fb31 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -1,3 +1,4 @@ +// -*- c++ -*- /************************************************************************ * Headers for the ABF and histogram biases * ************************************************************************/ @@ -100,7 +101,7 @@ public: cvm::real update(); -private: +protected: /// n-dim histogram colvar_grid_count *grid; @@ -108,8 +109,12 @@ private: 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(); - std::ofstream grid_os; /// Stream for writing grid to disk + cvm::ofstream grid_os; /// Stream for writing grid to disk std::istream& read_restart(std::istream&); std::ostream& write_restart(std::ostream&); diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index e8e104366b..0276359f0d 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -86,7 +86,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key) colvars[i]->enable(colvar::task_grid); } - hills_energy = new colvar_grid_scalar (colvars); + hills_energy = new colvar_grid_scalar(colvars); hills_energy_gradients = new colvar_grid_gradient(colvars); } else { rebin_grids = false; @@ -430,11 +430,11 @@ cvm::real colvarbias_meta::update() new_hills_energy->lower_boundaries = new_lower_boundaries; new_hills_energy->upper_boundaries = new_upper_boundaries; - new_hills_energy->create(new_sizes, 0.0, 1); + new_hills_energy->setup(new_sizes, 0.0, 1); new_hills_energy_gradients->lower_boundaries = new_lower_boundaries; new_hills_energy_gradients->upper_boundaries = new_upper_boundaries; - new_hills_energy_gradients->create(new_sizes, 0.0, colvars.size()); + new_hills_energy_gradients->setup(new_sizes, 0.0, colvars.size()); new_hills_energy->map_grid(*hills_energy); new_hills_energy_gradients->map_grid(*hills_energy_gradients); @@ -746,6 +746,8 @@ void colvarbias_meta::calc_hills_force(size_t const &i, break; case colvarvalue::type_notset: + case colvarvalue::type_all: + default: break; } } @@ -959,7 +961,7 @@ void colvarbias_meta::update_replicas_registry() (replicas.back())->comm = multiple_replicas; if (use_grids) { - (replicas.back())->hills_energy = new colvar_grid_scalar (colvars); + (replicas.back())->hills_energy = new colvar_grid_scalar(colvars); (replicas.back())->hills_energy_gradients = new colvar_grid_gradient(colvars); } } @@ -1552,7 +1554,7 @@ void colvarbias_meta::write_pmf() { // allocate a new grid to store the pmf colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy); - pmf->create(); + pmf->setup(); std::string fes_file_name_prefix(cvm::output_prefix); @@ -1582,7 +1584,7 @@ void colvarbias_meta::write_pmf() "."+cvm::to_str(cvm::step_absolute()) : "") + ".pmf"); cvm::backup_file(fes_file_name.c_str()); - std::ofstream fes_os(fes_file_name.c_str()); + cvm::ofstream fes_os(fes_file_name.c_str()); pmf->write_multicol(fes_os); fes_os.close(); } @@ -1606,7 +1608,7 @@ void colvarbias_meta::write_pmf() "."+cvm::to_str(cvm::step_absolute()) : "") + ".pmf"); cvm::backup_file(fes_file_name.c_str()); - std::ofstream fes_os(fes_file_name.c_str()); + cvm::ofstream fes_os(fes_file_name.c_str()); pmf->write_multicol(fes_os); fes_os.close(); } @@ -1621,7 +1623,7 @@ void colvarbias_meta::write_replica_state_file() // write down also the restart for the other replicas: TODO: this // is duplicated code, that could be cleaned up later cvm::backup_file(replica_state_file.c_str()); - std::ofstream rep_state_os(replica_state_file.c_str()); + cvm::ofstream rep_state_os(replica_state_file.c_str()); if (!rep_state_os.good()) cvm::fatal_error("Error: in opening file \""+ replica_state_file+"\" for writing.\n"); diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index de5b6d2622..e2f175ad92 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -61,7 +61,7 @@ protected: /// Write the hill logfile bool b_hills_traj; /// Logfile of hill management (creation and deletion) - std::ofstream hills_traj_os; + cvm::ofstream hills_traj_os; /// \brief List of hills used on this bias (total); if a grid is /// employed, these don't need to be updated at every time step @@ -223,7 +223,7 @@ protected: std::string replica_hills_file; /// \brief Output stream corresponding to replica_hills_file - std::ofstream replica_hills_os; + cvm::ofstream replica_hills_os; /// Position within replica_hills_file (when reading it) int replica_hills_file_pos; diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 69348da364..62bd8a274e 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -815,8 +815,6 @@ class colvar::h_bond : public colvar::cvc { protected: - /// Atoms involved in the component - cvm::atom acceptor, donor; /// \brief "Cutoff" distance between acceptor and donor cvm::real r0; /// Integer exponent of the function numerator @@ -1141,6 +1139,8 @@ class colvar::cartesian protected: /// Atom group cvm::atom_group atoms; + /// Which Cartesian coordinates to include + std::vector axes; public: cartesian(std::string const &conf); cartesian(); diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index b300417598..726b6befcc 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -54,10 +54,6 @@ colvar::angle::angle() void colvar::angle::calc_value() { - group1.read_positions(); - group2.read_positions(); - group3.read_positions(); - cvm::atom_pos const g1_pos = group1.center_of_mass(); cvm::atom_pos const g2_pos = group2.center_of_mass(); cvm::atom_pos const g3_pos = group3.center_of_mass(); @@ -215,11 +211,6 @@ colvar::dihedral::dihedral() void colvar::dihedral::calc_value() { - group1.read_positions(); - group2.read_positions(); - group3.read_positions(); - group4.read_positions(); - cvm::atom_pos const g1_pos = group1.center_of_mass(); cvm::atom_pos const g2_pos = group2.center_of_mass(); cvm::atom_pos const g3_pos = group3.center_of_mass(); diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index d6509aedc0..7541c8bb23 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -62,8 +62,8 @@ cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec, if (calculate_gradients) { cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0); cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x, - (2.0/(r0_vec.y*r0_vec.y))*diff.y, - (2.0/(r0_vec.z*r0_vec.z))*diff.z); + (2.0/(r0_vec.y*r0_vec.y))*diff.y, + (2.0/(r0_vec.z*r0_vec.z))*diff.z); A1.grad += (-1.0)*dFdl2*dl2dx; A2.grad += dFdl2*dl2dx; } @@ -134,10 +134,10 @@ void colvar::coordnum::calc_value() if (b_anisotropic) { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - x.real_value += switching_function (r0_vec, en, ed, *ai1, group2_com_atom); + x.real_value += switching_function(r0_vec, en, ed, *ai1, group2_com_atom); } else { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - x.real_value += switching_function (r0, en, ed, *ai1, group2_com_atom); + x.real_value += switching_function(r0, en, ed, *ai1, group2_com_atom); } } else { @@ -145,12 +145,12 @@ void colvar::coordnum::calc_value() if (b_anisotropic) { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { - x.real_value += switching_function (r0_vec, en, ed, *ai1, *ai2); + x.real_value += switching_function(r0_vec, en, ed, *ai1, *ai2); } } else { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { - x.real_value += switching_function (r0, en, ed, *ai1, *ai2); + x.real_value += switching_function(r0, en, ed, *ai1, *ai2); } } } @@ -168,10 +168,10 @@ void colvar::coordnum::calc_gradients() if (b_anisotropic) { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - switching_function (r0_vec, en, ed, *ai1, group2_com_atom); + switching_function(r0_vec, en, ed, *ai1, group2_com_atom); } else { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) - switching_function (r0, en, ed, *ai1, group2_com_atom); + switching_function(r0, en, ed, *ai1, group2_com_atom); } group2.set_weighted_gradient(group2_com_atom.grad); @@ -181,27 +181,15 @@ void colvar::coordnum::calc_gradients() if (b_anisotropic) { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { - switching_function (r0_vec, en, ed, *ai1, *ai2); + switching_function(r0_vec, en, ed, *ai1, *ai2); } } else { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { - switching_function (r0, en, ed, *ai1, *ai2); + switching_function(r0, en, ed, *ai1, *ai2); } } } - - // if (cvm::debug()) { - // for (size_t i = 0; i < group1.size(); i++) { - // cvm::log("atom["+cvm::to_str (group1[i].id+1)+"] gradient: "+ - // cvm::to_str (group1[i].grad)+"\n"); - // } - - // for (size_t i = 0; i < group2.size(); i++) { - // cvm::log("atom["+cvm::to_str (group2[i].id+1)+"] gradient: "+ - // cvm::to_str (group2[i].grad)+"\n"); - // } - // } } void colvar::coordnum::apply_force(colvarvalue const &force) @@ -234,8 +222,8 @@ colvar::h_bond::h_bond(std::string const &conf) cvm::fatal_error("Error: either acceptor or donor undefined.\n"); } - acceptor = cvm::atom(a_num); - donor = cvm::atom(d_num); + cvm::atom acceptor = cvm::atom(a_num); + cvm::atom donor = cvm::atom(d_num); atom_groups.push_back(new cvm::atom_group); atom_groups[0]->add_atom(acceptor); atom_groups[0]->add_atom(donor); @@ -253,12 +241,10 @@ colvar::h_bond::h_bond(std::string const &conf) } -colvar::h_bond::h_bond(cvm::atom const &acceptor_i, - cvm::atom const &donor_i, - cvm::real r0_i, int en_i, int ed_i) - : acceptor(acceptor_i), - donor(donor_i), - r0(r0_i), en(en_i), ed(ed_i) +colvar::h_bond::h_bond(cvm::atom const &acceptor, + cvm::atom const &donor, + cvm::real r0_i, int en_i, int ed_i) + : r0(r0_i), en(en_i), ed(ed_i) { function_type = "h_bond"; x.type(colvarvalue::type_scalar); @@ -278,30 +264,25 @@ colvar::h_bond::h_bond() colvar::h_bond::~h_bond() { - for (unsigned int i=0; i (r0, en, ed, acceptor, donor); + x.real_value = colvar::coordnum::switching_function(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]); } void colvar::h_bond::calc_gradients() { - colvar::coordnum::switching_function (r0, en, ed, acceptor, donor); - (*atom_groups[0])[0].grad = acceptor.grad; - (*atom_groups[0])[1].grad = donor.grad; + colvar::coordnum::switching_function(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]); } void colvar::h_bond::apply_force(colvarvalue const &force) { - acceptor.apply_force(force.real_value * acceptor.grad); - donor.apply_force (force.real_value * donor.grad); + (atom_groups[0])->apply_colvar_force(force); } @@ -339,23 +320,27 @@ colvar::selfcoordnum::selfcoordnum() void colvar::selfcoordnum::calc_value() { x.real_value = 0.0; - - for (size_t i = 0; i < group1.size() - 1; i++) - for (size_t j = i + 1; j < group1.size(); j++) - x.real_value += colvar::coordnum::switching_function (r0, en, ed, group1[i], group1[j]); + for (size_t i = 0; i < group1.size() - 1; i++) { + for (size_t j = i + 1; j < group1.size(); j++) { + x.real_value += colvar::coordnum::switching_function(r0, en, ed, group1[i], group1[j]); + } + } } void colvar::selfcoordnum::calc_gradients() { - for (size_t i = 0; i < group1.size() - 1; i++) - for (size_t j = i + 1; j < group1.size(); j++) - colvar::coordnum::switching_function (r0, en, ed, group1[i], group1[j]); + for (size_t i = 0; i < group1.size() - 1; i++) { + for (size_t j = i + 1; j < group1.size(); j++) { + colvar::coordnum::switching_function(r0, en, ed, group1[i], group1[j]); + } + } } void colvar::selfcoordnum::apply_force(colvarvalue const &force) { - if (!group1.noforce) + if (!group1.noforce) { group1.apply_colvar_force(force.real_value); + } } diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 6b9317ab2e..f253aeffc6 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -447,9 +447,10 @@ colvar::distance_inv::distance_inv(std::string const &conf) for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { - if (ai1->id == ai2->id) - cvm::error("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n"); - return; + if (ai1->id == ai2->id) { + cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n"); + return; + } } } @@ -853,10 +854,11 @@ colvar::rmsd::rmsd(std::string const &conf) 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) { cvm::error("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); + "if provided, must be non-zero.\n"); return; + } } else { // if not, rely on existing atom indices for the group atoms.create_sorted_ids(); @@ -1294,17 +1296,40 @@ colvar::cartesian::cartesian(std::string const &conf) parse_group(conf, "atoms", atoms); atom_groups.push_back(&atoms); + bool use_x, use_y, use_z; + get_keyval(conf, "useX", use_x, true); + get_keyval(conf, "useY", use_y, true); + get_keyval(conf, "useZ", use_z, true); + + axes.clear(); + if (use_x) axes.push_back(0); + if (use_y) axes.push_back(1); + if (use_z) axes.push_back(2); + + if (axes.size() == 0) { + cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n"); + } + x.type(colvarvalue::type_vector); - x.vector1d_value.resize(atoms.size() * 3); + x.vector1d_value.resize(atoms.size() * axes.size()); } void colvar::cartesian::calc_value() { - int ia, j; + size_t const dim = axes.size(); + size_t ia, j; for (ia = 0; ia < atoms.size(); ia++) { - for (j = 0; j < 3; j++) { - x.vector1d_value[3*ia + j] = atoms[ia].pos[j]; + for (j = 0; j < dim; j++) { + 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]; + } } } } @@ -1320,14 +1345,24 @@ void colvar::cartesian::calc_gradients() void colvar::cartesian::apply_force(colvarvalue const &force) { - int ia, j; + size_t const dim = axes.size(); + size_t ia, j; if (!atoms.noforce) { cvm::rvector f; - for (ia = 0; ia < atoms.size(); ia++) { - for (j = 0; j < 3; j++) { - f[j] = force.vector1d_value[3*ia + j]; + 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); } - atoms[ia].apply_force(f); } } } diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index ddd3113fc0..cd37629282 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -66,8 +66,11 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) for (size_t i = 0; i < residues.size()-2; i++) { theta.push_back(new colvar::angle(cvm::atom(r[i ], "CA", sid), - cvm::atom(r[i+1], "CA", sid), - cvm::atom(r[i+2], "CA", sid))); + cvm::atom(r[i+1], "CA", sid), + cvm::atom(r[i+2], "CA", sid))); + atom_groups.push_back(theta.back()->atom_groups[0]); + atom_groups.push_back(theta.back()->atom_groups[1]); + atom_groups.push_back(theta.back()->atom_groups[2]); } } else { @@ -85,8 +88,9 @@ colvar::alpha_angles::alpha_angles(std::string const &conf) for (size_t i = 0; i < residues.size()-4; i++) { hb.push_back(new colvar::h_bond(cvm::atom(r[i ], "O", sid), - cvm::atom(r[i+4], "N", sid), - r0, en, ed)); + cvm::atom(r[i+4], "N", sid), + r0, en, ed)); + atom_groups.push_back(hb.back()->atom_groups[0]); } } else { diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index 4c2c3965e3..76f7688b7c 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -10,16 +10,18 @@ colvar_grid_count::colvar_grid_count() : colvar_grid() -{} +{ + mult = 1; +} colvar_grid_count::colvar_grid_count(std::vector const &nx_i, size_t const &def_count) - : colvar_grid(nx_i, def_count) + : colvar_grid(nx_i, def_count, 1) {} colvar_grid_count::colvar_grid_count(std::vector &colvars, size_t const &def_count) - : colvar_grid(colvars, def_count) + : colvar_grid(colvars, def_count, 1) {} std::istream & colvar_grid_count::read_restart(std::istream &is) diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 17ff6a4fb2..be2d7cfc8f 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -15,7 +15,7 @@ /// \brief Grid of values of a function of several collective /// variables \param T The data type /// -/// Only scalar colvars supported so far +/// Only scalar colvars supported so far: vector colvars are treated as arrays template class colvar_grid : public colvarparse { protected: @@ -56,7 +56,7 @@ protected: addr += ix[i]*nxc[i]; if (cvm::debug()) { if (ix[i] >= nx[i]) { - cvm::error("Error: exceeding bounds in colvar_grid.\n"); + cvm::error("Error: exceeding bounds in colvar_grid.\n", BUG_ERROR); return 0; } } @@ -125,35 +125,51 @@ public: return mult; } - /// \brief Allocate data (allow initialization also after construction) - void create(std::vector const &nx_i, - T const &t = T(), - size_t const &mult_i = 1) + /// \brief Allocate data + int setup(std::vector const &nx_i, + T const &t = T(), + size_t const &mult_i = 1) { - mult = mult_i; - nd = nx_i.size(); - nxc.resize(nd); - nx = nx_i; + if (cvm::debug()) { + cvm::log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+ + ", dimensions = "+cvm::to_str(nx_i)+".\n"); + } + mult = mult_i; + + data.clear(); + + nx = nx_i; + nd = nx.size(); + + nxc.resize(nd); + + // setup dimensions nt = mult; for (int i = nd-1; i >= 0; i--) { - if (nx_i[i] <= 0) { - cvm::error("Error: providing an invalid number of points, "+ - cvm::to_str(nx_i[i])+".\n"); - return; + if (nx[i] <= 0) { + cvm::error("Error: providing an invalid number of grid points, "+ + cvm::to_str(nx[i])+".\n", BUG_ERROR); + return COLVARS_ERROR; } nxc[i] = nt; nt *= nx[i]; } + if (cvm::debug()) { + cvm::log("Total number of grid elements = "+cvm::to_str(nt)+".\n"); + } + data.reserve(nt); data.assign(nt, t); + + return COLVARS_OK; } /// \brief Allocate data (allow initialization also after construction) - void create() + int setup() { - create(this->nx, T(), this->mult); + return setup(this->nx, T(), this->mult); } /// \brief Reset data (in case the grid is being reused) @@ -168,6 +184,8 @@ public: { save_delimiters = false; nd = nt = 0; + mult = 1; + this->setup(); } /// Destructor @@ -176,20 +194,20 @@ public: /// \brief "Almost copy-constructor": only copies configuration /// parameters from another grid, but doesn't reallocate stuff; - /// create() must be called after that; + /// setup() must be called after that; colvar_grid(colvar_grid const &g) : nd(g.nd), - nx(g.nx), - mult(g.mult), - data(), - cv(g.cv), - actual_value(g.actual_value), - lower_boundaries(g.lower_boundaries), - upper_boundaries(g.upper_boundaries), - periodic(g.periodic), - hard_lower_boundaries(g.hard_lower_boundaries), - hard_upper_boundaries(g.hard_upper_boundaries), - widths(g.widths), - has_data(false) + nx(g.nx), + mult(g.mult), + data(), + cv(g.cv), + actual_value(g.actual_value), + lower_boundaries(g.lower_boundaries), + upper_boundaries(g.upper_boundaries), + periodic(g.periodic), + hard_lower_boundaries(g.hard_lower_boundaries), + hard_upper_boundaries(g.hard_upper_boundaries), + widths(g.widths), + has_data(false) { save_delimiters = false; } @@ -199,48 +217,68 @@ public: /// the function at each point (optional) \param mult_i Multiplicity /// of each value colvar_grid(std::vector const &nx_i, - T const &t = T(), - size_t const &mult_i = 1) : has_data(false) + T const &t = T(), + size_t mult_i = 1) + : has_data(false) { save_delimiters = false; - this->create(nx_i, t, mult_i); + this->setup(nx_i, t, mult_i); } /// \brief Constructor from a vector of colvars colvar_grid(std::vector const &colvars, - T const &t = T(), - size_t const &mult_i = 1, - bool margin = false) - : cv(colvars), has_data(false) + T const &t = T(), + size_t mult_i = 1, + bool margin = false) + : has_data(false) { save_delimiters = false; + this->init_from_colvars(colvars, t, mult_i, margin); + } - std::vector nx_i; + int init_from_colvars(std::vector const &colvars, + T const &t = T(), + size_t mult_i = 1, + bool margin = false) + { + if (cvm::debug()) { + cvm::log("Reading grid configuration from collective variables.\n"); + } - if (cvm::debug()) + cv = colvars; + nd = colvars.size(); + mult = mult_i; + + size_t i; + + for (i = 0; i < cv.size(); i++) { + if (!cv[i]->tasks[colvar::task_lower_boundary] || + !cv[i]->tasks[colvar::task_upper_boundary]) { + cvm::error("Tried to initialize a grid on a " + "variable with undefined boundaries.\n", INPUT_ERROR); + return COLVARS_ERROR; + } + } + + if (cvm::debug()) { cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+ - " collective variables.\n"); + " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n"); + } - for (size_t i = 0; i < cv.size(); i++) { + for (i = 0; i < cv.size(); i++) { if (cv[i]->value().type() != colvarvalue::type_scalar) { cvm::error("Colvar grids can only be automatically " - "constructed for scalar variables. " - "ABF and histogram can not be used; metadynamics " - "can be used with useGrids disabled.\n"); - return; + "constructed for scalar variables. " + "ABF and histogram can not be used; metadynamics " + "can be used with useGrids disabled.\n", INPUT_ERROR); + return COLVARS_ERROR; } if (cv[i]->width <= 0.0) { cvm::error("Tried to initialize a grid on a " - "variable with negative or zero width.\n"); - return; - } - - if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) { - cvm::error("Tried to initialize a grid on a " - "variable with undefined boundaries.\n"); - return; + "variable with negative or zero width.\n", INPUT_ERROR); + return COLVARS_ERROR; } widths.push_back(cv[i]->width); @@ -271,35 +309,49 @@ public: lower_boundaries.push_back(cv[i]->lower_boundary); upper_boundaries.push_back(cv[i]->upper_boundary); } - - - { - cvm::real nbins = ( upper_boundaries[i].real_value - - lower_boundaries[i].real_value ) / widths[i]; - int nbins_round = (int)(nbins+0.5); - - if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) { - cvm::log("Warning: grid interval("+ - cvm::to_str(lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+ - cvm::to_str(upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+ - ") is not commensurate to its bin width("+ - cvm::to_str(widths[i], cvm::cv_width, cvm::cv_prec)+").\n"); - upper_boundaries[i].real_value = lower_boundaries[i].real_value + - (nbins_round * widths[i]); - } - - if (cvm::debug()) - cvm::log("Number of points is "+cvm::to_str((int) nbins_round)+ - " for the colvar no. "+cvm::to_str(i+1)+".\n"); - - nx_i.push_back(nbins_round); - } - } - create(nx_i, t, mult_i); + this->init_from_boundaries(); + return this->setup(); } + int init_from_boundaries(T const &t = T(), + size_t const &mult_i = 1) + { + if (cvm::debug()) { + cvm::log("Configuring grid dimensions from colvars boundaries.\n"); + } + + // these will have to be updated + nx.clear(); + nxc.clear(); + nt = 0; + + for (size_t i = 0; i < lower_boundaries.size(); i++) { + + cvm::real nbins = ( upper_boundaries[i].real_value - + lower_boundaries[i].real_value ) / widths[i]; + int nbins_round = (int)(nbins+0.5); + + if (std::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) { + cvm::log("Warning: grid interval("+ + cvm::to_str(lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+ + cvm::to_str(upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+ + ") is not commensurate to its bin width("+ + cvm::to_str(widths[i], cvm::cv_width, cvm::cv_prec)+").\n"); + upper_boundaries[i].real_value = lower_boundaries[i].real_value + + (nbins_round * widths[i]); + } + + if (cvm::debug()) + cvm::log("Number of points is "+cvm::to_str((int) nbins_round)+ + " for the colvar no. "+cvm::to_str(i+1)+".\n"); + + nx.push_back(nbins_round); + } + + return COLVARS_OK; + } /// Wrap an index vector around periodic boundary conditions /// also checks validity of non-periodic indices @@ -309,14 +361,16 @@ public: if (periodic[i]) { ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result } else { - if (ix[i] < 0 || ix[i] >= nx[i]) - cvm::error("Trying to wrap illegal index vector(non-PBC): " - + cvm::to_str(ix)); + if (ix[i] < 0 || ix[i] >= nx[i]) { + cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: " + + cvm::to_str(ix), BUG_ERROR); return; + } } } } + /// \brief Report the bin corresponding to the current value of variable i inline int current_bin_scalar(int const i) const { @@ -373,7 +427,7 @@ public: cvm::error("Error: trying to subtract two grids with " "different multiplicity.\n"); return; - } + } if (other_grid.data.size() != this->data.size()) { cvm::error("Error: trying to subtract two grids with " @@ -395,8 +449,8 @@ public: if (other_grid.multiplicity() != this->multiplicity()) { cvm::error("Error: trying to copy two grids with " "different multiplicity.\n"); - return; - } + return; + } if (other_grid.data.size() != this->data.size()) { cvm::error("Error: trying to copy two grids with " @@ -673,9 +727,11 @@ public: return os; } - - bool parse_params(std::string const &conf) + /// Read a grid definition from a config string + int parse_params(std::string const &conf) { + if (cvm::debug()) cvm::log("Reading grid configuration from string.\n"); + std::vector old_nx = nx; std::vector old_lb = lower_boundaries; @@ -688,36 +744,52 @@ public: cvm::to_str(nd_in)+") than the grid defined " "in the configuration file("+cvm::to_str(nd)+ ").\n"); - return false; + return COLVARS_ERROR; } } colvarparse::get_keyval(conf, "lower_boundaries", lower_boundaries, lower_boundaries, colvarparse::parse_silent); - colvarparse::get_keyval(conf, "upper_boundaries", upper_boundaries, upper_boundaries, colvarparse::parse_silent); + // support also camel case + colvarparse::get_keyval(conf, "lowerBoundaries", + lower_boundaries, lower_boundaries, colvarparse::parse_silent); + colvarparse::get_keyval(conf, "upperBoundaries", + upper_boundaries, upper_boundaries, colvarparse::parse_silent); + colvarparse::get_keyval(conf, "widths", widths, widths, colvarparse::parse_silent); colvarparse::get_keyval(conf, "sizes", nx, nx, colvarparse::parse_silent); + if (nd < lower_boundaries.size()) nd = lower_boundaries.size(); + + if (! actual_value.size()) actual_value.assign(nd, false); + if (! periodic.size()) periodic.assign(nd, false); + if (! widths.size()) widths.assign(nd, 1.0); + bool new_params = false; - for (size_t i = 0; i < nd; i++) { - if ( (old_nx[i] != nx[i]) || - (std::sqrt(cv[i]->dist2(old_lb[i], + if (old_nx.size()) { + for (size_t i = 0; i < nd; i++) { + if ( (old_nx[i] != nx[i]) || + (std::sqrt(cv[i]->dist2(old_lb[i], lower_boundaries[i])) > 1.0E-10) ) { - new_params = true; + new_params = true; + } } + } else { + new_params = true; } // reallocate the array in case the grid params have just changed if (new_params) { - data.resize(0); - this->create(nx, T(), mult); + init_from_boundaries(); + // data.resize(0); // no longer needed: setup calls clear() + return this->setup(nx, T(), mult); } - return true; + return COLVARS_OK; } /// \brief Check that the grid information inside (boundaries, @@ -834,6 +906,7 @@ public: << periodic[i] << "\n"; } + for (std::vector ix = new_index(); index_ok(ix); incr(ix) ) { if (ix.back() == 0) { diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 2528088926..f737a14fc7 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -24,12 +24,12 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in) // TODO relax this error to handle multiple molecules in VMD // once the module is not static anymore cvm::error("Error: trying to allocate the collective " - "variable module twice.\n"); + "variable module twice.\n"); return; } cvm::log(cvm::line_marker); cvm::log("Initializing the collective variables module, version "+ - cvm::to_str(COLVARS_VERSION)+".\n"); + cvm::to_str(COLVARS_VERSION)+".\n"); // set initial default values @@ -56,14 +56,14 @@ int colvarmodule::config_file(char const *config_filename) { cvm::log(cvm::line_marker); cvm::log("Reading new configuration from file \""+ - std::string(config_filename)+"\":\n"); + std::string(config_filename)+"\":\n"); // open the configfile config_s.open(config_filename); - if (!config_s) { + if (!config_s.is_open()) { cvm::error("Error: in opening configuration file \""+ - std::string(config_filename)+"\".\n", - FILE_ERROR); + std::string(config_filename)+"\".\n", + FILE_ERROR); return COLVARS_ERROR; } @@ -134,8 +134,10 @@ int colvarmodule::config(std::string &conf) cvm::log("Collective variables module (re)initialized.\n"); cvm::log(cvm::line_marker); - // configuration might have changed, better redo the labels - write_traj_label(cv_traj_os); + if (cv_traj_os.is_open()) { + // configuration might have changed, better redo the labels + write_traj_label(cv_traj_os); + } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -151,24 +153,24 @@ int colvarmodule::parse_global_params(std::string const &conf) parse->get_keyval(conf, "analysis", b_analysis, b_analysis); parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size, - debug_gradients_step_size, - colvarparse::parse_silent); + debug_gradients_step_size, + colvarparse::parse_silent); parse->get_keyval(conf, "eigenvalueCrossingThreshold", - colvarmodule::rotation::crossing_threshold, - colvarmodule::rotation::crossing_threshold, - colvarparse::parse_silent); + colvarmodule::rotation::crossing_threshold, + colvarmodule::rotation::crossing_threshold, + colvarparse::parse_silent); parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq); parse->get_keyval(conf, "colvarsRestartFrequency", - restart_out_freq, restart_out_freq); + restart_out_freq, restart_out_freq); // if this is true when initializing, it means // we are continuing after a reset(): default to true parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append); parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false, - colvarparse::parse_silent); + colvarparse::parse_silent); if (use_scripted_forces && !proxy->force_script_defined) { cvm::fatal_error("User script for scripted colvar forces not found."); @@ -193,8 +195,8 @@ int colvarmodule::parse_colvars(std::string const &conf) colvars.push_back(new colvar(colvar_conf)); if (cvm::get_error() || ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) { - cvm::log("Error while constructing colvar number " + - cvm::to_str(colvars.size()) + " : deleting."); + cvm::log("Error while constructing colvar number " + + cvm::to_str(colvars.size()) + " : deleting."); delete colvars.back(); // the colvar destructor updates the colvars array return COLVARS_ERROR; } @@ -214,8 +216,8 @@ int colvarmodule::parse_colvars(std::string const &conf) if (colvars.size()) cvm::log(cvm::line_marker); cvm::log("Collective variables initialized, "+ - cvm::to_str(colvars.size())+ - " in total.\n"); + cvm::to_str(colvars.size())+ + " in total.\n"); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } @@ -223,9 +225,9 @@ int colvarmodule::parse_colvars(std::string const &conf) bool colvarmodule::check_new_bias(std::string &conf, char const *key) { if (cvm::get_error() || - (biases.back()->check_keywords(conf, key) != COLVARS_OK)) { + (biases.back()->check_keywords(conf, key) != COLVARS_OK)) { cvm::log("Error while constructing bias number " + - cvm::to_str(biases.size()) + " : deleting.\n"); + cvm::to_str(biases.size()) + " : deleting.\n"); delete biases.back(); // the bias destructor updates the biases array return true; } @@ -330,7 +332,7 @@ colvar *colvarmodule::colvar_by_name(std::string const &name) { int colvarmodule::change_configuration(std::string const &bias_name, - std::string const &conf) + std::string const &conf) { // This is deprecated; supported strategy is to delete the bias // and parse the new config @@ -357,7 +359,7 @@ std::string colvarmodule::read_colvar(std::string const &name) } cvm::real colvarmodule::energy_difference(std::string const &bias_name, - std::string const &conf) + std::string const &conf) { cvm::increase_depth(); colvarbias *b; @@ -449,7 +451,7 @@ int colvarmodule::calc() { if (cvm::debug()) { cvm::log(cvm::line_marker); cvm::log("Collective variables module, step no. "+ - cvm::to_str(cvm::step_absolute())+"\n"); + cvm::to_str(cvm::step_absolute())+"\n"); } // calculate collective variables and their gradients @@ -533,7 +535,7 @@ int colvarmodule::calc() { // equation of motion (extended system) if (cvm::debug()) cvm::log("Updating the internal degrees of freedom " - "of colvars (if they have any).\n"); + "of colvars (if they have any).\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { total_colvar_energy += (*cvi)->update(); @@ -564,10 +566,10 @@ int colvarmodule::calc() { if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log("Writing the state file \""+ - restart_out_name+"\".\n"); + restart_out_name+"\".\n"); proxy->backup_file(restart_out_name.c_str()); restart_out_os.open(restart_out_name.c_str()); - if (!write_restart(restart_out_os)) + if (!restart_out_os.is_open() || !write_restart(restart_out_os)) cvm::error("Error: in writing restart file.\n"); restart_out_os.close(); } @@ -594,7 +596,7 @@ int colvarmodule::calc() { if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+ - cv_traj_name+"\".\n"); + cv_traj_name+"\".\n"); cv_traj_os.flush(); } } @@ -694,8 +696,8 @@ int colvarmodule::setup_input() std::ifstream input_is(restart_in_name.c_str()); if (!input_is.good()) { cvm::error("Error: in opening restart file \""+ - std::string(restart_in_name)+"\".\n", - FILE_ERROR); + std::string(restart_in_name)+"\".\n", + FILE_ERROR); return COLVARS_ERROR; } else { cvm::log("Restarting from file \""+restart_in_name+"\".\n"); @@ -725,9 +727,9 @@ int colvarmodule::setup_output() output_prefix = proxy->output_prefix(); if (output_prefix.size()) { cvm::log("The final output state file will be \""+ - (output_prefix.size() ? - std::string(output_prefix+".colvars.state") : - std::string("colvars.state"))+"\".\n"); + (output_prefix.size() ? + std::string(output_prefix+".colvars.state") : + std::string("colvars.state"))+"\".\n"); // cvm::log (cvm::line_marker); } @@ -752,8 +754,8 @@ std::istream & colvarmodule::read_restart(std::istream &is) if (is >> colvarparse::read_block("configuration", restart_conf)) { if (it_restart_from_state_file) { parse->get_keyval(restart_conf, "step", - it_restart, (size_t) 0, - colvarparse::parse_silent); + it_restart, (size_t) 0, + colvarparse::parse_silent); it = it_restart; } } @@ -767,8 +769,8 @@ std::istream & colvarmodule::read_restart(std::istream &is) cvi++) { if ( !((*cvi)->read_restart(is)) ) { cvm::error("Error: in reading restart configuration for collective variable \""+ - (*cvi)->name+"\".\n", - INPUT_ERROR); + (*cvi)->name+"\".\n", + INPUT_ERROR); } } @@ -776,10 +778,11 @@ std::istream & colvarmodule::read_restart(std::istream &is) for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { - if (!((*bi)->read_restart(is))) + if (!((*bi)->read_restart(is))) { cvm::error("Error: in reading restart configuration for bias \""+ - (*bi)->name+"\".\n", - INPUT_ERROR); + (*bi)->name+"\".\n", + INPUT_ERROR); + } } cvm::decrease_depth(); @@ -802,11 +805,11 @@ int colvarmodule::write_output_files() std::string(output_prefix+".colvars.state") : std::string("colvars.state")); cvm::log("Saving collective variables state to \""+out_name+"\".\n"); - proxy->backup_file(out_name.c_str()); - std::ofstream out(out_name.c_str()); - out.setf(std::ios::scientific, std::ios::floatfield); - this->write_restart(out); - out.close(); + + std::ostream * os = proxy->output_stream(out_name); + os->setf(std::ios::scientific, std::ios::floatfield); + this->write_restart(*os); + proxy->close_output_stream(out_name); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); @@ -824,11 +827,11 @@ int colvarmodule::write_output_files() int colvarmodule::read_traj(char const *traj_filename, - size_t traj_read_begin, - size_t traj_read_end) + size_t traj_read_begin, + size_t traj_read_end) { cvm::log("Opening trajectory file \""+ - std::string(traj_filename)+"\".\n"); + std::string(traj_filename)+"\".\n"); std::ifstream traj_is(traj_filename); while (true) { @@ -839,7 +842,7 @@ int colvarmodule::read_traj(char const *traj_filename, do { if (!colvarparse::getline_nocomments(traj_is, line)) { cvm::log("End of file \""+std::string(traj_filename)+ - "\" reached, or corrupted file.\n"); + "\" reached, or corrupted file.\n"); traj_is.close(); return false; } @@ -867,8 +870,8 @@ int colvarmodule::read_traj(char const *traj_filename, (it > traj_read_end) ) { std::cerr << "\n"; cvm::error("Reached the end of the trajectory, " - "read_end = "+cvm::to_str(traj_read_end)+"\n", - FILE_ERROR); + "read_end = "+cvm::to_str(traj_read_end)+"\n", + FILE_ERROR); return COLVARS_ERROR; } @@ -877,9 +880,9 @@ int colvarmodule::read_traj(char const *traj_filename, cvi++) { if (!(*cvi)->read_traj(is)) { cvm::error("Error: in reading colvar \""+(*cvi)->name+ - "\" from trajectory file \""+ - std::string(traj_filename)+"\".\n", - FILE_ERROR); + "\" from trajectory file \""+ + std::string(traj_filename)+"\".\n", + FILE_ERROR); return COLVARS_ERROR; } } @@ -927,18 +930,18 @@ int colvarmodule::open_traj_file(std::string const &file_name) // (re)open trajectory file if (cv_traj_append) { cvm::log("Appending to colvar trajectory file \""+file_name+ - "\".\n"); + "\".\n"); cv_traj_os.open(file_name.c_str(), std::ios::app); } else { cvm::log("Writing to colvar trajectory file \""+file_name+ - "\".\n"); + "\".\n"); proxy->backup_file(file_name.c_str()); - cv_traj_os.open(file_name.c_str(), std::ios::out); + cv_traj_os.open(file_name.c_str()); } if (!cv_traj_os.is_open()) { cvm::error("Error: cannot write to file \""+file_name+"\".\n", - FILE_ERROR); + FILE_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -1046,10 +1049,11 @@ void cvm::exit(std::string const &message) int cvm::read_index_file(char const *filename) { std::ifstream is(filename, std::ios::binary); - if (!is.good()) + if (!is.good()) { cvm::error("Error: in opening index file \""+ - std::string(filename)+"\".\n", - FILE_ERROR); + std::string(filename)+"\".\n", + FILE_ERROR); + } while (is.good()) { char open, close; @@ -1062,8 +1066,8 @@ int cvm::read_index_file(char const *filename) names_i++) { if (*names_i == group_name) { cvm::error("Error: the group name \""+group_name+ - "\" appears in multiple index files.\n", - FILE_ERROR); + "\" appears in multiple index files.\n", + FILE_ERROR); } } cvm::index_group_names.push_back(group_name); @@ -1093,7 +1097,7 @@ int cvm::read_index_file(char const *filename) } cvm::log("The following index groups were read from the index file \""+ - std::string(filename)+"\":\n"); + std::string(filename)+"\":\n"); std::list::iterator names_i = index_group_names.begin(); std::list >::iterator lists_i = index_groups.begin(); for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) { @@ -1103,18 +1107,18 @@ int cvm::read_index_file(char const *filename) } int cvm::load_atoms(char const *file_name, - std::vector &atoms, - std::string const &pdb_field, - double const pdb_field_value) + std::vector &atoms, + std::string const &pdb_field, + double const pdb_field_value) { return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value); } int cvm::load_coords(char const *file_name, - std::vector &pos, - const std::vector &indices, - std::string const &pdb_field, - double const pdb_field_value) + std::vector &pos, + const std::vector &indices, + std::string const &pdb_field, + double const pdb_field_value) { // Differentiate between PDB and XYZ files // for XYZ files, use CVM internal parser @@ -1134,8 +1138,8 @@ int cvm::load_coords(char const *file_name, int cvm::load_coords_xyz(char const *filename, - std::vector &pos, - const std::vector &indices) + std::vector &pos, + const std::vector &indices) { std::ifstream xyz_is(filename); unsigned int natoms; @@ -1144,7 +1148,7 @@ int cvm::load_coords_xyz(char const *filename, if ( ! (xyz_is >> natoms) ) { cvm::error("Error: cannot parse XYZ file " - + std::string(filename) + ".\n", INPUT_ERROR); + + std::string(filename) + ".\n", INPUT_ERROR); } // skip comment line std::getline(xyz_is, line); diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 3c78d9c60b..d21b4a2885 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -4,7 +4,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2014-11-21" +#define COLVARS_VERSION "2015-02-04" #endif #ifndef COLVARS_DEBUG @@ -45,6 +45,10 @@ #include #include +#ifdef NAMD_VERSION +// use Lustre-friendly wrapper to POSIX write() +#include "fstream_namd.h" +#endif class colvarparse; class colvar; @@ -239,6 +243,12 @@ public: /// (Re)initialize the output trajectory and state file (does not write it yet) int setup_output(); +#ifdef NAMD_VERSION + typedef ofstream_namd ofstream; +#else + typedef std::ofstream ofstream; +#endif + /// Read the input restart file std::istream & read_restart(std::istream &is); /// Write the output restart file @@ -443,6 +453,7 @@ public: /// Pseudo-random number with Gaussian distribution static real rand_gaussian(void); + protected: /// Configuration file @@ -452,16 +463,16 @@ protected: colvarparse *parse; /// Name of the trajectory file - std::string cv_traj_name; + std::string cv_traj_name; /// Collective variables output trajectory file - std::ofstream cv_traj_os; + colvarmodule::ofstream cv_traj_os; /// Appending to the existing trajectory file? - bool cv_traj_append; + bool cv_traj_append; /// Output restart file - std::ofstream restart_out_os; + colvarmodule::ofstream restart_out_os; /// \brief Counter for the current depth in the object hierarchy (useg e.g. in output static size_t depth; @@ -470,12 +481,12 @@ protected: static bool use_scripted_forces; public: + /// \brief Pointer to the proxy object, used to retrieve atomic data /// from the hosting program; it is static in order to be accessible /// from static functions in the colvarmodule class static colvarproxy *proxy; - /// Increase the depth (number of indentations in the output) static void increase_depth(); diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 291c0108e4..84777149c0 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -47,20 +47,22 @@ size_t colvarparse::dummy_pos = 0; if (data.size()) { \ std::istringstream is(data); \ TYPE x(def_value); \ - if (is >> x) \ + if (is >> x) { \ value = x; \ - else \ + } else { \ cvm::error("Error: in parsing \""+ \ std::string(key)+"\".\n", INPUT_ERROR); \ + } \ if (parse_mode != parse_silent) { \ cvm::log("# "+std::string(key)+" = "+ \ cvm::to_str(value)+"\n"); \ } \ } else { \ \ - if (b_found_any) \ + if (b_found_any) { \ cvm::error("Error: improper or missing value " \ "for \""+std::string(key)+"\".\n", INPUT_ERROR); \ + } \ value = def_value; \ if (parse_mode != parse_silent) { \ cvm::log("# "+std::string(key)+" = \""+ \ @@ -110,19 +112,21 @@ size_t colvarparse::dummy_pos = 0; if (data_count == 0) \ cvm::fatal_error("Error: in parsing \""+ \ std::string(key)+"\".\n"); \ - if (data_count > 1) \ + if (data_count > 1) { \ cvm::error("Error: multiple values " \ "are not allowed for keyword \""+ \ std::string(key)+"\".\n", INPUT_ERROR); \ + } \ if (parse_mode != parse_silent) { \ cvm::log("# "+std::string(key)+" = "+ \ cvm::to_str(value)+"\n"); \ } \ } else { \ \ - if (b_found_any) \ + if (b_found_any) { \ cvm::error("Error: improper or missing value " \ "for \""+std::string(key)+"\".\n", INPUT_ERROR); \ + } \ value = def_value; \ if (parse_mode != parse_silent) { \ cvm::log("# "+std::string(key)+" = "+ \ @@ -185,11 +189,12 @@ size_t colvarparse::dummy_pos = 0; size_t i = 0; \ for ( ; i < values.size(); i++) { \ TYPE x(values[i]); \ - if (is >> x) \ + if (is >> x) { \ values[i] = x; \ - else \ + } else { \ cvm::error("Error: in parsing \""+ \ std::string(key)+"\".\n", INPUT_ERROR); \ + } \ } \ } \ \ @@ -200,9 +205,10 @@ size_t colvarparse::dummy_pos = 0; \ } else { \ \ - if (b_found_any) \ + if (b_found_any) { \ cvm::error("Error: improper or missing values for \""+ \ std::string(key)+"\".\n", INPUT_ERROR); \ + } \ \ for (size_t i = 0; i < values.size(); i++) \ values[i] = def_values[ (i > def_values.size()) ? 0 : i ]; \ @@ -476,16 +482,6 @@ bool colvarparse::key_lookup(std::string const &conf, } } - // check it is not between quotes - // if ( (conf.find_last_of ("\"", - // conf.find_last_of (white_space, pos)) != - // std::string::npos) && - // (conf.find_first_of ("\"", - // conf.find_first_of (white_space, pos)) != - // std::string::npos) ) - // return false; - - // save the pointer for a future call (when iterating over multiple // valid instances of the same keyword) save_pos = pos + key.size(); @@ -493,7 +489,7 @@ bool colvarparse::key_lookup(std::string const &conf, // get the remainder of the line size_t pl = conf.rfind("\n", pos); size_t line_begin = (pl == std::string::npos) ? 0 : pos; - size_t nl = conf.find ("\n", pos); + size_t nl = conf.find("\n", pos); size_t line_end = (nl == std::string::npos) ? conf.size() : nl; std::string line(conf, line_begin, (line_end-line_begin)); @@ -563,7 +559,7 @@ bool colvarparse::key_lookup(std::string const &conf, if (data.size() && save_delimiters) { data_begin_pos.push_back(conf.find(data, pos+key.size())); - data_end_pos.push_back (data_begin_pos.back()+data.size()); + data_end_pos.push_back(data_begin_pos.back()+data.size()); // std::cerr << "key = " << key << ", data = \"" // << data << "\", data_begin, data_end = " // << data_begin_pos.back() << ", " << data_end_pos.back() diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index af348464ea..2c7b8e000e 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -92,36 +92,36 @@ public: /// wrapper class (colvarvalue.h). #define _get_keyval_scalar_proto_(_type_,_def_value_) \ - bool get_keyval(std::string const &conf, \ - char const *key, \ - _type_ &value, \ - _type_ const &def_value = _def_value_, \ - Parse_Mode const parse_mode = parse_normal) + bool get_keyval(std::string const &conf, \ + char const *key, \ + _type_ &value, \ + _type_ const &def_value = _def_value_, \ + Parse_Mode const parse_mode = parse_normal) - _get_keyval_scalar_proto_(int, (int)0); - _get_keyval_scalar_proto_(size_t, (size_t)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()); - _get_keyval_scalar_proto_(cvm::quaternion, cvm::quaternion()); - _get_keyval_scalar_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset)); - _get_keyval_scalar_proto_(bool, false); + _get_keyval_scalar_proto_(int, (int)0); + _get_keyval_scalar_proto_(size_t, (size_t)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()); + _get_keyval_scalar_proto_(cvm::quaternion, cvm::quaternion()); + _get_keyval_scalar_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset)); + _get_keyval_scalar_proto_(bool, false); #define _get_keyval_vector_proto_(_type_,_def_value_) \ - bool get_keyval(std::string const &conf, \ - char const *key, \ - std::vector<_type_> &values, \ - std::vector<_type_> const &def_values = \ - std::vector<_type_> (0, static_cast<_type_>(_def_value_)), \ - Parse_Mode const parse_mode = parse_normal) + bool get_keyval(std::string const &conf, \ + char const *key, \ + std::vector<_type_> &values, \ + std::vector<_type_> const &def_values = \ + std::vector<_type_> (0, static_cast<_type_>(_def_value_)), \ + Parse_Mode const parse_mode = parse_normal) - _get_keyval_vector_proto_(int, 0); - _get_keyval_vector_proto_(size_t, 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()); - _get_keyval_vector_proto_(cvm::quaternion, cvm::quaternion()); - _get_keyval_vector_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset)); + _get_keyval_vector_proto_(int, 0); + _get_keyval_vector_proto_(size_t, 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()); + _get_keyval_vector_proto_(cvm::quaternion, cvm::quaternion()); + _get_keyval_vector_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset)); /// \brief Check that all the keywords within "conf" are in the list diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index e9388f68fb..47806990b7 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -3,6 +3,8 @@ #ifndef COLVARPROXY_H #define COLVARPROXY_H +#include +#include #include "colvarmodule.h" #include "colvarvalue.h" @@ -128,25 +130,25 @@ public: /// \brief Get the PBC-aware distance vector between two positions virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) = 0; + cvm::atom_pos const &pos2) = 0; /// \brief Get the PBC-aware square distance between two positions; /// may be implemented independently from position_distance() for optimization purposes virtual cvm::real position_dist2(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2); + cvm::atom_pos const &pos2); /// \brief Get the closest periodic image to a reference position /// \param pos The position to look for the closest periodic image /// \param ref_pos The reference position virtual void select_closest_image(cvm::atom_pos &pos, - cvm::atom_pos const &ref_pos) = 0; + cvm::atom_pos const &ref_pos) = 0; /// \brief Perform select_closest_image() on a set of atomic positions /// /// After that, distance vectors can then be calculated directly, /// without using position_distance() void select_closest_images(std::vector &pos, - cvm::atom_pos const &ref_pos); + cvm::atom_pos const &ref_pos); // **************** SCRIPTING INTERFACE **************** @@ -164,13 +166,13 @@ public: virtual int run_force_callback() { return COLVARS_NOT_IMPLEMENTED; } virtual int run_colvar_callback(std::string const &name, - std::vector const &cvcs, - colvarvalue &value) + std::vector const &cvcs, + colvarvalue &value) { return COLVARS_NOT_IMPLEMENTED; } virtual int run_colvar_gradient_callback(std::string const &name, - std::vector const &cvcs, - std::vector > &gradient) + std::vector const &cvcs, + std::vector > &gradient) { return COLVARS_NOT_IMPLEMENTED; } // **************** INPUT/OUTPUT **************** @@ -193,27 +195,78 @@ public: /// "filename" is a PDB file, use this field to determine which are /// the atoms to be set virtual int load_atoms(char const *filename, - std::vector &atoms, - std::string const &pdb_field, - double const pdb_field_value = 0.0) = 0; + std::vector &atoms, + std::string const &pdb_field, + double const pdb_field_value = 0.0) = 0; /// \brief Load the coordinates for a group of atoms from a file /// (usually a PDB); if "pos" is already allocated, the number of its /// elements must match the number of atoms in "filename" virtual int load_coords(char const *filename, - std::vector &pos, - const std::vector &indices, - std::string const &pdb_field, - double const pdb_field_value = 0.0) = 0; + std::vector &pos, + const std::vector &indices, + std::string const &pdb_field, + double const pdb_field_value = 0.0) = 0; + +protected: + + /// \brief Open output files: by default, these are ofstream objects. + /// Allows redefinition to implement different output mechanisms + std::list output_files; + /// \brief Identifiers for output_stream objects: by default, these are the names of the files + std::list output_stream_names; + +public: + + // TODO the following definitions may be moved to a .cpp file + + /// \brief Returns a reference to the given output channel; + /// if this is not open already, then open it + virtual std::ostream * output_stream(std::string const &output_name) + { + std::list::iterator osi = output_files.begin(); + std::list::iterator osni = output_stream_names.begin(); + for ( ; osi != output_files.end(); osi++, osni++) { + if (*osni == output_name) { + return *osi; + } + } + output_stream_names.push_back(output_name); + std::ofstream * os = new std::ofstream(output_name.c_str()); + if (!os->is_open()) { + cvm::error("Error: cannot write to file \""+output_name+"\".\n", + FILE_ERROR); + } + output_files.push_back(os); + return os; + } + + /// \brief Closes the given output channel + virtual int close_output_stream(std::string const &output_name) + { + std::list::iterator osi = output_files.begin(); + std::list::iterator osni = output_stream_names.begin(); + for ( ; osi != output_files.end(); osi++, osni++) { + if (*osni == output_name) { + ((std::ofstream *) (*osi))->close(); + output_files.erase(osi); + output_stream_names.erase(osni); + return COLVARS_OK; + } + } + return COLVARS_ERROR; + } /// \brief Rename the given file, before overwriting it virtual int backup_file(char const *filename) - { return COLVARS_NOT_IMPLEMENTED; } + { + return COLVARS_NOT_IMPLEMENTED; + } }; inline void colvarproxy::select_closest_images(std::vector &pos, - cvm::atom_pos const &ref_pos) + cvm::atom_pos const &ref_pos) { for (std::vector::iterator pi = pos.begin(); pi != pos.end(); ++pi) { @@ -222,7 +275,7 @@ inline void colvarproxy::select_closest_images(std::vector &pos, } inline cvm::real colvarproxy::position_dist2(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) + cvm::atom_pos const &pos2) { return (position_distance(pos1, pos2)).norm2(); } diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index 61cfdf2695..aeb2a9facc 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -324,7 +324,6 @@ protected: public: T * data; size_t length; - friend class matrix2d; inline row(T * const row_data, size_t const inner_length) : data(row_data), length(inner_length) {}