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

This commit is contained in:
sjplimp 2016-04-28 14:48:56 +00:00
parent 107e28c77a
commit 17fd5898df
15 changed files with 300 additions and 303 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
lib/colvars/colvar.cpp
View File

@ -649,7 +649,7 @@ int colvar::parse_analysis(std::string const &conf)
} else {
cvm::log("Unknown type of correlation function, \""+
acf_type_str+"\".\n");
cvm::set_error_bits(INPUT_ERROR);
cvm::set_error_bit(INPUT_ERROR);
}
get_keyval(conf, "corrFuncOffset", acf_offset, 0);
@ -720,11 +720,11 @@ int colvar::calc()
// Note: if anything is added here, it should be added also in the SMP block of calc_colvars()
int error_code = COLVARS_OK;
if (is_enabled(f_cv_active)) {
error_code |= update_cvc_flags();
error_code |= calc_cvcs();
error_code |= collect_cvc_data();
cvm::combine_errors(error_code, update_cvc_flags());
cvm::combine_errors(error_code, calc_cvcs());
cvm::combine_errors(error_code, collect_cvc_data());
}
return COLVARS_OK;
return error_code;
}
@ -735,15 +735,15 @@ int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
cvm::log("Calculating colvar \""+this->name+"\", components "+
cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
error_code |= check_cvc_range(first_cvc, num_cvcs);
cvm::combine_errors(error_code, check_cvc_range(first_cvc, num_cvcs));
if (error_code != COLVARS_OK) {
return error_code;
}
error_code |= calc_cvc_values(first_cvc, num_cvcs);
error_code |= calc_cvc_gradients(first_cvc, num_cvcs);
error_code |= calc_cvc_sys_forces(first_cvc, num_cvcs);
error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs);
cvm::combine_errors(error_code, calc_cvc_values(first_cvc, num_cvcs));
cvm::combine_errors(error_code, calc_cvc_gradients(first_cvc, num_cvcs));
cvm::combine_errors(error_code, calc_cvc_sys_forces(first_cvc, num_cvcs));
cvm::combine_errors(error_code, calc_cvc_Jacobians(first_cvc, num_cvcs));
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\".\n");
@ -759,12 +759,12 @@ int colvar::collect_cvc_data()
int error_code = COLVARS_OK;
error_code |= collect_cvc_values();
error_code |= collect_cvc_gradients();
error_code |= collect_cvc_sys_forces();
error_code |= collect_cvc_Jacobians();
cvm::combine_errors(error_code, collect_cvc_values());
cvm::combine_errors(error_code, collect_cvc_gradients());
cvm::combine_errors(error_code, collect_cvc_sys_forces());
cvm::combine_errors(error_code, collect_cvc_Jacobians());
error_code |= calc_colvar_properties();
cvm::combine_errors(error_code, calc_colvar_properties());
if (cvm::debug())
cvm::log("Done calculating colvar \""+this->name+"\"'s properties.\n");
@ -879,6 +879,11 @@ int colvar::calc_cvc_gradients(int first_cvc, size_t num_cvcs)
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();
if (cvcs[i]->is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients for " + cvcs[i]->description);
cvcs[i]->debug_gradients(cvcs[i]->atom_groups[ig]);
}
}
}
cvm::decrease_depth();
@ -1297,7 +1302,7 @@ bool colvar::periodic_boundaries(colvarvalue const &lb, colvarvalue const &ub) c
if ( (!is_enabled(f_cv_lower_boundary)) || (!is_enabled(f_cv_upper_boundary)) ) {
cvm::log("Error: checking periodicity for collective variable \""+this->name+"\" "
"requires lower and upper boundaries to be defined.\n");
cvm::set_error_bits(INPUT_ERROR);
cvm::set_error_bit(INPUT_ERROR);
}
if (period > 0.0) {

77
lib/colvars/colvaratoms.cpp
View File

@ -297,7 +297,7 @@ int cvm::atom_group::parse(std::string const &conf)
std::string numbers_conf = "";
size_t pos = 0;
while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) {
parse_error |= add_atom_numbers(numbers_conf);
cvm::combine_errors(parse_error, add_atom_numbers(numbers_conf));
numbers_conf = "";
}
}
@ -306,7 +306,7 @@ int cvm::atom_group::parse(std::string const &conf)
std::string index_group_name;
if (get_keyval(group_conf, "indexGroup", index_group_name)) {
// use an index group from the index file read globally
parse_error |= add_index_group(index_group_name);
cvm::combine_errors(parse_error, add_index_group(index_group_name));
}
}
@ -315,7 +315,7 @@ int cvm::atom_group::parse(std::string const &conf)
size_t pos = 0;
while (key_lookup(group_conf, "atomNumbersRange",
range_conf, pos)) {
parse_error |= add_atom_numbers_range(range_conf);
cvm::combine_errors(parse_error, add_atom_numbers_range(range_conf));
range_conf = "";
}
}
@ -342,8 +342,8 @@ int cvm::atom_group::parse(std::string const &conf)
cvm::error("Error: more instances of \"atomNameResidueRange\" than "
"values of \"psfSegID\".\n", INPUT_ERROR);
} else {
parse_error |= add_atom_name_residue_range(psf_segids.size() ? *psii : std::string(""),
range_conf);
cvm::combine_errors(parse_error, add_atom_name_residue_range(psf_segids.size() ?
*psii : std::string(""), range_conf));
if (psf_segids.size()) psii++;
}
range_conf = "";
@ -407,7 +407,7 @@ int cvm::atom_group::parse(std::string const &conf)
index = (cvm::proxy)->init_atom_group(atoms_ids);
}
parse_error |= parse_fitting_options(group_conf);
cvm::combine_errors(parse_error, parse_fitting_options(group_conf));
// TODO move this to colvarparse object
check_keywords(group_conf, key.c_str());
@ -612,7 +612,6 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
// regardless of the configuration, fit gradients must be calculated by refPositionsGroup
ref_pos_group->b_fit_gradients = this->b_fit_gradients;
this->b_fit_gradients = false;
}
atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this;
@ -682,7 +681,7 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
"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());
rot.request_group1_gradients(group_for_fit->size());
}
}
@ -778,13 +777,19 @@ int cvm::atom_group::calc_required_properties()
void cvm::atom_group::calc_apply_roto_translation()
{
// store the laborarory-frame COGs for when they are needed later
cog_orig = this->center_of_geometry();
if (ref_pos_group) {
ref_pos_group->cog_orig = ref_pos_group->center_of_geometry();
}
if (b_center) {
// center on the origin first
cvm::atom_pos const cog = ref_pos_group ?
cvm::atom_pos const rpg_cog = ref_pos_group ?
ref_pos_group->center_of_geometry() : this->center_of_geometry();
apply_translation(-1.0 * cog);
apply_translation(-1.0 * rpg_cog);
if (ref_pos_group) {
ref_pos_group->apply_translation(-1.0 * cog);
ref_pos_group->apply_translation(-1.0 * rpg_cog);
}
}
@ -945,26 +950,23 @@ void cvm::atom_group::calc_fit_gradients()
{
if (b_dummy) return;
if ((!b_center) && (!b_rotate)) return; // no fit
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
cvm::rvector atom_grad;
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(atoms[i].grad) :
atoms[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;
}
atom_grad += atoms[i].grad;
}
if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad);
atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
for (size_t j = 0; j < group_for_fit->size(); j++) {
group_for_fit->fit_gradients[j] = atom_grad;
}
}
@ -972,23 +974,26 @@ void cvm::atom_group::calc_fit_gradients()
// add the rotation matrix contribution to the gradients
cvm::rotation const rot_inv = rot.inverse();
cvm::atom_pos const cog = center_of_geometry();
for (size_t i = 0; i < this->size(); i++) {
cvm::atom_pos const pos_orig = rot_inv.rotate((b_center ? (atoms[i].pos - cog) : (atoms[i].pos)));
// compute centered, unrotated position
cvm::atom_pos const pos_orig =
rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
// 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, atoms[i].grad);
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, atoms[i].grad);
// multiply by \cdot {\partial q}/\partial\vec{x}_j and add it to the fit gradients
// multiply by {\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");
}
@ -1138,17 +1143,9 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force)
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 * group_for_fit->fit_gradients[j]));
}
} else {
for (size_t j = 0; j < group_for_fit->size(); j++) {
(*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
}
// Fit gradients are already calculated in "laboratory" frame
for (size_t j = 0; j < group_for_fit->size(); j++) {
(*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
}
}

9
lib/colvars/colvaratoms.h
View File

@ -286,7 +286,7 @@ public:
bool b_user_defined_fit;
/// \brief Whether or not the derivatives of the roto-translation
/// should be included when calculating the colvar's gradients (default: no)
/// should be included when calculating the colvar's gradients (default: yes)
bool b_fit_gradients;
/// \brief use reference coordinates for b_center or b_rotate
@ -359,10 +359,17 @@ public:
/// \brief Calculate the center of geometry of the atomic positions, assuming
/// that they are already pbc-wrapped
int calc_center_of_geometry();
private:
/// \brief Center of geometry
cvm::atom_pos cog;
/// \brief Center of geometry before any fitting
cvm::atom_pos cog_orig;
public:
/// \brief Return the center of geometry of the atomic positions
inline cvm::atom_pos center_of_geometry() const
{

5
lib/colvars/colvarbias_histogram.cpp
View File

@ -95,7 +95,7 @@ int colvarbias_histogram::update()
{
int error_code = COLVARS_OK;
// update base class
error_code |= colvarbias::update();
cvm::combine_errors(error_code, colvarbias::update());
if (cvm::debug()) {
cvm::log("Updating histogram bias " + this->name);
@ -148,7 +148,8 @@ int colvarbias_histogram::update()
write_output_files();
}
return error_code | cvm::get_error();
cvm::combine_errors(error_code, cvm::get_error());
return error_code;
}

136
lib/colvars/colvarcomp.cpp
View File

@ -26,7 +26,12 @@ colvar::cvc::cvc(std::string const &conf)
init_cvc_requires();
get_keyval(conf, "name", this->name, std::string(""), parse_silent);
if (get_keyval(conf, "name", this->name, std::string(""), parse_silent)) {
// Temporary description until child object is initialized
description = "cvc " + name;
} else {
description = "uninitialized cvc";
}
get_keyval(conf, "componentCoeff", sup_coeff, 1.0);
get_keyval(conf, "componentExp", sup_np, 1);
@ -34,6 +39,8 @@ colvar::cvc::cvc(std::string const &conf)
get_keyval(conf, "period", period, 0.0);
get_keyval(conf, "wrapAround", wrap_center, 0.0);
// All cvcs implement this
provide(f_cvc_debug_gradient);
{
bool b_debug_gradient;
get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent);
@ -157,6 +164,8 @@ 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")
// NOTE: this assumes that groups for this cvc are non-overlapping,
// since atom coordinates are modified only within the current group
if (group->b_dummy) return;
@ -168,31 +177,34 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
// 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_rotate || group->b_center)
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]);
}
cvm::atom_group *group_for_fit = group->ref_pos_group ? group->ref_pos_group : group;
cvm::atom_pos fit_gradient_sum, gradient_sum;
// print the values of the fit gradients
if (group->b_rotate || group->b_center) {
if (group->b_fit_gradients) {
size_t j;
// fit_gradients are in the simulation frame: we should print them in the rotated frame
cvm::log("Fit gradients:\n");
for (j = 0; j < group_for_fit->fit_gradients.size(); j++) {
cvm::log((group->ref_pos_group ? std::string("refPosGroup") : group->key) +
"[" + cvm::to_str(j) + "] = " +
(group->b_rotate ?
cvm::to_str(rot_0.rotate(group_for_fit->fit_gradients[j])) :
cvm::to_str(group_for_fit->fit_gradients[j])));
}
}
}
// debug the gradients
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;
cvm::rvector const atom_grad = (group->b_rotate ?
rot_inv.rotate((*group)[ia].grad) :
(*group)[ia].grad);
gradient_sum += atom_grad;
for (size_t id = 0; id < 3; id++) {
// (re)read original positions
@ -205,55 +217,59 @@ void colvar::cvc::debug_gradients(cvm::atom_group *group)
if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0];
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->fit_gradients.size() && (group->ref_pos_group == NULL)) ?
21, 14)+"\n");
cvm::real const dx_pred = (group->fit_gradients.size()) ?
(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");
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");
std::fabs(x_1 - x_0), 12, 5)+"\n");
}
}
/*
* The code below is WIP
*/
// if (group->ref_pos_group != NULL) {
// cvm::atom_group &ref = *group->ref_pos_group;
// group->calc_fit_gradients();
//
// for (size_t ia = 0; ia < ref.size(); ia++) {
//
// for (size_t id = 0; id < 3; id++) {
// // (re)read original positions
// group->read_positions();
// ref.read_positions();
// // change one coordinate
// ref[ia].pos[id] += cvm::debug_gradients_step_size;
// group->update();
// calc_value();
// cvm::real const x_1 = x.real_value;
// cvm::log("refPosGroup 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->fit_gradients.size() && (group->ref_pos_group == NULL)) ?
// // cvm::real const dx_pred = (group->fit_gradients.size()) ?
// // (cvm::debug_gradients_step_size * (atom_grad[id] + group->fit_gradients[ia][id])) :
// // (cvm::debug_gradients_step_size * atom_grad[id]);
// cvm::real const dx_pred = cvm::debug_gradients_step_size * ref.fit_gradients[ia][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");
// }
// }
// }
if ((group->b_fit_gradients) && (group->ref_pos_group != NULL)) {
cvm::atom_group *ref_group = group->ref_pos_group;
group->read_positions();
group->calc_required_properties();
for (size_t ia = 0; ia < ref_group->size(); ia++) {
// fit gradients are in the unrotated (simulation) frame
cvm::rvector const atom_grad = ref_group->fit_gradients[ia];
fit_gradient_sum += atom_grad;
for (size_t id = 0; id < 3; id++) {
// (re)read original positions
group->read_positions();
ref_group->read_positions();
// change one coordinate
(*ref_group)[ia].pos[id] += cvm::debug_gradients_step_size;
group->calc_required_properties();
calc_value();
cvm::real const x_1 = x.real_value;
cvm::log("refPosGroup 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 = 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");
}
}
}
cvm::log("Gradient sum: " + cvm::to_str(gradient_sum) +
" Fit gradient sum: " + cvm::to_str(fit_gradient_sum) +
" Total " + cvm::to_str(gradient_sum + fit_gradient_sum));
return;
}

6
lib/colvars/colvarcomp_angles.cpp
View File

@ -84,12 +84,6 @@ void colvar::angle::calc_gradients()
group1->set_weighted_gradient(dxdr1);
group2->set_weighted_gradient((dxdr1 + dxdr3) * (-1.0));
group3->set_weighted_gradient(dxdr3);
if (is_enabled(f_cvc_debug_gradient)) {
debug_gradients(group1);
debug_gradients(group2);
debug_gradients(group3);
}
}
void colvar::angle::calc_force_invgrads()

69
lib/colvars/colvarcomp_distances.cpp
View File

@ -167,7 +167,7 @@ colvar::distance_z::distance_z(std::string const &conf)
// this group is optional
ref2 = parse_group(conf, "ref2", true);
if (ref2->size()) {
if (ref2 && ref2->size()) {
cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
fixed_axis = false;
if (key_lookup(conf, "axis"))
@ -249,15 +249,6 @@ void colvar::distance_z::calc_gradients()
cvm::position_distance(main->center_of_mass(), ref1->center_of_mass()) + x.real_value * axis ));
}
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients for group main:\n");
debug_gradients(main);
cvm::log("Debugging gradients for group ref1:\n");
debug_gradients(ref1);
cvm::log("Debugging gradients for group ref2:\n");
debug_gradients(ref2);
}
}
void colvar::distance_z::calc_force_invgrads()
@ -282,7 +273,7 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
if (!ref1->noforce)
ref1->apply_colvar_force(force.real_value);
if (ref2->size() && !ref2->noforce)
if (ref2 && ref2->size() && !ref2->noforce)
ref2->apply_colvar_force(force.real_value);
if (!main->noforce)
@ -382,7 +373,7 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
if (!ref1->noforce)
ref1->apply_colvar_force(force.real_value);
if (ref2->size() && !ref2->noforce)
if (ref2 && ref2->size() && !ref2->noforce)
ref2->apply_colvar_force(force.real_value);
if (!main->noforce)
@ -593,10 +584,6 @@ void colvar::distance_pairs::calc_value()
void colvar::distance_pairs::calc_gradients()
{
// will be calculated on the fly in apply_force()
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients:\n");
debug_gradients(group1);
}
}
void colvar::distance_pairs::apply_force(colvarvalue const &force)
@ -667,11 +654,6 @@ void colvar::gyration::calc_gradients()
for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
ai->grad = drdx * ai->pos;
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients:\n");
debug_gradients(atoms);
}
}
@ -731,11 +713,6 @@ void colvar::inertia::calc_gradients()
for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
ai->grad = 2.0 * ai->pos;
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients:\n");
debug_gradients(atoms);
}
}
@ -786,11 +763,6 @@ void colvar::inertia_z::calc_gradients()
for (cvm::atom_iter ai = atoms->begin(); ai != atoms->end(); ai++) {
ai->grad = 2.0 * (ai->pos * axis) * axis;
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients:\n");
debug_gradients(atoms);
}
}
@ -811,7 +783,7 @@ colvar::rmsd::rmsd(std::string const &conf)
atoms = parse_group(conf, "atoms");
if (atoms->size() == 0) {
if (!atoms || atoms->size() == 0) {
cvm::error("Error: \"atoms\" must contain at least 1 atom to compute RMSD.");
return;
}
@ -891,18 +863,11 @@ colvar::rmsd::rmsd(std::string const &conf)
cvm::log("This is a standard minimum RMSD, derivatives of the optimal rotation "
"will not be computed as they cancel out in the gradients.");
atoms->b_fit_gradients = false;
}
if (atoms->b_rotate) {
// TODO: finer-grained control of this would require exposing a
// "request_Jacobian_derivative()" method to the colvar, and the same
// from the colvar to biases
// TODO: this should not be enabled here anyway, as it is not specific of the
// component - instead it should be decided in a generic way by the atom group
// 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());
}
}
@ -930,11 +895,6 @@ 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 (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients:\n");
debug_gradients(atoms);
}
}
@ -1031,7 +991,7 @@ colvar::eigenvector::eigenvector(std::string const &conf)
cvm::log("Using reference positions from input file.\n");
if (ref_pos.size() != atoms->size()) {
cvm::error("Error: reference positions do not "
"match the number of requested atoms->\n");
"match the number of requested atoms.\n");
return;
}
}
@ -1064,9 +1024,14 @@ colvar::eigenvector::eigenvector(std::string const &conf)
}
}
if (ref_pos.size() == 0) {
cvm::error("Error: reference positions were not provided.\n", INPUT_ERROR);
return;
}
if (ref_pos.size() != atoms->size()) {
cvm::error("Error: reference positions were not provided, or do not "
"match the number of requested atoms->\n");
cvm::error("Error: reference positions do not "
"match the number of requested atoms.\n", INPUT_ERROR);
return;
}
@ -1087,6 +1052,8 @@ colvar::eigenvector::eigenvector(std::string const &conf)
atoms->b_rotate = true;
atoms->ref_pos = ref_pos;
atoms->center_ref_pos();
atoms->b_fit_gradients = false; // cancel out if group is fitted on itself
// and cvc is translationally invariant
// request the calculation of the derivatives of the rotation defined by the atom group
atoms->rot.request_group1_gradients(atoms->size());
@ -1216,11 +1183,6 @@ void colvar::eigenvector::calc_gradients()
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = eigenvec[ia];
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging gradients:\n");
debug_gradients(atoms);
}
}
@ -1313,6 +1275,7 @@ colvar::cartesian::cartesian(std::string const &conf)
if (axes.size() == 0) {
cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");
return;
}
x.type(colvarvalue::type_vector);

22
lib/colvars/colvarcomp_rotations.cpp
View File

@ -57,14 +57,14 @@ colvar::orientation::orientation(std::string const &conf)
cvm::log("Centering the reference coordinates: it is "
"assumed that each atom is the closest "
"periodic image to the center of geometry.\n");
cvm::rvector cog(0.0, 0.0, 0.0);
cvm::rvector ref_cog(0.0, 0.0, 0.0);
size_t i;
for (i = 0; i < ref_pos.size(); i++) {
cog += ref_pos[i];
ref_cog += ref_pos[i];
}
cog /= cvm::real(ref_pos.size());
ref_cog /= cvm::real(ref_pos.size());
for (i = 0; i < ref_pos.size(); i++) {
ref_pos[i] -= cog;
ref_pos[i] -= ref_cog;
}
get_keyval(conf, "closestToQuaternion", ref_quat, cvm::quaternion(1.0, 0.0, 0.0, 0.0));
@ -87,6 +87,7 @@ colvar::orientation::orientation()
void colvar::orientation::calc_value()
{
rot.b_debug_gradients = is_enabled(f_cvc_debug_gradient);
atoms_cog = atoms->center_of_geometry();
rot.calc_optimal_rotation(ref_pos, atoms->positions_shifted(-1.0 * atoms_cog));
@ -163,10 +164,6 @@ void colvar::orientation_angle::calc_gradients()
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging orientationAngle component gradients:\n");
debug_gradients(atoms);
}
}
@ -210,10 +207,6 @@ void colvar::orientation_proj::calc_gradients()
for (size_t ia = 0; ia < atoms->size(); ia++) {
(*atoms)[ia].grad = (dxdq0 * (rot.dQ0_2[ia])[0]);
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging orientationProj component gradients:\n");
debug_gradients(atoms);
}
}
@ -272,11 +265,6 @@ void colvar::tilt::calc_gradients()
(*atoms)[ia].grad += (dxdq[iq] * (rot.dQ0_2[ia])[iq]);
}
}
if (is_enabled(f_cvc_debug_gradient)) {
cvm::log("Debugging tilt component gradients:\n");
debug_gradients(atoms);
}
}

3
lib/colvars/colvardeps.cpp
View File

@ -355,6 +355,9 @@ void cvm::deps::init_cvc_requires() {
f_description(f_cvc_inv_gradient, "inverse gradient");
f_req_self(f_cvc_inv_gradient, f_cvc_gradient);
f_description(f_cvc_debug_gradient, "debug gradient");
f_req_self(f_cvc_debug_gradient, f_cvc_gradient);
f_description(f_cvc_Jacobian, "Jacobian");
f_req_self(f_cvc_Jacobian, f_cvc_inv_gradient);

78
lib/colvars/colvarmodule.cpp
View File

@ -316,7 +316,8 @@ int colvarmodule::parse_biases(std::string const &conf)
int colvarmodule::catch_input_errors(int result)
{
if (result != COLVARS_OK || get_error()) {
set_error_bits(result | INPUT_ERROR);
set_error_bit(result);
set_error_bit(INPUT_ERROR);
parse->reset();
return get_error();
}
@ -468,27 +469,27 @@ int colvarmodule::calc()
cvm::to_str(cvm::step_absolute())+"\n");
}
error_code |= calc_colvars();
combine_errors(error_code, calc_colvars());
// set biasing forces to zero before biases are calculated and summed over
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end(); cvi++) {
(*cvi)->reset_bias_force();
}
error_code |= calc_biases();
error_code |= update_colvar_forces();
combine_errors(error_code, calc_biases());
combine_errors(error_code, update_colvar_forces());
if (cvm::b_analysis) {
error_code |= analyze();
combine_errors(error_code, analyze());
}
// write trajectory files, if needed
if (cv_traj_freq && cv_traj_name.size()) {
error_code |= write_traj_files();
combine_errors(error_code, write_traj_files());
}
// write restart files, if needed
if (restart_out_freq && restart_out_name.size()) {
error_code |= write_restart_files();
combine_errors(error_code, write_restart_files());
}
return error_code;
@ -520,7 +521,7 @@ int colvarmodule::calc_colvars()
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
error_code |= (*cvi)->update_cvc_flags();
combine_errors(error_code, (*cvi)->update_cvc_flags());
size_t num_items = (*cvi)->num_active_cvcs();
colvars_smp.reserve(colvars_smp.size() + num_items);
@ -535,11 +536,11 @@ int colvarmodule::calc_colvars()
cvm::decrease_depth();
// calculate colvar components in parallel
error_code |= proxy->smp_colvars_loop();
combine_errors(error_code, proxy->smp_colvars_loop());
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
error_code |= (*cvi)->collect_cvc_data();
combine_errors(error_code, (*cvi)->collect_cvc_data());
}
cvm::decrease_depth();
@ -548,7 +549,7 @@ int colvarmodule::calc_colvars()
// calculate colvars one at a time
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
error_code |= (*cvi)->calc();
combine_errors(error_code, (*cvi)->calc());
if (cvm::get_error()) {
return COLVARS_ERROR;
}
@ -556,7 +557,8 @@ int colvarmodule::calc_colvars()
cvm::decrease_depth();
}
return error_code | cvm::get_error();
combine_errors(error_code, cvm::get_error());
return error_code;
}
@ -575,21 +577,21 @@ int colvarmodule::calc_biases()
if (use_scripted_forces && !scripting_after_biases) {
// calculate biases and scripted forces in parallel
error_code |= proxy->smp_biases_script_loop();
combine_errors(error_code, proxy->smp_biases_script_loop());
} else {
// calculate biases in parallel
error_code |= proxy->smp_biases_loop();
combine_errors(error_code, proxy->smp_biases_loop());
}
} else {
if (use_scripted_forces && !scripting_after_biases) {
error_code |= calc_scripted_forces();
combine_errors(error_code, calc_scripted_forces());
}
cvm::increase_depth();
for (bi = biases.begin(); bi != biases.end(); bi++) {
error_code |= (*bi)->update();
combine_errors(error_code, (*bi)->update());
if (cvm::get_error()) {
return COLVARS_ERROR;
}
@ -627,7 +629,7 @@ int colvarmodule::update_colvar_forces()
cvm::decrease_depth();
if (use_scripted_forces && scripting_after_biases) {
error_code |= calc_scripted_forces();
combine_errors(error_code, calc_scripted_forces());
}
cvm::real total_colvar_energy = 0.0;
@ -875,17 +877,17 @@ int colvarmodule::setup_output()
std::string(""));
if (cv_traj_freq && cv_traj_name.size()) {
error_code |= open_traj_file(cv_traj_name);
combine_errors(error_code, open_traj_file(cv_traj_name));
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
error_code |= (*bi)->setup_output();
combine_errors(error_code, (*bi)->setup_output());
}
if (error_code != COLVARS_OK || cvm::get_error()) {
set_error_bits(FILE_ERROR);
set_error_bit(FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -904,6 +906,13 @@ std::istream & colvarmodule::read_restart(std::istream &is)
colvarparse::parse_silent);
it = it_restart;
}
std::string restart_version;
parse->get_keyval(restart_conf, "version",
restart_version, std::string(""),
colvarparse::parse_silent);
if (restart_version.size() && (restart_version != std::string(COLVARS_VERSION))) {
cvm::log("This state file was generated with version "+restart_version+"\n");
}
}
is.clear();
parse->clear_keyword_registry();
@ -1060,6 +1069,7 @@ std::ostream & colvarmodule::write_restart(std::ostream &os)
<< " step " << std::setw(it_width)
<< it << "\n"
<< " dt " << dt() << "\n"
<< " version " << std::string(COLVARS_VERSION) << "\n"
<< "}\n\n";
cvm::increase_depth();
@ -1213,14 +1223,30 @@ size_t & cvm::depth()
}
void colvarmodule::set_error_bits(int code)
void colvarmodule::set_error_bit(int code)
{
if (code > 0) {
cvm::fatal_error("Error: set_error_bit() received positive error code.\n");
return;
}
proxy->smp_lock();
errorCode |= code;
errorCode |= COLVARS_ERROR;
errorCode = -1 * ((-1 * errorCode) | (-1 * code) | (-1 * COLVARS_ERROR));
proxy->smp_unlock();
}
bool colvarmodule::get_error_bit(int code)
{
return bool((-1 * errorCode) & (-1 * code));
}
void colvarmodule::combine_errors(int &target, int const code)
{
if (code > 0 || target > 0) {
cvm::fatal_error("Error: combine_errors() received positive error code.\n");
return;
}
target = -1 * ((-1 * target) | (-1 * code));
}
void colvarmodule::clear_error()
{
@ -1232,7 +1258,7 @@ void colvarmodule::clear_error()
void cvm::error(std::string const &message, int code)
{
set_error_bits(code);
set_error_bit(code);
proxy->error(message);
}
@ -1241,7 +1267,7 @@ void cvm::fatal_error(std::string const &message)
{
// TODO once all non-fatal errors have been set to be handled by error(),
// set DELETE_COLVARS here for VMD to handle it
set_error_bits(FATAL_ERROR);
set_error_bit(FATAL_ERROR);
proxy->fatal_error(message);
}
@ -1394,7 +1420,7 @@ colvarproxy *colvarmodule::proxy = NULL;
// static runtime data
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03;
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07;
int colvarmodule::errorCode = 0;
long colvarmodule::it = 0;
long colvarmodule::it_restart = 0;

32
lib/colvars/colvarmodule.h
View File

@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2016-04-14"
#define COLVARS_VERSION "2016-04-27"
#endif
#ifndef COLVARS_DEBUG
@ -20,19 +20,19 @@
/// shared between all object instances) to be accessed from other
/// objects.
// Error codes
// Error codes are the negative of powers-of-two
// as a result, error codes should NOT be bitwise-ORed but only
// accessed through set_error_bit() and get_error_bit()
#define COLVARS_OK 0
#define COLVARS_ERROR -1
#define GENERAL_ERROR -1 // TODO this can be simply merged with COLVARS_ERROR
#define COLVARS_NOT_IMPLEMENTED (-1<<1)
#define INPUT_ERROR (-1<<2) // out of bounds or inconsistent input
#define BUG_ERROR (-1<<3) // Inconsistent state indicating bug
#define FILE_ERROR (-1<<4)
#define MEMORY_ERROR (-1<<5)
#define FATAL_ERROR (-1<<6) // Should be set, or not, together with other bits
#define DELETE_COLVARS (-1<<7) // Instruct the caller to delete cvm
#define COLVARS_NO_SUCH_FRAME (-1<<8) // Cannot load the requested frame
#define COLVARS_NOT_IMPLEMENTED (-1*(1<<1))
#define INPUT_ERROR (-1*(1<<2)) // out of bounds or inconsistent input
#define BUG_ERROR (-1*(1<<3)) // Inconsistent state indicating bug
#define FILE_ERROR (-1*(1<<4))
#define MEMORY_ERROR (-1*(1<<5))
#define FATAL_ERROR (-1*(1<<6)) // Should be set, or not, together with other bits
#define DELETE_COLVARS (-1*(1<<7)) // Instruct the caller to delete cvm
#define COLVARS_NO_SUCH_FRAME (-1*(1<<8)) // Cannot load the requested frame
#include <iostream>
#include <iomanip>
@ -115,7 +115,11 @@ protected:
public:
static void set_error_bits(int code);
static void set_error_bit(int code);
static bool get_error_bit(int code);
static void combine_errors(int &target, int const code);
static inline int get_error()
{
@ -395,7 +399,7 @@ public:
static void fatal_error(std::string const &message);
/// Print a message to the main log and set global error code
static void error(std::string const &message, int code = GENERAL_ERROR);
static void error(std::string const &message, int code = COLVARS_ERROR);
/// Print a message to the main log and exit normally
static void exit(std::string const &message);

3
lib/colvars/colvarproxy.h
View File

@ -69,10 +69,11 @@ public:
virtual cvm::real rand_gaussian(void) = 0;
/// \brief Get the current frame number
// Negative return values indicate error
virtual int frame() { return COLVARS_NOT_IMPLEMENTED; }
/// \brief Set the current frame number
// return 0 on success, -1 on failure
// return frame number on success, COLVARS_NO_SUCH_FRAME otherwise
virtual int frame(int) { return COLVARS_NOT_IMPLEMENTED; }
/// \brief Prefix to be used for input files (restarts, not

15
lib/colvars/colvarscript.cpp
View File

@ -53,17 +53,16 @@ int colvarscript::run(int argc, char const *argv[]) {
}
if (cmd == "delete") {
colvars->reset();
// Note: the delete bit may be ignored by some backends
// it is mostly useful in VMD
colvars->set_error_bits(DELETE_COLVARS);
colvars->set_error_bit(DELETE_COLVARS);
return COLVARS_OK;
}
if (cmd == "update") {
error_code |= proxy->update_input();
error_code |= colvars->calc();
error_code |= proxy->update_output();
cvm::combine_errors(error_code, proxy->update_input());
cvm::combine_errors(error_code, colvars->calc());
cvm::combine_errors(error_code, proxy->update_output());
if (error_code) {
result += "Error updating the colvars module.\n";
}
@ -143,8 +142,8 @@ int colvarscript::run(int argc, char const *argv[]) {
}
proxy->output_prefix_str = argv[2];
int error = 0;
error |= colvars->setup_output();
error |= colvars->write_output_files();
cvm::combine_errors(error, colvars->setup_output());
cvm::combine_errors(error, colvars->write_output_files());
return error ? COLVARSCRIPT_ERROR : COLVARS_OK;
}
@ -406,7 +405,7 @@ Input and output:\n\
list biases -- return a list of all biases\n\
load <file name> -- load a state file (requires configuration)\n\
save <file name> -- save a state file (requires configuration)\n\
update -- recalculate colvars and biases based\n\
update -- recalculate colvars and biases\n\
printframe -- return a summary of the current frame\n\
printframelabels -- return labels to annotate printframe's output\n";

2
lib/colvars/colvarscript.h
View File

@ -9,7 +9,7 @@
#include "colvarbias.h"
#include "colvarproxy.h"
// TODO merge these into colvarmodule.h
// Only these error values are part of the scripting interface
#define COLVARSCRIPT_ERROR -1
#define COLVARSCRIPT_OK 0

109
lib/colvars/colvartypes.cpp
View File

@ -312,10 +312,8 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
S_backup.resize(4,4);
S_backup = S;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log("S = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
}
if (b_debug_gradients) {
cvm::log("S = "+cvm::to_str(cvm::to_str(S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
}
diagonalize_matrix(S, S_eigval, S_eigvec);
@ -344,25 +342,23 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
q_old = q;
}
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
"\n");
}
if (b_debug_gradients) {
cvm::log("L0 = "+cvm::to_str(L0, cvm::cv_width, cvm::cv_prec)+
", Q0 = "+cvm::to_str(Q0, cvm::cv_width, cvm::cv_prec)+
", Q0*Q0 = "+cvm::to_str(Q0.inner(Q0), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L1 = "+cvm::to_str(L1, cvm::cv_width, cvm::cv_prec)+
", Q1 = "+cvm::to_str(Q1, cvm::cv_width, cvm::cv_prec)+
", Q0*Q1 = "+cvm::to_str(Q0.inner(Q1), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L2 = "+cvm::to_str(L2, cvm::cv_width, cvm::cv_prec)+
", Q2 = "+cvm::to_str(Q2, cvm::cv_width, cvm::cv_prec)+
", Q0*Q2 = "+cvm::to_str(Q0.inner(Q2), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log("L3 = "+cvm::to_str(L3, cvm::cv_width, cvm::cv_prec)+
", Q3 = "+cvm::to_str(Q3, cvm::cv_width, cvm::cv_prec)+
", Q0*Q3 = "+cvm::to_str(Q0.inner(Q3), cvm::cv_width, cvm::cv_prec)+
"\n");
}
// calculate derivatives of L0 and Q0 with respect to each atom in
@ -472,46 +468,43 @@ void colvarmodule::rotation::calc_optimal_rotation(std::vector<cvm::atom_pos> co
}
}
if (cvm::debug()) {
if (b_debug_gradients) {
if (b_debug_gradients) {
cvm::matrix2d<cvm::real> S_new(4, 4);
cvm::vector1d<cvm::real> S_new_eigval(4);
cvm::matrix2d<cvm::real> S_new_eigvec(4, 4);
cvm::matrix2d<cvm::real> S_new(4, 4);
cvm::vector1d<cvm::real> S_new_eigval(4);
cvm::matrix2d<cvm::real> S_new_eigvec(4, 4);
// make an infitesimal move along each cartesian coordinate of
// this atom, and solve again the eigenvector problem
for (size_t comp = 0; comp < 3; comp++) {
// make an infitesimal move along each cartesian coordinate of
// this atom, and solve again the eigenvector problem
for (size_t comp = 0; comp < 3; comp++) {
S_new = S_backup;
// diagonalize the new overlap matrix
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
S_new[i][j] +=
colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
}
S_new = S_backup;
// diagonalize the new overlap matrix
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
S_new[i][j] +=
colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
}
// cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
cvm::real const &L0_new = S_new_eigval[0];
cvm::quaternion const Q0_new(S_new_eigvec[0]);
cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
cvm::log( "|(l_0+dl_0) - l_0^new|/l_0 = "+
cvm::to_str(std::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
", |(q_0+dq_0) - q_0^new| = "+
cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
"\n");
}
// cvm::log("S_new = "+cvm::to_str(cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
diagonalize_matrix(S_new, S_new_eigval, S_new_eigvec);
cvm::real const &L0_new = S_new_eigval[0];
cvm::quaternion const Q0_new(S_new_eigvec[0]);
cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
cvm::quaternion const DQ0(dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
cvm::log( "|(l_0+dl_0) - l_0^new|/l_0 = "+
cvm::to_str(std::fabs(L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
", |(q_0+dq_0) - q_0^new| = "+
cvm::to_str((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
}