From 5f0578f95a42545cfb2428b3206b3cd3077cde3c Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 22 Apr 2016 20:50:51 -0400 Subject: [PATCH] update colvars to version of 2016-04-20 --- lib/colvars/colvar.cpp | 36 +++---- lib/colvars/colvaratoms.cpp | 60 +++++++----- lib/colvars/colvaratoms.h | 9 +- lib/colvars/colvarbias_histogram.cpp | 5 +- lib/colvars/colvarcomp.cpp | 132 +++++++++++++++----------- lib/colvars/colvarcomp.h | 3 + lib/colvars/colvarcomp_angles.cpp | 6 -- lib/colvars/colvarcomp_distances.cpp | 50 ++-------- lib/colvars/colvarcomp_rotations.cpp | 13 --- lib/colvars/colvardeps.cpp | 3 + lib/colvars/colvarmodule.cpp | 76 ++++++++++----- lib/colvars/colvarmodule.h | 32 ++++--- lib/colvars/colvarproxy.h | 3 +- lib/colvars/colvarscript.cpp | 15 ++- lib/colvars/colvarscript.h | 2 +- src/USER-COLVARS/colvarproxy_lammps.h | 1 - src/USER-COLVARS/fix_colvars.cpp | 1 + 17 files changed, 236 insertions(+), 211 deletions(-) diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index ebc3bf44f1..1ee35864bb 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -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,10 @@ 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)) { + cvcs[i]->debug_gradients(cvcs[i]->atom_groups[ig]); + } } } cvm::decrease_depth(); @@ -1297,7 +1301,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) { diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index f13f6993b6..2ab20ca4ec 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -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; @@ -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 = cog; + if (ref_pos_group) { + ref_pos_group->cog_orig = ref_pos_group->cog; + } + 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); } } @@ -951,20 +956,21 @@ void cvm::atom_group::calc_fit_gradients() 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 += (-1.0)/(cvm::real(group_for_fit->size())) * atoms[i].grad; + } + + // need to use the gradients in original frame + // we only rotate the sum for efficiency + if (b_rotate) atom_grad = (rot.inverse()).rotate(atom_grad); + + for (size_t j = 0; j < group_for_fit->size(); j++) { + group_for_fit->fit_gradients[j] = atom_grad; } } @@ -972,16 +978,19 @@ 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))); + // restore original position for this atom + cvm::atom_pos const pos_orig = + rot_inv.rotate((b_center ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos))) + + (ref_pos_group ? ref_pos_group->cog_orig : cog_orig); + + // 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 for (size_t iq = 0; iq < 4; iq++) { group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq]; @@ -989,6 +998,7 @@ void cvm::atom_group::calc_fit_gradients() } } } + if (cvm::debug()) cvm::log("Done calculating fit gradients.\n"); } diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index 827601fc4e..518ae53479 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -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 { diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index bb55669614..ea91a5a55b 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -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; } diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index 2f9ea78398..9e4d75ecaa 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -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,15 +39,14 @@ colvar::cvc::cvc(std::string const &conf) get_keyval(conf, "period", period, 0.0); get_keyval(conf, "wrapAround", wrap_center, 0.0); - { - bool b_debug_gradient; - get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent); - if (b_debug_gradient) enable(f_cvc_debug_gradient); - } + get_keyval(conf, "debugGradients", b_debug_gradient, false, parse_silent); // Attempt scalable calculations when in parallel? (By default yes, if available) get_keyval(conf, "scalable", b_try_scalable, true); + // All cvcs implement this + provide(f_cvc_debug_gradient); + if (cvm::debug()) cvm::log("Done initializing cvc base object.\n"); } @@ -103,6 +107,7 @@ int colvar::cvc::setup() if (b_try_scalable && is_available(f_cvc_scalable)) { enable(f_cvc_scalable); } + if (b_debug_gradient) enable(f_cvc_debug_gradient); return COLVARS_OK; } @@ -157,6 +162,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,25 +175,36 @@ 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(); + cvm::atom_group *group_for_fit = group->ref_pos_group ? group->ref_pos_group : group; + + // 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 original frame: we should print them in the rotated frame 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]); + for (j = 0; j < group_for_fit->fit_gradients.size(); j++) { + group_for_fit->fit_gradients[j] = rot_0.rotate(group_for_fit->fit_gradients[j]); } } - cvm::log("fit_gradients = "+cvm::to_str(group->fit_gradients)+"\n"); + + 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) + "] = " + cvm::to_str(group_for_fit->fit_gradients[j])); + } + 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 (j = 0; j < group_for_fit->fit_gradients.size(); j++) { + group_for_fit->fit_gradients[j] = rot_inv.rotate(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 @@ -205,55 +223,55 @@ 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++) { + + // tests are best conducted in the unrotated (simulation) frame + cvm::rvector const atom_grad = ref_group->fit_gradients[ia]; + + 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"); + } + } + } return; } diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 923b2c6bd1..027cd3499f 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -224,6 +224,9 @@ public: /// \brief Whether or not this CVC will be computed in parallel whenever possible bool b_try_scalable; +private: + bool b_debug_gradient; + protected: /// \brief Cached value diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 519288a01f..1f825a17c4 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -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() diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 8d15c5d67a..a2a5102b8b 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -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() @@ -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); - } } @@ -930,11 +902,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 +998,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 +1031,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; } @@ -1216,11 +1188,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 +1280,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); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 505f95e072..3bfab768ab 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -163,10 +163,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 +206,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 +264,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); - } } diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp index f7d829968b..5821c8b24a 100644 --- a/lib/colvars/colvardeps.cpp +++ b/lib/colvars/colvardeps.cpp @@ -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); diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 9991513aca..0eef090793 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -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::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::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); } diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index eac071805c..9bc5b4c886 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -4,7 +4,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2016-04-14" +#define COLVARS_VERSION "2016-04-20" #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 #include @@ -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); diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 585a40cdfc..b80b575729 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -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 diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index bd1c8878fb..cabb520248 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -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 -- load a state file (requires configuration)\n\ save -- 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"; diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h index 25138387ca..ccd3ff82c5 100644 --- a/lib/colvars/colvarscript.h +++ b/lib/colvars/colvarscript.h @@ -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 diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index 6368e965f7..0ce7e29a0c 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -10,7 +10,6 @@ #include "lammps.h" #include "domain.h" #include "force.h" -#include "random_park.h" #include "update.h" #include diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index c46b301119..3e86cd13a4 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -29,6 +29,7 @@ #include "group.h" #include "memory.h" #include "modify.h" +#include "random_park.h" #include "respa.h" #include "universe.h" #include "update.h"