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@9395 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2013-02-01 22:39:47 +00:00
parent f0c4f5680b
commit df9c70770d
10 changed files with 626 additions and 373 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

61
lib/colvars/colvar.cpp
View File

@ -292,6 +292,11 @@ colvar::colvar (std::string const &conf)
cvm::fatal_error ("Error: trying to expand boundaries that already "
"cover a whole period of a periodic colvar.\n");
}
if (expand_boundaries && hard_lower_boundary && hard_upper_boundary) {
cvm::fatal_error ("Error: inconsistent configuration "
"(trying to expand boundaries with both "
"hardLowerBoundary and hardUpperBoundary enabled).\n");
}
{
bool b_extended_lagrangian;
@ -676,12 +681,39 @@ void colvar::calc()
if (cvm::debug())
cvm::log ("Calculating colvar \""+this->name+"\".\n");
// prepare atom groups for calculation
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]);
atoms.reset_atoms_data();
atoms.read_positions();
if (atoms.b_center || atoms.b_rotate) {
atoms.calc_apply_roto_translation();
}
// each atom group will take care of its own ref_pos_group, if defined
}
}
if (tasks[task_output_velocity]) {
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvcs[i]->atom_groups[ig]->read_velocities();
}
}
}
if (tasks[task_system_force]) {
for (size_t i = 0; i < cvcs.size(); i++) {
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvcs[i]->atom_groups[ig]->read_system_forces();
}
}
}
// calculate the value of the colvar
x.reset();
if (x.type() == colvarvalue::type_scalar) {
// polynomial combination allowed
// scalar variable, polynomial combination allowed
for (size_t i = 0; i < cvcs.size(); i++) {
cvm::increase_depth();
(cvcs[i])->calc_value();
@ -697,6 +729,7 @@ void colvar::calc()
(cvcs[i])->value().real_value );
}
} else {
// only linear combination allowed
for (size_t i = 0; i < cvcs.size(); i++) {
cvm::increase_depth();
@ -716,12 +749,29 @@ void colvar::calc()
cvm::to_str (x, cvm::cv_width, cvm::cv_prec)+".\n");
if (tasks[task_gradients]) {
// calculate the gradients
if (cvm::debug())
cvm::log ("Calculating gradients of colvar \""+this->name+"\".\n");
for (size_t i = 0; i < cvcs.size(); i++) {
// calculate the gradients of each component
cvm::increase_depth();
(cvcs[i])->calc_gradients();
// if requested, propagate (via chain rule) the gradients above
// to the atoms used to define the roto-translation
for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
if (cvcs[i]->atom_groups[ig]->b_fit_gradients)
cvcs[i]->atom_groups[ig]->calc_fit_gradients();
}
cvm::decrease_depth();
}
if (cvm::debug())
cvm::log ("Done calculating gradients of colvar \""+this->name+"\".\n");
if (tasks[task_collect_gradients]) {
// Collect the atomic gradients inside colvar object
for (int a = 0; a < atomic_gradients.size(); a++) {
@ -760,6 +810,10 @@ void colvar::calc()
}
if (tasks[task_system_force]) {
if (cvm::debug())
cvm::log ("Calculating system force of colvar \""+this->name+"\".\n");
ft.reset();
if(!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) {
@ -778,6 +832,9 @@ void colvar::calc()
// correction internally: biases such as colvarbias_abf will handle it
ft += fj;
}
if (cvm::debug())
cvm::log ("Done calculating system force of colvar \""+this->name+"\".\n");
}
if (tasks[task_fdiff_velocity]) {

272
lib/colvars/colvaratoms.cpp
View File

@ -3,29 +3,76 @@
#include "colvaratoms.h"
// atom member functions depend tightly on the MD interface, and are
// member functions of the "atom" class depend tightly on the MD interface, and are
// thus defined in colvarproxy_xxx.cpp
// in this file only atom_group functions are defined
// atom_group member functions
// Note: "conf" is the configuration of the cvc who is using this atom group;
// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...)
cvm::atom_group::atom_group (std::string const &conf,
char const *key,
atom_group *ref_pos_group_in)
char const *key)
: 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
ref_pos_group (NULL),
b_fit_gradients (false),
noforce (false)
{
cvm::log ("Defining atom group \""+
std::string (key)+"\".\n");
parse (conf, key, ref_pos_group_in);
// real work is done by parse
parse (conf, key);
}
cvm::atom_group::atom_group (std::vector<cvm::atom> const &atoms)
: b_dummy (false), b_center (false), b_rotate (false),
ref_pos_group (NULL), b_fit_gradients (false),
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), b_fit_gradients (false),
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::parse (std::string const &conf,
char const *key,
atom_group *ref_pos_group_in)
char const *key)
{
std::string group_conf;
@ -226,6 +273,8 @@ void cvm::atom_group::parse (std::string const &conf,
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);
b_fit_gradients = get_keyval (group_conf, "fitGradients", b_fit_gradients, 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;
@ -236,9 +285,9 @@ void cvm::atom_group::parse (std::string const &conf,
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")) {
// instead of this group, define another group (refPositionsGroup)
// to be the one used to fit the coordinates
if (ref_pos_group) {
cvm::fatal_error ("Error: the atom group \""+
std::string (key)+"\" has already a reference group "
@ -252,15 +301,7 @@ void cvm::atom_group::parse (std::string const &conf,
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");
}
}
get_keyval (group_conf, "refPositions", ref_pos, ref_pos, mode);
std::string ref_pos_file;
if (get_keyval (group_conf, "refPositionsFile", ref_pos_file, std::string (""), mode)) {
@ -310,16 +351,22 @@ void cvm::atom_group::parse (std::string const &conf,
#endif
}
if (b_fit_gradients) {
group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::atom_pos (0.0, 0.0, 0.0));
rot.request_group1_gradients (group_for_fit->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");
"\" will be aligned to a fixed orientation given by the reference positions provided. "
"If the internal structure of the group changes too much (i.e. its RMSD is comparable "
"to its radius of gyration), the optimal rotation and its gradients may become discontinuous. "
"If that happens, use refPositionsGroup (or a different definition for it if already defined) "
"to align the coordinates.\n");
// initialize rot member data
rot.request_group1_gradients (this->size());
}
}
if (cvm::debug())
@ -336,50 +383,6 @@ void cvm::atom_group::parse (std::string const &conf,
}
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
@ -425,11 +428,6 @@ 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();
@ -437,13 +435,14 @@ void cvm::atom_group::read_positions()
if (ref_pos_group)
ref_pos_group->read_positions();
}
void cvm::atom_group::calc_apply_roto_translation()
{
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
// center on the origin first
cvm::atom_pos const cog = fit_group->center_of_geometry();
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
@ -453,9 +452,7 @@ void cvm::atom_group::read_positions()
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
// true, around the origin otherwise)
rot.calc_optimal_rotation (fit_group->positions(), ref_pos);
for (cvm::atom_iter ai = this->begin();
@ -465,7 +462,7 @@ void cvm::atom_group::read_positions()
}
if (b_center) {
// use the center of geometry of ref_pos
// align with the center of geometry of ref_pos
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
ai->pos += ref_pos_cog;
@ -473,7 +470,6 @@ void cvm::atom_group::read_positions()
}
}
void cvm::atom_group::apply_translation (cvm::rvector const &t)
{
if (b_dummy) return;
@ -484,7 +480,6 @@ void cvm::atom_group::apply_translation (cvm::rvector const &t)
}
}
void cvm::atom_group::apply_rotation (cvm::rotation const &rot)
{
if (b_dummy) return;
@ -538,24 +533,6 @@ void cvm::atom_group::read_system_forces()
}
}
/* 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)
@ -570,24 +547,6 @@ cvm::atom_pos cvm::atom_group::center_of_geometry() const
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)
@ -614,6 +573,57 @@ void cvm::atom_group::set_weighted_gradient (cvm::rvector const &grad)
}
void cvm::atom_group::calc_fit_gradients()
{
if (b_dummy) return;
if (cvm::debug())
cvm::log ("Calculating fit gradients.\n");
atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
group_for_fit->fit_gradients.assign (group_for_fit->size(), cvm::rvector (0.0, 0.0, 0.0));
if (b_center) {
// add the center of geometry contribution to the gradients
for (size_t i = 0; i < this->size(); i++) {
// need to bring the gradients in original frame first
cvm::rvector const atom_grad = b_rotate ?
(rot.inverse()).rotate ((*this)[i].grad) :
(*this)[i].grad;
for (size_t j = 0; j < group_for_fit->size(); j++) {
group_for_fit->fit_gradients[j] +=
(-1.0)/(cvm::real (group_for_fit->size())) *
atom_grad;
}
}
}
if (b_rotate) {
// add the rotation matrix contribution to the gradients
cvm::rotation const rot_inv = rot.inverse();
cvm::atom_pos const cog = this->center_of_geometry();
for (size_t i = 0; i < this->size(); i++) {
cvm::atom_pos const pos_orig = rot_inv.rotate ((b_center ? ((*this)[i].pos - cog) : ((*this)[i].pos)));
for (size_t j = 0; j < group_for_fit->size(); j++) {
// calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
cvm::quaternion const dxdq =
rot.q.position_derivative_inner (pos_orig, (*this)[i].grad);
// multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
for (size_t iq = 0; iq < 4; iq++) {
group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
}
}
}
}
if (cvm::debug())
cvm::log ("Done calculating fit gradients.\n");
}
std::vector<cvm::atom_pos> cvm::atom_group::positions() const
{
if (b_dummy)
@ -699,7 +709,7 @@ void cvm::atom_group::apply_colvar_force (cvm::real const &force)
if (b_rotate) {
// get the forces back from the rotated frame
// rotate forces back to the original frame
cvm::rotation const rot_inv = rot.inverse();
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
@ -708,13 +718,30 @@ void cvm::atom_group::apply_colvar_force (cvm::real const &force)
} 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);
}
}
if (b_fit_gradients) {
atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
// add the contribution from the roto-translational fit to the gradients
if (b_rotate) {
// rotate forces back to the original frame
cvm::rotation const rot_inv = rot.inverse();
for (size_t j = 0; j < group_for_fit->size(); j++) {
(*group_for_fit)[j].apply_force (rot_inv.rotate (force * fit_gradients[j]));
}
} else {
for (size_t j = 0; j < group_for_fit->size(); j++) {
(*group_for_fit)[j].apply_force (force * fit_gradients[j]);
}
}
}
}
@ -778,4 +805,3 @@ void cvm::atom_group::apply_forces (std::vector<cvm::rvector> const &forces)
}
}

65
lib/colvars/colvaratoms.h
View File

@ -140,39 +140,42 @@ public:
void create_sorted_ids (void);
/// \brief Before calculating colvars, move the group to overlap the
/// center of mass of reference coordinates
/// \brief When updating atomic coordinates, translate them to align with the
/// center of mass of the 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.
/// \brief When updating atom coordinates (and after
/// centering them if b_center is set), rotate the group to
/// align with the reference coordinates.
///
/// 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
/// The rotation calculated automatically if b_rotate is defined
cvm::rotation rot;
/// \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
/// \brief Indicates that the user has explicitly set centerReference or
/// rotateReference, and the corresponding reference:
/// cvc's (eg rmsd, eigenvector) will not override the user's choice
bool b_user_defined_fit;
/// \brief use these reference coordinates, for b_center, b_rotate, or to compute
/// certain colvar components (orientation, rmsd, etc)
/// \brief Whether or not the derivatives of the roto-translation
/// should be included when calculating the colvar's gradients (default: no)
bool b_fit_gradients;
/// \brief use reference coordinates for b_center or b_rotate
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)
/// \brief If b_center or b_rotate is true, use this group to
/// define the transformation (default: this group itself)
atom_group *ref_pos_group;
/// Total mass of the atom group
cvm::real total_mass;
@ -186,14 +189,12 @@ public:
/// 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);
char const *key);
/// \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);
char const *key);
/// \brief Initialize the group after a temporary vector of atoms
atom_group (std::vector<cvm::atom> const &atoms);
@ -208,10 +209,12 @@ public:
~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
/// true, calc_apply_roto_translation() will be called too
void read_positions();
/// \brief (Re)calculate the optimal roto-translation
void calc_apply_roto_translation();
/// \brief Save center of geometry fo ref positions,
/// then subtract it
void center_ref_pos();
@ -238,6 +241,8 @@ public:
{
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++)
ai->reset_data();
if (ref_pos_group)
ref_pos_group->reset_atoms_data();
}
/// \brief Return a copy of the current atom positions
@ -246,21 +251,15 @@ public:
/// \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
/// \brief Atom positions at the previous step
std::vector<cvm::atom_pos> old_pos;
@ -279,6 +278,12 @@ public:
/// \link center_of_mass() \endlink)
void set_weighted_gradient (cvm::rvector const &grad);
/// \brief Calculate the derivatives of the fitting transformation
void calc_fit_gradients();
/// \brief Derivatives of the fitting transformation
std::vector<cvm::atom_pos> fit_gradients;
/// \brief Used by a (scalar) colvar to apply its force on its \link
/// atom_group \endlink members
///
@ -290,7 +295,7 @@ public:
/// 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.
/// gradients are temporarily rotated 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

94
lib/colvars/colvarbias_meta.cpp
View File

@ -49,6 +49,15 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
get_keyval (conf, "hillWidth", hill_width, std::sqrt (2.0 * PI) / 2.0);
{
bool b_replicas = false;
get_keyval (conf, "multipleReplicas", b_replicas, false);
if (b_replicas)
comm = multiple_replicas;
else
comm = single_replica;
}
get_keyval (conf, "useGrids", use_grids, true);
if (use_grids) {
@ -84,15 +93,6 @@ colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
dump_replica_fes = false;
}
{
bool b_replicas = false;
get_keyval (conf, "multipleReplicas", b_replicas, false);
if (b_replicas)
comm = multiple_replicas;
else
comm = single_replica;
}
if (comm != single_replica) {
if (expand_grids)
@ -373,34 +373,36 @@ cvm::real colvarbias_meta::update()
int &new_size = new_sizes[i];
bool changed_lb = false, changed_ub = false;
if (curr_bin[i] < min_buffer) {
int const extra_points = (min_buffer - curr_bin[i]);
new_lb -= extra_points * colvars[i]->width;
new_size += extra_points;
// changed offset in this direction => the pointer needs to
// be changed, too
curr_bin[i] += extra_points;
if (!colvars[i]->hard_lower_boundary)
if (curr_bin[i] < min_buffer) {
int const extra_points = (min_buffer - curr_bin[i]);
new_lb -= extra_points * colvars[i]->width;
new_size += extra_points;
// changed offset in this direction => the pointer needs to
// be changed, too
curr_bin[i] += extra_points;
changed_lb = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new lower boundary for colvar \""+
colvars[i]->name+"\", at "+
cvm::to_str (new_lower_boundaries[i])+".\n");
}
changed_lb = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new lower boundary for colvar \""+
colvars[i]->name+"\", at "+
cvm::to_str (new_lower_boundaries[i])+".\n");
}
if (curr_bin[i] > new_size - min_buffer - 1) {
int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
new_ub += extra_points * colvars[i]->width;
new_size += extra_points;
if (!colvars[i]->hard_upper_boundary)
if (curr_bin[i] > new_size - min_buffer - 1) {
int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
new_ub += extra_points * colvars[i]->width;
new_size += extra_points;
changed_ub = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new upper boundary for colvar \""+
colvars[i]->name+"\", at "+
cvm::to_str (new_upper_boundaries[i])+".\n");
}
changed_ub = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new upper boundary for colvar \""+
colvars[i]->name+"\", at "+
cvm::to_str (new_upper_boundaries[i])+".\n");
}
if (changed_lb || changed_ub)
changed_grids = true;
@ -1521,7 +1523,8 @@ std::ostream & colvarbias_meta::write_restart (std::ostream& os)
void colvarbias_meta::write_pmf()
{
// allocate a new grid to store the pmf
colvar_grid_scalar *pmf = new colvar_grid_scalar (colvars);
colvar_grid_scalar *pmf = new colvar_grid_scalar (*hills_energy);
pmf->create();
std::string fes_file_name_prefix (cvm::output_prefix);
@ -1544,17 +1547,18 @@ void colvarbias_meta::write_pmf()
cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
pmf->multiply_constant (well_temper_scale);
}
std::string const fes_file_name (fes_file_name_prefix +
((comm != single_replica) ? ".partial" : "") +
(dump_fes_save ?
"."+cvm::to_str (cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file (fes_file_name.c_str());
std::ofstream fes_os (fes_file_name.c_str());
pmf->write_multicol (fes_os);
fes_os.close();
{
std::string const fes_file_name (fes_file_name_prefix +
((comm != single_replica) ? ".partial" : "") +
(dump_fes_save ?
"."+cvm::to_str (cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file (fes_file_name.c_str());
std::ofstream fes_os (fes_file_name.c_str());
pmf->write_multicol (fes_os);
fes_os.close();
}
}
if (comm != single_replica) {
// output the combined PMF from all replicas
pmf->reset();

73
lib/colvars/colvarcomp.cpp
View File

@ -71,3 +71,76 @@ void colvar::cvc::calc_Jacobian_derivative()
cvm::fatal_error ("Error: calculation of inverse gradients is not implemented "
"for colvar components of type \""+function_type+"\".\n");
}
void colvar::cvc::debug_gradients (cvm::atom_group &group)
{
// this function should work for any scalar variable:
// the only difference will be the name of the atom group (here, "group")
// collect into a vector for convenience
std::vector<cvm::rvector> gradients (group.size());
for (size_t i = 0; i < group.size(); i++) {
gradients[i] = group[i].grad;
}
cvm::rotation const rot_0 = group.rot;
cvm::rotation const rot_inv = group.rot.inverse();
cvm::real const x_0 = x.real_value;
cvm::log ("gradients = "+cvm::to_str (gradients)+"\n");
// it only makes sense to debug the fit gradients
// when the fitting group is the same as this group
if (group.b_fit_gradients && (group.ref_pos_group == NULL)) {
group.calc_fit_gradients();
if (group.b_rotate) {
// fit_gradients are in the original frame, we should print them in the rotated frame
for (size_t j = 0; j < group.fit_gradients.size(); j++) {
group.fit_gradients[j] = rot_0.rotate (group.fit_gradients[j]);
}
}
cvm::log ("fit_gradients = "+cvm::to_str (group.fit_gradients)+"\n");
if (group.b_rotate) {
for (size_t j = 0; j < group.fit_gradients.size(); j++) {
group.fit_gradients[j] = rot_inv.rotate (group.fit_gradients[j]);
}
}
}
for (size_t ia = 0; ia < group.size(); ia++) {
// tests are best conducted in the unrotated (simulation) frame
cvm::rvector const atom_grad = group.b_rotate ?
rot_inv.rotate (group[ia].grad) :
group[ia].grad;
for (size_t id = 0; id < 3; id++) {
// (re)read original positions
group.read_positions();
// change one coordinate
group[ia].pos[id] += cvm::debug_gradients_step_size;
// (re)do the fit (if defined)
if (group.b_center || group.b_rotate) {
group.calc_apply_roto_translation();
}
calc_value();
cvm::real const x_1 = x.real_value;
cvm::log ("Atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n");
cvm::log ("dx(actual) = "+cvm::to_str (x_1 - x_0,
21, 14)+"\n");
cvm::real const dx_pred = (group.b_fit_gradients && (group.ref_pos_group == NULL)) ?
(cvm::debug_gradients_step_size * (atom_grad[id] + group.fit_gradients[ia][id])) :
(cvm::debug_gradients_step_size * atom_grad[id]);
cvm::log ("dx(interp) = "+cvm::to_str (dx_pred,
21, 14)+"\n");
cvm::log ("|dx(actual) - dx(interp)|/|dx(actual)| = "+
cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) /
std::fabs (x_1 - x_0),
12, 5)+
".\n");
}
}
}

25
lib/colvars/colvarcomp.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARDEF_H
#define COLVARDEF_H
#ifndef COLVARCOMP_H
#define COLVARCOMP_H
#include <fstream>
#include <cmath>
@ -106,10 +106,6 @@ public:
/// Destructor
virtual ~cvc();
/// \brief If true, calc_gradients() will calculate
/// finite-difference gradients alongside the analytical ones and
/// report their differences
bool b_debug_gradients;
/// \brief If this flag is false (default), inverse gradients
/// (derivatives of atom coordinates with respect to x) are
@ -130,6 +126,12 @@ public:
/// order to apply forces
virtual void calc_gradients() = 0;
/// \brief If true, calc_gradients() will call debug_gradients() for every group needed
bool b_debug_gradients;
/// \brief Calculate finite-difference gradients alongside the analytical ones, for each Cartesian component
virtual void debug_gradients (cvm::atom_group &group);
/// \brief Calculate the total force from the system using the
/// inverse atomic gradients
virtual void calc_force_invgrads();
@ -527,11 +529,8 @@ public:
/// \brief Colvar component: moment of inertia of an atom group
/// (colvarvalue::type_scalar type, range [0:*))
class colvar::inertia
: public colvar::cvc
: public colvar::gyration
{
protected:
/// Atoms involved
cvm::atom_group atoms;
public:
/// Constructor
inertia (std::string const &conf);
@ -591,7 +590,10 @@ protected:
/// Reference coordinates
std::vector<cvm::atom_pos> ref_pos;
/// Eigenvector (of a normal or essential mode)
/// Geometric center of the reference coordinates
cvm::rvector ref_pos_center;
/// Eigenvector (of a normal or essential mode): will always have zero center
std::vector<cvm::rvector> eigenvec;
/// Inverse square norm of the eigenvector
@ -1138,7 +1140,6 @@ public:
};
// metrics functions for cvc implementations with a periodicity
inline cvm::real colvar::dihedral::dist2 (colvarvalue const &x1,

402
lib/colvars/colvarcomp_distances.cpp
View File

@ -8,9 +8,6 @@
/// \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)
@ -48,12 +45,6 @@ colvar::distance::distance()
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 {
@ -113,12 +104,6 @@ colvar::distance_vec::distance_vec()
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 {
@ -205,12 +190,6 @@ colvar::distance_z::distance_z()
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();
@ -219,8 +198,6 @@ void colvar::distance_z::calc_value()
main.center_of_mass());
}
} else {
ref2.reset_atoms_data();
ref2.read_positions();
if (b_no_PBC) {
dist_v = main.center_of_mass() -
@ -310,12 +287,6 @@ colvar::distance_xy::distance_xy()
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 {
@ -323,9 +294,6 @@ void colvar::distance_xy::calc_value()
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 {
@ -415,12 +383,6 @@ colvar::distance_dir::distance_dir()
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 {
@ -490,12 +452,6 @@ colvar::distance_inv::distance_inv()
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++) {
@ -561,6 +517,17 @@ colvar::gyration::gyration (std::string const &conf)
b_Jacobian_derivative = true;
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
if (atoms.b_user_defined_fit) {
cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
} else {
// default: fit everything
cvm::log ("No explicit fitting parameters: enabling centerReference by default to align with the origin\n");
cvm::log (" (rotateReference will not affect results for gyration and inertia, but it will for inertiaZ).\n");
atoms.b_center = true;
atoms.ref_pos.assign (1, cvm::rvector (0.0, 0.0, 0.0));
}
x.type (colvarvalue::type_scalar);
}
@ -576,10 +543,6 @@ colvar::gyration::gyration()
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();
@ -594,6 +557,11 @@ void colvar::gyration::calc_gradients()
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ai->grad = drdx * ai->pos;
}
if (b_debug_gradients) {
cvm::log ("Debugging gradients:\n");
debug_gradients (atoms);
}
}
@ -625,11 +593,11 @@ void colvar::gyration::apply_force (colvarvalue const &force)
colvar::inertia::inertia (std::string const &conf)
: cvc (conf)
: gyration (conf)
{
function_type = "inertia";
parse_group (conf, "atoms", atoms);
atom_groups.push_back (&atoms);
b_inverse_gradients = false;
b_Jacobian_derivative = false;
x.type (colvarvalue::type_scalar);
}
@ -643,10 +611,6 @@ colvar::inertia::inertia()
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();
@ -659,6 +623,11 @@ void colvar::inertia::calc_gradients()
for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
ai->grad = 2.0 * ai->mass * ai->pos;
}
if (b_debug_gradients) {
cvm::log ("Debugging gradients:\n");
debug_gradients (atoms);
}
}
@ -694,10 +663,6 @@ colvar::inertia_z::inertia_z()
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;
@ -711,6 +676,11 @@ 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;
}
if (b_debug_gradients) {
cvm::log ("Debugging gradients:\n");
debug_gradients (atoms);
}
}
@ -736,32 +706,15 @@ colvar::rmsd::rmsd (std::string const &conf)
if (atoms.b_dummy)
cvm::fatal_error ("Error: \"atoms\" cannot be a dummy atom.");
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.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;
}
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
// the following is a simplified version of the corresponding atom group options;
// we need this because the reference coordinates defined inside the atom group
// may be used only for fitting, and even more so if ref_pos_group is used
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()) {
@ -800,41 +753,37 @@ colvar::rmsd::rmsd (std::string const &conf)
}
}
if (!atoms.b_user_defined_fit) {
// Default case: reference positions for calculating the rmsd are also those used
if (atoms.b_user_defined_fit) {
cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".");
} else {
// Default: fit everything
cvm::log ("No explicit fitting parameters: enabling both centerReference and rotateReference to compute standard minimum RMSD.");
atoms.b_center = true;
atoms.b_rotate = true;
// default case: reference positions for calculating the rmsd are also those used
// for fitting
atoms.ref_pos = ref_pos;
atoms.center_ref_pos();
// 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 reference coordinates for Jacobian:
// this is only required for ABF, but we do both groups here for better caching
atoms.rot.request_group2_gradients (atoms.size());
}
}
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;
// rotational-translational fit is handled by the atom group
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;
x.real_value = std::sqrt (x.real_value);
}
@ -847,6 +796,11 @@ void colvar::rmsd::calc_gradients()
for (size_t ia = 0; ia < atoms.size(); ia++) {
atoms[ia].grad = (drmsddx2 * 2.0 * (atoms[ia].pos - ref_pos[ia]));
}
if (b_debug_gradients) {
cvm::log ("Debugging gradients:\n");
debug_gradients (atoms);
}
}
@ -908,11 +862,14 @@ void colvar::rmsd::calc_Jacobian_derivative()
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;
cvm::atom_pos &y = ref_pos[ia];
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];
for (size_t alpha = 0; alpha < 3; alpha++) {
for (size_t beta = 0; beta < 3; beta++) {
divergence += grad_rot_mat[beta][alpha][alpha] * y[beta];
// Note: equation was derived for inverse rotation (see colvars paper)
// so here the matrix is transposed
// (eq would give divergence += grad_rot_mat[alpha][beta][alpha] * y[beta];)
}
}
}
@ -935,23 +892,24 @@ colvar::eigenvector::eigenvector (std::string const &conf)
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");
}
}
{
bool const b_inline = get_keyval (conf, "refPositions", ref_pos, ref_pos);
if (b_inline) {
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)) {
if (b_inline)
cvm::fatal_error ("Error: refPositions and refPositionsFile cannot be specified at the same time.\n");
std::string file_col;
double file_col_value;
if (get_keyval (conf, "refPositionsCol", file_col, std::string (""))) {
@ -970,35 +928,59 @@ colvar::eigenvector::eigenvector (std::string const &conf)
}
}
// 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();
// save for later the geometric center of the provided positions (may not be the origin)
cvm::rvector ref_pos_center (0.0, 0.0, 0.0);
for (size_t i = 0; i < atoms.size(); i++) {
ref_pos_center += ref_pos[i];
}
ref_pos_center *= 1.0 / atoms.size();
// 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");
}
if (atoms.b_user_defined_fit) {
cvm::log ("WARNING: explicit fitting parameters were provided for atom group \"atoms\".\n");
} else {
// default: fit everything
cvm::log ("Enabling centerReference and rotateReference (minimize RMSD before calculating the projection).\n");
atoms.b_center = true;
atoms.b_rotate = true;
atoms.ref_pos = ref_pos;
atoms.center_ref_pos();
// 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 reference coordinates for Jacobian:
// this is only required for ABF, but we do both groups here for better caching
atoms.rot.request_group2_gradients (atoms.size());
}
{
bool const b_inline = get_keyval (conf, "vector", eigenvec, eigenvec);
// now load the eigenvector
if (b_inline) {
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");
}
if (b_inline)
cvm::fatal_error ("Error: vector and vectorFile cannot be specified at the same time.\n");
std::string file_col;
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");
if (get_keyval (conf, "vectorCol", file_col, std::string (""))) {
// use PDB flags if column is provided
bool found = get_keyval (conf, "vectorColValue", file_col_value, 0.0);
if (found && !file_col_value)
cvm::fatal_error ("Error: vectorColValue, "
"if provided, must be non-zero.\n");
} else {
// if not, use atom indices
atoms.create_sorted_ids();
}
eigenvec.resize (atoms.size());
cvm::load_coords (file_name.c_str(), eigenvec, atoms.sorted_ids, file_col, file_col_value);
@ -1010,34 +992,64 @@ colvar::eigenvector::eigenvector (std::string const &conf)
"and eigenvector must be defined.\n");
}
cvm::rvector center (0.0, 0.0, 0.0);
eigenvec_invnorm2 = 0.0;
cvm::rvector eig_center (0.0, 0.0, 0.0);
for (size_t i = 0; i < atoms.size(); i++) {
center += eigenvec[i];
eig_center += eigenvec[i];
}
center *= 1.0 / atoms.size();
eig_center *= 1.0 / atoms.size();
cvm::log ("Geometric center of the provided vector: "+cvm::to_str (eig_center)+"\n");
cvm::log ("Subtracting center of eigenvector components: " + cvm::to_str (center) + "\n");
bool b_vector_difference;
get_keyval (conf, "vectorDifference", b_vector_difference, false);
if (b_vector_difference) {
cvm::log ("Subtracting the reference positions from the provided vector.\n");
if (atoms.b_center) {
cvm::log ("centerReference is on: recentering the provided vector to the reference positions.\n");
// both sets should be centered on the origin for fitting
for (size_t i = 0; i < atoms.size(); i++) {
eigenvec[i] -= eig_center;
ref_pos[i] -= ref_pos_center;
}
}
if (atoms.b_rotate) {
cvm::log ("rotateReference is on: aligning the provided vector to the reference positions.\n");
atoms.rot.calc_optimal_rotation (eigenvec, ref_pos);
for (size_t i = 0; i < atoms.size(); i++) {
eigenvec[i] = atoms.rot.rotate (eigenvec[i]);
eigenvec[i] -= ref_pos[i];
}
}
if (atoms.b_center) {
// bring back the ref positions to where they were
for (size_t i = 0; i < atoms.size(); i++) {
ref_pos[i] -= ref_pos_center;
}
}
} else {
cvm::log ("Centering the provided vector to zero.\n");
for (size_t i = 0; i < atoms.size(); i++) {
eigenvec[i] -= eig_center;
}
}
cvm::log ("The first three components (v1x, v1y, v1z) of the resulting vector are: "+cvm::to_str (eigenvec[0])+".\n");
// for inverse gradients
eigenvec_invnorm2 = 0.0;
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];
@ -1060,7 +1072,14 @@ void colvar::eigenvector::calc_gradients()
atoms[ia].grad = eigenvec[ia];
}
/*
std::vector<cvm::rvector> gradients (atoms.size());
for (size_t i = 0; i < atoms.size(); i++) {
gradients[i] = atoms[i].grad;
}
std::vector<cvm::rvector> gradients_jh_correct (atoms.size());
/*
// Version B: complex, expensive, fragile. Right.
// gradient of the rotation matrix
@ -1080,7 +1099,6 @@ void colvar::eigenvector::calc_gradients()
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];
@ -1118,9 +1136,76 @@ void colvar::eigenvector::calc_gradients()
}
// Rotate correction term back to reference frame
atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
// atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
// atoms[ia].grad = eigenvec[ia] + quat0.rotate (rot_grad_term);
gradients_jh_correct[ia] = eigenvec[ia] + quat0.rotate (rot_grad_term);
}
*/
// TODO replace this with a call to debug_gradients()
if (b_debug_gradients) {
// the following code should work for any scalar variable;
// the only difference is the name of the atom group object (here, "atoms")
cvm::rotation const rot_0 = atoms.rot;
cvm::rotation const rot_inv = atoms.rot.inverse();
cvm::real const x_0 = x.real_value;
cvm::log ("gradients = "+cvm::to_str (gradients)+"\n");
/* cvm::log ("gradients(JH-correct) = "+cvm::to_str (gradients_jh_correct)+"\n"); */
if (atoms.b_fit_gradients) {
atoms.calc_fit_gradients();
if (atoms.b_rotate) {
// fit_gradients are in the original frame, output them in the rotated frame
for (size_t j = 0; j < atoms.fit_gradients.size(); j++) {
atoms.fit_gradients[j] = rot_0.rotate (atoms.fit_gradients[j]);
}
}
cvm::log ("fit_gradients = "+cvm::to_str (atoms.fit_gradients)+"\n");
if (atoms.b_rotate) {
for (size_t j = 0; j < atoms.fit_gradients.size(); j++) {
atoms.fit_gradients[j] = rot_inv.rotate (atoms.fit_gradients[j]);
}
}
}
for (size_t ia = 0; ia < atoms.size(); ia++) {
// tests are best conducted in the unrotated (simulation) frame
cvm::rvector const atom_grad = atoms.b_rotate ?
rot_inv.rotate (atoms[ia].grad) :
atoms[ia].grad;
for (size_t id = 0; id < 3; id++) {
// (re)read original positions
atoms.read_positions();
// change one coordinate
atoms[ia].pos[id] += cvm::debug_gradients_step_size;
// (re)do the fit (if defined)
if (atoms.b_center || atoms.b_rotate) {
atoms.calc_apply_roto_translation();
}
calc_value();
cvm::real const x_1 = x.real_value;
cvm::log ("Atom "+cvm::to_str (ia)+", component "+cvm::to_str (id)+":\n");
cvm::log ("dx(real) = "+cvm::to_str (x_1 - x_0)+"\n");
cvm::real const dx_pred = atoms.b_fit_gradients ?
(cvm::debug_gradients_step_size * (atom_grad[id] + atoms.fit_gradients[ia][id])) :
(cvm::debug_gradients_step_size * atom_grad[id]);
cvm::log ("dx(pred) = "+cvm::to_str (dx_pred)+"\n");
cvm::log ("|dx(real) - dx(pred)|/dx(real) = "+
cvm::to_str (std::fabs (x_1 - x_0 - dx_pred) /
std::fabs (x_1 - x_0),
cvm::cv_width, cvm::cv_prec)+
".\n");
}
}
}
*/
}
@ -1196,3 +1281,4 @@ void colvar::eigenvector::calc_Jacobian_derivative()

1
lib/colvars/colvargrid.h
View File

@ -175,6 +175,7 @@ public:
/// create() must be called after that;
colvar_grid (colvar_grid<T> const &g) : has_data (false),
nd (g.nd),
nx (g.nx),
mult (g.mult),
cv (g.cv),
lower_boundaries (g.lower_boundaries),

4
lib/colvars/colvarmodule.cpp
View File

@ -46,11 +46,11 @@ colvarmodule::colvarmodule (char const *config_filename,
parse->get_keyval (conf, "analysis", b_analysis, false);
parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03,
parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-07,
colvarparse::parse_silent);
parse->get_keyval (conf, "eigenvalueCrossingThreshold",
colvarmodule::rotation::crossing_threshold, 1.0e-04,
colvarmodule::rotation::crossing_threshold, 1.0e-02,
colvarparse::parse_silent);
parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, 100);

2
lib/colvars/colvarmodule.h
View File

@ -2,7 +2,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2012-11-20"
#define COLVARS_VERSION "2013-01-23"
#endif
#ifndef COLVARS_DEBUG