diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 10898f0fa9..34b171e088 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -129,6 +129,25 @@ void cvm::atom_group::parse (std::string const &conf, atom_indexes.clear(); } + + std::string index_group_name; + if (get_keyval (group_conf, "indexGroup", index_group_name)) { + // use an index group from the index file read globally + std::list::iterator names_i = cvm::index_group_names.begin(); + std::list >::iterator index_groups_i = cvm::index_groups.begin(); + for ( ; names_i != cvm::index_group_names.end() ; names_i++, index_groups_i++) { + if (*names_i == index_group_name) + break; + } + if (names_i == cvm::index_group_names.end()) { + cvm::fatal_error ("Error: could not find index group "+ + index_group_name+" among those provided by the index file.\n"); + } + this->reserve (index_groups_i->size()); + for (size_t i = 0; i < index_groups_i->size(); i++) { + this->push_back (cvm::atom ((*index_groups_i)[i])); + } + } } { diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 561bb4a688..dd6eb67297 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -41,10 +41,11 @@ colvarbias::colvarbias (std::string const &conf, char const *key) add_colvar (colvars_str[i]); } } - if (!colvars.size()) { cvm::fatal_error ("Error: no collective variables specified.\n"); } + + get_keyval (conf, "outputEnergy", b_output_energy, false); } @@ -82,6 +83,40 @@ void colvarbias::communicate_forces() } +void colvarbias::change_configuration(std::string const &conf) +{ + cvm::fatal_error ("Error: change_configuration() not implemented.\n"); +} + + +cvm::real colvarbias::energy_difference(std::string const &conf) +{ + cvm::fatal_error ("Error: energy_difference() not implemented.\n"); + return 0.; +} + + +std::ostream & colvarbias::write_traj_label (std::ostream &os) +{ + os << " "; + if (b_output_energy) + os << " E_" + << cvm::wrap_string (this->name, cvm::en_width-2); + return os; +} + + +std::ostream & colvarbias::write_traj (std::ostream &os) +{ + os << " "; + if (b_output_energy) + os << " " + << bias_energy; + return os; +} + + + colvarbias_harmonic::colvarbias_harmonic (std::string const &conf, char const *key) @@ -171,25 +206,16 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf, } } + get_keyval (conf, "outputCenters", b_output_centers, false); + get_keyval (conf, "outputAccumulatedWork", b_output_acc_work, false); + acc_work = 0.0; + if (cvm::debug()) cvm::log ("Done initializing a new harmonic restraint bias.\n"); } -void colvarbias::change_configuration(std::string const &conf) -{ - cvm::fatal_error ("Error: change_configuration() not implemented.\n"); -} - - -cvm::real colvarbias::energy_difference(std::string const &conf) -{ - cvm::fatal_error ("Error: energy_difference() not implemented.\n"); - return 0.; -} - - -void colvarbias_harmonic::change_configuration(std::string const &conf) +void colvarbias_harmonic::change_configuration (std::string const &conf) { get_keyval (conf, "forceConstant", force_k, force_k); if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) { @@ -379,6 +405,15 @@ cvm::real colvarbias_harmonic::update() colvar_centers[i]))+"\n"); } + if (b_output_acc_work) { + if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) { + for (size_t i = 0; i < colvars.size(); i++) { + // project forces on the calculated increments at this step + acc_work += colvar_forces[i] * centers_incr[i]; + } + } + } + if (cvm::debug()) cvm::log ("Current forces for the harmonic bias \""+ this->name+"\": "+cvm::to_str (colvar_forces)+".\n"); @@ -442,6 +477,11 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is) cvm::fatal_error ("Error: current stage is missing from the restart.\n"); } + if (b_output_acc_work) { + if (!get_keyval (conf, "accumulatedWork", acc_work)) + cvm::fatal_error ("Error: accumulatedWork is missing from the restart.\n"); + } + is >> brace; if (brace != "}") { cvm::fatal_error ("Error: corrupt restart information for harmonic bias \""+ @@ -484,9 +524,61 @@ std::ostream & colvarbias_harmonic::write_restart (std::ostream &os) << stage << "\n"; } + if (b_output_acc_work) { + os << " accumulatedWork " << acc_work << "\n"; + } + os << " }\n" << "}\n\n"; return os; } + +std::ostream & colvarbias_harmonic::write_traj_label (std::ostream &os) +{ + os << " "; + + if (b_output_energy) + os << " E_" + << cvm::wrap_string (this->name, cvm::en_width-2); + + if (b_output_centers) + for (size_t i = 0; i < colvars.size(); i++) { + size_t const this_cv_width = (colvars[i]->value()).output_width (cvm::cv_width); + os << " x0_" + << cvm::wrap_string (colvars[i]->name, this_cv_width-3); + } + + if (b_output_acc_work) + os << " W_" + << cvm::wrap_string (this->name, cvm::en_width-2); + + return os; +} + + +std::ostream & colvarbias_harmonic::write_traj (std::ostream &os) +{ + os << " "; + + if (b_output_energy) + os << " " + << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) + << bias_energy; + + if (b_output_centers) + for (size_t i = 0; i < colvars.size(); i++) { + os << " " + << std::setprecision (cvm::cv_prec) << std::setw (cvm::cv_width) + << colvar_centers[i]; + } + + if (b_output_acc_work) + os << " " + << std::setprecision (cvm::en_prec) << std::setw (cvm::en_width) + << acc_work; + + return os; +} + diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index 2130c4ba1c..f12f80b901 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -52,6 +52,13 @@ public: /// Write the bias configuration to a restart file virtual std::ostream & write_restart (std::ostream &os) = 0; + /// Write a label to the trajectory file (comment line) + virtual std::ostream & write_traj_label (std::ostream &os); + + /// Output quantities such as the bias energy to the trajectory file + virtual std::ostream & write_traj (std::ostream &os); + + protected: /// \brief Pointers to collective variables to which the bias is @@ -62,10 +69,12 @@ protected: /// \brief Current forces from this bias to the colvars std::vector colvar_forces; - /// \brief Current energy of this bias (colvar_forces should be - /// obtained by deriving this) + /// \brief Current energy of this bias (colvar_forces should be obtained by deriving this) cvm::real bias_energy; + /// Whether to write the current bias energy from this bias to the trajectory file + bool b_output_energy; + /// \brief Whether this bias has already accumulated information /// (when relevant) bool has_data; @@ -94,6 +103,12 @@ public: /// Write the bias configuration to a restart file virtual std::ostream & write_restart (std::ostream &os); + /// Write a label to the trajectory file (comment line) + virtual std::ostream & write_traj_label (std::ostream &os); + + /// Output quantities such as the bias energy to the trajectory file + virtual std::ostream & write_traj (std::ostream &os); + /// \brief Constructor colvarbias_harmonic (std::string const &conf, char const *key); @@ -109,27 +124,9 @@ protected: /// \brief Restraint centers without wrapping or constraints applied std::vector colvar_centers_raw; - /// \brief Restraint force constant - cvm::real force_k; - /// \brief Moving target? bool b_chg_centers; - /// \brief Changing force constant? - bool b_chg_force_k; - - /// \brief Restraint force constant (target value) - cvm::real target_force_k; - - /// \brief Equilibration steps for restraint FE calculation through TI - cvm::real target_equil_steps; - - /// \brief Restraint force constant (starting value) - cvm::real starting_force_k; - - /// \brief Lambda-schedule for custom varying force constant - std::vector lambda_schedule; - /// \brief New restraint centers std::vector target_centers; @@ -137,12 +134,41 @@ protected: /// (or stage) towards the new values (calculated from target_nsteps) std::vector centers_incr; + /// Whether to write the current restraint centers to the trajectory file + bool b_output_centers; + + /// Whether to write the current accumulated work to the trajectory file + bool b_output_acc_work; + + /// \brief Accumulated work + cvm::real acc_work; + + + /// \brief Restraint force constant + cvm::real force_k; + + /// \brief Changing force constant? + bool b_chg_force_k; + + /// \brief Restraint force constant (target value) + cvm::real target_force_k; + + /// \brief Restraint force constant (starting value) + cvm::real starting_force_k; + + /// \brief Lambda-schedule for custom varying force constant + std::vector lambda_schedule; + /// \brief Exponent for varying the force constant cvm::real force_k_exp; + + /// \brief Intermediate quantity to compute the restraint free energy + /// (in TI, would be the accumulating FE derivative) + cvm::real restraint_FE; - /// \brief Number of steps required to reach the target force constant - /// or restraint centers - size_t target_nsteps; + + /// \brief Equilibration steps for restraint FE calculation through TI + cvm::real target_equil_steps; /// \brief Number of stages over which to perform the change /// If zero, perform a continuous change @@ -150,10 +176,10 @@ protected: /// \brief Number of current stage of the perturbation int stage; - - /// \brief Intermediate quantity to compute the restraint free energy - /// (in TI, would be the accumulating FE derivative) - cvm::real restraint_FE; + + /// \brief Number of steps required to reach the target force constant + /// or restraint centers + size_t target_nsteps; }; diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 39a6d12954..5f344ea539 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -44,6 +44,11 @@ colvarmodule::colvarmodule (char const *config_filename, config_s.close(); } + std::string index_file_name; + if (parse->get_keyval (conf, "indexFile", index_file_name)) { + read_index_file (index_file_name.c_str()); + } + parse->get_keyval (conf, "analysis", b_analysis, false); parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07, @@ -478,6 +483,11 @@ void colvarmodule::calc() { cvi++) { (*cvi)->write_traj_label (cv_traj_os); } + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + (*bi)->write_traj_label (cv_traj_os); + } cv_traj_os << "\n"; if (cvm::debug()) cv_traj_os.flush(); @@ -494,6 +504,11 @@ void colvarmodule::calc() { cvi++) { (*cvi)->write_traj (cv_traj_os); } + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + (*bi)->write_traj (cv_traj_os); + } cv_traj_os << "\n"; if (cvm::debug()) cv_traj_os.flush(); @@ -564,7 +579,6 @@ colvarmodule::~colvarmodule() } delete parse; - proxy = NULL; } @@ -721,6 +735,55 @@ void cvm::exit (std::string const &message) } +void cvm::read_index_file (char const *filename) +{ + std::ifstream is (filename); + if (!is.good()) + fatal_error ("Error: in opening index file \""+ + std::string (filename)+"\".\n"); + // std::list::iterator names_i = cvm::index_group_names.begin(); + // std::list >::iterator lists_i = cvm::index_groups.begin(); + while (is.good()) { + char open, close; + std::string group_name; + if ( (is >> open) && (open == '[') && + (is >> group_name) && + (is >> close) && (close == ']') ) { + cvm::index_group_names.push_back (group_name); + cvm::index_groups.push_back (std::vector ()); + } else { + cvm::fatal_error ("Error: in parsing index file \""+ + std::string (filename)+"\".\n"); + } + + int atom_number = 1; + size_t pos = is.tellg(); + while ( (is >> atom_number) && (atom_number > 0) ) { + (cvm::index_groups.back()).push_back (atom_number); + pos = is.tellg(); + } + is.clear(); + is.seekg (pos, std::ios::beg); + std::string delim; + if ( (is >> delim) && (delim == "[") ) { + // new group + is.clear(); + is.seekg (pos, std::ios::beg); + } else { + break; + } + } + + cvm::log ("The following index groups were read from the index file \""+ + std::string (filename)+"\":\n"); + std::list::iterator names_i = cvm::index_group_names.begin(); + std::list >::iterator lists_i = cvm::index_groups.begin(); + for ( ; names_i != cvm::index_group_names.end() ; names_i++, lists_i++) { + cvm::log (" "+(*names_i)+" ("+cvm::to_str (lists_i->size())+" atoms).\n"); + } + +} + // static pointers std::vector colvarmodule::colvars; @@ -741,6 +804,8 @@ size_t colvarmodule::cv_traj_freq = 0; size_t colvarmodule::depth = 0; bool colvarmodule::b_analysis = false; cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-04; +std::list colvarmodule::index_group_names; +std::list > colvarmodule::index_groups; // file name prefixes @@ -932,15 +997,13 @@ cvm::quaternion::position_derivative_inner (cvm::rvector const &pos, + + // Calculate the optimal rotation between two groups, and implement it // as a quaternion. The method is the one documented in: Coutsias EA, // Seok C, Dill KA. Using quaternions to calculate RMSD. J Comput // Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254 - - - - void colvarmodule::rotation::build_matrix (std::vector const &pos1, std::vector const &pos2, matrix2d &S) diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 78c93b34e7..01f8d78cd2 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -2,7 +2,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2013-04-23" +#define COLVARS_VERSION "2013-05-14" #endif #ifndef COLVARS_DEBUG @@ -290,6 +290,16 @@ public: atom_pos const &ref_pos); + /// \brief Names of groups from a Gromacs .ndx file to be read at startup + static std::list index_group_names; + + /// \brief Groups from a Gromacs .ndx file read at startup + static std::list > index_groups; + + /// \brief Read a Gromacs .ndx file + static void read_index_file (char const *filename); + + /// \brief Create atoms from a file \param filename name of the file /// (usually a PDB) \param atoms array of the atoms to be allocated /// \param pdb_field (optiona) if "filename" is a PDB file, use this diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index f7d6aa6bf6..25f888c001 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -242,7 +242,7 @@ std::istream & operator >> (std::istream &is, colvarvalue &x) } -size_t colvarvalue::output_width (size_t const &real_width) +size_t colvarvalue::output_width (size_t const &real_width) const { switch (this->value_type) { case colvarvalue::type_scalar: diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index 134585af00..eb44593454 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -265,10 +265,10 @@ public: /// with a different type to this object void error_rside (Type const &vt) const; - ///�Give the number of characters required to output this + /// Give the number of characters required to output this /// colvarvalue, given the current type assigned and the number of /// characters for a real number - size_t output_width (size_t const &real_width); + size_t output_width (size_t const &real_width) const; // optimized routines for operations with an array; xv and inner as