jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@9130 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2012-11-30 18:03:00 +00:00
parent f3e6110453
commit 87215c925c
4 changed files with 101 additions and 78 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

140
lib/colvars/colvaratoms.cpp
View File

@ -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<cvm::atom_pos>::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<cvm::atom_pos>::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<cvm::atom_pos>::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<cvm::atom_pos>::iterator pi = ref_pos.begin();
pi != ref_pos.end(); pi++) {
(*pi) -= ref_pos_cog;
}
}
void cvm::atom_group::read_positions()
{

12
lib/colvars/colvaratoms.h
View File

@ -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);

25
lib/colvars/colvarcomp_distances.cpp
View File

@ -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;

2
lib/colvars/colvarmodule.h
View File

@ -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