From 87215c925c6fd213e76f6ef429bc6c4718d422e7 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Fri, 30 Nov 2012 18:03:00 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9130 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- lib/colvars/colvaratoms.cpp | 140 ++++++++++++++------------- lib/colvars/colvaratoms.h | 12 ++- lib/colvars/colvarcomp_distances.cpp | 25 +++-- lib/colvars/colvarmodule.h | 2 +- 4 files changed, 101 insertions(+), 78 deletions(-) diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 0e67cf3ba3..7a3ffa1903 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -12,7 +12,7 @@ cvm::atom_group::atom_group (std::string const &conf, char const *key, atom_group *ref_pos_group_in) - : b_center (false), b_rotate (false), b_prevent_fitting (false), + : b_center (false), b_rotate (false), b_user_defined_fit (false), ref_pos_group (NULL), // this is always set within parse(), // regardless of ref_pos_group_in noforce (false) @@ -223,43 +223,45 @@ void cvm::atom_group::parse (std::string const &conf, // FITTING OPTIONS - bool fit_defined_by_user = - ( get_keyval (group_conf, "centerReference", b_center, false, mode) || - get_keyval (group_conf, "rotateReference", b_rotate, false, mode) ); - if ((!b_rotate) && (!b_center) && fit_defined_by_user) - b_prevent_fitting = true; + bool b_defined_center = get_keyval (group_conf, "centerReference", b_center, false, mode); + bool b_defined_rotate = get_keyval (group_conf, "rotateReference", b_rotate, false, mode); + + // this cannot be shortened to one statement because lazy evaluation may prevent one + // function from being called! + b_user_defined_fit = b_defined_center || b_defined_rotate; // if ((b_center || b_rotate) && b_dummy) // cvm::fatal_error ("Error: cannot set \"centerReference\" or " // "\"rotateReference\" when \"dummyAtom\" is defined.\n"); - // instead of this group, define another group (refPositionsGroup) to be the one - // used to fit the coordinates - if (key_lookup (group_conf, "refPositionsGroup")) { - if (ref_pos_group) { - cvm::fatal_error ("Error: the atom group \""+ - std::string (key)+"\" has already a reference group " - "for the rototranslational fit, which was communicated by the " - "colvar component. You should not use refPositionsGroup " - "in this case.\n"); + if (b_center || b_rotate) { + + // instead of this group, define another group (refPositionsGroup) to be the one + // used to fit the coordinates + if (key_lookup (group_conf, "refPositionsGroup")) { + if (ref_pos_group) { + cvm::fatal_error ("Error: the atom group \""+ + std::string (key)+"\" has already a reference group " + "for the rototranslational fit, which was communicated by the " + "colvar component. You should not use refPositionsGroup " + "in this case.\n"); + } + cvm::log ("Within atom group \""+std::string (key)+"\":\n"); + ref_pos_group = new atom_group (group_conf, "refPositionsGroup"); } - cvm::log ("Within atom group \""+std::string (key)+"\":\n"); - ref_pos_group = new atom_group (group_conf, "refPositionsGroup"); - } - atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this; + atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this; - if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) { - if (ref_pos.size() != group_for_fit->size()) { - cvm::fatal_error ("Error: the number of reference positions provided ("+ - cvm::to_str (ref_pos.size())+ - ") does not match the number of atoms of group \""+ - std::string (key)+ - "\" ("+cvm::to_str (group_for_fit->size())+").\n"); + if (get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode)) { + if (ref_pos.size() != group_for_fit->size()) { + cvm::fatal_error ("Error: the number of reference positions provided ("+ + cvm::to_str (ref_pos.size())+ + ") does not match the number of atoms of group \""+ + std::string (key)+ + "\" ("+cvm::to_str (group_for_fit->size())+").\n"); + } } - } - { std::string ref_pos_file; if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) { @@ -284,50 +286,40 @@ void cvm::atom_group::parse (std::string const &conf, cvm::load_coords (ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids, ref_pos_col, ref_pos_col_value); } - } - if (ref_pos.size()) { + if (ref_pos.size()) { - if (b_rotate) { - if (ref_pos.size() != group_for_fit->size()) - cvm::fatal_error ("Error: the number of reference positions provided ("+ - cvm::to_str (ref_pos.size())+ - ") does not match the number of atoms within \""+ - std::string (key)+ - "\" ("+cvm::to_str (group_for_fit->size())+ - "): to perform a rotational fit, "+ - "these numbers should be equal.\n"); - } + if (b_rotate) { + if (ref_pos.size() != group_for_fit->size()) + cvm::fatal_error ("Error: the number of reference positions provided ("+ + cvm::to_str (ref_pos.size())+ + ") does not match the number of atoms within \""+ + std::string (key)+ + "\" ("+cvm::to_str (group_for_fit->size())+ + "): to perform a rotational fit, "+ + "these numbers should be equal.\n"); + } - // save the center of geometry of ref_pos and then subtract it from - // them; in this way it will be possible to use ref_pos also for - // the rotational fit - ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0); - std::vector::iterator pi = ref_pos.begin(); - for ( ; pi != ref_pos.end(); pi++) { - ref_pos_cog += *pi; - } - ref_pos_cog /= (cvm::real) ref_pos.size(); - for (std::vector::iterator pi = ref_pos.begin(); - pi != ref_pos.end(); pi++) { - (*pi) -= ref_pos_cog; - } - } else { + // save the center of geometry of ref_pos and subtract it + center_ref_pos(); + + } else { #if (! defined (COLVARS_STANDALONE)) - cvm::fatal_error ("Error: no reference positions provided.\n"); + cvm::fatal_error ("Error: no reference positions provided.\n"); #endif - } + } - if (b_rotate && !noforce) { - cvm::log ("Warning: atom group \""+std::string (key)+ - "\" will be fitted automatically onto a fixed orientation: " - "in few cases, torques applied on this group " - "may make the simulation unstable. " - "If this happens, set \"disableForces\" to yes " - "for this group.\n"); - // initialize rot member data - rot.request_group1_gradients (this->size()); + if (b_rotate && !noforce) { + cvm::log ("Warning: atom group \""+std::string (key)+ + "\" will be fitted automatically onto a fixed orientation: " + "in few cases, torques applied on this group " + "may make the simulation unstable. " + "If this happens, set \"disableForces\" to yes " + "for this group.\n"); + // initialize rot member data + rot.request_group1_gradients (this->size()); + } } if (cvm::debug()) @@ -410,6 +402,24 @@ void cvm::atom_group::create_sorted_ids (void) return; } +void cvm::atom_group::center_ref_pos() +{ + // save the center of geometry of ref_pos and then subtract it from + // them; in this way it will be possible to use ref_pos also for + // the rotational fit + // This is called either by atom_group::parse or by CVCs that set + // reference positions (eg. RMSD, eigenvector) + ref_pos_cog = cvm::atom_pos (0.0, 0.0, 0.0); + std::vector::iterator pi = ref_pos.begin(); + for ( ; pi != ref_pos.end(); pi++) { + ref_pos_cog += *pi; + } + ref_pos_cog /= (cvm::real) ref_pos.size(); + for (std::vector::iterator pi = ref_pos.begin(); + pi != ref_pos.end(); pi++) { + (*pi) -= ref_pos_cog; + } +} void cvm::atom_group::read_positions() { diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index 7f1f002956..6ffda6db56 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -156,10 +156,10 @@ public: /// Rotation between the group and its reference coordinates cvm::rotation rot; - /// \brief Indicates that b_center and b_rotate are not false by default, - /// but following the user's choice instead: cvc's will not override - /// the value of b_center or b_rotate - bool b_prevent_fitting; + /// \brief Indicates that the user has specified centerReference or + /// rotateReference (and the related reference coordinates): + /// cvc's (eg rmsd, eigenvector) should not override these settings + bool b_user_defined_fit; /// \brief use these reference coordinates, for b_center, b_rotate, or to compute /// certain colvar components (orientation, rmsd, etc) @@ -212,6 +212,10 @@ public: /// them void read_positions(); + /// \brief Save center of geometry fo ref positions, + /// then subtract it + void center_ref_pos(); + /// \brief Move all positions void apply_translation (cvm::rvector const &t); diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 82b7b8aa6d..bc77c5a4b3 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -736,10 +736,11 @@ colvar::rmsd::rmsd (std::string const &conf) if (atoms.b_dummy) cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom."); - if (!atoms.b_prevent_fitting) { - // fit everything, unless the user made an explicit choice - cvm::log ("No value was specified within \"atoms\" for \"centerReference\" or " - "\"rotateReference\": enabling both.\n"); + if (atoms.b_user_defined_fit) { + cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\"."); + cvm::log ("Not computing standard minimum RMSD (if that is the desired result, leave atom group parameters at default values)."); + } else { + // Default: fit everything // NOTE: this won't work for class V cvc's atoms.b_center = true; atoms.b_rotate = true; @@ -792,13 +793,19 @@ colvar::rmsd::rmsd (std::string const &conf) // if not, rely on existing atom indices for the group atoms.create_sorted_ids(); } + + ref_pos.resize (atoms.size()); cvm::load_coords (ref_pos_file.c_str(), ref_pos, atoms.sorted_ids, ref_pos_col, ref_pos_col_value); } } - - + if (!atoms.b_user_defined_fit) { + // Default case: reference positions for calculating the rmsd are also those used + // for fitting + atoms.ref_pos = ref_pos; + atoms.center_ref_pos(); + } } @@ -838,7 +845,7 @@ void colvar::rmsd::calc_gradients() 0.0; for (size_t ia = 0; ia < atoms.size(); ia++) { - atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - atoms.ref_pos[ia])); + atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - ref_pos[ia])); } } @@ -967,6 +974,7 @@ colvar::eigenvector::eigenvector (std::string const &conf) atoms.b_center = true; atoms.b_rotate = true; atoms.ref_pos = ref_pos; + atoms.center_ref_pos(); // now load the eigenvector if (get_keyval (conf, "vector", eigenvec, eigenvec)) { @@ -1008,8 +1016,9 @@ colvar::eigenvector::eigenvector (std::string const &conf) for (size_t i = 0; i < atoms.size(); i++) { center += eigenvec[i]; } + center *= 1.0 / atoms.size(); - cvm::log ("Subtracting sum of eigenvector components: " + cvm::to_str (center) + "\n"); + cvm::log ("Subtracting center of eigenvector components: " + cvm::to_str (center) + "\n"); for (size_t i = 0; i < atoms.size(); i++) { eigenvec[i] = eigenvec[i] - center; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index dfc0578ebd..2375fafa39 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -2,7 +2,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2012-10-03" +#define COLVARS_VERSION "2012-11-20" #endif #ifndef COLVARS_DEBUG