diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp
index 0e67cf3ba..7a3ffa190 100644
--- a/lib/colvars/colvaratoms.cpp
+++ b/lib/colvars/colvaratoms.cpp
@@ -1,771 +1,781 @@
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvaratoms.h"
// atom member functions depend tightly on the MD interface, and are
// thus defined in colvarproxy_xxx.cpp
// atom_group member functions
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)
{
cvm::log ("Defining atom group \""+
std::string (key)+"\".\n");
parse (conf, key, ref_pos_group_in);
}
void cvm::atom_group::parse (std::string const &conf,
char const *key,
atom_group *ref_pos_group_in)
{
std::string group_conf;
// save_delimiters is set to false for this call, because "conf" is
// not the config string of this group, but of its parent object
// (which has already taken care of the delimiters)
save_delimiters = false;
key_lookup (conf, key, group_conf, dummy_pos);
// restoring the normal value, because we do want keywords checked
// inside "group_conf"
save_delimiters = true;
if (group_conf.size() == 0) {
cvm::fatal_error ("Error: atom group \""+
std::string (key)+"\" is set, but "
"has no definition.\n");
}
cvm::increase_depth();
cvm::log ("Initializing atom group \""+std::string (key)+"\".\n");
// whether or not to include messages in the log
// colvarparse::Parse_Mode mode = parse_silent;
// {
// bool b_verbose;
// get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
// if (b_verbose) mode = parse_normal;
// }
colvarparse::Parse_Mode mode = parse_normal;
{
// get the atoms by numbers
std::vector<int> atom_indexes;
if (get_keyval (group_conf, "atomNumbers", atom_indexes, atom_indexes, mode)) {
if (atom_indexes.size()) {
this->reserve (this->size()+atom_indexes.size());
for (size_t i = 0; i < atom_indexes.size(); i++) {
this->push_back (cvm::atom (atom_indexes[i]));
}
} else
cvm::fatal_error ("Error: no numbers provided for \""
"atomNumbers\".\n");
}
}
{
std::string range_conf = "";
size_t pos = 0;
while (key_lookup (group_conf, "atomNumbersRange",
range_conf, pos)) {
if (range_conf.size()) {
std::istringstream is (range_conf);
int initial, final;
char dash;
if ( (is >> initial) && (initial > 0) &&
(is >> dash) && (dash == '-') &&
(is >> final) && (final > 0) ) {
for (int anum = initial; anum <= final; anum++) {
this->push_back (cvm::atom (anum));
}
range_conf = "";
continue;
}
}
cvm::fatal_error ("Error: no valid definition for \""
"atomNumbersRange\", \""+
range_conf+"\".\n");
}
}
std::vector<std::string> psf_segids;
get_keyval (group_conf, "psfSegID", psf_segids, std::vector<std::string> (), mode);
for (std::vector<std::string>::iterator psii = psf_segids.begin();
psii < psf_segids.end(); psii++) {
if ( (psii->size() == 0) || (psii->size() > 4) ) {
cvm::fatal_error ("Error: invalid segmend identifier provided, \""+
(*psii)+"\".\n");
}
}
{
std::string range_conf = "";
size_t pos = 0;
size_t range_count = 0;
std::vector<std::string>::iterator psii = psf_segids.begin();
while (key_lookup (group_conf, "atomNameResidueRange",
range_conf, pos)) {
range_count++;
if (range_count > psf_segids.size()) {
cvm::fatal_error ("Error: more instances of \"atomNameResidueRange\" than "
"values of \"psfSegID\".\n");
}
std::string const &psf_segid = psf_segids.size() ? *psii : std::string ("");
if (range_conf.size()) {
std::istringstream is (range_conf);
std::string atom_name;
int initial, final;
char dash;
if ( (is >> atom_name) && (atom_name.size()) &&
(is >> initial) && (initial > 0) &&
(is >> dash) && (dash == '-') &&
(is >> final) && (final > 0) ) {
for (int resid = initial; resid <= final; resid++) {
this->push_back (cvm::atom (resid, atom_name, psf_segid));
}
range_conf = "";
} else {
cvm::fatal_error ("Error: cannot parse definition for \""
"atomNameResidueRange\", \""+
range_conf+"\".\n");
}
} else {
cvm::fatal_error ("Error: atomNameResidueRange with empty definition.\n");
}
if (psf_segid.size())
psii++;
}
}
{
// read the atoms from a file
std::string atoms_file_name;
if (get_keyval (group_conf, "atomsFile", atoms_file_name, std::string (""), mode)) {
std::string atoms_col;
if (!get_keyval (group_conf, "atomsCol", atoms_col, std::string (""), mode)) {
cvm::fatal_error ("Error: parameter atomsCol is required if atomsFile is set.\n");
}
double atoms_col_value;
bool const atoms_col_value_defined = get_keyval (group_conf, "atomsColValue", atoms_col_value, 0.0, mode);
if (atoms_col_value_defined && (!atoms_col_value))
cvm::fatal_error ("Error: atomsColValue, "
"if provided, must be non-zero.\n");
cvm::load_atoms (atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
}
}
for (std::vector<cvm::atom>::iterator a1 = this->begin();
a1 != this->end(); a1++) {
std::vector<cvm::atom>::iterator a2 = a1;
++a2;
for ( ; a2 != this->end(); a2++) {
if (a1->id == a2->id) {
if (cvm::debug())
cvm::log ("Discarding doubly counted atom with number "+
cvm::to_str (a1->id+1)+".\n");
a2 = this->erase (a2);
if (a2 == this->end())
break;
}
}
}
if (get_keyval (group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) {
b_dummy = true;
this->total_mass = 1.0;
} else
b_dummy = false;
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");
#if (! defined (COLVARS_STANDALONE))
if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) {
cvm::fatal_error ("Error: no atoms defined for atom group \""+
std::string (key)+"\".\n");
}
#endif
if (!b_dummy) {
this->total_mass = 0.0;
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
this->total_mass += ai->mass;
}
}
if (!b_dummy)
get_keyval (group_conf, "disableForces", noforce, false, mode);
// 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)) {
if (ref_pos.size()) {
cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
"\"refPositions\" at the same time.\n");
}
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);
if (found && !ref_pos_col_value)
cvm::fatal_error ("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
} else {
// if not, rely on existing atom indices for the group
group_for_fit->create_sorted_ids();
}
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;
- }
+ // save the center of geometry of ref_pos and subtract it
+ center_ref_pos();
- } else {
+ } 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())
cvm::log ("Done initializing atom group with name \""+
std::string (key)+"\".\n");
this->check_keywords (group_conf, key);
cvm::log ("Atom group \""+std::string (key)+"\" defined, "+
cvm::to_str (this->size())+" initialized: total mass = "+
cvm::to_str (this->total_mass)+".\n");
cvm::decrease_depth();
}
cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
: b_dummy (false), b_center (false), b_rotate (false),
ref_pos_group (NULL), noforce (false)
{
this->reserve (atoms.size());
for (size_t i = 0; i < atoms.size(); i++) {
this->push_back (atoms[i]);
}
total_mass = 0.0;
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
total_mass += ai->mass;
}
}
cvm::atom_group::atom_group()
: b_dummy (false), b_center (false), b_rotate (false),
ref_pos_group (NULL), noforce (false)
{
total_mass = 0.0;
}
cvm::atom_group::~atom_group()
{
if (ref_pos_group) {
delete ref_pos_group;
ref_pos_group = NULL;
}
}
void cvm::atom_group::add_atom (cvm::atom const &a)
{
if (b_dummy) {
cvm::fatal_error ("Error: cannot add atoms to a dummy group.\n");
} else {
this->push_back (a);
total_mass += a.mass;
}
}
void cvm::atom_group::create_sorted_ids (void)
{
// Only do the work if the vector is not yet populated
if (sorted_ids.size())
return;
std::list<int> temp_id_list;
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
temp_id_list.push_back (ai->id);
}
temp_id_list.sort();
temp_id_list.unique();
if (temp_id_list.size() != this->size()) {
cvm::fatal_error ("Error: duplicate atom IDs in atom group? (found " +
cvm::to_str(temp_id_list.size()) +
" unique atom IDs instead of" +
cvm::to_str(this->size()) + ").\n");
}
sorted_ids = std::vector<int> (temp_id_list.begin(), temp_id_list.end());
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()
{
if (b_dummy) return;
#if (! defined (COLVARS_STANDALONE))
if (!this->size())
cvm::fatal_error ("Error: no atoms defined in the requested group.\n");
#endif
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->read_position();
}
if (ref_pos_group)
ref_pos_group->read_positions();
atom_group *fit_group = ref_pos_group ? ref_pos_group : this;
if (b_center) {
// store aside the current center of geometry (all positions will be
// set to the closest images to the first one) and then center on
// the origin
cvm::atom_pos const cog = fit_group->center_of_geometry();
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->pos -= cog;
}
}
if (b_rotate) {
// rotate the group (around the center of geometry if b_center is
// true, around the origin otherwise); store the rotation, in
// order to bring back the forces to the original frame before
// applying them
rot.calc_optimal_rotation (fit_group->positions(), ref_pos);
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->pos = rot.rotate (ai->pos);
}
}
if (b_center) {
// use the center of geometry of ref_pos
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->pos += ref_pos_cog;
}
}
}
void cvm::atom_group::apply_translation (cvm::rvector const &t)
{
if (b_dummy) return;
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->pos += t;
}
}
void cvm::atom_group::apply_rotation (cvm::rotation const &rot)
{
if (b_dummy) return;
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->pos = rot.rotate (ai->pos);
}
}
void cvm::atom_group::read_velocities()
{
if (b_dummy) return;
if (b_rotate) {
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->read_velocity();
ai->vel = rot.rotate (ai->vel);
}
} else {
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->read_velocity();
}
}
}
void cvm::atom_group::read_system_forces()
{
if (b_dummy) return;
if (b_rotate) {
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->read_system_force();
ai->system_force = rot.rotate (ai->system_force);
}
} else {
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->read_system_force();
}
}
}
/* This is deprecated.
cvm::atom_pos cvm::atom_group::center_of_geometry (cvm::atom_pos const &ref_pos)
{
if (b_dummy) {
cvm::select_closest_image (dummy_atom_pos, ref_pos);
return dummy_atom_pos;
}
cvm::atom_pos cog (0.0, 0.0, 0.0);
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
cvm::select_closest_image (ai->pos, ref_pos);
cog += ai->pos;
}
cog /= this->size();
return cog;
} */
cvm::atom_pos cvm::atom_group::center_of_geometry() const
{
if (b_dummy)
return dummy_atom_pos;
cvm::atom_pos cog (0.0, 0.0, 0.0);
for (cvm::atom_const_iter ai = this->begin();
ai != this->end(); ai++) {
cog += ai->pos;
}
cog /= this->size();
return cog;
}
/* This is deprecated.
cvm::atom_pos cvm::atom_group::center_of_mass (cvm::atom_pos const &ref_pos)
{
if (b_dummy) {
cvm::select_closest_image (dummy_atom_pos, ref_pos);
return dummy_atom_pos;
}
cvm::atom_pos com (0.0, 0.0, 0.0);
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
cvm::select_closest_image (ai->pos, ref_pos);
com += ai->mass * ai->pos;
}
com /= this->total_mass;
return com;
}*/
cvm::atom_pos cvm::atom_group::center_of_mass() const
{
if (b_dummy)
return dummy_atom_pos;
cvm::atom_pos com (0.0, 0.0, 0.0);
for (cvm::atom_const_iter ai = this->begin();
ai != this->end(); ai++) {
com += ai->mass * ai->pos;
}
com /= this->total_mass;
return com;
}
void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad)
{
if (b_dummy) return;
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->grad = (ai->mass/this->total_mass) * grad;
}
}
std::vector<cvm::atom_pos> cvm::atom_group::positions() const
{
if (b_dummy)
cvm::fatal_error ("Error: positions are not available "
"from a dummy atom group.\n");
std::vector<cvm::atom_pos> x (this->size(), 0.0);
cvm::atom_const_iter ai = this->begin();
std::vector<cvm::atom_pos>::iterator xi = x.begin();
for ( ; ai != this->end(); xi++, ai++) {
*xi = ai->pos;
}
return x;
}
std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted (cvm::rvector const &shift) const
{
if (b_dummy)
cvm::fatal_error ("Error: positions are not available "
"from a dummy atom group.\n");
std::vector<cvm::atom_pos> x (this->size(), 0.0);
cvm::atom_const_iter ai = this->begin();
std::vector<cvm::atom_pos>::iterator xi = x.begin();
for ( ; ai != this->end(); xi++, ai++) {
*xi = (ai->pos + shift);
}
return x;
}
std::vector<cvm::rvector> cvm::atom_group::velocities() const
{
if (b_dummy)
cvm::fatal_error ("Error: velocities are not available "
"from a dummy atom group.\n");
std::vector<cvm::rvector> v (this->size(), 0.0);
cvm::atom_const_iter ai = this->begin();
std::vector<cvm::atom_pos>::iterator vi = v.begin();
for ( ; ai != this->end(); vi++, ai++) {
*vi = ai->vel;
}
return v;
}
std::vector<cvm::rvector> cvm::atom_group::system_forces() const
{
if (b_dummy)
cvm::fatal_error ("Error: system forces are not available "
"from a dummy atom group.\n");
std::vector<cvm::rvector> f (this->size(), 0.0);
cvm::atom_const_iter ai = this->begin();
std::vector<cvm::atom_pos>::iterator fi = f.begin();
for ( ; ai != this->end(); fi++, ai++) {
*fi = ai->system_force;
}
return f;
}
cvm::rvector cvm::atom_group::system_force() const
{
if (b_dummy)
cvm::fatal_error ("Error: system forces are not available "
"from a dummy atom group.\n");
cvm::rvector f (0.0);
for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
f += ai->system_force;
}
return f;
}
void cvm::atom_group::apply_colvar_force (cvm::real const &force)
{
if (b_dummy)
return;
if (noforce)
cvm::fatal_error ("Error: sending a force to a group that has "
"\"disableForces\" defined.\n");
if (b_rotate) {
// get the forces back from the rotated frame
cvm::rotation const rot_inv = rot.inverse();
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->apply_force (rot_inv.rotate (force * ai->grad));
}
} else {
// no need to manipulate gradients, they are still in the original
// frame
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->apply_force (force * ai->grad);
}
}
}
void cvm::atom_group::apply_force (cvm::rvector const &force)
{
if (b_dummy)
return;
if (noforce)
cvm::fatal_error ("Error: sending a force to a group that has "
"\"disableForces\" defined.\n");
if (b_rotate) {
cvm::rotation const rot_inv = rot.inverse();
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->apply_force (rot_inv.rotate ((ai->mass/this->total_mass) * force));
}
} else {
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->apply_force ((ai->mass/this->total_mass) * force);
}
}
}
void cvm::atom_group::apply_forces (std::vector<cvm::rvector> const &forces)
{
if (b_dummy)
return;
if (noforce)
cvm::fatal_error ("Error: sending a force to a group that has "
"\"disableForces\" defined.\n");
if (forces.size() != this->size()) {
cvm::fatal_error ("Error: trying to apply an array of forces to an atom "
"group which does not have the same length.\n");
}
if (b_rotate) {
cvm::rotation const rot_inv = rot.inverse();
cvm::atom_iter ai = this->begin();
std::vector<cvm::rvector>::const_iterator fi = forces.begin();
for ( ; ai != this->end(); fi++, ai++) {
ai->apply_force (rot_inv.rotate (*fi));
}
} else {
cvm::atom_iter ai = this->begin();
std::vector<cvm::rvector>::const_iterator fi = forces.begin();
for ( ; ai != this->end(); fi++, ai++) {
ai->apply_force (*fi);
}
}
}
diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h
index 7f1f00295..6ffda6db5 100644
--- a/lib/colvars/colvaratoms.h
+++ b/lib/colvars/colvaratoms.h
@@ -1,318 +1,322 @@
// -*- c++ -*-
#ifndef COLVARATOMS_H
#define COLVARATOMS_H
#include "colvarmodule.h"
#include "colvarparse.h"
/// \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
/// numeric id, all acting independently: forces communicated through
/// these instances will be summed together.
///
/// Read/write operations depend on the underlying code: hence, some
/// member functions are defined in colvarproxy_xxx.h.
class colvarmodule::atom {
protected:
/// \brief Index in the list of atoms involved by the colvars (\b
/// NOT in the global topology!)
size_t index;
public:
/// Internal identifier (zero-based)
int id;
/// Mass
cvm::real mass;
/// \brief Current position (copied from the program, can be
/// manipulated)
cvm::atom_pos pos;
/// \brief Current velocity (copied from the program, can be
/// manipulated)
cvm::rvector vel;
/// \brief System force at the previous step (copied from the
/// program, can be manipulated)
cvm::rvector system_force;
/// \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
/// most frequent case. For more complex types of \link
/// colvarvalue \endlink objects, atomic gradients should be
/// defined within the specific \link cvc \endlink
/// implementation
cvm::rvector grad;
/// \brief Default constructor, setting id to a non-valid value
inline atom() {}
/// \brief Initialize an atom for collective variable calculation
/// and get its internal identifier \param atom_number Atom index in
/// the system topology (starting from 1)
atom (int const &atom_number);
/// \brief Initialize an atom for collective variable calculation
/// and get its internal identifier \param residue Residue number
/// \param atom_name Name of the atom in the residue \param
/// segment_id For PSF topologies, the segment identifier; for other
/// type of topologies, may not be required
atom (cvm::residue_id const &residue,
std::string const &atom_name,
std::string const &segment_id = std::string (""));
/// Copy constructor
atom (atom const &a);
/// Destructor
~atom();
/// Set non-constant data (everything except id and mass) to zero
inline void reset_data() {
pos = atom_pos (0.0);
vel = grad = system_force = rvector (0.0);
}
/// Get the current position
void read_position();
/// Get the current velocity
void read_velocity();
/// Get the system force
void read_system_force();
/// \brief Apply a force to the atom
///
/// The force will be used later by the MD integrator, the
/// collective variables module does not integrate equations of
/// motion. Multiple calls to this function by either the same
/// \link atom \endlink object or different objects with identical
/// \link id \endlink, will all add to the existing MD force.
void apply_force (cvm::rvector const &new_force);
};
/// \brief Group of \link atom \endlink objects, mostly used by a
/// \link cvc \endlink
///
/// This class inherits from \link colvarparse \endlink and from
/// std::vector<colvarmodule::atom>, and hence all functions and
/// operators (including the bracket operator, group[i]) can be used
/// on an \link atom_group \endlink object. It can be initialized as
/// a vector, or by parsing a keyword in the configuration.
class colvarmodule::atom_group
: public std::vector<cvm::atom>,
public colvarparse
{
public:
// Note: all members here are kept public, to make possible to any
// object accessing and manipulating them
/// \brief If this option is on, this group merely acts as a wrapper
/// for a fixed position; any calls to atoms within or to
/// functions that return disaggregated data will fail
bool b_dummy;
/// \brief dummy atom position
cvm::atom_pos dummy_atom_pos;
/// Sorted list of zero-based (internal) atom ids
/// (populated on-demand by create_sorted_ids)
std::vector<int> sorted_ids;
/// Allocates and populates the sorted list of atom ids
void create_sorted_ids (void);
/// \brief Before calculating colvars, move the group to overlap the
/// center of mass of reference coordinates
bool b_center;
/// \brief Right after updating atom coordinates (and after
/// centering coordinates, if b_center is true), rotate the group to
/// overlap the reference coordinates. You should not manipulate
/// atoms individually if you turn on this flag.
///
/// Note: gradients will be calculated in the rotated frame: when
/// forces will be applied, they will rotated back to the original
/// frame
bool b_rotate;
/// 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)
std::vector<cvm::atom_pos> ref_pos;
/// \brief Center of geometry of the reference coordinates; regardless
/// of whether b_center is true, ref_pos is centered to zero at
/// initialization, and ref_pos_cog serves to center the positions
cvm::atom_pos ref_pos_cog;
/// \brief In case b_center or b_rotate is true, fit this group to
/// the reference positions (default: the parent group itself)
atom_group *ref_pos_group;
/// Total mass of the atom group
cvm::real total_mass;
/// \brief Don't apply any force on this group (use its coordinates
/// only to calculate a colvar)
bool noforce;
/// \brief Initialize the group by looking up its configuration
/// string in conf and parsing it; this is actually done by parse(),
/// which is a member function so that a group can be initialized
/// also after construction
atom_group (std::string const &conf,
char const *key,
atom_group *ref_pos_group = NULL);
/// \brief Initialize the group by looking up its configuration
/// string in conf and parsing it
void parse (std::string const &conf,
char const *key,
atom_group *ref_pos_group = NULL);
/// \brief Initialize the group after a temporary vector of atoms
atom_group (std::vector<cvm::atom> const &atoms);
/// \brief Add an atom to this group
void add_atom (cvm::atom const &a);
/// \brief Default constructor
atom_group();
/// \brief Destructor
~atom_group();
/// \brief Get the current positions; if b_center or b_rotate are
/// true, center and/or rotate the coordinates right after reading
/// 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);
/// \brief Rotate all positions
void apply_rotation (cvm::rotation const &q);
/// \brief Get the current velocities; this must be called always
/// *after* read_positions(); if b_rotate is defined, the same
/// rotation applied to the coordinates will be used
void read_velocities();
/// \brief Get the current system_forces; this must be called always
/// *after* read_positions(); if b_rotate is defined, the same
/// rotation applied to the coordinates will be used
void read_system_forces();
/// Call reset_data() for each atom
inline void reset_atoms_data()
{
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++)
ai->reset_data();
}
/// \brief Return a copy of the current atom positions
std::vector<cvm::atom_pos> positions() const;
/// \brief Return a copy of the current atom positions, shifted by a constant vector
std::vector<cvm::atom_pos> positions_shifted (cvm::rvector const &shift) const;
/// \brief Return the center of geometry of the positions \param ref_pos
/// Use the closest periodic images to this position
cvm::atom_pos center_of_geometry (cvm::atom_pos const &ref_pos);
/// \brief Return the center of geometry of the positions, assuming
/// that coordinates are already pbc-wrapped
cvm::atom_pos center_of_geometry() const;
/// \brief Return the center of mass of the positions \param ref_pos
/// Use the closest periodic images to this position
cvm::atom_pos center_of_mass (cvm::atom_pos const &ref_pos);
/// \brief Return the center of mass of the positions, assuming that
/// coordinates are already pbc-wrapped
cvm::atom_pos center_of_mass() const;
/// \brief Store atom positions from the previous step
std::vector<cvm::atom_pos> old_pos;
/// \brief Return a copy of the current atom velocities
std::vector<cvm::rvector> velocities() const;
/// \brief Return a copy of the system forces
std::vector<cvm::rvector> system_forces() const;
/// \brief Return a copy of the aggregated total force on the group
cvm::rvector system_force() const;
/// \brief Shorthand: save the specified gradient on each atom,
/// weighting with the atom mass (mostly used in combination with
/// \link center_of_mass() \endlink)
void set_weighted_gradient (cvm::rvector const &grad);
/// \brief Used by a (scalar) colvar to apply its force on its \link
/// atom_group \endlink members
///
/// The (scalar) force is multiplied by the colvar gradient for each
/// atom; this should be used when a colvar with scalar \link
/// colvarvalue \endlink type is used (this is the most frequent
/// case: for colvars with a non-scalar type, the most convenient
/// solution is to sum together the Cartesian forces from all the
/// colvar components, and use apply_force() or apply_forces()). If
/// the group is being rotated to a reference frame (e.g. to express
/// the colvar independently from the solute rotation), the
/// gradients are temporarily to the original frame.
void apply_colvar_force (cvm::real const &force);
/// \brief Apply a force "to the center of mass", i.e. the force is
/// distributed on each atom according to its mass
///
/// If the group is being rotated to a reference frame (e.g. to
/// express the colvar independently from the solute rotation), the
/// force is rotated back to the original frame. Colvar gradients
/// are not used, either because they were not defined (e.g because
/// the colvar has not a scalar value) or the biases require to
/// micromanage the force.
void apply_force (cvm::rvector const &force);
/// \brief Apply an array of forces directly on the individual
/// atoms; the length of the specified vector must be the same of
/// this \link atom_group \endlink.
///
/// If the group is being rotated to a reference frame (e.g. to
/// express the colvar independently from the solute rotation), the
/// forces are rotated back to the original frame. Colvar gradients
/// are not used, either because they were not defined (e.g because
/// the colvar has not a scalar value) or the biases require to
/// micromanage the forces.
void apply_forces (std::vector<cvm::rvector> const &forces);
};
#endif
diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp
index 82b7b8aa6..bc77c5a4b 100644
--- a/lib/colvars/colvarcomp_distances.cpp
+++ b/lib/colvars/colvarcomp_distances.cpp
@@ -1,1189 +1,1198 @@
#include <cmath>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
#include "colvar.h"
#include "colvarcomp.h"
/// \file cvc_distance.cpp \brief Collective variables
/// determining various type of distances between two groups
// "twogroup" flag defaults to true; set to false by selfCoordNum
// (only distance-derived component based on only one group)
colvar::distance::distance (std::string const &conf, bool twogroups)
: cvc (conf)
{
function_type = "distance";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) {
cvm::log ("Computing distance using absolute positions (not minimal-image)");
}
if (twogroups && get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
cvm::log ("Computing system force on group 1 only");
}
parse_group (conf, "group1", group1);
atom_groups.push_back (&group1);
if (twogroups) {
parse_group (conf, "group2", group2);
atom_groups.push_back (&group2);
}
x.type (colvarvalue::type_scalar);
}
colvar::distance::distance()
: cvc ()
{
function_type = "distance";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
b_1site_force = false;
x.type (colvarvalue::type_scalar);
}
void colvar::distance::calc_value()
{
group1.reset_atoms_data();
group2.reset_atoms_data();
group1.read_positions();
group2.read_positions();
if (b_no_PBC) {
dist_v = group2.center_of_mass() - group1.center_of_mass();
} else {
dist_v = cvm::position_distance (group1.center_of_mass(),
group2.center_of_mass());
}
x.real_value = dist_v.norm();
}
void colvar::distance::calc_gradients()
{
cvm::rvector const u = dist_v.unit();
group1.set_weighted_gradient (-1.0 * u);
group2.set_weighted_gradient ( u);
}
void colvar::distance::calc_force_invgrads()
{
group1.read_system_forces();
if ( b_1site_force ) {
ft.real_value = -1.0 * (group1.system_force() * dist_v.unit());
} else {
group2.read_system_forces();
ft.real_value = 0.5 * ((group2.system_force() - group1.system_force()) * dist_v.unit());
}
}
void colvar::distance::calc_Jacobian_derivative()
{
jd.real_value = x.real_value ? (2.0 / x.real_value) : 0.0;
}
void colvar::distance::apply_force (colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_colvar_force (force.real_value);
if (!group2.noforce)
group2.apply_colvar_force (force.real_value);
}
colvar::distance_vec::distance_vec (std::string const &conf)
: distance (conf)
{
function_type = "distance_vec";
x.type (colvarvalue::type_vector);
}
colvar::distance_vec::distance_vec()
: distance()
{
function_type = "distance_vec";
x.type (colvarvalue::type_vector);
}
void colvar::distance_vec::calc_value()
{
group1.reset_atoms_data();
group2.reset_atoms_data();
group1.read_positions();
group2.read_positions();
if (b_no_PBC) {
x.rvector_value = group2.center_of_mass() - group1.center_of_mass();
} else {
x.rvector_value = cvm::position_distance (group1.center_of_mass(),
group2.center_of_mass());
}
}
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
}
void colvar::distance_vec::apply_force (colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_force (-1.0 * force.rvector_value);
if (!group2.noforce)
group2.apply_force ( force.rvector_value);
}
colvar::distance_z::distance_z (std::string const &conf)
: cvc (conf)
{
function_type = "distance_z";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
x.type (colvarvalue::type_scalar);
// TODO detect PBC from MD engine (in simple cases)
// and then update period in real time
if (period != 0.0)
b_periodic = true;
if ((wrap_center != 0.0) && (period == 0.0)) {
cvm::fatal_error ("Error: wrapAround was defined in a distanceZ component,"
" but its period has not been set.\n");
}
parse_group (conf, "main", main);
parse_group (conf, "ref", ref1);
atom_groups.push_back (&main);
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\"");
fixed_axis = false;
if (key_lookup (conf, "axis"))
cvm::log ("Warning: explicit axis definition will be ignored!");
} else {
if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
if (axis.norm2() == 0.0)
cvm::fatal_error ("Axis vector is zero!");
if (axis.norm2() != 1.0) {
axis = axis.unit();
cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
}
}
fixed_axis = true;
}
if (get_keyval (conf, "forceNoPBC", b_no_PBC, false)) {
cvm::log ("Computing distance using absolute positions (not minimal-image)");
}
if (get_keyval (conf, "oneSiteSystemForce", b_1site_force, false)) {
cvm::log ("Computing system force on group \"main\" only");
}
}
colvar::distance_z::distance_z()
{
function_type = "distance_z";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
x.type (colvarvalue::type_scalar);
}
void colvar::distance_z::calc_value()
{
main.reset_atoms_data();
ref1.reset_atoms_data();
main.read_positions();
ref1.read_positions();
if (fixed_axis) {
if (b_no_PBC) {
dist_v = main.center_of_mass() - ref1.center_of_mass();
} else {
dist_v = cvm::position_distance (ref1.center_of_mass(),
main.center_of_mass());
}
} else {
ref2.reset_atoms_data();
ref2.read_positions();
if (b_no_PBC) {
dist_v = main.center_of_mass() -
(0.5 * (ref1.center_of_mass() + ref2.center_of_mass()));
axis = ref2.center_of_mass() - ref1.center_of_mass();
} else {
dist_v = cvm::position_distance (0.5 * (ref1.center_of_mass() +
ref2.center_of_mass()), main.center_of_mass());
axis = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
}
axis_norm = axis.norm();
axis = axis.unit();
}
x.real_value = axis * dist_v;
this->wrap (x);
}
void colvar::distance_z::calc_gradients()
{
main.set_weighted_gradient ( axis );
if (fixed_axis) {
ref1.set_weighted_gradient (-1.0 * axis);
} else {
if (b_no_PBC) {
ref1.set_weighted_gradient ( 1.0 / axis_norm * (main.center_of_mass() - ref2.center_of_mass() -
x.real_value * axis ));
ref2.set_weighted_gradient ( 1.0 / axis_norm * (ref1.center_of_mass() - main.center_of_mass() +
x.real_value * axis ));
} else {
ref1.set_weighted_gradient ( 1.0 / axis_norm * (
cvm::position_distance (ref2.center_of_mass(), main.center_of_mass()) - x.real_value * axis ));
ref2.set_weighted_gradient ( 1.0 / axis_norm * (
cvm::position_distance (main.center_of_mass(), ref1.center_of_mass()) + x.real_value * axis ));
}
}
}
void colvar::distance_z::calc_force_invgrads()
{
main.read_system_forces();
if (fixed_axis && !b_1site_force) {
ref1.read_system_forces();
ft.real_value = 0.5 * ((main.system_force() - ref1.system_force()) * axis);
} else {
ft.real_value = main.system_force() * axis;
}
}
void colvar::distance_z::calc_Jacobian_derivative()
{
jd.real_value = 0.0;
}
void colvar::distance_z::apply_force (colvarvalue const &force)
{
if (!ref1.noforce)
ref1.apply_colvar_force (force.real_value);
if (ref2.size() && !ref2.noforce)
ref2.apply_colvar_force (force.real_value);
if (!main.noforce)
main.apply_colvar_force (force.real_value);
}
colvar::distance_xy::distance_xy (std::string const &conf)
: distance_z (conf)
{
function_type = "distance_xy";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
x.type (colvarvalue::type_scalar);
}
colvar::distance_xy::distance_xy()
: distance_z()
{
function_type = "distance_xy";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
x.type (colvarvalue::type_scalar);
}
void colvar::distance_xy::calc_value()
{
ref1.reset_atoms_data();
main.reset_atoms_data();
ref1.read_positions();
main.read_positions();
if (b_no_PBC) {
dist_v = main.center_of_mass() - ref1.center_of_mass();
} else {
dist_v = cvm::position_distance (ref1.center_of_mass(),
main.center_of_mass());
}
if (!fixed_axis) {
ref2.reset_atoms_data();
ref2.read_positions();
if (b_no_PBC) {
v12 = ref2.center_of_mass() - ref1.center_of_mass();
} else {
v12 = cvm::position_distance (ref1.center_of_mass(), ref2.center_of_mass());
}
axis_norm = v12.norm();
axis = v12.unit();
}
dist_v_ortho = dist_v - (dist_v * axis) * axis;
x.real_value = dist_v_ortho.norm();
}
void colvar::distance_xy::calc_gradients()
{
// Intermediate quantity (r_P3 / r_12 where P is the projection
// of 3 (main) on the plane orthogonal to 12, containing 1 (ref1))
cvm::real A;
cvm::real x_inv;
if (x.real_value == 0.0) return;
x_inv = 1.0 / x.real_value;
if (fixed_axis) {
ref1.set_weighted_gradient (-1.0 * x_inv * dist_v_ortho);
main.set_weighted_gradient ( x_inv * dist_v_ortho);
} else {
if (b_no_PBC) {
v13 = main.center_of_mass() - ref1.center_of_mass();
} else {
v13 = cvm::position_distance (ref1.center_of_mass(), main.center_of_mass());
}
A = (dist_v * axis) / axis_norm;
ref1.set_weighted_gradient ( (A - 1.0) * x_inv * dist_v_ortho);
ref2.set_weighted_gradient ( -A * x_inv * dist_v_ortho);
main.set_weighted_gradient ( 1.0 * x_inv * dist_v_ortho);
}
}
void colvar::distance_xy::calc_force_invgrads()
{
main.read_system_forces();
if (fixed_axis && !b_1site_force) {
ref1.read_system_forces();
ft.real_value = 0.5 / x.real_value * ((main.system_force() - ref1.system_force()) * dist_v_ortho);
} else {
ft.real_value = 1.0 / x.real_value * main.system_force() * dist_v_ortho;
}
}
void colvar::distance_xy::calc_Jacobian_derivative()
{
jd.real_value = x.real_value ? (1.0 / x.real_value) : 0.0;
}
void colvar::distance_xy::apply_force (colvarvalue const &force)
{
if (!ref1.noforce)
ref1.apply_colvar_force (force.real_value);
if (ref2.size() && !ref2.noforce)
ref2.apply_colvar_force (force.real_value);
if (!main.noforce)
main.apply_colvar_force (force.real_value);
}
colvar::distance_dir::distance_dir (std::string const &conf)
: distance (conf)
{
function_type = "distance_dir";
x.type (colvarvalue::type_unitvector);
}
colvar::distance_dir::distance_dir()
: distance()
{
function_type = "distance_dir";
x.type (colvarvalue::type_unitvector);
}
void colvar::distance_dir::calc_value()
{
group1.reset_atoms_data();
group2.reset_atoms_data();
group1.read_positions();
group2.read_positions();
if (b_no_PBC) {
dist_v = group2.center_of_mass() - group1.center_of_mass();
} else {
dist_v = cvm::position_distance (group1.center_of_mass(),
group2.center_of_mass());
}
x.rvector_value = dist_v.unit();
}
void colvar::distance_dir::calc_gradients()
{
// gradients are computed on the fly within apply_force()
// Note: could be a problem if a future bias relies on gradient
// calculations...
}
void colvar::distance_dir::apply_force (colvarvalue const &force)
{
// remove the radial force component
cvm::real const iprod = force.rvector_value * x.rvector_value;
cvm::rvector const force_tang = force.rvector_value - iprod * x.rvector_value;
if (!group1.noforce)
group1.apply_force (-1.0 * force_tang);
if (!group2.noforce)
group2.apply_force ( force_tang);
}
colvar::distance_inv::distance_inv (std::string const &conf)
: distance (conf)
{
function_type = "distance_inv";
get_keyval (conf, "exponent", exponent, 6);
if (exponent%2) {
cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n");
}
if (exponent <= 0) {
cvm::fatal_error ("Error: negative or zero exponent provided.\n");
}
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
if (ai1->id == ai2->id)
cvm::fatal_error ("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n");
}
}
b_inverse_gradients = false;
b_Jacobian_derivative = false;
x.type (colvarvalue::type_scalar);
}
colvar::distance_inv::distance_inv()
{
function_type = "distance_inv";
exponent = 6;
b_inverse_gradients = false;
b_Jacobian_derivative = false;
b_1site_force = false;
x.type (colvarvalue::type_scalar);
}
void colvar::distance_inv::calc_value()
{
group1.reset_atoms_data();
group2.reset_atoms_data();
group1.read_positions();
group2.read_positions();
x.real_value = 0.0;
if (b_no_PBC) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
cvm::rvector const dv = ai2->pos - ai1->pos;
cvm::real const d2 = dv.norm2();
cvm::real dinv = 1.0;
for (int ne = 0; ne < exponent/2; ne++)
dinv *= 1.0/d2;
x.real_value += dinv;
cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
ai1->grad += -1.0 * dsumddv;
ai2->grad += dsumddv;
}
}
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos);
cvm::real const d2 = dv.norm2();
cvm::real dinv = 1.0;
for (int ne = 0; ne < exponent/2; ne++)
dinv *= 1.0/d2;
x.real_value += dinv;
cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv;
ai1->grad += -1.0 * dsumddv;
ai2->grad += dsumddv;
}
}
}
x.real_value *= 1.0 / cvm::real (group1.size() * group2.size());
x.real_value = std::pow (x.real_value, -1.0/(cvm::real (exponent)));
}
void colvar::distance_inv::calc_gradients()
{
cvm::real const dxdsum = (-1.0/(cvm::real (exponent))) * std::pow (x.real_value, exponent+1) / cvm::real (group1.size() * group2.size());
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
ai1->grad *= dxdsum;
}
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
ai2->grad *= dxdsum;
}
}
void colvar::distance_inv::apply_force (colvarvalue const &force)
{
if (!group1.noforce)
group1.apply_colvar_force (force.real_value);
if (!group2.noforce)
group2.apply_colvar_force (force.real_value);
}
colvar::gyration::gyration (std::string const &conf)
: cvc (conf)
{
function_type = "gyration";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
x.type (colvarvalue::type_scalar);
}
colvar::gyration::gyration()
{
function_type = "gyration";
b_inverse_gradients = true;
b_Jacobian_derivative = true;
x.type (colvarvalue::type_scalar);
}
void colvar::gyration::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions();
atoms.apply_translation ((-1.0) * atoms.center_of_geometry());
x.real_value = 0.0;
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
x.real_value += (ai->pos).norm2();
}
x.real_value = std::sqrt (x.real_value / cvm::real (atoms.size()));
}
void colvar::gyration::calc_gradients()
{
cvm::real const drdx = 1.0/(cvm::real (atoms.size()) * x.real_value);
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ai->grad = drdx * ai->pos;
}
}
void colvar::gyration::calc_force_invgrads()
{
atoms.read_system_forces();
cvm::real const dxdr = 1.0/x.real_value;
ft.real_value = 0.0;
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ft.real_value += dxdr * ai->pos * ai->system_force;
}
}
void colvar::gyration::calc_Jacobian_derivative()
{
jd = x.real_value ? (3.0 * cvm::real (atoms.size()) - 4.0) / x.real_value : 0.0;
}
void colvar::gyration::apply_force (colvarvalue const &force)
{
if (!atoms.noforce)
atoms.apply_colvar_force (force.real_value);
}
colvar::inertia::inertia (std::string const &conf)
: cvc (conf)
{
function_type = "inertia";
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
x.type (colvarvalue::type_scalar);
}
colvar::inertia::inertia()
{
function_type = "inertia";
x.type (colvarvalue::type_scalar);
}
void colvar::inertia::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions();
atoms.apply_translation ((-1.0) * atoms.center_of_mass());
x.real_value = 0.0;
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
x.real_value += ai->mass * (ai->pos).norm2();
}
}
void colvar::inertia::calc_gradients()
{
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ai->grad = 2.0 * ai->mass * ai->pos;
}
}
void colvar::inertia::apply_force (colvarvalue const &force)
{
if (!atoms.noforce)
atoms.apply_colvar_force (force.real_value);
}
colvar::inertia_z::inertia_z (std::string const &conf)
: inertia (conf)
{
function_type = "inertia_z";
if (get_keyval (conf, "axis", axis, cvm::rvector (0.0, 0.0, 1.0))) {
if (axis.norm2() == 0.0)
cvm::fatal_error ("Axis vector is zero!");
if (axis.norm2() != 1.0) {
axis = axis.unit();
cvm::log ("The normalized axis is: "+cvm::to_str (axis)+".\n");
}
}
x.type (colvarvalue::type_scalar);
}
colvar::inertia_z::inertia_z()
{
function_type = "inertia_z";
x.type (colvarvalue::type_scalar);
}
void colvar::inertia_z::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions();
atoms.apply_translation ((-1.0) * atoms.center_of_mass());
x.real_value = 0.0;
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
cvm::real const iprod = ai->pos * axis;
x.real_value += ai->mass * iprod * iprod;
}
}
void colvar::inertia_z::calc_gradients()
{
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ai->grad = 2.0 * ai->mass * (ai->pos * axis) * axis;
}
}
void colvar::inertia_z::apply_force (colvarvalue const &force)
{
if (!atoms.noforce)
atoms.apply_colvar_force (force.real_value);
}
colvar::rmsd::rmsd (std::string const &conf)
: cvc (conf)
{
b_inverse_gradients = true;
b_Jacobian_derivative = true;
function_type = "rmsd";
x.type (colvarvalue::type_scalar);
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
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;
}
if (atoms.b_center || atoms.b_rotate) {
// request the calculation of the derivatives of the rotation defined by the atom group
atoms.rot.request_group1_gradients (atoms.size());
// request derivatives of optimal rotation wrt 2nd group for Jacobian: this is only
// required for ABF, but we do both groups here for better caching
atoms.rot.request_group2_gradients (atoms.size());
if (atoms.ref_pos_group != NULL) {
cvm::log ("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
"Jacobian derivatives of the RMSD will not be calculated.\n");
b_Jacobian_derivative = false;
}
}
// the following is a simplified version of the corresponding atom group options:
// to define the reference coordinates to compute this variable
if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
cvm::log ("Using reference positions from configuration file to calculate the variable.\n");
if (ref_pos.size() != atoms.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 \"atoms\" ("+
cvm::to_str (atoms.size())+").\n");
}
}
{
std::string ref_pos_file;
if (get_keyval (conf, "refPositionsFile", ref_pos_file, std::string (""))) {
if (ref_pos.size()) {
cvm::fatal_error ("Error: cannot specify \"refPositionsFile\" and "
"\"refPositions\" at the same time.\n");
}
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);
if (found && !ref_pos_col_value)
cvm::fatal_error ("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
} else {
// 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();
+ }
}
void colvar::rmsd::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions();
// rotational fit is done internally
// atoms_cog = atoms.center_of_geometry();
// rot.calc_optimal_rotation (ref_pos, atoms.positions_shifted (-1.0 * atoms_cog));
// cvm::real group_pos_sum2 = 0.0;
// for (size_t i = 0; i < atoms.size(); i++) {
// group_pos_sum2 += (atoms[i].pos - atoms_cog).norm2();
// }
// // value of the RMSD (Coutsias et al)
// cvm::real const MSD = 1.0/(cvm::real (atoms.size())) *
// ( group_pos_sum2 + ref_pos_sum2 - 2.0 * rot.lambda );
// x.real_value = (MSD > 0.0) ? std::sqrt (MSD) : 0.0;
x.real_value = 0.0;
for (size_t ia = 0; ia < atoms.size(); ia++) {
x.real_value += (atoms[ia].pos - ref_pos[ia]).norm2();
}
x.real_value /= cvm::real (atoms.size()); // MSD
x.real_value = (x.real_value > 0.0) ? std::sqrt (x.real_value) : 0.0;
}
void colvar::rmsd::calc_gradients()
{
cvm::real const drmsddx2 = (x.real_value > 0.0) ?
0.5 / (x.real_value * cvm::real (atoms.size())) :
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]));
}
}
void colvar::rmsd::apply_force (colvarvalue const &force)
{
if (!atoms.noforce)
atoms.apply_colvar_force (force.real_value);
}
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;
}
ft.real_value *= atoms.size();
}
void colvar::rmsd::calc_Jacobian_derivative()
{
// divergence of the rotated coordinates (including only derivatives of the rotation matrix)
cvm::real divergence = 0.0;
if (atoms.b_rotate) {
// gradient of the rotation matrix
cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
// 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++) {
// Gradient of optimal quaternion wrt current Cartesian position
cvm::vector1d< cvm::rvector, 4 > &dq = atoms.rot.dQ0_1[ia];
g11 = 2.0 * (atoms.rot.q)[1]*dq[1];
g22 = 2.0 * (atoms.rot.q)[2]*dq[2];
g33 = 2.0 * (atoms.rot.q)[3]*dq[3];
g01 = (atoms.rot.q)[0]*dq[1] + (atoms.rot.q)[1]*dq[0];
g02 = (atoms.rot.q)[0]*dq[2] + (atoms.rot.q)[2]*dq[0];
g03 = (atoms.rot.q)[0]*dq[3] + (atoms.rot.q)[3]*dq[0];
g12 = (atoms.rot.q)[1]*dq[2] + (atoms.rot.q)[2]*dq[1];
g13 = (atoms.rot.q)[1]*dq[3] + (atoms.rot.q)[3]*dq[1];
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);
cvm::atom_pos &x = atoms[ia].pos;
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
divergence += grad_rot_mat[i][j][i] * x[j];
}
}
}
}
jd.real_value = x.real_value > 0.0 ? (3.0 * atoms.size() - 4.0 - divergence) / x.real_value : 0.0;
}
colvar::eigenvector::eigenvector (std::string const &conf)
: cvc (conf)
{
b_inverse_gradients = true;
b_Jacobian_derivative = true;
function_type = "eigenvector";
x.type (colvarvalue::type_scalar);
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
if (atoms.b_rotate) {
cvm::fatal_error ("Error: rotateReference should be disabled:"
"eigenvector component will set it internally.");
}
if (get_keyval (conf, "refPositions", ref_pos, ref_pos)) {
cvm::log ("Using reference positions from input file.\n");
if (ref_pos.size() != atoms.size()) {
cvm::fatal_error ("Error: reference positions do not "
"match the number of requested atoms.\n");
}
}
{
std::string file_name;
if (get_keyval (conf, "refPositionsFile", file_name)) {
std::string file_col;
double file_col_value;
if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) {
// use PDB flags if column is provided
bool found = get_keyval (conf, "refPositionsColValue", file_col_value, 0.0);
if (found && !file_col_value)
cvm::fatal_error ("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
} else {
// if not, use atom indices
atoms.create_sorted_ids();
}
ref_pos.resize (atoms.size());
cvm::load_coords (file_name.c_str(), ref_pos, atoms.sorted_ids, file_col, file_col_value);
}
}
// Set mobile frame of reference for atom group
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)) {
cvm::log ("Using vector components from input file.\n");
if (eigenvec.size() != atoms.size()) {
cvm::fatal_error ("Error: vector components do not "
"match the number of requested atoms.\n");
}
}
{
std::string file_name;
if (get_keyval (conf, "vectorFile", file_name)) {
std::string file_col;
if (!get_keyval (conf, "vectorCol", file_col, std::string (""))) {
cvm::fatal_error ("Error: parameter vectorCol is required if vectorFile is set.\n");
}
double file_col_value;
bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
if (found && !file_col_value)
cvm::fatal_error ("Error: eigenvectorColValue, "
"if provided, must be non-zero.\n");
eigenvec.resize (atoms.size());
cvm::load_coords (file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value);
}
}
if (!ref_pos.size() || !eigenvec.size()) {
cvm::fatal_error ("Error: both reference coordinates"
"and eigenvector must be defined.\n");
}
cvm::rvector center (0.0, 0.0, 0.0);
eigenvec_invnorm2 = 0.0;
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;
eigenvec_invnorm2 += eigenvec[i].norm2();
}
eigenvec_invnorm2 = 1.0 / eigenvec_invnorm2;
// request derivatives of optimal rotation wrt 2nd group
// for Jacobian
atoms.rot.request_group1_gradients(atoms.size());
atoms.rot.request_group2_gradients(atoms.size());
}
void colvar::eigenvector::calc_value()
{
atoms.reset_atoms_data();
atoms.read_positions(); // this will also update atoms.rot
x.real_value = 0.0;
for (size_t i = 0; i < atoms.size(); i++) {
x.real_value += (atoms[i].pos - ref_pos[i]) * eigenvec[i];
}
}
void colvar::eigenvector::calc_gradients()
{
// There are two versions of this code
// The simple version is not formally exact, but its
// results are numerically indistinguishable from the
// exact one. The exact one is more expensive and possibly
// less stable in cases where the optimal rotation
// becomes ill-defined.
// Version A: simple, intuitive, cheap, robust. Wrong.
// Works just fine in practice.
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = eigenvec[ia];
}
/*
// Version B: complex, expensive, fragile. Right.
// gradient of the rotation matrix
cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
cvm::quaternion &quat0 = atoms.rot.q;
// gradients of products of 2 quaternion components
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
// a term that involves the rotation gradients
cvm::rvector rot_grad_term;
cvm::atom_pos x_relative;
cvm::atom_pos atoms_cog = atoms.center_of_geometry();
for (size_t ia = 0; ia < atoms.size(); ia++) {
// Gradient of optimal quaternion wrt current Cartesian position
// WARNING: we want derivatives wrt the FIRST group here (unlike RMSD)
cvm::vector1d< cvm::rvector, 4 > &dq = atoms.rot.dQ0_1[ia];
g11 = 2.0 * quat0[1]*dq[1];
g22 = 2.0 * quat0[2]*dq[2];
g33 = 2.0 * quat0[3]*dq[3];
g01 = quat0[0]*dq[1] + quat0[1]*dq[0];
g02 = quat0[0]*dq[2] + quat0[2]*dq[0];
g03 = quat0[0]*dq[3] + quat0[3]*dq[0];
g12 = quat0[1]*dq[2] + quat0[2]*dq[1];
g13 = quat0[1]*dq[3] + quat0[3]*dq[1];
g23 = quat0[2]*dq[3] + quat0[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);
// this term needs to be rotated back, so we sum it separately
rot_grad_term.reset();
for (size_t ja = 0; ja < atoms.size(); ja++) {
x_relative = atoms[ja].pos - atoms_cog;
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
rot_grad_term += eigenvec[ja][i] * grad_rot_mat[i][j] * x_relative[j];
}
}
}
// Rotate correction term back to reference frame
atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
}
*/
}
void colvar::eigenvector::apply_force (colvarvalue const &force)
{
if (!atoms.noforce)
atoms.apply_colvar_force (force.real_value);
}
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;
}
}
void colvar::eigenvector::calc_Jacobian_derivative()
{
// gradient of the rotation matrix
cvm::matrix2d <cvm::rvector, 3, 3> grad_rot_mat;
cvm::quaternion &quat0 = atoms.rot.q;
// gradients of products of 2 quaternion components
cvm::rvector g11, g22, g33, g01, g02, g03, g12, g13, g23;
cvm::atom_pos x_relative;
cvm::real sum = 0.0;
for (size_t ia = 0; ia < atoms.size(); ia++) {
// Gradient of optimal quaternion wrt current Cartesian position
// trick: d(R^-1)/dx = d(R^t)/dx = (dR/dx)^t
// we can just transpose the derivatives of the direct matrix
cvm::vector1d< cvm::rvector, 4 > &dq_1 = atoms.rot.dQ0_1[ia];
g11 = 2.0 * quat0[1]*dq_1[1];
g22 = 2.0 * quat0[2]*dq_1[2];
g33 = 2.0 * quat0[3]*dq_1[3];
g01 = quat0[0]*dq_1[1] + quat0[1]*dq_1[0];
g02 = quat0[0]*dq_1[2] + quat0[2]*dq_1[0];
g03 = quat0[0]*dq_1[3] + quat0[3]*dq_1[0];
g12 = quat0[1]*dq_1[2] + quat0[2]*dq_1[1];
g13 = quat0[1]*dq_1[3] + quat0[3]*dq_1[1];
g23 = quat0[2]*dq_1[3] + quat0[3]*dq_1[2];
// 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);
for (size_t i = 0; i < 3; i++) {
for (size_t j = 0; j < 3; j++) {
sum += grad_rot_mat[i][j][i] * eigenvec[ia][j];
}
}
}
jd.real_value = sum * std::sqrt (eigenvec_invnorm2);
}
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index dfc0578eb..2375fafa3 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -1,493 +1,493 @@
#ifndef COLVARMODULE_H
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2012-10-03"
+#define COLVARS_VERSION "2012-11-20"
#endif
#ifndef COLVARS_DEBUG
#define COLVARS_DEBUG false
#endif
/// \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
/// shared between all object instances) to be accessed from other
/// objects.
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <cmath>
#include <vector>
#include <list>
class colvarparse;
class colvar;
class colvarbias;
class colvarproxy;
/// \brief Collective variables module (main class)
///
/// 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
class colvarmodule {
private:
/// Impossible to initialize the main object without arguments
colvarmodule();
public:
friend class colvarproxy;
/// Defining an abstract real number allows to switch precision
typedef double real;
/// Residue identifier
typedef int residue_id;
class rvector;
template <class T,
size_t const length> class vector1d;
template <class T,
size_t const outer_length,
size_t const inner_length> class matrix2d;
class quaternion;
class rotation;
/// \brief Atom position (different type name from rvector, to make
/// possible future PBC-transparent implementations)
typedef rvector atom_pos;
/// \brief 3x3 matrix of real numbers
class rmatrix;
// allow these classes to access protected data
class atom;
class atom_group;
friend class atom;
friend class atom_group;
typedef std::vector<atom>::iterator atom_iter;
typedef std::vector<atom>::const_iterator atom_const_iter;
/// Current step number
static size_t it;
/// Starting step number for this run
static size_t it_restart;
/// Return the current step number from the beginning of this run
static inline size_t step_relative()
{
return it - it_restart;
}
/// Return the current step number from the beginning of the whole
/// calculation
static inline size_t step_absolute()
{
return it;
}
/// If true, get it_restart from the state file; if set to false,
/// the MD program is providing it
bool it_restart_from_state_file;
/// \brief Finite difference step size (if there is no dynamics, or
/// if gradients need to be tested independently from the size of
/// dt)
static real debug_gradients_step_size;
/// Prefix for all output files for this run
static std::string output_prefix;
/// Prefix for files from a previous run (including restart/output)
static std::string input_prefix;
/// input restart file name (determined from input_prefix)
static std::string restart_in_name;
/// Array of collective variables
static std::vector<colvar *> colvars;
/* TODO: implement named CVCs
/// Array of named (reusable) collective variable components
static std::vector<cvc *> cvcs;
/// Named cvcs register themselves at initialization time
inline void register_cvc (cvc *p) {
cvcs.push_back(p);
}
*/
/// Array of collective variable biases
static std::vector<colvarbias *> biases;
/// \brief Number of ABF biases initialized (in normal conditions
/// should be 1)
static size_t n_abf_biases;
/// \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
/// number)
static size_t n_harm_biases;
/// \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()
{
return COLVARS_DEBUG;
}
/// \brief Constructor \param config_name Configuration file name
/// \param restart_name (optional) Restart file name
colvarmodule (char const *config_name,
colvarproxy *proxy_in);
/// Destructor
~colvarmodule();
/// Initialize collective variables
void init_colvars (std::string const &conf);
/// Initialize collective variable biases
void init_biases (std::string const &conf);
/// Load new configuration - force constant and/or centers only
void change_configuration(std::string const &name, std::string const &conf);
/// Calculate change in energy from using alternate configuration
real energy_difference(std::string const &name, std::string const &conf);
/// Calculate collective variables and biases
void calc();
/// Read the input restart file
std::istream & read_restart (std::istream &is);
/// Write the output restart file
std::ostream & write_restart (std::ostream &os);
/// Write all output files (called by the proxy)
void write_output_files();
/// \brief Call colvarproxy::backup_file()
static void backup_file (char const *filename);
/// Perform analysis
void analyze();
/// \brief Read a collective variable trajectory (post-processing
/// only, not called at runtime)
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<typename T> 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<typename T> static std::string to_str (std::vector<T> const &x,
size_t const &width = 0,
size_t const &prec = 0);
/// Reduce the number of characters in a string
static inline std::string wrap_string (std::string const &s,
size_t const &nchars)
{
if (!s.size())
return std::string (nchars, ' ');
else
return ( (s.size() <= size_t (nchars)) ?
(s+std::string (nchars-s.size(), ' ')) :
(std::string (s, 0, nchars)) );
}
/// Number of characters to represent a time step
static size_t const it_width;
/// Number of digits to represent a collective variables value(s)
static size_t const cv_prec;
/// Number of characters to represent a collective variables value(s)
static size_t const cv_width;
/// Number of digits to represent the collective variables energy
static size_t const en_prec;
/// Number of characters to represent the collective variables energy
static size_t const en_width;
/// Line separator in the log output
static std::string const line_marker;
// proxy functions
/// \brief Value of the unit for atomic coordinates with respect to
/// angstroms (used by some variables for hard-coded default values)
static real unit_angstrom();
/// \brief Boltmann constant
static real boltzmann();
/// \brief Temperature of the simulation (K)
static real temperature();
/// \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);
/// Print a message to the main log and exit with error code
static void fatal_error (std::string const &message);
/// Print a message to the main log and exit normally
static void exit (std::string const &message);
/// \brief Get the distance between two atomic positions with pbcs handled
/// correctly
static rvector position_distance (atom_pos const &pos1,
atom_pos const &pos2);
/// \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)
static real position_dist2 (atom_pos const &pos1,
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 (optional) The reference position
static void select_closest_image (atom_pos &pos,
atom_pos const &ref_pos);
/// \brief Perform select_closest_image() on a set of atomic positions
///
/// After that, distance vectors can then be calculated directly,
/// without using position_distance()
static void select_closest_images (std::vector<atom_pos> &pos,
atom_pos const &ref_pos);
/// \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
/// field to determine which are the atoms to be set
static void load_atoms (char const *filename,
std::vector<atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from a file
/// (usually a PDB); the number of atoms in "filename" must match
/// the number of elements in "pos"
static void load_coords (char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// Frequency for collective variables trajectory output
static size_t cv_traj_freq;
/// \brief True if only analysis is performed and not a run
static bool b_analysis;
/// Frequency for saving output restarts
static size_t restart_out_freq;
/// Output restart file name
std::string restart_out_name;
/// Pseudo-random number with Gaussian distribution
static real rand_gaussian (void);
protected:
/// Configuration file
std::ifstream config_s;
/// Configuration file parser object
colvarparse *parse;
/// Name of the trajectory file
std::string cv_traj_name;
/// Collective variables output trajectory file
std::ofstream cv_traj_os;
/// Appending to the existing trajectory file?
bool cv_traj_append;
/// Output restart file
std::ofstream restart_out_os;
/// \brief Pointer to the proxy object, used to retrieve atomic data
/// from the hosting program; it is static in order to be accessible
/// from static functions in the colvarmodule class
static colvarproxy *proxy;
/// \brief Counter for the current depth in the object hierarchy (useg e.g. in outpu
static size_t depth;
public:
/// Increase the depth (number of indentations in the output)
static void increase_depth();
/// Decrease the depth (number of indentations in the output)
static void decrease_depth();
};
/// Shorthand for the frequently used type prefix
typedef colvarmodule cvm;
#include "colvartypes.h"
std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
std::istream & operator >> (std::istream &is, cvm::rvector &v);
template<typename T> std::string cvm::to_str (T const &x,
size_t const &width,
size_t const &prec) {
std::ostringstream os;
if (width) os.width (width);
if (prec) {
os.setf (std::ios::scientific, std::ios::floatfield);
os.precision (prec);
}
os << x;
return os.str();
}
template<typename T> std::string cvm::to_str (std::vector<T> const &x,
size_t const &width,
size_t const &prec) {
if (!x.size()) return std::string ("");
std::ostringstream os;
if (prec) {
os.setf (std::ios::scientific, std::ios::floatfield);
}
os << "{ ";
if (width) os.width (width);
if (prec) os.precision (prec);
os << x[0];
for (size_t i = 1; i < x.size(); i++) {
os << ", ";
if (width) os.width (width);
if (prec) os.precision (prec);
os << x[i];
}
os << " }";
return os.str();
}
#include "colvarproxy.h"
inline cvm::real cvm::unit_angstrom()
{
return proxy->unit_angstrom();
}
inline cvm::real cvm::boltzmann()
{
return proxy->boltzmann();
}
inline cvm::real cvm::temperature()
{
return proxy->temperature();
}
inline cvm::real cvm::dt()
{
return proxy->dt();
}
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)
{
proxy->select_closest_image (pos, ref_pos);
}
inline void cvm::select_closest_images (std::vector<atom_pos> &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_images (pos, ref_pos);
}
inline cvm::rvector cvm::position_distance (atom_pos const &pos1,
atom_pos const &pos2)
{
return proxy->position_distance (pos1, pos2);
}
inline cvm::real cvm::position_dist2 (cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return proxy->position_dist2 (pos1, pos2);
}
inline void cvm::load_atoms (char const *file_name,
std::vector<cvm::atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value)
{
proxy->load_atoms (file_name, atoms, pdb_field, pdb_field_value);
}
inline void cvm::load_coords (char const *file_name,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value)
{
proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value);
}
inline void cvm::backup_file (char const *filename)
{
proxy->backup_file (filename);
}
inline cvm::real cvm::rand_gaussian (void)
{
return proxy->rand_gaussian();
}
#endif
// Emacs
// Local Variables:
// mode: C++
// End: