From 4da20dce992f7366bb6b0e29b0055f61f2b96720 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 27 Jun 2013 22:48:27 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10118 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/colvars/Makefile.g++ | 14 +-- lib/colvars/Makefile.mingw32-cross | 70 ++++++------ lib/colvars/README | 17 ++- lib/colvars/colvar.cpp | 81 ++++++++------ lib/colvars/colvar.h | 7 +- lib/colvars/colvaratoms.cpp | 57 +++++++--- lib/colvars/colvaratoms.h | 25 +++-- lib/colvars/colvarbias.cpp | 136 ++++++++++++++++++++---- lib/colvars/colvarbias.h | 86 +++++++++------ lib/colvars/colvarbias_abf.cpp | 60 ++++++++--- lib/colvars/colvarbias_abf.h | 6 +- lib/colvars/colvarbias_meta.cpp | 72 +++++++------ lib/colvars/colvarbias_meta.h | 14 +-- lib/colvars/colvarcomp.cpp | 4 +- lib/colvars/colvarcomp.h | 153 ++++++++++++++------------- lib/colvars/colvarcomp_angles.cpp | 4 +- lib/colvars/colvarcomp_coordnums.cpp | 67 ++---------- lib/colvars/colvarcomp_distances.cpp | 68 ++++++------ lib/colvars/colvarcomp_protein.cpp | 16 +-- lib/colvars/colvarcomp_rotations.cpp | 4 +- lib/colvars/colvargrid.cpp | 4 +- lib/colvars/colvargrid.h | 54 +++++----- lib/colvars/colvarmodule.cpp | 153 ++++++++++++++++++++------- lib/colvars/colvarmodule.h | 48 ++++++--- lib/colvars/colvarparse.cpp | 34 +++--- lib/colvars/colvarparse.h | 10 +- lib/colvars/colvarproxy.h | 78 ++++++++------ lib/colvars/colvartypes.h | 74 ++++++------- lib/colvars/colvarvalue.cpp | 14 +-- lib/colvars/colvarvalue.h | 24 ++--- 30 files changed, 883 insertions(+), 571 deletions(-) diff --git a/lib/colvars/Makefile.g++ b/lib/colvars/Makefile.g++ index 825e24a862..272d93b801 100644 --- a/lib/colvars/Makefile.g++ +++ b/lib/colvars/Makefile.g++ @@ -1,4 +1,4 @@ -# library build makefile for colvars module +# library build -*- makefile -*- for colvars module # which file will be copied to Makefile.lammps @@ -14,9 +14,9 @@ SHELL = /bin/sh # ------ DEFINITIONS ------ -SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ - colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ - colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ +SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ + colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ + colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp LIB = libcolvars.a @@ -25,11 +25,13 @@ EXE = #colvars_standalone # ------ MAKE PROCEDURE ------ -default: $(LIB) $(EXE) +default: $(LIB) $(EXE) Makefile.lammps + +Makefile.lammps: + @cp $(EXTRAMAKE) Makefile.lammps $(LIB): $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) - @cp $(EXTRAMAKE) Makefile.lammps colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) $(CXX) -o $@ $(CXXFLAGS) $^ diff --git a/lib/colvars/Makefile.mingw32-cross b/lib/colvars/Makefile.mingw32-cross index 42435e9f4c..05bdbf32a5 100644 --- a/lib/colvars/Makefile.mingw32-cross +++ b/lib/colvars/Makefile.mingw32-cross @@ -1,4 +1,4 @@ -# library build makefile for colvars module +# library build -*- makefile -*- for colvars module # which file will be copied to Makefile.lammps @@ -6,35 +6,42 @@ EXTRAMAKE = Makefile.lammps.empty # ------ SETTINGS ------ -CXX = i686-pc-mingw32-g++ +CXX = i686-w64-mingw32-g++ CXXFLAGS = -O2 -march=i686 -mtune=generic -mfpmath=387 -mpc64 \ -fno-rtti -fno-exceptions -finline-functions \ -ffast-math -funroll-loops -fstrict-aliasing \ -Wall -W -Wno-uninitialized -ARCHIVE = i686-pc-mingw32-ar +ARCHIVE = i686-w64-mingw32-ar ARCHFLAG = -rscv SHELL = /bin/sh # ------ DEFINITIONS ------ -SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ - colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ - colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ +SRC = colvaratoms.cpp colvarbias_abf.cpp colvarbias.cpp colvarbias_meta.cpp \ + colvar.cpp colvarcomp_angles.cpp colvarcomp.cpp colvarcomp_coordnums.cpp \ + colvarcomp_distances.cpp colvarcomp_protein.cpp colvarcomp_rotations.cpp \ colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp -LIB = libcolvars.a -OBJ = $(SRC:.cpp=.o) +DIR = Obj_mingw32/ +LIB = $(DIR)libcolvars.a +OBJ = $(SRC:%.cpp=$(DIR)%.o) EXE = #colvars_standalone # ------ MAKE PROCEDURE ------ -default: $(LIB) $(EXE) +default: $(DIR) $(LIB) $(EXE) Makefile.lammps -$(LIB): $(OBJ) +$(DIR): + mkdir $(DIR) + +Makefile.lammps: + @cp $(EXTRAMAKE) Makefile.lammps + +$(LIB): $(DIR) $(OBJ) $(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ) @cp $(EXTRAMAKE) Makefile.lammps -colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) +$(DIR)colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) $(CXX) -o $@ $(CXXFLAGS) $^ # ------ MAKE FLAGS ------ @@ -46,58 +53,59 @@ colvars_standalone: colvars_main.o colvarproxy_standalone.o $(LIB) # ------ COMPILE RULES ------ -.cpp.o: - $(CXX) $(CXXFLAGS) -c $< +$(DIR)%.o: %.cpp + $(CXX) $(CXXFLAGS) -c $< -o $@ # ------ DEPENDENCIES ------ # -colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvars_main.o: colvars_main.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarproxy_standalone.h colvaratoms.h colvarparse.h colvarvalue.h -colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \ +$(DIR)colvarproxy_standalone.o: colvarproxy_standalone.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvaratoms.h colvarparse.h colvarvalue.h \ colvarproxy_standalone.h -colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarparse.h colvarvalue.h colvaratoms.h -colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ +$(DIR)colvarbias_abf.o: colvarbias_abf.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarbias_abf.h \ colvarbias.h colvargrid.h -colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvarbias.o: colvarbias.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarbias.h colvar.h colvarparse.h -colvarbias_meta.o: colvarbias_meta.cpp colvar.h colvarmodule.h \ +$(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 -colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvar.o: colvar.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h colvar.h colvarcomp.h colvaratoms.h -colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \ +$(DIR)colvarcomp_angles.o: colvarcomp_angles.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvar.h colvarvalue.h colvarparse.h colvarcomp.h \ colvaratoms.h -colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvar.h colvarparse.h colvarcomp.h colvaratoms.h -colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \ +$(DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarparse.h colvarvalue.h colvaratoms.h \ colvar.h colvarcomp.h -colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \ +$(DIR)colvarcomp_distances.o: colvarcomp_distances.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h -colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \ +$(DIR)colvarcomp_protein.o: colvarcomp_protein.cpp colvarmodule.h colvartypes.h \ colvarproxy.h colvarvalue.h colvarparse.h colvar.h colvarcomp.h \ colvaratoms.h -colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \ +$(DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp colvarmodule.h \ colvartypes.h colvarproxy.h colvarvalue.h colvarparse.h colvar.h \ colvarcomp.h colvaratoms.h -colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)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 \ +$(DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarparse.h colvarvalue.h colvar.h colvarbias.h colvarbias_meta.h \ colvargrid.h colvarbias_abf.h -colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h -colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ +$(DIR)colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h # ------ CLEAN ------ clean: - -rm *.o *~ $(LIB) + -rm $(DIR)*.o *~ $(LIB) + -rmdir $(DIR) diff --git a/lib/colvars/README b/lib/colvars/README index 9c6b1f6c09..d6efc333a5 100644 --- a/lib/colvars/README +++ b/lib/colvars/README @@ -4,9 +4,24 @@ that allows enhanced sampling in molecular dynamics simulations. The module is written to maximize performance, portability, flexibility of usage for the user, and extensibility for the developer. -The following publication describes the principles of +The development of the colvars library is now hosted on github at: +http://colvars.github.io/ +You can use this site to get access to the latest development sources +and the up-to-date documentation. + +Copy of the specific documentation is also in + doc/PDF/colvars-refman-lammps.pdf + +Please report bugs and request new features at: +https://github.com/colvars/colvars/issues + +The following publications describe the principles of the implementation of this library: + Using collective variables to drive molecular dynamics simulations, + Giacomo Fiorin , Michael L. Klein & Jérôme Hénin (2013): + Molecular Physics DOI:10.1080/00268976.2013.813594 + Exploring Multidimensional Free Energy Landscapes Using Time-Dependent Biases on Collective Variables, J. Hénin, G. Fiorin, C. Chipot, and M. L. Klein, diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 1c3d9c8215..f6652852eb 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -14,7 +14,7 @@ colvar::colvar (std::string const &conf) { cvm::log ("Initializing a new collective variable.\n"); - + get_keyval (conf, "name", this->name, (std::string ("colvar")+cvm::to_str (cvm::colvars.size()+1))); @@ -26,7 +26,7 @@ colvar::colvar (std::string const &conf) "\", as another colvar.\n"); } - // all tasks disabled by default + // all tasks disabled by default for (size_t i = 0; i < task_ntot; i++) { tasks[i] = false; } @@ -152,7 +152,7 @@ colvar::colvar (std::string const &conf) // Decide whether the colvar is periodic // Used to wrap extended DOF if extendedLagrangian is on - if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1 + if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1 && (cvcs[0])->sup_coeff == 1.0 ) { this->b_periodic = true; this->period = (cvcs[0])->period; @@ -167,6 +167,9 @@ colvar::colvar (std::string const &conf) // check the available features of each cvc for (size_t i = 0; i < cvcs.size(); i++) { + if ((cvcs[i])->b_debug_gradients) + enable (task_gradients); + if ((cvcs[i])->sup_np != 1) { if (cvm::debug() && b_linear) cvm::log ("Warning: You are using a non-linear polynomial " @@ -307,7 +310,7 @@ colvar::colvar (std::string const &conf) cvm::log ("Enabling the extended Lagrangian term for colvar \""+ this->name+"\".\n"); - + enable (task_extended_lagrangian); xr.type (this->type()); @@ -525,10 +528,10 @@ void colvar::enable (colvar::task const &t) if ( !tasks[task_extended_lagrangian] ) { enable (task_gradients); - if (!b_Jacobian_force) + if (!b_Jacobian_force) cvm::fatal_error ("Error: colvar \""+this->name+ "\" does not have Jacobian forces implemented.\n"); - if (!b_linear) + if (!b_linear) cvm::fatal_error ("Error: colvar \""+this->name+ "\" must be defined as a linear combination " "to calculate the Jacobian force.\n"); @@ -664,12 +667,24 @@ void colvar::disable (colvar::task const &t) } +void colvar::setup() { + // loop over all components to reset masses of all groups + for (size_t i = 0; i < cvcs.size(); i++) { + for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { + cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); + atoms.read_positions(); + atoms.reset_mass(name,i,ig); + } + } +} + + colvar::~colvar() { for (size_t i = 0; i < cvcs.size(); i++) { delete cvcs[i]; } -} +} @@ -682,6 +697,8 @@ void colvar::calc() cvm::log ("Calculating colvar \""+this->name+"\".\n"); // prepare atom groups for calculation + if (cvm::debug()) + cvm::log ("Collecting data from atom groups.\n"); for (size_t i = 0; i < cvcs.size(); i++) { for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); @@ -697,19 +714,21 @@ void colvar::calc() for (size_t i = 0; i < cvcs.size(); i++) { for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvcs[i]->atom_groups[ig]->read_velocities(); - } + } } } if (tasks[task_system_force]) { for (size_t i = 0; i < cvcs.size(); i++) { for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvcs[i]->atom_groups[ig]->read_system_forces(); - } + } } } // calculate the value of the colvar + if (cvm::debug()) + cvm::log ("Calculating colvar components.\n"); x.reset(); if (x.type() == colvarvalue::type_scalar) { // polynomial combination allowed @@ -727,7 +746,7 @@ void colvar::calc() ( ((cvcs[i])->sup_np != 1) ? std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np) : (cvcs[i])->value().real_value ); - } + } } else { // only linear combination allowed @@ -741,7 +760,7 @@ void colvar::calc() cvm::to_str ((cvcs[i])->value(), cvm::cv_width, cvm::cv_prec)+".\n"); x += (cvcs[i])->sup_coeff * (cvcs[i])->value(); - } + } } if (cvm::debug()) @@ -762,9 +781,9 @@ void colvar::calc() // if requested, propagate (via chain rule) the gradients above // to the atoms used to define the roto-translation for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { - if (cvcs[i]->atom_groups[ig]->b_fit_gradients) + if (cvcs[i]->atom_groups[ig]->b_fit_gradients) cvcs[i]->atom_groups[ig]->calc_fit_gradients(); - } + } cvm::decrease_depth(); } @@ -801,7 +820,7 @@ void colvar::calc() for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { int a = std::lower_bound (atom_ids.begin(), atom_ids.end(), cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin(); - atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad; + atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad; } } } @@ -925,7 +944,7 @@ cvm::real colvar::update() } if (tasks[task_Jacobian_force]) { - + cvm::increase_depth(); for (size_t i = 0; i < cvcs.size(); i++) { (cvcs[i])->calc_Jacobian_derivative(); @@ -941,8 +960,8 @@ cvm::real colvar::update() fj *= cvm::boltzmann() * cvm::temperature(); // the instantaneous Jacobian force was not included in the reported system force; - // instead, it is subtracted from the applied force (silent Jacobian correction) - if (! tasks[task_report_Jacobian_force]) + // instead, it is subtracted from the applied force (silent Jacobian correction) + if (! tasks[task_report_Jacobian_force]) f -= fj; } @@ -966,9 +985,9 @@ cvm::real colvar::update() // leap to v_(i+1/2) if (tasks[task_langevin]) { vr -= dt * ext_gamma * vr.real_value; - vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass; - } - vr += (0.5 * dt) * fr / ext_mass; + vr += dt * ext_sigma * cvm::rand_gaussian() / ext_mass; + } + vr += (0.5 * dt) * fr / ext_mass; xr += dt * vr; xr.apply_constraints(); if (this->b_periodic) this->wrap (xr); @@ -989,13 +1008,13 @@ cvm::real colvar::update() void colvar::communicate_forces() { if (cvm::debug()) - cvm::log ("Communicating forces from colvar \""+this->name+"\".\n"); + cvm::log ("Communicating forces from colvar \""+this->name+"\".\n"); if (x.type() == colvarvalue::type_scalar) { for (size_t i = 0; i < cvcs.size(); i++) { cvm::increase_depth(); - (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff * + (cvcs[i])->apply_force (f * (cvcs[i])->sup_coeff * cvm::real ((cvcs[i])->sup_np) * (std::pow ((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1)) ); @@ -1012,7 +1031,7 @@ void colvar::communicate_forces() } if (cvm::debug()) - cvm::log ("Done communicating forces from colvar \""+this->name+"\".\n"); + cvm::log ("Done communicating forces from colvar \""+this->name+"\".\n"); } @@ -1200,7 +1219,7 @@ std::istream & colvar::read_traj (std::istream &is) is >> ft; - if (tasks[task_extended_lagrangian]) { + if (tasks[task_extended_lagrangian]) { is >> fr; ft_reported = fr; } else { @@ -1248,7 +1267,7 @@ std::ostream & colvar::write_restart (std::ostream &os) { os << "}\n\n"; return os; -} +} std::ostream & colvar::write_traj_label (std::ostream & os) @@ -1418,7 +1437,7 @@ inline void history_add_value (size_t const &history_length, inline void history_incr (std::list< std::list > &history, std::list< std::list >::iterator &history_p) { - if ((++history_p) == history.end()) + if ((++history_p) == history.end()) history_p = history.begin(); } @@ -1435,7 +1454,7 @@ void colvar::calc_acf() // first-step operations - colvar *cfcv = (acf_colvar_name.size() ? + colvar *cfcv = (acf_colvar_name.size() ? cvm::colvar_p (acf_colvar_name) : this); if (cfcv->type() != this->type()) @@ -1474,7 +1493,7 @@ void colvar::calc_acf() } else { - colvar *cfcv = (acf_colvar_name.size() ? + colvar *cfcv = (acf_colvar_name.size() ? cvm::colvar_p (acf_colvar_name) : this); @@ -1539,7 +1558,7 @@ void colvar::calc_vel_acf (std::list &v_list, // inner products of previous velocities with current (acf_i and // vs_i are updated) - colvarvalue::inner_opt (v, vs_i, v_list.end(), acf_i); + colvarvalue::inner_opt (v, vs_i, v_list.end(), acf_i); acf_nframes++; } @@ -1559,7 +1578,7 @@ void colvar::calc_coor_acf (std::list &x_list, *(acf_i++) += x.norm2(); - colvarvalue::inner_opt (x, xs_i, x_list.end(), acf_i); + colvarvalue::inner_opt (x, xs_i, x_list.end(), acf_i); acf_nframes++; } @@ -1581,7 +1600,7 @@ void colvar::calc_p2coor_acf (std::list &x_list, // value of P2(0) = 1 *(acf_i++) += 1.0; - colvarvalue::p2leg_opt (x, xs_i, x_list.end(), acf_i); + colvarvalue::p2leg_opt (x, xs_i, x_list.end(), acf_i); acf_nframes++; } diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index c16355dd16..6f6692e904 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -215,7 +215,7 @@ protected: cvm::real ext_gamma; /// Amplitude of Gaussian white noise for Langevin extended dynamics cvm::real ext_sigma; - + /// \brief Harmonic restraint force colvarvalue fr; @@ -288,6 +288,9 @@ public: /// Disable the specified task void disable (colvar::task const &t); + /// Get ready for a run and possibly re-initialize internal data + void setup(); + /// Destructor ~colvar(); @@ -380,7 +383,7 @@ protected: colvarvalue x_old; /// Time series of values and velocities used in correlation - /// functions + /// functions std::list< std::list > acf_x_history, acf_v_history; /// Time series of values and velocities used in correlation /// functions (pointers)x diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 10898f0fa9..c211168ff2 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -26,8 +26,8 @@ cvm::atom_group::atom_group (std::string const &conf, cvm::atom_group::atom_group (std::vector const &atoms) - : b_dummy (false), b_center (false), b_rotate (false), - ref_pos_group (NULL), b_fit_gradients (false), + : b_dummy (false), b_center (false), b_rotate (false), + ref_pos_group (NULL), b_fit_gradients (false), noforce (false) { this->reserve (atoms.size()); @@ -43,8 +43,8 @@ cvm::atom_group::atom_group (std::vector const &atoms) cvm::atom_group::atom_group() - : b_dummy (false), b_center (false), b_rotate (false), - ref_pos_group (NULL), b_fit_gradients (false), + : b_dummy (false), b_center (false), b_rotate (false), + ref_pos_group (NULL), b_fit_gradients (false), noforce (false) { total_mass = 0.0; @@ -71,6 +71,18 @@ void cvm::atom_group::add_atom (cvm::atom const &a) } +void cvm::atom_group::reset_mass(std::string &name, int i, int j) +{ + total_mass = 0.0; + for (cvm::atom_iter ai = this->begin(); + ai != this->end(); ai++) { + total_mass += ai->mass; + } + cvm::log ("Re-initialized atom group "+name+":"+cvm::to_str (i)+"/"+ + cvm::to_str (j)+". "+ cvm::to_str (this->size())+ + " atoms: total mass = "+cvm::to_str (this->total_mass)+".\n"); +} + void cvm::atom_group::parse (std::string const &conf, char const *key) { @@ -129,6 +141,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])); + } + } } { @@ -234,7 +265,7 @@ void cvm::atom_group::parse (std::string const &conf, } } - for (std::vector::iterator a1 = this->begin(); + for (std::vector::iterator a1 = this->begin(); a1 != this->end(); a1++) { std::vector::iterator a2 = a1; ++a2; @@ -253,10 +284,10 @@ void cvm::atom_group::parse (std::string const &conf, if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) { b_dummy = true; this->total_mass = 1.0; - } else + } else b_dummy = false; - if (b_dummy && (this->size())) + if (b_dummy && (this->size())) cvm::fatal_error ("Error: cannot set up group \""+ std::string (key)+"\" as a dummy atom " "and provide it with atom definitions.\n"); @@ -328,7 +359,7 @@ void cvm::atom_group::parse (std::string const &conf, std::string ref_pos_col; double ref_pos_col_value; - + if (get_keyval (group_conf, "refPositionsCol", ref_pos_col, std::string (""), mode)) { // if provided, use PDB column to select coordinates bool found = get_keyval (group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode); @@ -486,7 +517,7 @@ void cvm::atom_group::calc_apply_roto_translation() } } -void cvm::atom_group::apply_translation (cvm::rvector const &t) +void cvm::atom_group::apply_translation (cvm::rvector const &t) { if (b_dummy) return; @@ -496,7 +527,7 @@ void cvm::atom_group::apply_translation (cvm::rvector const &t) } } -void cvm::atom_group::apply_rotation (cvm::rotation const &rot) +void cvm::atom_group::apply_rotation (cvm::rotation const &rot) { if (b_dummy) return; @@ -617,11 +648,11 @@ void cvm::atom_group::calc_fit_gradients() } if (b_rotate) { - + // add the rotation matrix contribution to the gradients cvm::rotation const rot_inv = rot.inverse(); cvm::atom_pos const cog = this->center_of_geometry(); - + for (size_t i = 0; i < this->size(); i++) { cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos))); @@ -785,7 +816,7 @@ void cvm::atom_group::apply_force (cvm::rvector const &force) for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force ((ai->mass/this->total_mass) * force); - } + } } } diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index d2be1dfe6e..e4d95537f8 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -8,7 +8,7 @@ /// \brief Stores numeric id, mass and all mutable data for an atom, /// mostly used by a \link cvc \endlink -/// +/// /// This class may be used (although not necessarily) to keep atomic /// data (id, mass, position and collective variable derivatives) /// altogether. There may be multiple instances with identical @@ -23,7 +23,7 @@ protected: /// \brief Index in the list of atoms involved by the colvars (\b /// NOT in the global topology!) - size_t index; + int index; public: @@ -47,7 +47,7 @@ public: /// \brief Gradient of a scalar collective variable with respect /// to this atom - /// + /// /// This can only handle a scalar collective variable (i.e. when /// the \link colvarvalue::real_value \endlink member is used /// from the \link colvarvalue \endlink class), which is also the @@ -57,8 +57,8 @@ public: /// implementation cvm::rvector grad; - /// \brief Default constructor, setting id to a non-valid value - inline atom() {} + /// \brief Default constructor, setting id and index to invalid numbers + atom() : id (-1), index (-1) { reset_data(); } /// \brief Initialize an atom for collective variable calculation /// and get its internal identifier \param atom_number Atom index in @@ -121,7 +121,7 @@ class colvarmodule::atom_group public colvarparse { public: - // Note: all members here are kept public, to make possible to any + // Note: all members here are kept public, to allow any // object accessing and manipulating them @@ -138,7 +138,7 @@ public: /// Allocates and populates the sorted list of atom ids void create_sorted_ids (void); - + /// \brief When updating atomic coordinates, translate them to align with the /// center of mass of the reference coordinates @@ -175,7 +175,7 @@ public: /// \brief If b_center or b_rotate is true, use this group to /// define the transformation (default: this group itself) atom_group *ref_pos_group; - + /// Total mass of the atom group cvm::real total_mass; @@ -199,9 +199,14 @@ public: /// \brief Initialize the group after a temporary vector of atoms atom_group (std::vector const &atoms); - /// \brief Add an atom to this group + /// \brief Add an atom to this group void add_atom (cvm::atom const &a); + /// \brief Re-initialize the total mass of a group. + /// This is needed in case the hosting MD code has an option to + /// change atom masses after their initialization. + void reset_mass (std::string &name, int i, int j); + /// \brief Default constructor atom_group(); @@ -216,7 +221,7 @@ public: void calc_apply_roto_translation(); /// \brief Save center of geometry fo ref positions, - /// then subtract it + /// then subtract it void center_ref_pos(); /// \brief Move all positions diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 561bb4a688..58ec049b3c 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); } @@ -59,7 +60,7 @@ void colvarbias::add_colvar (std::string const &cv_name) cvp->enable (colvar::task_gradients); if (cvm::debug()) cvm::log ("Applying this bias to collective variable \""+ - cvp->name+"\".\n"); + cvp->name+"\".\n"); colvars.push_back (cvp); colvar_forces.push_back (colvarvalue (cvp->type())); } else { @@ -79,13 +80,47 @@ void colvarbias::communicate_forces() } colvars[i]->add_bias_force (colvar_forces[i]); } -} +} + + +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) - : colvarbias (conf, key), + : colvarbias (conf, key), target_nsteps (0), target_nstages (0) { @@ -171,25 +206,22 @@ 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) +colvarbias_harmonic::~colvarbias_harmonic () { - cvm::fatal_error ("Error: change_configuration() not implemented.\n"); + if (cvm::n_harm_biases > 0) + cvm::n_harm_biases -= 1; } -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)) { @@ -317,7 +349,7 @@ cvm::real colvarbias_harmonic::update() if (target_equil_steps == 0 || cvm::step_absolute() % target_nsteps >= target_equil_steps) { // Start averaging after equilibration period, if requested - + // Square distance normalized by square colvar width cvm::real dist_sq = 0.0; for (size_t i = 0; i < colvars.size(); i++) { @@ -335,7 +367,7 @@ cvm::real colvarbias_harmonic::update() cvm::log ("Lambda= " + cvm::to_str (lambda) + " dA/dLambda= " + cvm::to_str (restraint_FE / cvm::real(target_nsteps - target_equil_steps))); - + // ...and move on to the next one if (stage < target_nstages) { @@ -363,7 +395,7 @@ cvm::real colvarbias_harmonic::update() if (cvm::debug()) cvm::log ("Done updating the harmonic bias \""+this->name+"\".\n"); - + // Force and energy calculation for (size_t i = 0; i < colvars.size(); i++) { colvar_forces[i] = @@ -379,6 +411,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"); @@ -408,7 +449,7 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is) return is; } -// int id = -1; +// int id = -1; std::string name = ""; // if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) && // (id != this->id) ) || @@ -442,6 +483,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 +530,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..5a2c2d8c5f 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -14,7 +14,7 @@ public: /// Name of this bias std::string name; - + /// Add a new collective variable to this bias void add_colvar (std::string const &cv_name); @@ -35,7 +35,7 @@ public: 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); @@ -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,11 +103,17 @@ 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); /// Destructor - virtual inline ~colvarbias_harmonic() {} + virtual ~colvarbias_harmonic(); protected: @@ -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 Number of steps required to reach the target force constant - /// or restraint centers - size_t target_nsteps; + /// \brief Intermediate quantity to compute the restraint free energy + /// (in TI, would be the accumulating FE derivative) + cvm::real restraint_FE; + + + /// \brief Equilibration steps for restraint FE calculation through TI + cvm::real target_equil_steps; /// \brief Number of stages over which to perform the change /// If zero, perform a continuous change @@ -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/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 677bedffbd..dcbfba6acd 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -15,7 +15,7 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key) { if (cvm::temperature() == 0.0) cvm::log ("WARNING: ABF should not be run without a thermostat or at 0 Kelvin!\n"); - + // ************* parsing general ABF options *********************** get_keyval (conf, "applyBias", apply_bias, true); @@ -32,7 +32,7 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key) } get_keyval (conf, "fullSamples", full_samples, 200); - if ( full_samples <= 1 ) full_samples = 1; + if ( full_samples <= 1 ) full_samples = 1; min_samples = full_samples / 2; // full_samples - min_samples >= 1 is guaranteed @@ -79,7 +79,21 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key) } // Here we could check for orthogonality of the Cartesian coordinates - // and make it just a warning if some parameter is set? + // and make it just a warning if some parameter is set? + } + + if (get_keyval (conf, "maxForce", max_force)) { + if (max_force.size() != colvars.size()) { + cvm::fatal_error ("Error: Number of parameters to maxForce does not match number of colvars."); + } + for (size_t i=0; i 0 ) { read_gradients_samples (); } - + cvm::log ("Finished ABF setup.\n"); } @@ -114,6 +128,9 @@ colvarbias_abf::~colvarbias_abf() } delete [] force; + + if (cvm::n_abf_biases > 0) + cvm::n_abf_biases -= 1; } @@ -125,7 +142,7 @@ cvm::real colvarbias_abf::update() if (cvm::debug()) cvm::log ("Updating ABF bias " + this->name); if (cvm::step_relative() == 0) { - + // At first timestep, do only: // initialization stuff (file operations relying on n_abf_biases // compute current value of colvars @@ -158,7 +175,7 @@ cvm::real colvarbias_abf::update() } // save bin for next timestep - force_bin = bin; + force_bin = bin; // Reset biasing forces from previous timestep for (size_t i=0; ivalue (bin)); if ( fact != 0.0 ) { - if ( (colvars.size() == 1) && colvars[0]->periodic_boundaries() ) { // Enforce a zero-mean bias on periodic, 1D coordinates - colvar_forces[0].real_value += fact * (grad[0] / cvm::real (count) - gradients->average ()); + // in other words: boundary condition is that the biasing potential is periodic + colvar_forces[0].real_value = fact * (grad[0] / cvm::real (count) - gradients->average ()); } else { for (size_t i=0; i max_force[i] * max_force[i] ) { + colvar_forces[i].real_value = (colvar_forces[i].real_value > 0 ? max_force[i] : -1.0 * max_force[i]); + } } } } @@ -203,7 +226,7 @@ cvm::real colvarbias_abf::update() // otherwise, backup and replace write_gradients_samples (output_prefix + ".hist", (cvm::step_absolute() > 0)); } - return 0.0; // TODO compute bias energy whenever possible (i.e. 1D with updateBias off) + return 0.0; } @@ -250,7 +273,7 @@ void colvarbias_abf::read_gradients_samples () samples_in_name = input_prefix[i] + ".count"; gradients_in_name = input_prefix[i] + ".grad"; // For user-provided files, the per-bias naming scheme may not apply - + std::ifstream is; cvm::log ("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name); @@ -272,7 +295,7 @@ void colvarbias_abf::read_gradients_samples () std::ostream & colvarbias_abf::write_restart (std::ostream& os) { - std::ios::fmtflags flags (os.flags ()); + std::ios::fmtflags flags (os.flags ()); os.setf(std::ios::fmtflags (0), std::ios::floatfield); // default floating-point format os << "abf {\n" @@ -325,7 +348,7 @@ std::istream & colvarbias_abf::read_restart (std::istream& is) if ( name == "" ) { cvm::fatal_error ("Error: \"abf\" block in the restart file has no name.\n"); } - + if ( !(is >> key) || !(key == "samples")) { cvm::log ("Error: in reading restart configuration for ABF bias \""+ this->name+"\" at position "+ @@ -381,7 +404,7 @@ colvarbias_histogram::colvarbias_histogram (std::string const &conf, char const bin.assign (colvars.size(), 0); out_name = cvm::output_prefix + "." + this->name + ".dat"; - cvm::log ("Histogram will be written to file " + out_name); + cvm::log ("Histogram will be written to file " + out_name); cvm::log ("Finished histogram setup.\n"); } @@ -395,6 +418,9 @@ colvarbias_histogram::~colvarbias_histogram() delete grid; grid = NULL; } + + if (cvm::n_histo_biases > 0) + cvm::n_histo_biases -= 1; } /// Update the grid @@ -452,7 +478,7 @@ std::istream & colvarbias_histogram::read_restart (std::istream& is) cvm::fatal_error ("Error: \"histogram\" block in the restart file " "has no name.\n"); } - + if ( !(is >> key) || !(key == "grid")) { cvm::log ("Error: in reading restart configuration for histogram \""+ this->name+"\" at position "+ diff --git a/lib/colvars/colvarbias_abf.h b/lib/colvars/colvarbias_abf.h index f107b78248..e82d2f08ac 100644 --- a/lib/colvars/colvarbias_abf.h +++ b/lib/colvars/colvarbias_abf.h @@ -46,6 +46,10 @@ private: bool b_history_files; size_t history_freq; + /// Cap applied biasing force? + bool cap_force; + std::vector max_force; + // Internal data and methods std::vector bin, force_bin; @@ -84,7 +88,7 @@ private: std::vector bin; std::string out_name; - int output_freq; + int output_freq; void write_grid (); std::ofstream grid_os; /// Stream for writing grid to disk diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 2d03b5c8d8..3ddc6af6c5 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -54,7 +54,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) get_keyval (conf, "multipleReplicas", b_replicas, false); if (b_replicas) comm = multiple_replicas; - else + else comm = single_replica; } @@ -76,7 +76,8 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) } get_keyval (conf, "keepHills", keep_hills, false); - get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true); + if (! get_keyval (conf, "writeFreeEnergyFile", dump_fes, true)) + get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent); get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false); for (size_t i = 0; i < colvars.size(); i++) { @@ -113,7 +114,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) "when using more than one replica.\n"); get_keyval (conf, "replicasRegistry", - replicas_registry_file, + replicas_registry_file, (this->name+".replicas.registry.txt")); get_keyval (conf, "replicaUpdateFrequency", @@ -171,7 +172,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) // 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()) @@ -188,7 +189,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) // 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(), + 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 \""+ @@ -232,7 +233,7 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key) cvm::log("Well-tempered metadynamics is used.\n"); cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\n"); } - + if (cvm::debug()) cvm::log ("Done initializing the metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n"); @@ -258,6 +259,9 @@ colvarbias_meta::~colvarbias_meta() if (hills_traj_os.good()) hills_traj_os.close(); + + if (cvm::n_meta_biases > 0) + cvm::n_meta_biases -= 1; } @@ -266,7 +270,7 @@ colvarbias_meta::~colvarbias_meta() // Hill management member functions // ********************************************************************** -std::list::const_iterator +std::list::const_iterator colvarbias_meta::create_hill (colvarbias_meta::hill const &h) { hill_iter const hills_end = hills.end(); @@ -324,7 +328,7 @@ colvarbias_meta::delete_hill (hill_iter &h) } if (hills_traj_os.good()) { - // output to the trajectory + // output to the trajectory hills_traj_os << "# DELETED this hill: " << (hills.back()).output_traj() << "\n"; @@ -459,9 +463,9 @@ cvm::real colvarbias_meta::update() cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin); cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); create_hill (hill ((hill_weight*exp_weight), colvars, hill_width)); - } else { + } else { create_hill (hill (hill_weight, colvars, hill_width)); - } + } break; case multiple_replicas: @@ -470,9 +474,9 @@ cvm::real colvarbias_meta::update() cvm::real const hills_energy_sum_here = hills_energy->value(curr_bin); cvm::real const exp_weight = std::exp(-hills_energy_sum_here/(bias_temperature*cvm::boltzmann())); create_hill (hill ((hill_weight*exp_weight), colvars, hill_width, replica_id)); - } else { + } else { create_hill (hill (hill_weight, colvars, hill_width, replica_id)); - } + } if (replica_hills_os.good()) { replica_hills_os << hills.back(); } else { @@ -591,11 +595,11 @@ cvm::real colvarbias_meta::update() calc_hills_force (ic, new_hills_begin, hills.end(), colvar_forces); } - if (cvm::debug()) + if (cvm::debug()) cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+ ", hills forces = "+cvm::to_str (colvar_forces)+".\n"); - if (cvm::debug()) + if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ ": adding the forces from the other replicas.\n"); @@ -611,7 +615,7 @@ cvm::real colvarbias_meta::update() replicas[ir]->hills.end(), colvar_forces); } - if (cvm::debug()) + if (cvm::debug()) cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+ ", hills forces = "+cvm::to_str (colvar_forces)+".\n"); } @@ -641,7 +645,7 @@ void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter h_first, } for (hill_iter h = h_first; h != h_last; h++) { - + // compute the gaussian exponent cvm::real cv_sqdev = 0.0; for (size_t i = 0; i < colvars.size(); i++) { @@ -684,8 +688,8 @@ void colvarbias_meta::calc_hills_force (size_t const &i, if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; cvm::real const half_width = 0.5 * h->widths[i]; - forces[i].real_value += - ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * + forces[i].real_value += + ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * (colvars[i]->dist2_lgrad (x, center)).real_value ); } break; @@ -993,7 +997,7 @@ void colvarbias_meta::read_replica_files() } // (re)read the state file if necessary - if ( (! (replicas[ir])->has_data) || + if ( (! (replicas[ir])->has_data) || (! (replicas[ir])->replica_state_file_in_sync) ) { cvm::log ("Metadynamics bias \""+this->name+"\""+ @@ -1013,11 +1017,11 @@ void colvarbias_meta::read_replica_files() } is.close(); } - + // now read the hills added after writing the state file if ((replicas[ir])->replica_hills_file.size()) { - if (cvm::debug()) + if (cvm::debug()) cvm::log ("Metadynamics bias \""+this->name+"\""+ ": checking for new hills from replica \""+ (replicas[ir])->replica_id+"\" in the file \""+ @@ -1064,7 +1068,7 @@ void colvarbias_meta::read_replica_files() "\" at position "+ cvm::to_str ((replicas[ir])->replica_hills_file_pos)+".\n"); - // test whether this is the end of the file + // test whether this is the end of the file is.seekg (0, std::ios::end); if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) { (replicas[ir])->update_status++; @@ -1108,7 +1112,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) size_t const start_pos = is.tellg(); if (comm == single_replica) { - // if using a multiple replicas scheme, output messages + // if using a multiple replicas scheme, output messages // are printed before and after calling this function cvm::log ("Restarting metadynamics bias \""+this->name+"\""+ ".\n"); @@ -1118,7 +1122,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) !(is >> brace) || !(brace == "{") || !(is >> colvarparse::read_block ("configuration", conf)) ) { - if (comm == single_replica) + if (comm == single_replica) cvm::log ("Error: in reading restart configuration for metadynamics bias \""+ this->name+"\""+ ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+ @@ -1239,7 +1243,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) is.clear(); is.seekg (hills_energy_gradients_pos, std::ios::beg); grids_from_restart_file = false; - if (!rebin_grids) { + 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+"\""+ @@ -1283,7 +1287,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) // read the hills explicitly written (if there are any) while (read_hill (is)) { - if (cvm::debug()) + if (cvm::debug()) cvm::log ("Read a previously saved hill under the " "metadynamics bias \""+ this->name+"\", created at step "+ @@ -1364,7 +1368,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) hills.erase (hills.begin(), old_hills_end); hills_off_grid.erase (hills_off_grid.begin(), old_hills_off_grid_end); } - + has_data = true; if (comm != single_replica) { @@ -1372,7 +1376,7 @@ std::istream & colvarbias_meta::read_restart (std::istream& is) } return is; -} +} std::istream & colvarbias_meta::read_hill (std::istream &is) @@ -1404,7 +1408,7 @@ std::istream & colvarbias_meta::read_hill (std::istream &is) std::vector h_centers (colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { - h_centers[i].type ((colvars[i]->value()).type()); + h_centers[i].type ((colvars[i]->value()).type()); } { // it is safer to read colvarvalue objects one at a time; @@ -1421,7 +1425,7 @@ std::istream & colvarbias_meta::read_hill (std::istream &is) get_keyval (data, "widths", h_widths, 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); @@ -1517,7 +1521,7 @@ std::ostream & colvarbias_meta::write_restart (std::ostream& os) } return os; -} +} void colvarbias_meta::write_pmf() @@ -1527,7 +1531,7 @@ void colvarbias_meta::write_pmf() pmf->create(); std::string fes_file_name_prefix (cvm::output_prefix); - + if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) { // if this is not the only free energy integrator, append // this bias's name, to distinguish it from the output of the other @@ -1598,7 +1602,7 @@ void colvarbias_meta::write_replica_state_file() cvm::fatal_error ("Error: in opening file \""+ replica_state_file+"\" for writing.\n"); - rep_state_os.setf (std::ios::scientific, std::ios::floatfield); + rep_state_os.setf (std::ios::scientific, std::ios::floatfield); rep_state_os << "\n" << "metadynamics {\n" << " configuration {\n" @@ -1672,7 +1676,7 @@ std::string colvarbias_meta::hill::output_traj() << std::setw (cvm::en_width) << W << "\n"; return os.str(); -} +} std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h) diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 1f51e5469c..8c3eb91003 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -33,7 +33,7 @@ public: /// Destructor virtual ~colvarbias_meta(); - + virtual cvm::real update(); virtual std::istream & read_restart (std::istream &is); @@ -84,7 +84,7 @@ protected: std::istream & read_hill (std::istream &is); /// \brief step present in a state file - /// + /// /// When using grids and reading state files containing them /// (multiple replicas), this is used to check whether a hill is /// newer or older than the grids @@ -147,10 +147,10 @@ protected: /// time steps, appending the step number to each file bool dump_fes_save; - /// \brief Whether to use well-tempered metadynamics - bool well_tempered; + /// \brief Whether to use well-tempered metadynamics + bool well_tempered; - /// \brief Biasing temperature in well-tempered metadynamics + /// \brief Biasing temperature in well-tempered metadynamics cvm::real bias_temperature; /// \brief Try to read the restart information by allocating new @@ -282,7 +282,7 @@ public: centers[i] = cv[i]->value(); widths[i] = cv[i]->width * hill_width; } - if (cvm::debug()) + if (cvm::debug()) cvm::log ("New hill, applied to "+cvm::to_str (cv.size())+ " collective variables, with centers "+ cvm::to_str (centers)+", widths "+ @@ -322,7 +322,7 @@ public: /// Destructor inline ~hill() {} - + /// Get the energy inline cvm::real energy() { diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index 25da537da8..316454e8ef 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -39,7 +39,7 @@ colvar::cvc::cvc (std::string const &conf) } -void colvar::cvc::parse_group (std::string const &conf, +void colvar::cvc::parse_group (std::string const &conf, char const *group_key, cvm::atom_group &group, bool optional) @@ -120,7 +120,7 @@ void colvar::cvc::debug_gradients (cvm::atom_group &group) for (size_t id = 0; id < 3; id++) { // (re)read original positions group.read_positions(); - // change one coordinate + // change one coordinate group[ia].pos[id] += cvm::debug_gradients_step_size; // (re)do the fit (if defined) if (group.b_center || group.b_rotate) { diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index a8a06814f4..964a733964 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -1,6 +1,16 @@ #ifndef COLVARCOMP_H #define COLVARCOMP_H +// Declaration of colvar::cvc base class and derived ones. +// +// Future cvc's could be declared on additional header files. +// After the declaration of a new derived class, its metric +// functions must be reimplemented as well. +// If the new cvc has no symmetry or periodicity, +// this can be done straightforwardly by using the macro: +// simple_scalar_dist_functions (derived_class) + + #include #include @@ -40,7 +50,7 @@ /// 2. add a call to the parser in colvar.C, within the function colvar::colvar() \par /// 3. declare the class in colvarcomp.h \par /// 4. implement the class in one of the files colvarcomp_*.C -/// +/// /// /// The cvm::atom and cvm::atom_group classes are available to /// transparently communicate with the simulation program. However, @@ -58,7 +68,7 @@ public: std::string name; /// \brief Description of the type of collective variable - /// + /// /// Normally this string is set by the parent \link colvar \endlink /// object within its constructor, when all \link cvc \endlink /// objects are initialized; therefore the main "config string" @@ -94,7 +104,7 @@ public: /// \brief Within the constructor, make a group parse its own /// options from the provided configuration string - void parse_group (std::string const &conf, + void parse_group (std::string const &conf, char const *group_key, cvm::atom_group &group, bool optional = false); @@ -164,7 +174,7 @@ public: /// \brief Square distance between x1 and x2 (can be redefined to /// transparently implement constraints, symmetries and /// periodicities) - /// + /// /// colvar::cvc::dist2() and the related functions are /// declared as "const" functions, but not "static", because /// additional parameters defining the metrics (e.g. the @@ -182,7 +192,7 @@ public: /// to provide a gradient which is compatible with the constraint, /// i.e. already deprived of its component normal to the constraint /// hypersurface. - /// + /// /// Finally, another useful application, if you are performing very /// many operations with these functions, could be to override the /// \link colvarvalue \endlink member functions and access directly @@ -591,7 +601,7 @@ protected: std::vector ref_pos; /// Geometric center of the reference coordinates - cvm::rvector ref_pos_center; + cvm::atom_pos ref_pos_center; /// Eigenvector (of a normal or essential mode): will always have zero center std::vector eigenvec; @@ -691,7 +701,7 @@ protected: /// \brief Compute system force on first site only to avoid unwanted /// coupling to other colvars (see e.g. Ciccotti et al., 2005) bool b_1site_force; - + public: /// Initialize by parsing the configuration @@ -759,7 +769,7 @@ public: static cvm::real switching_function (cvm::real const &r0, int const &exp_num, int const &exp_den, cvm::atom &A1, cvm::atom &A2); - + template /// \brief Calculate a coordination number through the function /// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec @@ -808,7 +818,7 @@ public: static cvm::real switching_function (cvm::real const &r0, int const &exp_num, int const &exp_den, cvm::atom &A1, cvm::atom &A2); - + virtual cvm::real dist2 (colvarvalue const &x1, colvarvalue const &x2) const; virtual colvarvalue dist2_lgrad (colvarvalue const &x1, @@ -822,7 +832,7 @@ public: /// \brief Colvar component: hydrogen bond, defined as the product of /// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol)) /// (colvarvalue::type_scalar type, range [0:1]) -class colvar::h_bond +class colvar::h_bond : public colvar::cvc { protected: @@ -887,7 +897,7 @@ public: // alpha_dihedrals(); // virtual inline ~alpha_dihedrals() {} // virtual void calc_value(); -// virtual void calc_gradients(); +// virtual void calc_gradients(); // virtual void apply_force (colvarvalue const &force); // virtual cvm::real dist2 (colvarvalue const &x1, // colvarvalue const &x2) const; @@ -921,7 +931,7 @@ protected: /// List of hydrogen bonds std::vector hb; - /// Contribution of the hb terms + /// Contribution of the hb terms cvm::real hb_coeff; public: @@ -1140,6 +1150,62 @@ public: }; + +// metrics functions for cvc implementations + +// simple definitions of the distance functions; these are useful only +// for optimization (the type check performed in the default +// colvarcomp functions is skipped) + +// definitions assuming the scalar type + +#define simple_scalar_dist_functions(TYPE) \ + \ + inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1, \ + colvarvalue const &x2) const \ + { \ + return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \ + } \ + \ + inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1, \ + colvarvalue const &x2) const \ + { \ + return 2.0 * (x1.real_value - x2.real_value); \ + } \ + \ + inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1, \ + colvarvalue const &x2) const \ + { \ + return this->dist2_lgrad (x2, x1); \ + } \ + \ + inline cvm::real colvar::TYPE::compare (colvarvalue const &x1, \ + colvarvalue const &x2) const \ + { \ + return this->dist2_lgrad (x1, x2); \ + } \ + \ + + simple_scalar_dist_functions (distance) + // NOTE: distance_z has explicit functions, see below + simple_scalar_dist_functions (distance_xy) + simple_scalar_dist_functions (distance_inv) + simple_scalar_dist_functions (angle) + simple_scalar_dist_functions (coordnum) + simple_scalar_dist_functions (selfcoordnum) + simple_scalar_dist_functions (h_bond) + simple_scalar_dist_functions (gyration) + simple_scalar_dist_functions (inertia) + simple_scalar_dist_functions (inertia_z) + simple_scalar_dist_functions (rmsd) + simple_scalar_dist_functions (orientation_angle) + simple_scalar_dist_functions (tilt) + simple_scalar_dist_functions (eigenvector) + // simple_scalar_dist_functions (alpha_dihedrals) + simple_scalar_dist_functions (alpha_angles) + simple_scalar_dist_functions (dihedPC) + + // metrics functions for cvc implementations with a periodicity inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1, @@ -1177,12 +1243,12 @@ inline void colvar::dihedral::wrap (colvarvalue &x) const if ((x.real_value - wrap_center) >= 180.0) { x.real_value -= 360.0; return; - } + } if ((x.real_value - wrap_center) < -180.0) { x.real_value += 360.0; return; - } + } return; } @@ -1222,68 +1288,16 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const if ((x.real_value - wrap_center) >= 180.0) { x.real_value -= 360.0; return; - } + } if ((x.real_value - wrap_center) < -180.0) { x.real_value += 360.0; return; - } + } return; } -// simple definitions of the distance functions; these are useful only -// for optimization (the type check performed in the default -// colvarcomp functions is skipped) - -// definitions assuming the scalar type - -#define simple_scalar_dist_functions(TYPE) \ - \ - inline cvm::real colvar::TYPE::dist2 (colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return (x1.real_value - x2.real_value)*(x1.real_value - x2.real_value); \ - } \ - \ - inline colvarvalue colvar::TYPE::dist2_lgrad (colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return 2.0 * (x1.real_value - x2.real_value); \ - } \ - \ - inline colvarvalue colvar::TYPE::dist2_rgrad (colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return this->dist2_lgrad (x2, x1); \ - } \ - \ - inline cvm::real colvar::TYPE::compare (colvarvalue const &x1, \ - colvarvalue const &x2) const \ - { \ - return this->dist2_lgrad (x1, x2); \ - } \ - \ - - simple_scalar_dist_functions (distance) - // NOTE: distance_z has explicit functions, see below - simple_scalar_dist_functions (distance_xy) - simple_scalar_dist_functions (distance_inv) - simple_scalar_dist_functions (angle) - simple_scalar_dist_functions (coordnum) - simple_scalar_dist_functions (selfcoordnum) - simple_scalar_dist_functions (h_bond) - simple_scalar_dist_functions (gyration) - simple_scalar_dist_functions (inertia) - simple_scalar_dist_functions (inertia_z) - simple_scalar_dist_functions (rmsd) - simple_scalar_dist_functions (orientation_angle) - simple_scalar_dist_functions (tilt) - simple_scalar_dist_functions (eigenvector) - // simple_scalar_dist_functions (alpha_dihedrals) - simple_scalar_dist_functions (alpha_angles) - simple_scalar_dist_functions (dihedPC) - // Projected distance // Differences should always be wrapped around 0 (ignoring wrap_center) @@ -1422,9 +1436,6 @@ inline cvm::real colvar::orientation::compare (colvarvalue const &x1, } - - - #endif diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 414c01df7e..bc2aee0598 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -77,7 +77,7 @@ void colvar::angle::calc_gradients() { cvm::real const cos_theta = (r21*r23)/(r21l*r23l); cvm::real const dxdcos = -1.0 / std::sqrt (1.0 - cos_theta*cos_theta); - + dxdr1 = (180.0/PI) * dxdcos * (1.0/r21l) * ( r23/r23l + (-1.0) * cos_theta * r21/r21l ); @@ -246,7 +246,7 @@ void colvar::dihedral::calc_gradients() cvm::real rA = A.norm(); cvm::rvector B = cvm::rvector::outer (r23, r34); cvm::real rB = B.norm(); - cvm::rvector C = cvm::rvector::outer (r23, A); + cvm::rvector C = cvm::rvector::outer (r23, A); cvm::real rC = C.norm(); cvm::real const cos_phi = (A*B)/(rA*rB); diff --git a/lib/colvars/colvarcomp_coordnums.cpp b/lib/colvars/colvarcomp_coordnums.cpp index a33f95e072..d08a5e505c 100644 --- a/lib/colvars/colvarcomp_coordnums.cpp +++ b/lib/colvars/colvarcomp_coordnums.cpp @@ -72,11 +72,14 @@ cvm::real colvar::coordnum::switching_function (cvm::rvector const &r0_vec, colvar::coordnum::coordnum (std::string const &conf) : distance (conf), b_anisotropic (false), b_group2_center_only (false) -{ +{ function_type = "coordnum"; x.type (colvarvalue::type_scalar); // group1 and group2 are already initialized by distance() + if (group1.b_dummy) + cvm::fatal_error ("Error: only group2 is allowed to be a dummy atom\n"); + // need to specify this explicitly because the distance() constructor // has set it to true @@ -105,7 +108,7 @@ colvar::coordnum::coordnum (std::string const &conf) cvm::fatal_error ("Error: odd exponents provided, can only use even ones.\n"); } - get_keyval (conf, "group2CenterOnly", b_group2_center_only, false); + get_keyval (conf, "group2CenterOnly", b_group2_center_only, group2.b_dummy); } @@ -121,18 +124,10 @@ void colvar::coordnum::calc_value() { x.real_value = 0.0; - // these are necessary: for each atom, gradients are summed together - // by multiple calls to switching_function() - group1.reset_atoms_data(); - group2.reset_atoms_data(); - - group1.read_positions(); - group2.read_positions(); - if (b_group2_center_only) { // create a fake atom to hold the group2 com coordinates - cvm::atom group2_com_atom (group2[0]); + cvm::atom group2_com_atom; group2_com_atom.pos = group2.center_of_mass(); if (b_anisotropic) { @@ -165,9 +160,10 @@ void colvar::coordnum::calc_gradients() if (b_group2_center_only) { // create a fake atom to hold the group2 com coordinates - cvm::atom group2_com_atom (group2[0]); + cvm::atom group2_com_atom; group2_com_atom.pos = group2.center_of_mass(); + if (b_anisotropic) { for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) switching_function (r0_vec, en, ed, *ai1, group2_com_atom); @@ -288,14 +284,6 @@ colvar::h_bond::~h_bond() void colvar::h_bond::calc_value() { - // this is necessary, because switching_function() will sum the new - // gradient to the current one - acceptor.reset_data(); - donor.reset_data(); - - acceptor.read_position(); - donor.read_position(); - x.real_value = colvar::coordnum::switching_function (r0, en, ed, acceptor, donor); } @@ -315,40 +303,11 @@ void colvar::h_bond::apply_force (colvarvalue const &force) } -// Self-coordination number for a group - -template -cvm::real colvar::selfcoordnum::switching_function (cvm::real const &r0, - int const &en, - int const &ed, - cvm::atom &A1, - cvm::atom &A2) -{ - cvm::rvector const diff = cvm::position_distance (A1.pos, A2.pos); - cvm::real const l2 = diff.norm2()/(r0*r0); - - // Assume en and ed are even integers, and avoid sqrt in the following - int const en2 = en/2; - int const ed2 = ed/2; - - cvm::real const xn = std::pow (l2, en2); - cvm::real const xd = std::pow (l2, ed2); - cvm::real const func = (1.0-xn)/(1.0-xd); - - 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*r0))*diff; - A1.grad += (-1.0)*dFdl2*dl2dx; - A2.grad += dFdl2*dl2dx; - } - - return func; -} colvar::selfcoordnum::selfcoordnum (std::string const &conf) : distance (conf, false) -{ +{ function_type = "selfcoordnum"; x.type (colvarvalue::type_scalar); @@ -379,13 +338,9 @@ void colvar::selfcoordnum::calc_value() { x.real_value = 0.0; - // for each atom, gradients are summed by multiple calls to switching_function() - group1.reset_atoms_data(); - group1.read_positions(); - for (size_t i = 0; i < group1.size() - 1; i++) for (size_t j = i + 1; j < group1.size(); j++) - x.real_value += switching_function (r0, en, ed, group1[i], group1[j]); + x.real_value += colvar::coordnum::switching_function (r0, en, ed, group1[i], group1[j]); } @@ -393,7 +348,7 @@ 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++) - switching_function (r0, en, ed, group1[i], group1[j]); + colvar::coordnum::switching_function (r0, en, ed, group1[i], group1[j]); } void colvar::selfcoordnum::apply_force (colvarvalue const &force) diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index a371cff605..339ba95f7d 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -113,7 +113,7 @@ void colvar::distance_vec::calc_value() } void colvar::distance_vec::calc_gradients() -{ +{ // gradients are not stored: a 3x3 matrix for each atom would be // needed to store just the identity matrix } @@ -140,7 +140,7 @@ colvar::distance_z::distance_z (std::string const &conf) // TODO detect PBC from MD engine (in simple cases) // and then update period in real time if (period != 0.0) - b_periodic = true; + b_periodic = true; if ((wrap_center != 0.0) && (period == 0.0)) { cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component," @@ -153,7 +153,7 @@ colvar::distance_z::distance_z (std::string const &conf) atom_groups.push_back (&ref1); // this group is optional parse_group (conf, "ref2", ref2, true); - + if (ref2.size()) { atom_groups.push_back (&ref2); cvm::log ("Using axis joining the centers of mass of groups \"ref\" and \"ref2\""); @@ -522,7 +522,7 @@ colvar::gyration::gyration (std::string const &conf) cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); } else { atoms.b_center = true; - atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0)); + atoms.ref_pos.assign (1, cvm::atom_pos (0.0, 0.0, 0.0)); } x.type (colvarvalue::type_scalar); @@ -700,7 +700,7 @@ colvar::rmsd::rmsd (std::string const &conf) parse_group (conf, "atoms", atoms); atom_groups.push_back (&atoms); - if (atoms.b_dummy) + if (atoms.b_dummy) cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom."); if (atoms.ref_pos_group != NULL) { @@ -732,7 +732,7 @@ colvar::rmsd::rmsd (std::string const &conf) std::string ref_pos_col; double ref_pos_col_value; - + 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); @@ -771,7 +771,7 @@ colvar::rmsd::rmsd (std::string const &conf) } } - + void colvar::rmsd::calc_value() { // rotational-translational fit is handled by the atom group @@ -813,9 +813,9 @@ void colvar::rmsd::calc_force_invgrads() { atoms.read_system_forces(); ft.real_value = 0.0; - + // Note: gradient square norm is 1/N_atoms - + for (size_t ia = 0; ia < atoms.size(); ia++) { ft.real_value += atoms[ia].grad * atoms[ia].system_force; } @@ -832,7 +832,7 @@ void colvar::rmsd::calc_Jacobian_derivative() // gradient of the rotation matrix cvm::matrix2d grad_rot_mat; - // gradients of products of 2 quaternion components + // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; for (size_t ia = 0; ia < atoms.size(); ia++) { @@ -850,15 +850,15 @@ void colvar::rmsd::calc_Jacobian_derivative() g23 = (atoms.rot.q)[2]*dq[3] + (atoms.rot.q)[3]*dq[2]; // Gradient of the rotation matrix wrt current Cartesian position - grad_rot_mat[0][0] = -2.0 * (g22 + g33); - grad_rot_mat[1][0] = 2.0 * (g12 + g03); - grad_rot_mat[2][0] = 2.0 * (g13 - g02); - grad_rot_mat[0][1] = 2.0 * (g12 - g03); - grad_rot_mat[1][1] = -2.0 * (g11 + g33); - grad_rot_mat[2][1] = 2.0 * (g01 + g23); - grad_rot_mat[0][2] = 2.0 * (g02 + g13); - grad_rot_mat[1][2] = 2.0 * (g23 - g01); - grad_rot_mat[2][2] = -2.0 * (g11 + g22); + grad_rot_mat[0][0] = -2.0 * (g22 + g33); + grad_rot_mat[1][0] = 2.0 * (g12 + g03); + grad_rot_mat[2][0] = 2.0 * (g13 - g02); + grad_rot_mat[0][1] = 2.0 * (g12 - g03); + grad_rot_mat[1][1] = -2.0 * (g11 + g33); + grad_rot_mat[2][1] = 2.0 * (g01 + g23); + grad_rot_mat[0][2] = 2.0 * (g02 + g13); + grad_rot_mat[1][2] = 2.0 * (g23 - g01); + grad_rot_mat[2][2] = -2.0 * (g11 + g22); cvm::atom_pos &y = ref_pos[ia]; @@ -991,7 +991,7 @@ colvar::eigenvector::eigenvector (std::string const &conf) "and eigenvector must be defined.\n"); } - cvm::rvector eig_center (0.0, 0.0, 0.0); + cvm::atom_pos eig_center (0.0, 0.0, 0.0); for (size_t i = 0; i < atoms.size(); i++) { eig_center += eigenvec[i]; } @@ -1053,7 +1053,7 @@ colvar::eigenvector::eigenvector (std::string const &conf) } } - + void colvar::eigenvector::calc_value() { x.real_value = 0.0; @@ -1087,7 +1087,7 @@ void colvar::eigenvector::calc_force_invgrads() { atoms.read_system_forces(); ft.real_value = 0.0; - + for (size_t ia = 0; ia < atoms.size(); ia++) { ft.real_value += eigenvec_invnorm2 * atoms[ia].grad * atoms[ia].system_force; @@ -1101,10 +1101,10 @@ void colvar::eigenvector::calc_Jacobian_derivative() cvm::matrix2d grad_rot_mat; cvm::quaternion &quat0 = atoms.rot.q; - // gradients of products of 2 quaternion components + // gradients of products of 2 quaternion components cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23; - cvm::atom_pos x_relative; + cvm::atom_pos x_relative; cvm::real sum = 0.0; for (size_t ia = 0; ia < atoms.size(); ia++) { @@ -1126,15 +1126,15 @@ void colvar::eigenvector::calc_Jacobian_derivative() // Gradient of the inverse rotation matrix wrt current Cartesian position // (transpose of the gradient of the direct rotation) - grad_rot_mat[0][0] = -2.0 * (g22 + g33); - grad_rot_mat[0][1] = 2.0 * (g12 + g03); - grad_rot_mat[0][2] = 2.0 * (g13 - g02); - grad_rot_mat[1][0] = 2.0 * (g12 - g03); - grad_rot_mat[1][1] = -2.0 * (g11 + g33); - grad_rot_mat[1][2] = 2.0 * (g01 + g23); - grad_rot_mat[2][0] = 2.0 * (g02 + g13); - grad_rot_mat[2][1] = 2.0 * (g23 - g01); - grad_rot_mat[2][2] = -2.0 * (g11 + g22); + grad_rot_mat[0][0] = -2.0 * (g22 + g33); + grad_rot_mat[0][1] = 2.0 * (g12 + g03); + grad_rot_mat[0][2] = 2.0 * (g13 - g02); + grad_rot_mat[1][0] = 2.0 * (g12 - g03); + grad_rot_mat[1][1] = -2.0 * (g11 + g33); + grad_rot_mat[1][2] = 2.0 * (g01 + g23); + grad_rot_mat[2][0] = 2.0 * (g02 + g13); + grad_rot_mat[2][1] = 2.0 * (g23 - g01); + grad_rot_mat[2][2] = -2.0 * (g11 + g22); for (size_t i = 0; i < 3; i++) { for (size_t j = 0; j < 3; j++) { @@ -1143,7 +1143,7 @@ void colvar::eigenvector::calc_Jacobian_derivative() } } - jd.real_value = sum * std::sqrt (eigenvec_invnorm2); + jd.real_value = sum * std::sqrt (eigenvec_invnorm2); } diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index b30d50e1c1..1750782b7a 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -122,7 +122,7 @@ void colvar::alpha_angles::calc_value() if (theta.size()) { - cvm::real const theta_norm = + cvm::real const theta_norm = (1.0-hb_coeff) / cvm::real (theta.size()); for (size_t i = 0; i < theta.size(); i++) { @@ -162,7 +162,7 @@ void colvar::alpha_angles::calc_value() void colvar::alpha_angles::calc_gradients() { - for (size_t i = 0; i < theta.size(); i++) + for (size_t i = 0; i < theta.size(); i++) (theta[i])->calc_gradients(); for (size_t i = 0; i < hb.size(); i++) @@ -175,9 +175,9 @@ void colvar::alpha_angles::apply_force (colvarvalue const &force) if (theta.size()) { - cvm::real const theta_norm = + cvm::real const theta_norm = (1.0-hb_coeff) / cvm::real (theta.size()); - + for (size_t i = 0; i < theta.size(); i++) { cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol; @@ -185,10 +185,10 @@ void colvar::alpha_angles::apply_force (colvarvalue const &force) (1.0 - std::pow (t, (int) 4)) ); cvm::real const dfdt = - 1.0/(1.0 - std::pow (t, (int) 4)) * + 1.0/(1.0 - std::pow (t, (int) 4)) * ( (-2.0 * t) + (-1.0*f)*(-4.0 * std::pow (t, (int) 3)) ); - (theta[i])->apply_force (theta_norm * + (theta[i])->apply_force (theta_norm * dfdt * (1.0/theta_tol) * force.real_value ); } @@ -216,7 +216,7 @@ void colvar::alpha_angles::apply_force (colvarvalue const &force) // of this cvc (alpha) // This is true of all cvcs with sub-cvcs, and those // that do not calculate explicit gradients - // SO: we need a flag giving the availability of + // SO: we need a flag giving the availability of // atomic gradients colvar::dihedPC::dihedPC (std::string const &conf) : cvc (conf) @@ -260,7 +260,7 @@ colvar::dihedPC::dihedPC (std::string const &conf) std::string vecFileName; int vecNumber; if (get_keyval (conf, "vectorFile", vecFileName, vecFileName)) { - get_keyval (conf, "vectorNumber", vecNumber, 0); + get_keyval (conf, "vectorNumber", vecNumber, 0); if (vecNumber < 1) cvm::fatal_error ("A positive value of vectorNumber is required."); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 50c9e0da2e..b77c563878 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -74,7 +74,7 @@ colvar::orientation::orientation (std::string const &conf) } - + colvar::orientation::orientation() : cvc () { @@ -160,7 +160,7 @@ void colvar::orientation_angle::calc_value() void colvar::orientation_angle::calc_gradients() { cvm::real const dxdq0 = - ( ((rot.q).q0 * (rot.q).q0 < 1.0) ? + ( ((rot.q).q0 * (rot.q).q0 < 1.0) ? ((180.0 / PI) * (-2.0) / std::sqrt (1.0 - ((rot.q).q0 * (rot.q).q0))) : 0.0 ); diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index c6099bcb8f..a300b91340 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -141,7 +141,7 @@ void colvar_grid_gradient::write_1D_integral (std::ostream &os) { cvm::real bin, min, integral; std::vector int_vals; - + os << "# xi A(xi)\n"; if ( cv.size() != 1 ) { @@ -162,7 +162,7 @@ void colvar_grid_gradient::write_1D_integral (std::ostream &os) } for (std::vector ix = new_index(); index_ok (ix); incr (ix), bin += 1.0 ) { - + if (samples) { size_t const samples_here = samples->value (ix); if (samples_here) diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index 80a7b70b16..ccaf046039 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -266,11 +266,11 @@ public: } - { + { 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)+" - "+ @@ -284,7 +284,7 @@ public: 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); } @@ -336,7 +336,7 @@ public: { return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin); } - + /// \brief Same as the standard version, but uses different parameters inline colvarvalue bin_to_value_scalar (int const &i_bin, colvarvalue const &new_offset, @@ -349,7 +349,7 @@ public: inline void set_value (std::vector const &ix, T const &t, size_t const &imult = 0) - { + { data[this->address (ix)+imult] = t; has_data = true; } @@ -367,7 +367,7 @@ public: /// \brief Add a constant to all elements (fast loop) inline void add_constant (T const &t) { - for (size_t i = 0; i < nt; i++) + for (size_t i = 0; i < nt; i++) data[i] += t; has_data = true; } @@ -375,11 +375,11 @@ public: /// \brief Multiply all elements by a scalar constant (fast loop) inline void multiply_constant (cvm::real const &a) { - for (size_t i = 0; i < nt; i++) + for (size_t i = 0; i < nt; i++) data[i] *= a; } - + /// \brief Get the bin indices corresponding to the provided values of /// the colvars inline std::vector const get_colvars_index (std::vector const &values) const @@ -432,7 +432,7 @@ public: /// \brief Add data from another grid of the same type - /// + /// /// Note: this function maps other_grid inside this one regardless /// of whether it fits or not. void map_grid (colvar_grid const &other_grid) @@ -482,11 +482,11 @@ public: if (other_grid.multiplicity() != this->multiplicity()) cvm::fatal_error ("Error: trying to sum togetehr two grids with values of " "different multiplicity.\n"); - if (scale_factor != 1.0) + if (scale_factor != 1.0) for (size_t i = 0; i < data.size(); i++) { data[i] += scale_factor * other_grid.data[i]; } - else + else // skip multiplication if possible for (size_t i = 0; i < data.size(); i++) { data[i] += other_grid.data[i]; @@ -571,30 +571,30 @@ public: std::ostream & write_params (std::ostream &os) { os << "grid_parameters {\n n_colvars " << nd << "\n"; - + os << " lower_boundaries "; - for (size_t i = 0; i < nd; i++) + for (size_t i = 0; i < nd; i++) os << " " << lower_boundaries[i]; os << "\n"; os << " upper_boundaries "; - for (size_t i = 0; i < nd; i++) + for (size_t i = 0; i < nd; i++) os << " " << upper_boundaries[i]; os << "\n"; os << " widths "; - for (size_t i = 0; i < nd; i++) + for (size_t i = 0; i < nd; i++) os << " " << widths[i]; os << "\n"; os << " sizes "; - for (size_t i = 0; i < nd; i++) + for (size_t i = 0; i < nd; i++) os << " " << nx[i]; os << "\n"; os << "}\n"; return os; - } + } bool parse_params (std::string const &conf) @@ -648,9 +648,9 @@ public: { for (size_t i = 0; i < nd; i++) { if ( (std::sqrt (cv[i]->dist2 (cv[i]->lower_boundary, - lower_boundaries[i])) > 1.0E-10) || + lower_boundaries[i])) > 1.0E-10) || (std::sqrt (cv[i]->dist2 (cv[i]->upper_boundary, - upper_boundaries[i])) > 1.0E-10) || + upper_boundaries[i])) > 1.0E-10) || (std::sqrt (cv[i]->dist2 (cv[i]->width, widths[i])) > 1.0E-10) ) { cvm::fatal_error ("Error: restart information for a grid is " @@ -669,11 +669,11 @@ public: // matter: boundaries should be EXACTLY the same (otherwise, // map_grid() should be used) if ( (std::fabs (other_grid.lower_boundaries[i] - - lower_boundaries[i]) > 1.0E-10) || + lower_boundaries[i]) > 1.0E-10) || (std::fabs (other_grid.upper_boundaries[i] - - upper_boundaries[i]) > 1.0E-10) || + upper_boundaries[i]) > 1.0E-10) || (std::fabs (other_grid.widths[i] - - widths[i]) > 1.0E-10) || + widths[i]) > 1.0E-10) || (data.size() != other_grid.data.size()) ) { cvm::fatal_error ("Error: inconsistency between " "two grids that are supposed to be equal, " @@ -681,7 +681,7 @@ public: } } } - + /// \brief Write the grid data without labels, as they are /// represented in memory @@ -714,7 +714,7 @@ public: std::istream & read_raw (std::istream &is) { size_t const start_pos = is.tellg(); - + for (std::vector ix = new_index(); index_ok (ix); incr (ix)) { for (size_t imult = 0; imult < mult; imult++) { T new_value; @@ -845,7 +845,7 @@ std::istream & read_multicol (std::istream &is, bool add = false) if ( !(is >> x) ) end_of_file = true; bin[i] = value_to_bin_scalar (x, i); } - if (end_of_file) break; + if (end_of_file) break; for (size_t imult = 0; imult < mult; imult++) { is >> new_value[imult]; @@ -904,7 +904,7 @@ public: ++(data[this->address (ix)]); } - /// \brief Get the binned count indexed by ix from the newly read data + /// \brief Get the binned count indexed by ix from the newly read data inline size_t const & new_count (std::vector const &ix, size_t const &imult = 0) { @@ -993,7 +993,7 @@ public: A0 += data[address (ix)]; grad[n] = 0.5 * (A1 - A0) / widths[n]; } - return grad; + return grad; } /// \brief Return the value of the function at ix divided by its diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 98e05ca00d..03fe18da94 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -9,14 +9,14 @@ colvarmodule::colvarmodule (char const *config_filename, colvarproxy *proxy_in) -{ +{ // pointer to the proxy object if (proxy == NULL) { proxy = proxy_in; parse = new colvarparse(); } else { - cvm::fatal_error ("Error: trying to allocate twice the collective " - "variable module.\n"); + cvm::fatal_error ("Error: trying to allocate the collective " + "variable module twice.\n"); } cvm::log (cvm::line_marker); @@ -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, @@ -80,12 +85,10 @@ colvarmodule::colvarmodule (char const *config_filename, std::string (output_prefix+".colvars.state") : std::string ("colvars.state"))+"\".\n"); - cv_traj_name = + cv_traj_name = (output_prefix.size() ? std::string (output_prefix+".colvars.traj") : std::string ("colvars.traj")); - cvm::log ("The trajectory file will be \""+ - cv_traj_name+"\".\n"); if (cv_traj_freq) { // open trajectory file @@ -94,6 +97,8 @@ colvarmodule::colvarmodule (char const *config_filename, "\".\n"); cv_traj_os.open (cv_traj_name.c_str(), std::ios::app); } else { + cvm::log ("Writing to colvar trajectory file \""+cv_traj_name+ + "\".\n"); proxy->backup_file (cv_traj_name.c_str()); cv_traj_os.open (cv_traj_name.c_str(), std::ios::out); } @@ -135,7 +140,7 @@ colvarmodule::colvarmodule (char const *config_filename, } -std::istream & colvarmodule::read_restart (std::istream &is) +std::istream & colvarmodule::read_restart (std::istream &is) { { // read global restart information @@ -192,7 +197,7 @@ void colvarmodule::init_colvars (std::string const &conf) (colvars.back())->check_keywords (colvar_conf, "colvar"); cvm::decrease_depth(); } else { - cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n"); + cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n"); } colvar_conf = ""; } @@ -347,13 +352,13 @@ void colvarmodule::calc() { } // calculate collective variables and their gradients - if (cvm::debug()) + if (cvm::debug()) cvm::log ("Calculating collective variables.\n"); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->calc(); + (*cvi)->calc(); } cvm::decrease_depth(); @@ -365,7 +370,7 @@ void colvarmodule::calc() { for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { - total_bias_energy += (*bi)->update(); + total_bias_energy += (*bi)->update(); } cvm::decrease_depth(); @@ -376,12 +381,12 @@ void colvarmodule::calc() { for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->reset_bias_force(); + (*cvi)->reset_bias_force(); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { - (*bi)->communicate_forces(); + (*bi)->communicate_forces(); } cvm::decrease_depth(); @@ -393,12 +398,12 @@ void colvarmodule::calc() { for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->analyse(); + (*cvi)->analyse(); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { - (*bi)->analyse(); + (*bi)->analyse(); } cvm::decrease_depth(); } @@ -426,7 +431,7 @@ void colvarmodule::calc() { cvi != colvars.end(); cvi++) { if ((*cvi)->tasks[colvar::task_gradients]) - (*cvi)->communicate_forces(); + (*cvi)->communicate_forces(); } cvm::decrease_depth(); @@ -442,7 +447,7 @@ void colvarmodule::calc() { if (!write_restart (restart_out_os)) cvm::fatal_error ("Error: in writing restart file.\n"); restart_out_os.close(); - } + } } // write trajectory file, if needed @@ -471,15 +476,20 @@ void colvarmodule::calc() { if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) { cv_traj_os << "# " << cvm::wrap_string ("step", cvm::it_width-2) << " "; - if (cvm::debug()) + if (cvm::debug()) cv_traj_os.flush(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj_label (cv_traj_os); } + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + (*bi)->write_traj_label (cv_traj_os); + } cv_traj_os << "\n"; - if (cvm::debug()) + if (cvm::debug()) cv_traj_os.flush(); } cvm::decrease_depth(); @@ -494,13 +504,18 @@ 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(); } cvm::decrease_depth(); - if (restart_out_freq) { + if (restart_out_freq) { // flush the trajectory file if we are at the restart frequency if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { @@ -527,7 +542,7 @@ void colvarmodule::analyze() cvi != colvars.end(); cvi++) { cvm::increase_depth(); - (*cvi)->analyse(); + (*cvi)->analyse(); cvm::decrease_depth(); } @@ -536,12 +551,20 @@ void colvarmodule::analyze() bi != biases.end(); bi++) { cvm::increase_depth(); - (*bi)->analyse(); + (*bi)->analyse(); cvm::decrease_depth(); } } +void colvarmodule::setup() +{ + // loop over all components of all colvars to reset masses of all groups + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); cvi++) { + (*cvi)->setup(); + } +} colvarmodule::~colvarmodule() { @@ -565,7 +588,7 @@ colvarmodule::~colvarmodule() delete parse; proxy = NULL; -} +} void colvarmodule::write_output_files() @@ -590,7 +613,7 @@ void colvarmodule::write_output_files() (*cvi)->write_output_files(); } cvm::decrease_depth(); - + // do not close to avoid problems with multiple NAMD runs cv_traj_os.flush(); } @@ -631,7 +654,7 @@ bool colvarmodule::read_traj (char const *traj_filename, continue; - } else { + } else { if ((it % 1000) == 0) std::cerr << "Reading from trajectory, step = " << it @@ -721,6 +744,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 +813,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 @@ -827,7 +901,7 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q) std::string ("euler")) ) { // parse the Euler angles - + char sep; cvm::real phi, theta, psi; if ( !(is >> sep) || !(sep == '(') || @@ -846,7 +920,7 @@ std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q) // parse the quaternion components - is.seekg (start_pos, std::ios::beg); + is.seekg (start_pos, std::ios::beg); char sep; if ( !(is >> sep) || !(sep == '(') || !(is >> q.q0) || !(is >> sep) || !(sep == ',') || @@ -932,15 +1006,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) @@ -959,7 +1031,7 @@ void colvarmodule::rotation::build_matrix (std::vector const &pos C.zx() += pos1[i].z * pos2[i].x; C.zy() += pos1[i].z * pos2[i].y; C.zz() += pos1[i].z * pos2[i].z; - } + } // build the "overlap" matrix, whose eigenvectors are stationary // points of the RMSD in the space of rotations @@ -1068,11 +1140,12 @@ void colvarmodule::rotation::calc_optimal_rotation if (q_old.norm2() > 0.0) { q.match (q_old); if (q_old.inner (q) < (1.0 - crossing_threshold)) { - cvm::log ("Warning: discontinuous rotation!\n"); + cvm::log ("Warning: one molecular orientation has changed by more than "+ + cvm::to_str (crossing_threshold)+": discontinuous rotation ?\n"); } } q_old = q; - + if (cvm::debug()) { if (b_debug_gradients) { cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+ @@ -1143,8 +1216,8 @@ void colvarmodule::rotation::calc_optimal_rotation for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { dq0_1[p] += - (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] + - (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] + + (Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] + + (Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] + (Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p]; } } @@ -1193,15 +1266,15 @@ void colvarmodule::rotation::calc_optimal_rotation for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { dq0_2[p] += - (Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] + - (Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] + + (Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] + + (Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] + (Q3[i] * ds_2[i][j] * Q0[j]) / (L0-L3) * Q3[p]; } } } if (cvm::debug()) { - + if (b_debug_gradients) { matrix2d S_new; @@ -1216,7 +1289,7 @@ void colvarmodule::rotation::calc_optimal_rotation // diagonalize the new overlap matrix for (size_t i = 0; i < 4; i++) { for (size_t j = 0; j < 4; j++) { - S_new[i][j] += + S_new[i][j] += colvarmodule::debug_gradients_step_size * ds_2[i][j][comp]; } } diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 7449ae23f7..2587f95d88 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -2,16 +2,16 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2013-04-17" +#define COLVARS_VERSION "2013-06-19" #endif #ifndef COLVARS_DEBUG #define COLVARS_DEBUG false #endif -/// \file colvarmodule.h +/// \file colvarmodule.h /// \brief Collective variables main module -/// +/// /// This file declares the main class for defining and manipulating /// collective variables: there should be only one instance of this /// class, because several variables are made static (i.e. they are @@ -40,7 +40,7 @@ class colvarproxy; /// Class to control the collective variables calculation. An object /// (usually one) of this class is spawned from the MD program, /// containing all i/o routines and general interface. -/// +/// /// At initialization, the colvarmodule object creates a proxy object /// to provide a transparent interface between the MD program and the /// child objects @@ -54,7 +54,7 @@ private: public: friend class colvarproxy; - + /// Defining an abstract real number allows to switch precision typedef double real; /// Residue identifier @@ -147,7 +147,7 @@ public: /// \brief Number of histograms initialized (no limit on the /// number) static size_t n_histo_biases; - + /// \brief Whether debug output should be enabled (compile-time option) static inline bool debug() { @@ -157,8 +157,9 @@ public: /// \brief Constructor \param config_name Configuration file name /// \param restart_name (optional) Restart file name - colvarmodule (char const *config_name, + colvarmodule (char const *config_name, colvarproxy *proxy_in); + /// Destructor ~colvarmodule(); @@ -168,6 +169,11 @@ public: /// Initialize collective variable biases void init_biases (std::string const &conf); + /// Re-initialize data at the beginning of a run. For use with + /// MD codes that can change system parameters like atom masses + /// between run commands. + void setup(); + /// Load new configuration for the given bias - /// currently works for harmonic (force constant and/or centers) void change_configuration (std::string const &bias_name, std::string const &conf); @@ -194,16 +200,16 @@ public: bool read_traj (char const *traj_filename, size_t traj_read_begin, size_t traj_read_end); - + /// Get the pointer of a colvar from its name (returns NULL if not found) static colvar * colvar_p (std::string const &name); /// Quick conversion of an object to a string - template static std::string to_str (T const &x, + template static std::string to_str (T const &x, size_t const &width = 0, size_t const &prec = 0); /// Quick conversion of a vector of objects to a string - template static std::string to_str (std::vector const &x, + template static std::string to_str (std::vector const &x, size_t const &width = 0, size_t const &prec = 0); @@ -247,10 +253,10 @@ public: /// \brief Time step of MD integrator (fs) static real dt(); - + /// Request calculation of system force from MD engine static void request_system_force(); - + /// Print a message to the main log static void log (std::string const &message); @@ -278,7 +284,7 @@ public: /// \brief Get the closest periodic image to a reference position /// \param pos The position to look for the closest periodic image - /// \param ref_pos (optional) The reference position + /// \param ref_pos (optional) The reference position static void select_closest_image (atom_pos &pos, atom_pos const &ref_pos); @@ -290,6 +296,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 @@ -371,7 +387,7 @@ std::ostream & operator << (std::ostream &os, cvm::rvector const &v); std::istream & operator >> (std::istream &is, cvm::rvector &v); -template std::string cvm::to_str (T const &x, +template std::string cvm::to_str (T const &x, size_t const &width, size_t const &prec) { std::ostringstream os; @@ -384,7 +400,7 @@ template std::string cvm::to_str (T const &x, return os.str(); } -template std::string cvm::to_str (std::vector const &x, +template std::string cvm::to_str (std::vector const &x, size_t const &width, size_t const &prec) { if (!x.size()) return std::string (""); @@ -434,7 +450,7 @@ inline void cvm::request_system_force() { proxy->request_system_force (true); } - + inline void cvm::select_closest_image (atom_pos &pos, atom_pos const &ref_pos) { diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index e08746c547..2e62072f34 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -241,11 +241,11 @@ _get_keyval_vector_ (colvarvalue); -bool colvarparse::get_keyval (std::string const &conf, - char const *key, - bool &value, - bool const &def_value, - Parse_Mode const parse_mode) +bool colvarparse::get_keyval (std::string const &conf, + char const *key, + bool &value, + bool const &def_value, + Parse_Mode const parse_mode) { std::string data; bool b_found = false, b_found_any = false; @@ -277,7 +277,7 @@ bool colvarparse::get_keyval (std::string const &conf, (data == std::string ("false")) ) { value = false; } else - cvm::fatal_error ("Error: boolean values only are allowed " + cvm::fatal_error ("Error: boolean values only are allowed " "for \""+std::string (key)+"\".\n"); if (parse_mode != parse_silent) { cvm::log ("# "+std::string (key)+" = "+ @@ -371,8 +371,8 @@ void colvarparse::check_keywords (std::string &conf, char const *key) std::string uk; std::istringstream line_is (line); line_is >> uk; - if (cvm::debug()) - cvm::log ("Checking the validity of \""+uk+"\" from line:\n" + line); + // if (cvm::debug()) + // cvm::log ("Checking the validity of \""+uk+"\" from line:\n" + line); uk = to_lower_cppstr (uk); bool found_keyword = false; @@ -426,8 +426,8 @@ bool colvarparse::key_lookup (std::string const &conf, colvarparse::dummy_pos = 0; // start from the first occurrence of key - size_t pos = conf_lower.find (key, save_pos); - + size_t pos = conf_lower.find (key, save_pos); + // iterate over all instances until it finds the isolated keyword while (true) { @@ -442,7 +442,7 @@ bool colvarparse::key_lookup (std::string const &conf, if ( std::string ("\n"+white_space+ "}").find (conf[pos-1]) == std::string::npos ) { - // none of the valid delimiting characters is on the left of key + // none of the valid delimiting characters is on the left of key b_isolated_left = false; } } @@ -451,7 +451,7 @@ bool colvarparse::key_lookup (std::string const &conf, if ( std::string ("\n"+white_space+ "{").find (conf[pos+key.size()]) == std::string::npos ) { - // none of the valid delimiting characters is on the right of key + // none of the valid delimiting characters is on the right of key b_isolated_right = false; } } @@ -461,7 +461,7 @@ bool colvarparse::key_lookup (std::string const &conf, bool const b_isolated = (b_isolated_left && b_isolated_right && b_not_within_block); - + if (b_isolated) { // found it break; @@ -490,7 +490,7 @@ bool colvarparse::key_lookup (std::string const &conf, size_t line_begin = (pl == std::string::npos) ? 0 : 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)); + std::string line (conf, line_begin, (line_end-line_begin)); size_t data_begin = (to_lower_cppstr (line)).find (key) + key.size(); data_begin = line.find_first_not_of (white_space, data_begin+1); @@ -527,9 +527,9 @@ bool colvarparse::key_lookup (std::string const &conf, line_begin = line_end; nl = conf.find ('\n', line_begin+1); - if (nl == std::string::npos) + if (nl == std::string::npos) line_end = conf.size(); - else + else line_end = nl; line.append (conf, line_begin, (line_end-line_begin)); @@ -565,7 +565,7 @@ bool colvarparse::key_lookup (std::string const &conf, // << "\n"; } } - + save_pos = line_end; return true; diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index e7e1d96275..d3b3968a79 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -53,15 +53,15 @@ public: {} /// How a keyword is parsed in a string - enum Parse_Mode { + enum Parse_Mode { /// \brief (default) Read the first instance of a keyword (if /// any), report its value, and print a warning when there is more /// than one - parse_normal, + parse_normal, /// \brief Like parse_normal, but don't send any message to the log /// (useful e.g. in restart files when such messages are very /// numerous and redundant) - parse_silent + parse_silent }; /// \fn get_keyval bool const get_keyval (std::string const &conf, @@ -112,7 +112,7 @@ public: 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 ("")); @@ -169,7 +169,7 @@ public: /// data (optional) holds the string provided after "key", if any /// \param save_pos (optional) stores the position of the keyword /// within "conf", useful when doing multiple calls \param - /// save_delimiters (optional) + /// save_delimiters (optional) bool key_lookup (std::string const &conf, char const *key, std::string &data = dummy_string, diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 0af4826c96..c1416f13b7 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -16,6 +16,12 @@ public: /// Pointer to the instance of colvarmodule colvarmodule *colvars; + /// Default destructor + virtual inline ~colvarproxy() {} + + + // **************** SYSTEM-WIDE PHYSICAL QUANTITIES **************** + /// \brief Value of the unit for atomic coordinates with respect to /// angstroms (used by some variables for hard-coded default values) virtual cvm::real unit_angstrom() = 0; @@ -29,21 +35,11 @@ public: /// \brief Time step of the simulation (fs) virtual cvm::real dt() = 0; - /// Pass restraint energy value for current timestep to MD engine - virtual void add_energy (cvm::real energy) = 0; + /// \brief Pseudo-random number with Gaussian distribution + virtual cvm::real rand_gaussian (void) = 0; - /// Tell the proxy whether system forces are needed - virtual void request_system_force (bool yesno) = 0; - - /// Print a message to the main log - virtual void log (std::string const &message) = 0; - - /// Print a message to the main log and exit with error code - virtual void fatal_error (std::string const &message) = 0; - - /// Print a message to the main log and exit normally - virtual void exit (std::string const &message) = 0; + // **************** SIMULATION PARAMETERS **************** /// \brief Prefix to be used for input files (restarts, not /// configuration) @@ -59,25 +55,32 @@ public: /// \brief Restarts will be fritten each time this number of steps has passed virtual size_t restart_frequency() = 0; + + + // **************** ACCESS ATOMIC DATA **************** + + /// Pass restraint energy value for current timestep to MD engine + virtual void add_energy (cvm::real energy) = 0; + + /// Tell the proxy whether system forces are needed (may not always be available) + virtual void request_system_force (bool yesno) = 0; + + + // **************** PERIODIC BOUNDARY CONDITIONS **************** - /// \brief Get the simple distance vector between two positions - /// (with periodic boundary conditions handled transparently) + /// \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; - /// \brief Get the square distance between two positions (with - /// periodic boundary conditions handled transparently) - /// - /// Note: in the case of periodic boundary conditions, this provides - /// an analytical square distance (while taking the square of - /// position_distance() would produce leads to a cusp) + /// \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) = 0; + 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 + /// \param ref_pos The reference position virtual void select_closest_image (cvm::atom_pos &pos, cvm::atom_pos const &ref_pos) = 0; @@ -89,6 +92,18 @@ public: cvm::atom_pos const &ref_pos); + + // **************** INPUT/OUTPUT **************** + + /// Print a message to the main log + virtual void log (std::string const &message) = 0; + + /// Print a message to the main log and exit with error code + virtual void fatal_error (std::string const &message) = 0; + + /// Print a message to the main log and exit normally + virtual void exit (std::string const &message) = 0; + /// \brief Read atom identifiers from a file \param filename name of /// the file (usually a PDB) \param atoms array to which atoms read /// from "filename" will be appended \param pdb_field (optiona) if @@ -97,7 +112,7 @@ public: virtual void load_atoms (char const *filename, std::vector &atoms, std::string const pdb_field, - double const pdb_field_value = 0.0) = 0; + double const pdb_field_value = 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 @@ -108,14 +123,9 @@ public: std::string const pdb_field, double const pdb_field_value = 0.0) = 0; - /// \brief Rename the given file, under the convention provided by - /// the MD program - virtual void backup_file (char const *filename) = 0; + /// \brief Rename the given file, before overwriting it + virtual void backup_file (char const *filename) {} - /// \brief Pseudo-random number with Gaussian distribution - virtual cvm::real rand_gaussian (void) = 0; - - virtual inline ~colvarproxy() {} }; @@ -128,6 +138,12 @@ 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) +{ + return (position_distance (pos1, pos2)).norm2(); +} + #endif diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index fa18310ac4..4c285631a8 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -19,7 +19,7 @@ class colvarmodule::rvector { public: cvm::real x, y, z; - + inline rvector() : x (0.0), y (0.0), z (0.0) {} @@ -55,7 +55,7 @@ public: } - inline cvm::rvector & operator = (cvm::real const &v) + inline cvm::rvector & operator = (cvm::real const &v) { x = v; y = v; @@ -63,28 +63,28 @@ public: return *this; } - inline void operator += (cvm::rvector const &v) + inline void operator += (cvm::rvector const &v) { x += v.x; y += v.y; z += v.z; } - inline void operator -= (cvm::rvector const &v) + inline void operator -= (cvm::rvector const &v) { x -= v.x; y -= v.y; z -= v.z; } - inline void operator *= (cvm::real const &v) + inline void operator *= (cvm::real const &v) { x *= v; y *= v; z *= v; } - inline void operator /= (cvm::real const& v) + inline void operator /= (cvm::real const& v) { x /= v; y /= v; @@ -113,57 +113,57 @@ public: } - static inline cvm::rvector outer (cvm::rvector const &v1, cvm::rvector const &v2) + static inline cvm::rvector outer (cvm::rvector const &v1, cvm::rvector const &v2) { return cvm::rvector ( v1.y*v2.z - v2.y*v1.z, -v1.x*v2.z + v2.x*v1.z, v1.x*v2.y - v2.x*v1.y); } - friend inline cvm::rvector operator - (cvm::rvector const &v) + friend inline cvm::rvector operator - (cvm::rvector const &v) { return cvm::rvector (-v.x, -v.y, -v.z); } - friend inline int operator == (cvm::rvector const &v1, cvm::rvector const &v2) + friend inline int operator == (cvm::rvector const &v1, cvm::rvector const &v2) { return (v1.x == v2.x) && (v1.y == v2.y) && (v1.z == v2.z); } - friend inline int operator != (cvm::rvector const &v1, cvm::rvector const &v2) + friend inline int operator != (cvm::rvector const &v1, cvm::rvector const &v2) { return (v1.x != v2.x) || (v1.y != v2.y) || (v1.z != v2.z); } - friend inline cvm::rvector operator + (cvm::rvector const &v1, cvm::rvector const &v2) + friend inline cvm::rvector operator + (cvm::rvector const &v1, cvm::rvector const &v2) { return cvm::rvector (v1.x + v2.x, v1.y + v2.y, v1.z + v2.z); } - friend inline cvm::rvector operator - (cvm::rvector const &v1, cvm::rvector const &v2) + friend inline cvm::rvector operator - (cvm::rvector const &v1, cvm::rvector const &v2) { return cvm::rvector (v1.x - v2.x, v1.y - v2.y, v1.z - v2.z); } - friend inline cvm::real operator * (cvm::rvector const &v1, cvm::rvector const &v2) + friend inline cvm::real operator * (cvm::rvector const &v1, cvm::rvector const &v2) { return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z; } - friend inline cvm::rvector operator * (cvm::real const &a, cvm::rvector const &v) + friend inline cvm::rvector operator * (cvm::real const &a, cvm::rvector const &v) { return cvm::rvector (a*v.x, a*v.y, a*v.z); } - friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real const &a) + friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real const &a) { return cvm::rvector (a*v.x, a*v.y, a*v.z); } - friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real const &a) + friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real const &a) { return cvm::rvector (v.x/a, v.y/a, v.z/a); } - + }; @@ -186,7 +186,7 @@ public: { return length; } - + /// Default constructor inline vector1d (T const &t = T()) { @@ -251,7 +251,7 @@ public: return prod; } - /// Formatted output + /// Formatted output friend std::ostream & operator << (std::ostream &os, vector1d const &v) { @@ -310,7 +310,7 @@ public: this->alloc(); reset(); } - + /// Constructor from a 2-d C array inline matrix2d (T const **m) { @@ -380,7 +380,7 @@ public: // } // } - /// Formatted output + /// Formatted output friend std::ostream & operator << (std::ostream &os, matrix2d const &m) { @@ -454,16 +454,16 @@ public: inline cvm::real zz() const { return array[2][2]; } /// Constructor from a 2-d C array - inline rmatrix (cvm::real const **m) - : cvm::matrix2d (m) + inline rmatrix (cvm::real const **m) + : cvm::matrix2d (m) {} /// Default constructor - inline rmatrix() + inline rmatrix() : cvm::matrix2d() {} - /// Constructor component by component + /// Constructor component by component inline rmatrix (cvm::real const &xxi, cvm::real const &xyi, cvm::real const &xzi, @@ -472,7 +472,7 @@ public: cvm::real const &yzi, cvm::real const &zxi, cvm::real const &zyi, - cvm::real const &zzi) + cvm::real const &zzi) : cvm::matrix2d() { this->xx() = xxi; @@ -488,7 +488,7 @@ public: /// Destructor inline ~rmatrix() - {} + {} /// Return the determinant inline cvm::real determinant() const @@ -528,7 +528,7 @@ public: }; - + inline cvm::rvector operator * (cvm::rmatrix const &m, cvm::rvector const &r) { @@ -680,7 +680,7 @@ public: return std::sqrt (this->norm2()); } - /// Return the conjugate quaternion + /// Return the conjugate quaternion inline cvm::quaternion conjugate() const { return cvm::quaternion (q0, -q1, -q2, -q3); @@ -720,7 +720,7 @@ public: return cvm::quaternion (0.0, v.x, v.y, v.z); } /// Return the vector component - inline cvm::rvector get_vector() const + inline cvm::rvector get_vector() const { return cvm::rvector (q1, q2, q3); } @@ -1002,11 +1002,11 @@ public: if (q.q0 != 0.0) { // cvm::real const x = iprod/q.q0; - + cvm::real const dspindx = (180.0/PI) * 2.0 * (1.0 / (1.0 + (iprod*iprod)/(q.q0*q.q0))); - return - cvm::quaternion ( dspindx * (iprod * (-1.0) / (q.q0*q.q0)), + return + cvm::quaternion ( dspindx * (iprod * (-1.0) / (q.q0*q.q0)), dspindx * ((1.0/q.q0) * axis.x), dspindx * ((1.0/q.q0) * axis.y), dspindx * ((1.0/q.q0) * axis.z)); @@ -1023,11 +1023,11 @@ public: inline cvm::real cos_theta (cvm::rvector const &axis) const { cvm::rvector const q_vec = q.get_vector(); - cvm::real const alpha = + cvm::real const alpha = (180.0/PI) * 2.0 * std::atan2 (axis * q_vec, q.q0); cvm::real const cos_spin_2 = std::cos (alpha * (PI/180.0) * 0.5); - cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ? + cvm::real const cos_theta_2 = ( (cos_spin_2 != 0.0) ? (q.q0 / cos_spin_2) : (0.0) ); // cos(2t) = 2*cos(t)^2 - 1 @@ -1044,7 +1044,7 @@ public: if (q.q0 != 0.0) { - cvm::real const d_cos_theta_dq0 = + cvm::real const d_cos_theta_dq0 = (4.0 * q.q0 / (cos_spin_2*cos_spin_2)) * (1.0 - (iprod*iprod)/(q.q0*q.q0) / (1.0 + (iprod*iprod)/(q.q0*q.q0))); @@ -1052,7 +1052,7 @@ public: (4.0 * q.q0 / (cos_spin_2*cos_spin_2) * (iprod/q.q0) / (1.0 + (iprod*iprod)/(q.q0*q.q0))); - return cvm::quaternion (d_cos_theta_dq0, + return cvm::quaternion (d_cos_theta_dq0, d_cos_theta_dqn * axis.x, d_cos_theta_dqn * axis.y, d_cos_theta_dqn * axis.z); diff --git a/lib/colvars/colvarvalue.cpp b/lib/colvars/colvarvalue.cpp index f7d6aa6bf6..18a3bd00c5 100644 --- a/lib/colvars/colvarvalue.cpp +++ b/lib/colvars/colvarvalue.cpp @@ -28,7 +28,7 @@ void colvarvalue::error_rside cvm::fatal_error ("Trying to assign a colvar value with type \""+ type_desc[this->value_type]+"\" to one with type \""+ type_desc[vt]+"\".\n"); -} +} void colvarvalue::error_lside (colvarvalue::Type const &vt) const @@ -36,7 +36,7 @@ void colvarvalue::error_lside cvm::fatal_error ("Trying to use a colvar value with type \""+ type_desc[vt]+"\" as one of type \""+ type_desc[this->value_type]+"\".\n"); -} +} @@ -125,7 +125,7 @@ void colvarvalue::p2leg_opt (colvarvalue const &x, break; case colvarvalue::type_vector: while (xvi != xv_end) { - cvm::real const cosine = + cvm::real const cosine = ((xvi)->rvector_value * x.rvector_value) / ((xvi)->rvector_value.norm() * x.rvector_value.norm()); xvi++; @@ -167,7 +167,7 @@ void colvarvalue::p2leg_opt (colvarvalue const &x, break; case colvarvalue::type_vector: while (xvi != xv_end) { - cvm::real const cosine = + cvm::real const cosine = ((xvi)->rvector_value * x.rvector_value) / ((xvi)->rvector_value.norm() * x.rvector_value.norm()); xvi++; @@ -225,11 +225,11 @@ std::istream & operator >> (std::istream &is, colvarvalue &x) switch (x.type()) { case colvarvalue::type_scalar: - is >> x.real_value; + is >> x.real_value; break; case colvarvalue::type_vector: case colvarvalue::type_unitvector: - is >> x.rvector_value; + is >> x.rvector_value; break; case colvarvalue::type_quaternion: is >> x.quaternion_value; @@ -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..e262e7febe 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -37,7 +37,7 @@ /// the problem, because \link colvarvalue \endlink objects are first /// initialized in the configuration, and the stream operation will be /// performed only when reading restart files. -/// +/// /// No problem of course with the output streams: \code os << x; /// \endcode will print a different output according to the value of /// colvarvalue::value_type, and the width of such output is returned @@ -92,7 +92,7 @@ public: cvm::quaternion quaternion_value; /// Current type of this colvarvalue object - Type value_type; + Type value_type; static inline bool type_checking() { @@ -146,7 +146,7 @@ public: case type_scalar: real_value = x.real_value; break; - case type_vector: + case type_vector: case type_unitvector: rvector_value = x.rvector_value; break; @@ -254,7 +254,7 @@ public: /// Ensure that the two types are the same within a binary operator void static check_types (colvarvalue const &x1, colvarvalue const &x2); - /// Undefined operation + /// Undefined operation void undef_op() const; /// Trying to assign this \link colvarvalue \endlink object to @@ -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 @@ -324,7 +324,7 @@ inline void colvarvalue::reset() case colvarvalue::type_notset: default: break; - } + } } @@ -344,7 +344,7 @@ inline void colvarvalue::apply_constraints() case colvarvalue::type_notset: default: break; - } + } } @@ -362,7 +362,7 @@ inline cvm::real colvarvalue::norm2() const case colvarvalue::type_notset: default: return 0.0; - } + } } @@ -383,7 +383,7 @@ inline colvarvalue colvarvalue::inverse() const case colvarvalue::type_notset: default: undef_op(); - } + } return colvarvalue(); } @@ -391,11 +391,11 @@ inline colvarvalue colvarvalue::inverse() const inline colvarvalue & colvarvalue::operator = (colvarvalue const &x) { if (this->value_type != type_notset) - if (this->value_type != x.value_type) + if (this->value_type != x.value_type) error_lside (x.value_type); this->value_type = x.value_type; - + switch (this->value_type) { case colvarvalue::type_scalar: this->real_value = x.real_value;