diff --git a/lib/colvars/Makefile.g++ b/lib/colvars/Makefile.g++ index 272d93b801..80422be3a3 100644 --- a/lib/colvars/Makefile.g++ +++ b/lib/colvars/Makefile.g++ @@ -17,7 +17,7 @@ SHELL = /bin/sh 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 + colvargrid.cpp colvarmodule.cpp colvarparse.cpp colvarvalue.cpp colvarbias_alb.cpp LIB = libcolvars.a OBJ = $(SRC:.cpp=.o) @@ -60,6 +60,9 @@ colvaratoms.o: colvaratoms.cpp colvarmodule.h colvartypes.h colvarproxy.h \ 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_alb.o: colvarbias_alb.cpp colvarmodule.h colvartypes.h \ + colvar.h colvarvalue.h colvarparse.h colvarbias_alb.h \ + colvarbias.h 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 \ @@ -89,7 +92,7 @@ colvargrid.o: colvargrid.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvargrid.h 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 + colvargrid.h colvarbias_abf.h colvarbias_alb.h colvarparse.o: colvarparse.cpp colvarmodule.h colvartypes.h colvarproxy.h \ colvarvalue.h colvarparse.h colvarvalue.o: colvarvalue.cpp colvarmodule.h colvartypes.h colvarproxy.h \ diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 18d9f5b964..328bc6a662 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -14,8 +14,9 @@ colvarbias::colvarbias (std::string const &conf, char const *key) if (to_lower_cppstr (key_str) == std::string ("abf")) { rank = cvm::n_abf_biases+1; } - if (to_lower_cppstr (key_str) == std::string ("harmonic")) { - rank = cvm::n_harm_biases+1; + if (to_lower_cppstr (key_str) == std::string ("harmonic") || + to_lower_cppstr (key_str) == std::string ("linear")) { + rank = cvm::n_rest_biases+1; } if (to_lower_cppstr (key_str) == std::string ("histogram")) { rank = cvm::n_histo_biases+1; @@ -118,20 +119,13 @@ std::ostream & colvarbias::write_traj (std::ostream &os) -colvarbias_harmonic::colvarbias_harmonic (std::string const &conf, +colvarbias_restraint::colvarbias_restraint (std::string const &conf, char const *key) : colvarbias (conf, key), target_nstages (0), target_nsteps (0) { get_keyval (conf, "forceConstant", force_k, 1.0); - for (size_t i = 0; i < colvars.size(); i++) { - if (colvars[i]->width != 1.0) - cvm::log ("The force constant for colvar \""+colvars[i]->name+ - "\" will be rescaled to "+ - cvm::to_str (force_k/(colvars[i]->width*colvars[i]->width))+ - " according to the specified width.\n"); - } // get the initial restraint centers colvar_centers.resize (colvars.size()); @@ -151,7 +145,7 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf, } if (colvar_centers.size() != colvars.size()) - cvm::fatal_error ("Error: number of harmonic centers does not match " + cvm::fatal_error ("Error: number of centers does not match " "that of collective variables.\n"); if (get_keyval (conf, "targetCenters", target_centers, colvar_centers)) { @@ -211,17 +205,17 @@ colvarbias_harmonic::colvarbias_harmonic (std::string const &conf, acc_work = 0.0; if (cvm::debug()) - cvm::log ("Done initializing a new harmonic restraint bias.\n"); + cvm::log ("Done initializing a new restraint bias.\n"); } -colvarbias_harmonic::~colvarbias_harmonic () +colvarbias_restraint::~colvarbias_restraint () { - if (cvm::n_harm_biases > 0) - cvm::n_harm_biases -= 1; + if (cvm::n_rest_biases > 0) + cvm::n_rest_biases -= 1; } -void colvarbias_harmonic::change_configuration (std::string const &conf) +void colvarbias_restraint::change_configuration (std::string const &conf) { get_keyval (conf, "forceConstant", force_k, force_k); if (get_keyval (conf, "centers", colvar_centers, colvar_centers)) { @@ -233,7 +227,7 @@ void colvarbias_harmonic::change_configuration (std::string const &conf) } -cvm::real colvarbias_harmonic::energy_difference (std::string const &conf) +cvm::real colvarbias_restraint::energy_difference (std::string const &conf) { std::vector alt_colvar_centers; cvm::real alt_force_k; @@ -252,20 +246,21 @@ cvm::real colvarbias_harmonic::energy_difference (std::string const &conf) } for (size_t i = 0; i < colvars.size(); i++) { - alt_bias_energy += 0.5 * alt_force_k / (colvars[i]->width * colvars[i]->width) * - colvars[i]->dist2 (colvars[i]->value(), alt_colvar_centers[i]); + alt_bias_energy += restraint_potential(restraint_convert_k(alt_force_k, colvars[i]->width), + colvars[i], + alt_colvar_centers[i]); } return alt_bias_energy - bias_energy; } -cvm::real colvarbias_harmonic::update() +cvm::real colvarbias_restraint::update() { bias_energy = 0.0; if (cvm::debug()) - cvm::log ("Updating the harmonic bias \""+this->name+"\".\n"); + cvm::log ("Updating the restraint bias \""+this->name+"\".\n"); // Setup first stage of staged variable force constant calculation if (b_chg_force_k && target_nstages && cvm::step_absolute() == 0) { @@ -277,7 +272,7 @@ cvm::real colvarbias_harmonic::update() } force_k = starting_force_k + (target_force_k - starting_force_k) * std::pow (lambda, force_k_exp); - cvm::log ("Harmonic restraint " + this->name + ", stage " + + cvm::log ("Restraint " + this->name + ", stage " + cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); cvm::log ("Setting force constant to " + cvm::to_str (force_k)); } @@ -297,7 +292,7 @@ cvm::real colvarbias_harmonic::update() (target_nsteps - cvm::step_absolute())); } if (cvm::debug()) - cvm::log ("Center increment for the harmonic bias \""+ + cvm::log ("Center increment for the restraint bias \""+ this->name+"\": "+cvm::to_str (centers_incr)+" at stage "+cvm::to_str (stage)+ ".\n"); } @@ -330,7 +325,7 @@ cvm::real colvarbias_harmonic::update() } if (cvm::debug()) - cvm::log ("Current centers for the harmonic bias \""+ + cvm::log ("Current centers for the restraint bias \""+ this->name+"\": "+cvm::to_str (colvar_centers)+".\n"); } @@ -380,7 +375,7 @@ cvm::real colvarbias_harmonic::update() } force_k = starting_force_k + (target_force_k - starting_force_k) * std::pow (lambda, force_k_exp); - cvm::log ("Harmonic restraint " + this->name + ", stage " + + cvm::log ("Restraint " + this->name + ", stage " + cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); cvm::log ("Setting force constant to " + cvm::to_str (force_k)); } @@ -394,21 +389,16 @@ cvm::real colvarbias_harmonic::update() } if (cvm::debug()) - cvm::log ("Done updating the harmonic bias \""+this->name+"\".\n"); + cvm::log ("Done updating the restraint bias \""+this->name+"\".\n"); // Force and energy calculation for (size_t i = 0; i < colvars.size(); i++) { - colvar_forces[i] = - (-0.5) * force_k / - (colvars[i]->width * colvars[i]->width) * - colvars[i]->dist2_lgrad (colvars[i]->value(), - colvar_centers[i]); - bias_energy += 0.5 * force_k / (colvars[i]->width * colvars[i]->width) * - colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]); - if (cvm::debug()) - cvm::log ("dist_grad["+cvm::to_str (i)+ - "] = "+cvm::to_str (colvars[i]->dist2_lgrad (colvars[i]->value(), - colvar_centers[i]))+"\n"); + colvar_forces[i] = -restraint_force(restraint_convert_k(force_k, colvars[i]->width), + colvars[i], + colvar_centers[i]); + bias_energy += restraint_potential(restraint_convert_k(force_k, colvars[i]->width), + colvars[i], + colvar_centers[i]); } if (b_output_acc_work) { @@ -421,26 +411,26 @@ cvm::real colvarbias_harmonic::update() } if (cvm::debug()) - cvm::log ("Current forces for the harmonic bias \""+ + cvm::log ("Current forces for the restraint bias \""+ this->name+"\": "+cvm::to_str (colvar_forces)+".\n"); return bias_energy; } -std::istream & colvarbias_harmonic::read_restart (std::istream &is) +std::istream & colvarbias_restraint::read_restart (std::istream &is) { size_t const start_pos = is.tellg(); - cvm::log ("Restarting harmonic bias \""+ + cvm::log ("Restarting restraint bias \""+ this->name+"\".\n"); std::string key, brace, conf; - if ( !(is >> key) || !(key == "harmonic") || + if ( !(is >> key) || !(key == "restraint" || key == "harmonic") || !(is >> brace) || !(brace == "{") || !(is >> colvarparse::read_block ("configuration", conf)) ) { - cvm::log ("Error: in reading restart configuration for harmonic bias \""+ + cvm::log ("Error: in reading restart configuration for restraint bias \""+ this->name+"\" at position "+ cvm::to_str (is.tellg())+" in stream.\n"); is.clear(); @@ -456,10 +446,10 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is) if ( (colvarparse::get_keyval (conf, "name", name, std::string (""), colvarparse::parse_silent)) && (name != this->name) ) cvm::fatal_error ("Error: in the restart file, the " - "\"harmonic\" block has a wrong name\n"); + "\"restraint\" block has a wrong name\n"); // if ( (id == -1) && (name == "") ) { if (name.size() == 0) { - cvm::fatal_error ("Error: \"harmonic\" block in the restart file " + cvm::fatal_error ("Error: \"restraint\" block in the restart file " "has no identifiers.\n"); } @@ -490,7 +480,7 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is) is >> brace; if (brace != "}") { - cvm::fatal_error ("Error: corrupt restart information for harmonic bias \""+ + cvm::fatal_error ("Error: corrupt restart information for restraint bias \""+ this->name+"\": no matching brace at position "+ cvm::to_str (is.tellg())+" in the restart file.\n"); is.setstate (std::ios::failbit); @@ -499,9 +489,9 @@ std::istream & colvarbias_harmonic::read_restart (std::istream &is) } -std::ostream & colvarbias_harmonic::write_restart (std::ostream &os) +std::ostream & colvarbias_restraint::write_restart (std::ostream &os) { - os << "harmonic {\n" + os << "restraint {\n" << " configuration {\n" // << " id " << this->id << "\n" << " name " << this->name << "\n"; @@ -541,7 +531,7 @@ std::ostream & colvarbias_harmonic::write_restart (std::ostream &os) } -std::ostream & colvarbias_harmonic::write_traj_label (std::ostream &os) +std::ostream & colvarbias_restraint::write_traj_label (std::ostream &os) { os << " "; @@ -564,7 +554,7 @@ std::ostream & colvarbias_harmonic::write_traj_label (std::ostream &os) } -std::ostream & colvarbias_harmonic::write_traj (std::ostream &os) +std::ostream & colvarbias_restraint::write_traj (std::ostream &os) { os << " "; @@ -588,3 +578,55 @@ std::ostream & colvarbias_harmonic::write_traj (std::ostream &os) return os; } +colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(std::string const &conf, char const *key) : + colvarbias_restraint(conf, key) { + for (size_t i = 0; i < colvars.size(); i++) { + if (colvars[i]->width != 1.0) + cvm::log ("The force constant for colvar \""+colvars[i]->name+ + "\" will be rescaled to "+ + cvm::to_str (restraint_convert_k(force_k, colvars[i]->width))+ + " according to the specified width.\n"); + } +} + +cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k, colvar* x, const colvarvalue &xcenter) const +{ + return 0.5 * k * x->dist2(x->value(), xcenter); +} + +colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k, colvar* x, const colvarvalue &xcenter) const +{ + return 0.5 * k * x->dist2_lgrad(x->value(), xcenter); +} + +cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k, cvm::real dist_measure) const +{ + return k / (dist_measure * dist_measure); +} + + +colvarbias_restraint_linear::colvarbias_restraint_linear(std::string const &conf, char const *key) : + colvarbias_restraint(conf, key) { + for (size_t i = 0; i < colvars.size(); i++) { + if (colvars[i]->width != 1.0) + cvm::log ("The force constant for colvar \""+colvars[i]->name+ + "\" will be rescaled to "+ + cvm::to_str (restraint_convert_k(force_k, colvars[i]->width))+ + " according to the specified width.\n"); + } +} + +cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k, colvar* x, const colvarvalue &xcenter) const +{ + return k * (x->value() - xcenter); +} + +colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k, colvar* x, const colvarvalue &xcenter) const +{ + return k; +} + +cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k, cvm::real dist_measure) const +{ + return k / dist_measure; +} diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index 5a2c2d8c5f..c4e6493151 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -82,9 +82,9 @@ protected: }; -/// \brief Harmonic restraint, optionally moving towards a target +/// \brief Bias restraint, optionally moving towards a target /// (implementation of \link colvarbias \endlink) -class colvarbias_harmonic : public colvarbias { +class colvarbias_restraint : public colvarbias { public: @@ -110,14 +110,23 @@ public: virtual std::ostream & write_traj (std::ostream &os); /// \brief Constructor - colvarbias_harmonic (std::string const &conf, char const *key); + colvarbias_restraint (std::string const &conf, char const *key); /// Destructor - virtual ~colvarbias_harmonic(); + virtual ~colvarbias_restraint(); protected: + /// \brief Potential function + virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0; + + /// \brief Force function + virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const = 0; + + ///\brief Unit scaling + virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const = 0; + /// \brief Restraint centers std::vector colvar_centers; @@ -143,7 +152,6 @@ protected: /// \brief Accumulated work cvm::real acc_work; - /// \brief Restraint force constant cvm::real force_k; @@ -182,6 +190,42 @@ protected: size_t target_nsteps; }; +/// \brief Harmonic bias restraint +/// (implementation of \link colvarbias_restraint \endlink) +class colvarbias_restraint_harmonic : public colvarbias_restraint { + +public: + colvarbias_restraint_harmonic(std::string const &conf, char const *key); + +protected: /// \brief Potential function + virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const; + + /// \brief Force function + virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const; + + ///\brief Unit scaling + virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const; + +}; + +/// \brief Linear bias restraint +/// (implementation of \link colvarbias_restraint \endlink) +class colvarbias_restraint_linear : public colvarbias_restraint { + +public: + colvarbias_restraint_linear(std::string const &conf, char const *key); + +protected: /// \brief Potential function + virtual cvm::real restraint_potential(cvm::real k, colvar* x, const colvarvalue& xcenter) const; + + /// \brief Force function + virtual colvarvalue restraint_force(cvm::real k, colvar* x, const colvarvalue& xcenter) const; + + ///\brief Unit scaling + virtual cvm::real restraint_convert_k(cvm::real k, cvm::real dist_measure) const; + +}; + #endif diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 964a733964..6a4315d9d1 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -601,7 +601,7 @@ protected: std::vector ref_pos; /// Geometric center of the reference coordinates - cvm::atom_pos ref_pos_center; + cvm::atom_pos ref_pos_center; /// Eigenvector (of a normal or essential mode): will always have zero center std::vector eigenvec; diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index ec92617259..3420618605 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -4,6 +4,7 @@ #include "colvarproxy.h" #include "colvar.h" #include "colvarbias.h" +#include "colvarbias_alb.h" #include "colvarbias_meta.h" #include "colvarbias_abf.h" @@ -247,10 +248,10 @@ void colvarmodule::init_biases (std::string const &conf) if (harm_conf.size()) { cvm::log (cvm::line_marker); cvm::increase_depth(); - biases.push_back (new colvarbias_harmonic (harm_conf, "harmonic")); + biases.push_back (new colvarbias_restraint_harmonic (harm_conf, "harmonic")); (biases.back())->check_keywords (harm_conf, "harmonic"); cvm::decrease_depth(); - n_harm_biases++; + n_rest_biases++; } else { cvm::log ("Warning: \"harmonic\" keyword found without configuration.\n"); } @@ -258,6 +259,45 @@ void colvarmodule::init_biases (std::string const &conf) } } + { + /// initialize linear restraints + std::string lin_conf = ""; + size_t lin_pos = 0; + while (parse->key_lookup (conf, "linear", lin_conf, lin_pos)) { + if (lin_conf.size()) { + cvm::log (cvm::line_marker); + cvm::increase_depth(); + biases.push_back (new colvarbias_restraint_linear (lin_conf, "linear")); + (biases.back())->check_keywords (lin_conf, "linear"); + cvm::decrease_depth(); + n_rest_biases++; + } else { + cvm::log ("Warning: \"linear\" keyword found without configuration.\n"); + } + lin_conf = ""; + } + } + + { + /// initialize adaptive linear biases + std::string alb_conf = ""; + size_t alb_pos = 0; + while (parse->key_lookup (conf, "ALB", alb_conf, alb_pos)) { + if (alb_conf.size()) { + cvm::log (cvm::line_marker); + cvm::increase_depth(); + biases.push_back (new colvarbias_alb (alb_conf, "ALB")); + (biases.back())->check_keywords (alb_conf, "ALB"); + cvm::decrease_depth(); + n_rest_biases++; + } else { + cvm::log ("Warning: \"ALB\" keyword found without configuration.\n"); + } + alb_conf = ""; + } + } + + { /// initialize histograms std::string histo_conf = ""; @@ -815,7 +855,6 @@ void cvm::read_index_file (char const *filename) } - void cvm::load_atoms (char const *file_name, std::vector &atoms, std::string const &pdb_field, @@ -893,7 +932,7 @@ void cvm::load_coords_xyz (char const *filename, std::vector colvarmodule::colvars; std::vector colvarmodule::biases; size_t colvarmodule::n_abf_biases = 0; -size_t colvarmodule::n_harm_biases = 0; +size_t colvarmodule::n_rest_biases = 0; size_t colvarmodule::n_histo_biases = 0; size_t colvarmodule::n_meta_biases = 0; colvarproxy *colvarmodule::proxy = NULL; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index fc9872e9bf..ca9a149c94 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -2,7 +2,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2014-05-08" +#define COLVARS_VERSION "2014-08-13" #endif #ifndef COLVARS_DEBUG @@ -142,9 +142,9 @@ public: /// \brief Number of metadynamics biases initialized (in normal /// conditions should be 1) static size_t n_meta_biases; - /// \brief Number of harmonic biases initialized (no limit on the + /// \brief Number of restraint biases initialized (no limit on the /// number) - static size_t n_harm_biases; + static size_t n_rest_biases; /// \brief Number of histograms initialized (no limit on the /// number) static size_t n_histo_biases;