diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index ef64989ad..cd0ffe8b6 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -1,280 +1,282 @@ /// -*- c++ -*- #include "colvarmodule.h" #include "colvarvalue.h" #include "colvar.h" #include "colvarcomp.h" colvar::cvc::cvc() : sup_coeff(1.0), sup_np(1), b_periodic(false), b_try_scalable(true) -{} +{ + init_cvc_requires(); +} colvar::cvc::cvc(std::string const &conf) : sup_coeff(1.0), sup_np(1), b_periodic(false), b_try_scalable(true) { if (cvm::debug()) cvm::log("Initializing cvc base object.\n"); init_cvc_requires(); 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); 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); if (b_debug_gradient) enable(f_cvc_debug_gradient); } // Attempt scalable calculations when in parallel? (By default yes, if available) get_keyval(conf, "scalable", b_try_scalable, true); if (cvm::debug()) cvm::log("Done initializing cvc base object.\n"); } cvm::atom_group *colvar::cvc::parse_group(std::string const &conf, char const *group_key, bool optional) { cvm::atom_group *group = NULL; if (key_lookup(conf, group_key)) { group = new cvm::atom_group; group->key = group_key; if (b_try_scalable) { if (is_available(f_cvc_scalable_com) && is_available(f_cvc_com_based)) { enable(f_cvc_scalable_com); enable(f_cvc_scalable); group->enable(f_ag_scalable_com); group->enable(f_ag_scalable); } // TODO check for other types of parallelism here if (is_enabled(f_cvc_scalable)) { cvm::log("Will enable scalable calculation for group \""+group->key+"\".\n"); } else { cvm::log("Scalable calculation is not available for group \""+group->key+"\" with the current configuration.\n"); } } if (group->parse(conf) == COLVARS_OK) { atom_groups.push_back(group); } else { cvm::error("Error parsing definition for atom group \""+ std::string(group_key)+"\".\n"); } } else { if (! optional) { cvm::error("Error: definition for atom group \""+ std::string(group_key)+"\" not found.\n"); } } return group; } int colvar::cvc::setup() { size_t i; description = "cvc " + name; for (i = 0; i < atom_groups.size(); i++) { add_child((cvm::deps *) atom_groups[i]); } if (b_try_scalable && is_available(f_cvc_scalable)) { enable(f_cvc_scalable); } return COLVARS_OK; } colvar::cvc::~cvc() { remove_all_children(); for (size_t i = 0; i < atom_groups.size(); i++) { if (atom_groups[i] != NULL) delete atom_groups[i]; } } void colvar::cvc::read_data() { size_t ig; for (ig = 0; ig < atom_groups.size(); ig++) { cvm::atom_group &atoms = *(atom_groups[ig]); atoms.reset_atoms_data(); atoms.read_positions(); atoms.calc_required_properties(); // each atom group will take care of its own ref_pos_group, if defined } //// Don't try to get atom velocities, as no back-end currently implements it // if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) { // for (i = 0; i < cvcs.size(); i++) { // for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { // cvcs[i]->atom_groups[ig]->read_velocities(); // } // } // } } void colvar::cvc::calc_force_invgrads() { cvm::fatal_error("Error: calculation of inverse gradients is not implemented " "for colvar components of type \""+function_type+"\".\n"); } 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") // 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; cvm::rotation const rot_0 = group->rot; cvm::rotation const rot_inv = group->rot.inverse(); cvm::real x_0 = x.real_value; if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0]; // cvm::log("gradients = "+cvm::to_str (gradients)+"\n"); 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); gradient_sum += atom_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; group->calc_required_properties(); calc_value(); cvm::real x_1 = x.real_value; 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()) ? (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"); } } 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; } // Static members std::vector colvar::cvc::cvc_features; diff --git a/lib/colvars/colvarcomp_protein.cpp b/lib/colvars/colvarcomp_protein.cpp index cd3762928..a12fab1c1 100644 --- a/lib/colvars/colvarcomp_protein.cpp +++ b/lib/colvars/colvarcomp_protein.cpp @@ -1,392 +1,402 @@ /// -*- c++ -*- #include #include "colvarmodule.h" #include "colvarvalue.h" #include "colvarparse.h" #include "colvar.h" #include "colvarcomp.h" ////////////////////////////////////////////////////////////////////// // alpha component ////////////////////////////////////////////////////////////////////// + // FIXME: this will not make collect_gradients work + // because gradients in individual atom groups + // are those of the sub-cvcs (angle, hb), not those + // of this cvc (alpha) + // This is true of all cvcs with sub-cvcs, and those + // that do not calculate explicit gradients + // SO: we need a flag giving the availability of + // atomic gradients + colvar::alpha_angles::alpha_angles(std::string const &conf) : cvc(conf) { if (cvm::debug()) cvm::log("Initializing alpha_angles object.\n"); function_type = "alpha_angles"; x.type(colvarvalue::type_scalar); std::string segment_id; get_keyval(conf, "psfSegID", segment_id, std::string("MAIN")); std::vector residues; { std::string residues_conf = ""; key_lookup(conf, "residueRange", residues_conf); if (residues_conf.size()) { std::istringstream is(residues_conf); int initial, final; char dash; if ( (is >> initial) && (initial > 0) && (is >> dash) && (dash == '-') && (is >> final) && (final > 0) ) { for (int rnum = initial; rnum <= final; rnum++) { residues.push_back(rnum); } } } else { cvm::fatal_error("Error: no residues defined in \"residueRange\".\n"); } } if (residues.size() < 5) { cvm::fatal_error("Error: not enough residues defined in \"residueRange\".\n"); } std::string const &sid = segment_id; std::vector const &r = residues; get_keyval(conf, "hBondCoeff", hb_coeff, 0.5); if ( (hb_coeff < 0.0) || (hb_coeff > 1.0) ) { cvm::fatal_error("Error: hBondCoeff must be defined between 0 and 1.\n"); } get_keyval(conf, "angleRef", theta_ref, 88.0); get_keyval(conf, "angleTol", theta_tol, 15.0); if (hb_coeff < 1.0) { for (size_t i = 0; i < residues.size()-2; i++) { theta.push_back(new colvar::angle(cvm::atom(r[i ], "CA", sid), cvm::atom(r[i+1], "CA", sid), cvm::atom(r[i+2], "CA", sid))); atom_groups.push_back(theta.back()->atom_groups[0]); atom_groups.push_back(theta.back()->atom_groups[1]); atom_groups.push_back(theta.back()->atom_groups[2]); } } else { cvm::log("The hBondCoeff specified will disable the Calpha-Calpha-Calpha angle terms.\n"); } { cvm::real r0; size_t en, ed; get_keyval(conf, "hBondCutoff", r0, (3.3 * cvm::unit_angstrom())); get_keyval(conf, "hBondExpNumer", en, 6); get_keyval(conf, "hBondExpDenom", ed, 8); if (hb_coeff > 0.0) { for (size_t i = 0; i < residues.size()-4; i++) { hb.push_back(new colvar::h_bond(cvm::atom(r[i ], "O", sid), cvm::atom(r[i+4], "N", sid), r0, en, ed)); atom_groups.push_back(hb.back()->atom_groups[0]); } } else { cvm::log("The hBondCoeff specified will disable the hydrogen bond terms.\n"); } } if (cvm::debug()) cvm::log("Done initializing alpha_angles object.\n"); } colvar::alpha_angles::alpha_angles() : cvc() { function_type = "alpha_angles"; x.type(colvarvalue::type_scalar); } colvar::alpha_angles::~alpha_angles() { while (theta.size() != 0) { delete theta.back(); theta.pop_back(); } while (hb.size() != 0) { delete hb.back(); hb.pop_back(); } } void colvar::alpha_angles::calc_value() { x.real_value = 0.0; if (theta.size()) { cvm::real const theta_norm = (1.0-hb_coeff) / cvm::real(theta.size()); for (size_t i = 0; i < theta.size(); i++) { (theta[i])->calc_value(); cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol; cvm::real const f = ( (1.0 - std::pow(t, (int) 2)) / (1.0 - std::pow(t, (int) 4)) ); x.real_value += theta_norm * f; if (cvm::debug()) cvm::log("Calpha-Calpha angle no. "+cvm::to_str(i+1)+" in \""+ this->name+"\" has a value of "+ (cvm::to_str((theta[i])->value().real_value))+ " degrees, f = "+cvm::to_str(f)+".\n"); } } if (hb.size()) { cvm::real const hb_norm = hb_coeff / cvm::real(hb.size()); for (size_t i = 0; i < hb.size(); i++) { (hb[i])->calc_value(); x.real_value += hb_norm * (hb[i])->value().real_value; if (cvm::debug()) cvm::log("Hydrogen bond no. "+cvm::to_str(i+1)+" in \""+ this->name+"\" has a value of "+ (cvm::to_str((hb[i])->value().real_value))+".\n"); } } } void colvar::alpha_angles::calc_gradients() { size_t i; for (i = 0; i < theta.size(); i++) (theta[i])->calc_gradients(); for (i = 0; i < hb.size(); i++) (hb[i])->calc_gradients(); } void colvar::alpha_angles::apply_force(colvarvalue const &force) { if (theta.size()) { cvm::real const theta_norm = (1.0-hb_coeff) / cvm::real(theta.size()); for (size_t i = 0; i < theta.size(); i++) { cvm::real const t = ((theta[i])->value().real_value-theta_ref)/theta_tol; cvm::real const f = ( (1.0 - std::pow(t, (int) 2)) / (1.0 - std::pow(t, (int) 4)) ); cvm::real const dfdt = 1.0/(1.0 - std::pow(t, (int) 4)) * ( (-2.0 * t) + (-1.0*f)*(-4.0 * std::pow(t, (int) 3)) ); (theta[i])->apply_force(theta_norm * dfdt * (1.0/theta_tol) * force.real_value ); } } if (hb.size()) { cvm::real const hb_norm = hb_coeff / cvm::real(hb.size()); for (size_t i = 0; i < hb.size(); i++) { (hb[i])->apply_force(0.5 * hb_norm * force.real_value); } } } ////////////////////////////////////////////////////////////////////// // dihedral principal component ////////////////////////////////////////////////////////////////////// // FIXME: this will not make collect_gradients work // because gradients in individual atom groups - // are those of the sub-cvcs (angle, hb), not those - // of this cvc (alpha) + // are those of the sub-cvcs (dihedral), not those + // of this cvc // This is true of all cvcs with sub-cvcs, and those // that do not calculate explicit gradients // SO: we need a flag giving the availability of // atomic gradients + colvar::dihedPC::dihedPC(std::string const &conf) : cvc(conf) { if (cvm::debug()) cvm::log("Initializing dihedral PC object.\n"); function_type = "dihedPC"; x.type(colvarvalue::type_scalar); std::string segment_id; get_keyval(conf, "psfSegID", segment_id, std::string("MAIN")); std::vector residues; { std::string residues_conf = ""; key_lookup(conf, "residueRange", residues_conf); if (residues_conf.size()) { std::istringstream is(residues_conf); int initial, final; char dash; if ( (is >> initial) && (initial > 0) && (is >> dash) && (dash == '-') && (is >> final) && (final > 0) ) { for (int rnum = initial; rnum <= final; rnum++) { residues.push_back(rnum); } } } else { cvm::fatal_error("Error: no residues defined in \"residueRange\".\n"); } } if (residues.size() < 2) { cvm::fatal_error("Error: dihedralPC requires at least two residues.\n"); } std::string const &sid = segment_id; std::vector const &r = residues; std::string vecFileName; int vecNumber; if (get_keyval(conf, "vectorFile", vecFileName, vecFileName)) { get_keyval(conf, "vectorNumber", vecNumber, 0); if (vecNumber < 1) cvm::fatal_error("A positive value of vectorNumber is required."); std::ifstream vecFile; vecFile.open(vecFileName.c_str()); if (!vecFile.good()) cvm::fatal_error("Error opening dihedral PCA vector file " + vecFileName + " for reading"); // TODO: adapt to different formats by setting this flag bool eigenvectors_as_columns = true; if (eigenvectors_as_columns) { // Carma-style dPCA file std::string line; cvm::real c; while (vecFile.good()) { getline(vecFile, line); if (line.length() < 2) break; std::istringstream ls(line); for (int i=0; i> c; coeffs.push_back(c); } } /* TODO Uncomment this when different formats are recognized else { // Eigenvectors as lines // Skip to the right line for (int i = 1; i> c; coeffs.push_back(c); } } */ vecFile.close(); } else { get_keyval(conf, "vector", coeffs, coeffs); } if ( coeffs.size() != 4 * (residues.size() - 1)) { cvm::fatal_error("Error: wrong number of coefficients: " + cvm::to_str(coeffs.size()) + ". Expected " + cvm::to_str(4 * (residues.size() - 1)) + " (4 coeffs per residue, minus one residue).\n"); } for (size_t i = 0; i < residues.size()-1; i++) { // Psi theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "N", sid), cvm::atom(r[i ], "CA", sid), cvm::atom(r[i ], "C", sid), cvm::atom(r[i+1], "N", sid))); // Phi (next res) theta.push_back(new colvar::dihedral(cvm::atom(r[i ], "C", sid), cvm::atom(r[i+1], "N", sid), cvm::atom(r[i+1], "CA", sid), cvm::atom(r[i+1], "C", sid))); } if (cvm::debug()) cvm::log("Done initializing dihedPC object.\n"); } colvar::dihedPC::dihedPC() : cvc() { function_type = "dihedPC"; x.type(colvarvalue::type_scalar); } colvar::dihedPC::~dihedPC() { while (theta.size() != 0) { delete theta.back(); theta.pop_back(); } } void colvar::dihedPC::calc_value() { x.real_value = 0.0; for (size_t i = 0; i < theta.size(); i++) { theta[i]->calc_value(); cvm::real const t = (PI / 180.) * theta[i]->value().real_value; x.real_value += coeffs[2*i ] * std::cos(t) + coeffs[2*i+1] * std::sin(t); } } void colvar::dihedPC::calc_gradients() { for (size_t i = 0; i < theta.size(); i++) { theta[i]->calc_gradients(); } } void colvar::dihedPC::apply_force(colvarvalue const &force) { for (size_t i = 0; i < theta.size(); i++) { cvm::real const t = (PI / 180.) * theta[i]->value().real_value; cvm::real const dcosdt = - (PI / 180.) * std::sin(t); cvm::real const dsindt = (PI / 180.) * std::cos(t); theta[i]->apply_force((coeffs[2*i ] * dcosdt + coeffs[2*i+1] * dsindt) * force); } } diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index b84f226ec..bdf6fcee0 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -1,1449 +1,1453 @@ /// -*- c++ -*- #include #include #include "colvarmodule.h" #include "colvarparse.h" #include "colvarproxy.h" #include "colvar.h" #include "colvarbias.h" #include "colvarbias_abf.h" #include "colvarbias_alb.h" #include "colvarbias_histogram.h" #include "colvarbias_meta.h" #include "colvarbias_restraint.h" #include "colvarscript.h" colvarmodule::colvarmodule(colvarproxy *proxy_in) { // pointer to the proxy object if (proxy == NULL) { proxy = proxy_in; parse = new colvarparse(); } else { // TODO relax this error to handle multiple molecules in VMD // once the module is not static anymore cvm::error("Error: trying to allocate the collective " "variable module twice.\n"); return; } cvm::log(cvm::line_marker); cvm::log("Initializing the collective variables module, version "+ cvm::to_str(COLVARS_VERSION)+".\n"); cvm::log("Please cite Fiorin et al, Mol Phys 2013 in any publication " "based on this calculation.\n"); + if (proxy->smp_enabled() == COLVARS_OK) { + cvm::log("SMP parallelism is available.\n"); + } + // set initial default values // "it_restart" will be set by the input state file, if any; // "it" should be updated by the proxy colvarmodule::it = colvarmodule::it_restart = 0; colvarmodule::it_restart_from_state_file = true; colvarmodule::use_scripted_forces = false; colvarmodule::b_analysis = false; colvarmodule::debug_gradients_step_size = 1.0e-07; colvarmodule::rotation::monitor_crossings = false; colvarmodule::rotation::crossing_threshold = 1.0e-02; colvarmodule::cv_traj_freq = 100; colvarmodule::restart_out_freq = proxy->restart_frequency(); // by default overwrite the existing trajectory file colvarmodule::cv_traj_append = false; } int colvarmodule::read_config_file(char const *config_filename) { cvm::log(cvm::line_marker); cvm::log("Reading new configuration from file \""+ std::string(config_filename)+"\":\n"); // open the configfile config_s.open(config_filename); if (!config_s.is_open()) { cvm::error("Error: in opening configuration file \""+ std::string(config_filename)+"\".\n", FILE_ERROR); return COLVARS_ERROR; } // read the config file into a string std::string conf = ""; std::string line; while (colvarparse::getline_nocomments(config_s, line)) { conf.append(line+"\n"); } config_s.close(); return parse_config(conf); } int colvarmodule::read_config_string(std::string const &config_str) { cvm::log(cvm::line_marker); cvm::log("Reading new configuration:\n"); std::istringstream config_s(config_str); // strip the comments away std::string conf = ""; std::string line; while (colvarparse::getline_nocomments(config_s, line)) { conf.append(line+"\n"); } return parse_config(conf); } int colvarmodule::parse_config(std::string &conf) { // parse global options if (catch_input_errors(parse_global_params(conf))) { return get_error(); } // parse the options for collective variables if (catch_input_errors(parse_colvars(conf))) { return get_error(); } // parse the options for biases if (catch_input_errors(parse_biases(conf))) { return get_error(); } // done parsing known keywords, check that all keywords found were valid ones if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) { return get_error(); } cvm::log(cvm::line_marker); cvm::log("Collective variables module (re)initialized.\n"); cvm::log(cvm::line_marker); if (cv_traj_os.is_open()) { // configuration might have changed, better redo the labels write_traj_label(cv_traj_os); } return get_error(); } int colvarmodule::parse_global_params(std::string const &conf) { std::string index_file_name; if (parse->get_keyval(conf, "indexFile", index_file_name)) { read_index_file(index_file_name.c_str()); } parse->get_keyval(conf, "analysis", b_analysis, b_analysis); parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size, debug_gradients_step_size, colvarparse::parse_silent); parse->get_keyval(conf, "monitorEigenvalueCrossing", colvarmodule::rotation::monitor_crossings, colvarmodule::rotation::monitor_crossings, colvarparse::parse_silent); parse->get_keyval(conf, "eigenvalueCrossingThreshold", colvarmodule::rotation::crossing_threshold, colvarmodule::rotation::crossing_threshold, colvarparse::parse_silent); parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq); parse->get_keyval(conf, "colvarsRestartFrequency", restart_out_freq, restart_out_freq); // if this is true when initializing, it means // we are continuing after a reset(): default to true parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append); parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false); parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true); if (use_scripted_forces && !proxy->force_script_defined) { cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::parse_colvars(std::string const &conf) { if (cvm::debug()) cvm::log("Initializing the collective variables.\n"); std::string colvar_conf = ""; size_t pos = 0; while (parse->key_lookup(conf, "colvar", colvar_conf, pos)) { if (colvar_conf.size()) { cvm::log(cvm::line_marker); cvm::increase_depth(); colvars.push_back(new colvar(colvar_conf)); if (cvm::get_error() || ((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) { cvm::log("Error while constructing colvar number " + cvm::to_str(colvars.size()) + " : deleting."); delete colvars.back(); // the colvar destructor updates the colvars array return COLVARS_ERROR; } cvm::decrease_depth(); } else { cvm::error("Error: \"colvar\" keyword found without any configuration.\n", INPUT_ERROR); return COLVARS_ERROR; } cvm::decrease_depth(); colvar_conf = ""; } if (!colvars.size()) { cvm::log("Warning: no collective variables defined.\n"); } if (colvars.size()) cvm::log(cvm::line_marker); cvm::log("Collective variables initialized, "+ cvm::to_str(colvars.size())+ " in total.\n"); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } bool colvarmodule::check_new_bias(std::string &conf, char const *key) { if (cvm::get_error() || (biases.back()->check_keywords(conf, key) != COLVARS_OK)) { cvm::log("Error while constructing bias number " + cvm::to_str(biases.size()) + " : deleting.\n"); delete biases.back(); // the bias destructor updates the biases array return true; } return false; } template int colvarmodule::parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count) { std::string bias_conf = ""; size_t conf_saved_pos = 0; while (parse->key_lookup(conf, keyword, bias_conf, conf_saved_pos)) { if (bias_conf.size()) { cvm::log(cvm::line_marker); cvm::increase_depth(); biases.push_back(new bias_type(bias_conf, keyword)); if (cvm::check_new_bias(bias_conf, keyword)) { return COLVARS_ERROR; } cvm::decrease_depth(); bias_count++; } else { cvm::error("Error: keyword \""+std::string(keyword)+"\" found without configuration.\n", INPUT_ERROR); return COLVARS_ERROR; } bias_conf = ""; } return COLVARS_OK; } int colvarmodule::parse_biases(std::string const &conf) { if (cvm::debug()) cvm::log("Initializing the collective variables biases.\n"); /// initialize ABF instances parse_biases_type(conf, "abf", n_abf_biases); /// initialize adaptive linear biases parse_biases_type(conf, "ALB", n_rest_biases); /// initialize harmonic restraints parse_biases_type(conf, "harmonic", n_rest_biases); /// initialize histograms parse_biases_type(conf, "histogram", n_histo_biases); /// initialize linear restraints parse_biases_type(conf, "linear", n_rest_biases); /// initialize metadynamics instances parse_biases_type(conf, "metadynamics", n_meta_biases); if (use_scripted_forces) { cvm::log(cvm::line_marker); cvm::increase_depth(); cvm::log("User forces script will be run at each bias update."); cvm::decrease_depth(); } for (size_t i = 0; i < biases.size(); i++) { biases[i]->enable(cvm::deps::f_cvb_active); if (cvm::debug()) biases[i]->print_state(); } if (biases.size() || use_scripted_forces) { cvm::log(cvm::line_marker); cvm::log("Collective variables biases initialized, "+ cvm::to_str(biases.size())+" in total.\n"); } else { if (!use_scripted_forces) { cvm::log("No collective variables biases were defined.\n"); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::catch_input_errors(int result) { if (result != COLVARS_OK || get_error()) { set_error_bit(result); set_error_bit(INPUT_ERROR); parse->reset(); return get_error(); } return COLVARS_OK; } colvarbias * colvarmodule::bias_by_name(std::string const &name) { for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if ((*bi)->name == name) { return (*bi); } } return NULL; } colvar *colvarmodule::colvar_by_name(std::string const &name) { for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ((*cvi)->name == name) { return (*cvi); } } return NULL; } int colvarmodule::change_configuration(std::string const &bias_name, std::string const &conf) { // This is deprecated; supported strategy is to delete the bias // and parse the new config cvm::increase_depth(); colvarbias *b; b = bias_by_name(bias_name); if (b == NULL) { cvm::error("Error: bias not found: " + bias_name); } b->change_configuration(conf); cvm::decrease_depth(); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::string colvarmodule::read_colvar(std::string const &name) { cvm::increase_depth(); colvar *c; std::stringstream ss; c = colvar_by_name(name); if (c == NULL) { cvm::fatal_error("Error: colvar not found: " + name); } ss << c->value(); cvm::decrease_depth(); return ss.str(); } cvm::real colvarmodule::energy_difference(std::string const &bias_name, std::string const &conf) { cvm::increase_depth(); colvarbias *b; cvm::real energy_diff = 0.; b = bias_by_name(bias_name); if (b == NULL) { cvm::fatal_error("Error: bias not found: " + bias_name); } energy_diff = b->energy_difference(conf); cvm::decrease_depth(); return energy_diff; } int colvarmodule::bias_current_bin(std::string const &bias_name) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { ret = b->current_bin(); } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::bias_bin_num(std::string const &bias_name) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { ret = b->bin_num(); } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::bias_bin_count(std::string const &bias_name, size_t bin_index) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { ret = b->bin_count(bin_index); } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::bias_share(std::string const &bias_name) { cvm::increase_depth(); int ret; colvarbias *b = bias_by_name(bias_name); if (b != NULL) { b->replica_share(); ret = COLVARS_OK; } else { cvm::error("Error: bias not found.\n"); ret = COLVARS_ERROR; } cvm::decrease_depth(); return ret; } int colvarmodule::calc() { int error_code = COLVARS_OK; if (cvm::debug()) { cvm::log(cvm::line_marker); cvm::log("Collective variables module, step no. "+ cvm::to_str(cvm::step_absolute())+"\n"); } 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(); } combine_errors(error_code, calc_biases()); combine_errors(error_code, update_colvar_forces()); if (cvm::b_analysis) { combine_errors(error_code, analyze()); } // write trajectory files, if needed if (cv_traj_freq && cv_traj_name.size()) { combine_errors(error_code, write_traj_files()); } // write restart files, if needed if (restart_out_freq && restart_out_name.size()) { combine_errors(error_code, write_restart_files()); } return error_code; } int colvarmodule::calc_colvars() { if (cvm::debug()) cvm::log("Calculating collective variables.\n"); // calculate collective variables and their gradients int error_code = COLVARS_OK; std::vector::iterator cvi; // if SMP support is available, split up the work if (proxy->smp_enabled() == COLVARS_OK) { // first, calculate how much work (currently, how many active CVCs) each colvar has colvars_smp.resize(0); colvars_smp_items.resize(0); colvars_smp.reserve(colvars.size()); colvars_smp_items.reserve(colvars.size()); // set up a vector containing all components size_t num_colvar_items = 0; cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { combine_errors(error_code, (*cvi)->update_cvc_flags()); size_t num_items = (*cvi)->num_active_cvcs(); colvars_smp.reserve(colvars_smp.size() + num_items); colvars_smp_items.reserve(colvars_smp_items.size() + num_items); for (size_t icvc = 0; icvc < num_items; icvc++) { colvars_smp.push_back(*cvi); colvars_smp_items.push_back(icvc); } num_colvar_items += num_items; } cvm::decrease_depth(); // calculate colvar components in parallel combine_errors(error_code, proxy->smp_colvars_loop()); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { combine_errors(error_code, (*cvi)->collect_cvc_data()); } cvm::decrease_depth(); } else { // calculate colvars one at a time cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { combine_errors(error_code, (*cvi)->calc()); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); } combine_errors(error_code, cvm::get_error()); return error_code; } int colvarmodule::calc_biases() { // update the biases and communicate their forces to the collective // variables if (cvm::debug() && biases.size()) cvm::log("Updating collective variable biases.\n"); std::vector::iterator bi; int error_code = COLVARS_OK; // if SMP support is available, split up the work if (proxy->smp_enabled() == COLVARS_OK) { if (use_scripted_forces && !scripting_after_biases) { // calculate biases and scripted forces in parallel combine_errors(error_code, proxy->smp_biases_script_loop()); } else { // calculate biases in parallel combine_errors(error_code, proxy->smp_biases_loop()); } } else { if (use_scripted_forces && !scripting_after_biases) { combine_errors(error_code, calc_scripted_forces()); } cvm::increase_depth(); for (bi = biases.begin(); bi != biases.end(); bi++) { combine_errors(error_code, (*bi)->update()); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); } cvm::real total_bias_energy = 0.0; for (bi = biases.begin(); bi != biases.end(); bi++) { total_bias_energy += (*bi)->get_energy(); } proxy->add_energy(total_bias_energy); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::update_colvar_forces() { int error_code = COLVARS_OK; std::vector::iterator cvi; std::vector::iterator bi; // sum the forces from all biases for each collective variable if (cvm::debug() && biases.size()) cvm::log("Collecting forces from all biases.\n"); cvm::increase_depth(); for (bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); if (use_scripted_forces && scripting_after_biases) { combine_errors(error_code, calc_scripted_forces()); } cvm::real total_colvar_energy = 0.0; // sum up the forces for each colvar, including wall forces // and integrate any internal // equation of motion (extended system) if (cvm::debug()) cvm::log("Updating the internal degrees of freedom " "of colvars (if they have any).\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { total_colvar_energy += (*cvi)->update_forces_energy(); if (cvm::get_error()) { return COLVARS_ERROR; } } cvm::decrease_depth(); proxy->add_energy(total_colvar_energy); // make collective variables communicate their forces to their // coupled degrees of freedom (i.e. atoms) if (cvm::debug()) cvm::log("Communicating forces from the colvars to the atoms.\n"); cvm::increase_depth(); for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ((*cvi)->is_enabled(cvm::deps::f_cv_gradient)) { (*cvi)->communicate_forces(); if (cvm::get_error()) { return COLVARS_ERROR; } } } cvm::decrease_depth(); return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::calc_scripted_forces() { // Run user force script, if provided, // potentially adding scripted forces to the colvars int res; res = proxy->run_force_callback(); if (res == COLVARS_NOT_IMPLEMENTED) { cvm::error("Colvar forces scripts are not implemented."); return COLVARS_NOT_IMPLEMENTED; } if (res != COLVARS_OK) { cvm::error("Error running user colvar forces script"); return COLVARS_ERROR; } return COLVARS_OK; } int colvarmodule::write_restart_files() { if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log("Writing the state file \""+ restart_out_name+"\".\n"); proxy->backup_file(restart_out_name.c_str()); restart_out_os.open(restart_out_name.c_str()); if (!restart_out_os.is_open() || !write_restart(restart_out_os)) cvm::error("Error: in writing restart file.\n"); restart_out_os.close(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::write_traj_files() { if (!cv_traj_os.is_open()) { open_traj_file(cv_traj_name); } // write labels in the traj file every 1000 lines and at first timestep if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) { write_traj_label(cv_traj_os); } if ((cvm::step_absolute() % cv_traj_freq) == 0) { write_traj(cv_traj_os); } if (restart_out_freq && cv_traj_os.is_open()) { // flush the trajectory file if we are at the restart frequency if ( (cvm::step_relative() > 0) && ((cvm::step_absolute() % restart_out_freq) == 0) ) { cvm::log("Synchronizing (emptying the buffer of) trajectory file \""+ cv_traj_name+"\".\n"); cv_traj_os.flush(); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::analyze() { if (cvm::debug()) { cvm::log("colvarmodule::analyze(), step = "+cvm::to_str(it)+".\n"); } if (cvm::step_relative() == 0) cvm::log("Performing analysis.\n"); // perform colvar-specific analysis for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { cvm::increase_depth(); (*cvi)->analyze(); cvm::decrease_depth(); } // perform bias-specific analysis for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { cvm::increase_depth(); (*bi)->analyze(); cvm::decrease_depth(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup() { // loop over all components of all colvars to reset masses of all groups for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->setup(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } colvarmodule::~colvarmodule() { if ((proxy->smp_thread_id() == COLVARS_NOT_IMPLEMENTED) || (proxy->smp_thread_id() == 0)) { reset(); delete parse; proxy = NULL; } } int colvarmodule::reset() { parse->reset(); cvm::log("Resetting the Collective Variables Module.\n"); // Iterate backwards because we are deleting the elements as we go for (std::vector::reverse_iterator bi = biases.rbegin(); bi != biases.rend(); bi++) { delete *bi; // the bias destructor updates the biases array } biases.clear(); // Iterate backwards because we are deleting the elements as we go for (std::vector::reverse_iterator cvi = colvars.rbegin(); cvi != colvars.rend(); cvi++) { delete *cvi; // the colvar destructor updates the colvars array } colvars.clear(); index_groups.clear(); index_group_names.clear(); if (cv_traj_os.is_open()) { // Do not close file here, as we might not be done with it yet. cv_traj_os.flush(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup_input() { // name of input state file restart_in_name = proxy->input_prefix().size() ? std::string(proxy->input_prefix()+".colvars.state") : std::string("") ; // read the restart configuration, if available if (restart_in_name.size()) { // read the restart file std::ifstream input_is(restart_in_name.c_str()); if (!input_is.good()) { cvm::error("Error: in opening restart file \""+ std::string(restart_in_name)+"\".\n", FILE_ERROR); return COLVARS_ERROR; } else { cvm::log("Restarting from file \""+restart_in_name+"\".\n"); read_restart(input_is); if (cvm::get_error() != COLVARS_OK) { return COLVARS_ERROR; } cvm::log(cvm::line_marker); } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::setup_output() { int error_code = 0; // output state file (restart) restart_out_name = proxy->restart_output_prefix().size() ? std::string(proxy->restart_output_prefix()+".colvars.state") : std::string(""); if (restart_out_name.size()) { cvm::log("The restart output state file will be \""+restart_out_name+"\".\n"); } output_prefix = proxy->output_prefix(); if (output_prefix.size()) { cvm::log("The final output state file will be \""+ (output_prefix.size() ? std::string(output_prefix+".colvars.state") : std::string("colvars.state"))+"\".\n"); // cvm::log (cvm::line_marker); } cv_traj_name = (output_prefix.size() ? std::string(output_prefix+".colvars.traj") : std::string("")); if (cv_traj_freq && cv_traj_name.size()) { combine_errors(error_code, open_traj_file(cv_traj_name)); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { combine_errors(error_code, (*bi)->setup_output()); } if (error_code != COLVARS_OK || cvm::get_error()) { set_error_bit(FILE_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::istream & colvarmodule::read_restart(std::istream &is) { { // read global restart information std::string restart_conf; if (is >> colvarparse::read_block("configuration", restart_conf)) { if (it_restart_from_state_file) { parse->get_keyval(restart_conf, "step", it_restart, (size_t) 0, 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(); } // colvars restart cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if ( !((*cvi)->read_restart(is)) ) { cvm::error("Error: in reading restart configuration for collective variable \""+ (*cvi)->name+"\".\n", INPUT_ERROR); } } // biases restart for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { if (!((*bi)->read_restart(is))) { cvm::error("Error: in reading restart configuration for bias \""+ (*bi)->name+"\".\n", INPUT_ERROR); } } cvm::decrease_depth(); return is; } int colvarmodule::backup_file(char const *filename) { return proxy->backup_file(filename); } int colvarmodule::write_output_files() { // if this is a simulation run (i.e. not a postprocessing), output data // must be written to be able to restart the simulation std::string const out_name = (output_prefix.size() ? std::string(output_prefix+".colvars.state") : std::string("colvars.state")); cvm::log("Saving collective variables state to \""+out_name+"\".\n"); std::ostream * os = proxy->output_stream(out_name); os->setf(std::ios::scientific, std::ios::floatfield); this->write_restart(*os); proxy->close_output_stream(out_name); cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_output_files(); } cvm::decrease_depth(); cvm::increase_depth(); for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_output_files(); } cvm::decrease_depth(); if (cv_traj_os.is_open()) { // do not close to avoid problems with multiple NAMD runs cv_traj_os.flush(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::read_traj(char const *traj_filename, long traj_read_begin, long traj_read_end) { cvm::log("Opening trajectory file \""+ std::string(traj_filename)+"\".\n"); std::ifstream traj_is(traj_filename); while (true) { while (true) { std::string line(""); do { if (!colvarparse::getline_nocomments(traj_is, line)) { cvm::log("End of file \""+std::string(traj_filename)+ "\" reached, or corrupted file.\n"); traj_is.close(); return false; } } while (line.find_first_not_of(colvarparse::white_space) == std::string::npos); std::istringstream is(line); if (!(is >> it)) return false; if ( (it < traj_read_begin) ) { if ((it % 1000) == 0) std::cerr << "Skipping trajectory step " << it << " \r"; continue; } else { if ((it % 1000) == 0) std::cerr << "Reading from trajectory, step = " << it << " \r"; if ( (traj_read_end > traj_read_begin) && (it > traj_read_end) ) { std::cerr << "\n"; cvm::error("Reached the end of the trajectory, " "read_end = "+cvm::to_str(traj_read_end)+"\n", FILE_ERROR); return COLVARS_ERROR; } for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { if (!(*cvi)->read_traj(is)) { cvm::error("Error: in reading colvar \""+(*cvi)->name+ "\" from trajectory file \""+ std::string(traj_filename)+"\".\n", FILE_ERROR); return COLVARS_ERROR; } } break; } } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::ostream & colvarmodule::write_restart(std::ostream &os) { os.setf(std::ios::scientific, std::ios::floatfield); os << "configuration {\n" << " step " << std::setw(it_width) << it << "\n" << " dt " << dt() << "\n" << " version " << std::string(COLVARS_VERSION) << "\n" << "}\n\n"; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_restart(os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_restart(os); } cvm::decrease_depth(); return os; } int colvarmodule::open_traj_file(std::string const &file_name) { if (cv_traj_os.is_open()) { return COLVARS_OK; } // (re)open trajectory file if (cv_traj_append) { cvm::log("Appending to colvar trajectory file \""+file_name+ "\".\n"); cv_traj_os.open(file_name.c_str(), std::ios::app); } else { cvm::log("Writing to colvar trajectory file \""+file_name+ "\".\n"); proxy->backup_file(file_name.c_str()); cv_traj_os.open(file_name.c_str()); } if (!cv_traj_os.is_open()) { cvm::error("Error: cannot write to file \""+file_name+"\".\n", FILE_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int colvarmodule::close_traj_file() { if (cv_traj_os.is_open()) { cv_traj_os.close(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } std::ostream & colvarmodule::write_traj_label(std::ostream &os) { if (!os.good()) { cvm::error("Cannot write to trajectory file."); return os; } os.setf(std::ios::scientific, std::ios::floatfield); os << "# " << cvm::wrap_string("step", cvm::it_width-2) << " "; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj_label(os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_traj_label(os); } os << "\n"; if (cvm::debug()) { os.flush(); } cvm::decrease_depth(); return os; } std::ostream & colvarmodule::write_traj(std::ostream &os) { os.setf(std::ios::scientific, std::ios::floatfield); os << std::setw(cvm::it_width) << it << " "; cvm::increase_depth(); for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { (*cvi)->write_traj(os); } for (std::vector::iterator bi = biases.begin(); bi != biases.end(); bi++) { (*bi)->write_traj(os); } os << "\n"; if (cvm::debug()) { os.flush(); } cvm::decrease_depth(); return os; } void cvm::log(std::string const &message) { size_t const d = depth(); if (d > 0) proxy->log((std::string(2*d, ' '))+message); else proxy->log(message); } void cvm::increase_depth() { (depth())++; } void cvm::decrease_depth() { if (depth() > 0) { (depth())--; } } size_t & cvm::depth() { // NOTE: do not call log() or error() here, to avoid recursion size_t const nt = proxy->smp_num_threads(); if (proxy->smp_enabled() == COLVARS_OK) { if (depth_v.size() != nt) { // update array of depths proxy->smp_lock(); if (depth_v.size() > 0) { depth_s = depth_v[0]; } depth_v.clear(); depth_v.assign(nt, depth_s); proxy->smp_unlock(); } return depth_v[proxy->smp_thread_id()]; } return depth_s; } 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 = -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() { proxy->smp_lock(); errorCode = COLVARS_OK; proxy->smp_unlock(); } void cvm::error(std::string const &message, int code) { set_error_bit(code); proxy->error(message); } 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_bit(FATAL_ERROR); proxy->fatal_error(message); } void cvm::exit(std::string const &message) { proxy->exit(message); } int cvm::read_index_file(char const *filename) { std::ifstream is(filename, std::ios::binary); if (!is.good()) { cvm::error("Error: in opening index file \""+ std::string(filename)+"\".\n", FILE_ERROR); } while (is.good()) { char open, close; std::string group_name; if ( (is >> open) && (open == '[') && (is >> group_name) && (is >> close) && (close == ']') ) { for (std::list::iterator names_i = index_group_names.begin(); names_i != index_group_names.end(); names_i++) { if (*names_i == group_name) { cvm::error("Error: the group name \""+group_name+ "\" appears in multiple index files.\n", FILE_ERROR); } } cvm::index_group_names.push_back(group_name); cvm::index_groups.push_back(std::vector ()); } else { cvm::error("Error: in parsing index file \""+ std::string(filename)+"\".\n", INPUT_ERROR); } int atom_number = 1; size_t pos = is.tellg(); while ( (is >> atom_number) && (atom_number > 0) ) { (cvm::index_groups.back()).push_back(atom_number); pos = is.tellg(); } is.clear(); is.seekg(pos, std::ios::beg); std::string delim; if ( (is >> delim) && (delim == "[") ) { // new group is.clear(); is.seekg(pos, std::ios::beg); } else { break; } } cvm::log("The following index groups were read from the index file \""+ std::string(filename)+"\":\n"); std::list::iterator names_i = index_group_names.begin(); std::list >::iterator lists_i = index_groups.begin(); for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) { cvm::log(" "+(*names_i)+" ("+cvm::to_str(lists_i->size())+" atoms).\n"); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } int cvm::load_atoms(char const *file_name, cvm::atom_group &atoms, std::string const &pdb_field, double const pdb_field_value) { return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value); } int cvm::load_coords(char const *file_name, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value) { // Differentiate between PDB and XYZ files // for XYZ files, use CVM internal parser // otherwise call proxy function for PDB std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name); if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) { if ( pdb_field.size() > 0 ) { cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR); return COLVARS_ERROR; } return cvm::load_coords_xyz(file_name, pos, indices); } else { return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value); } } int cvm::load_coords_xyz(char const *filename, std::vector &pos, const std::vector &indices) { std::ifstream xyz_is(filename); unsigned int natoms; char symbol[256]; std::string line; if ( ! (xyz_is >> natoms) ) { cvm::error("Error: cannot parse XYZ file " + std::string(filename) + ".\n", INPUT_ERROR); } // skip comment line std::getline(xyz_is, line); std::getline(xyz_is, line); xyz_is.width(255); std::vector::iterator pos_i = pos.begin(); if (pos.size() != natoms) { // Use specified indices int next = 0; // indices are zero-based std::vector::const_iterator index = indices.begin(); for ( ; pos_i != pos.end() ; pos_i++, index++) { while ( next < *index ) { std::getline(xyz_is, line); next++; } xyz_is >> symbol; xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; } } else { // Use all positions for ( ; pos_i != pos.end() ; pos_i++) { xyz_is >> symbol; xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2]; } } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } // static pointers std::vector colvarmodule::colvars; std::vector colvarmodule::biases; size_t colvarmodule::n_abf_biases = 0; size_t colvarmodule::n_rest_biases = 0; size_t colvarmodule::n_histo_biases = 0; size_t colvarmodule::n_meta_biases = 0; colvarproxy *colvarmodule::proxy = NULL; // static runtime data cvm::real colvarmodule::debug_gradients_step_size = 1.0e-07; int colvarmodule::errorCode = 0; long colvarmodule::it = 0; long colvarmodule::it_restart = 0; size_t colvarmodule::restart_out_freq = 0; size_t colvarmodule::cv_traj_freq = 0; size_t colvarmodule::depth_s = 0; std::vector colvarmodule::depth_v(0); bool colvarmodule::b_analysis = false; std::list colvarmodule::index_group_names; std::list > colvarmodule::index_groups; bool colvarmodule::use_scripted_forces = false; bool colvarmodule::scripting_after_biases = true; // file name prefixes std::string colvarmodule::output_prefix = ""; std::string colvarmodule::restart_in_name = ""; // i/o constants size_t const colvarmodule::it_width = 12; size_t const colvarmodule::cv_prec = 14; size_t const colvarmodule::cv_width = 21; size_t const colvarmodule::en_prec = 14; size_t const colvarmodule::en_width = 21; std::string const colvarmodule::line_marker = "----------------------------------------------------------------------\n"; diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index 3649d053d..785027884 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -1,676 +1,676 @@ /// -*- c++ -*- #ifndef COLVARMODULE_H #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2016-04-28" +#define COLVARS_VERSION "2016-05-10" #endif #ifndef COLVARS_DEBUG #define COLVARS_DEBUG false #endif /// \file colvarmodule.h /// \brief Collective variables main module /// /// This file declares the main class for defining and manipulating /// collective variables: there should be only one instance of this /// class, because several variables are made static (i.e. they are /// shared between all object instances) to be accessed from other /// objects. // 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 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 #include #include #include #include #include #include #include #ifdef NAMD_VERSION // use Lustre-friendly wrapper to POSIX write() #include "fstream_namd.h" #endif class colvarparse; class colvar; class colvarbias; class colvarproxy; class colvarscript; /// \brief Collective variables module (main class) /// /// Class to control the collective variables calculation. An object /// (usually one) of this class is spawned from the MD program, /// containing all i/o routines and general interface. /// /// At initialization, the colvarmodule object creates a proxy object /// to provide a transparent interface between the MD program and the /// child objects class colvarmodule { private: /// Impossible to initialize the main object without arguments colvarmodule(); public: /// Base class to handle mutual dependencies of most objects class deps; friend class colvarproxy; // TODO colvarscript should be unaware of colvarmodule's internals friend class colvarscript; /// Defining an abstract real number allows to switch precision typedef double real; /// Residue identifier typedef int residue_id; class rvector; template class vector1d; template class matrix2d; class quaternion; class rotation; /// \brief Atom position (different type name from rvector, to make /// possible future PBC-transparent implementations) typedef rvector atom_pos; /// \brief 3x3 matrix of real numbers class rmatrix; // allow these classes to access protected data class atom; class atom_group; friend class atom; friend class atom_group; typedef std::vector::iterator atom_iter; typedef std::vector::const_iterator atom_const_iter; /// Module-wide error state /// see constants at the top of this file protected: static int errorCode; public: 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() { return errorCode; } static void clear_error(); /// Current step number static long it; /// Starting step number for this run static long it_restart; /// Return the current step number from the beginning of this run static inline long step_relative() { return it - it_restart; } /// Return the current step number from the beginning of the whole /// calculation static inline long step_absolute() { return it; } /// If true, get it_restart from the state file; if set to false, /// the MD program is providing it bool it_restart_from_state_file; /// \brief Finite difference step size (if there is no dynamics, or /// if gradients need to be tested independently from the size of /// dt) static real debug_gradients_step_size; /// Prefix for all output files for this run static std::string output_prefix; /// input restart file name (determined from input_prefix) static std::string restart_in_name; /// Array of collective variables static std::vector colvars; /* TODO: implement named CVCs /// Array of named (reusable) collective variable components static std::vector cvcs; /// Named cvcs register themselves at initialization time inline void register_cvc(cvc *p) { cvcs.push_back(p); } */ /// Collective variables to be calculated on different threads; /// colvars with multple items (e.g. multiple active CVCs) are duplicated std::vector colvars_smp; /// Indexes of the items to calculate for each colvar std::vector colvars_smp_items; /// Array of collective variable biases static std::vector biases; /// \brief Number of ABF biases initialized (in normal conditions /// should be 1) static size_t n_abf_biases; /// \brief Number of metadynamics biases initialized (in normal /// conditions should be 1) static size_t n_meta_biases; /// \brief Number of restraint biases initialized (no limit on the /// number) static size_t n_rest_biases; /// \brief Number of histograms initialized (no limit on the /// number) static size_t n_histo_biases; /// \brief Whether debug output should be enabled (compile-time option) static inline bool debug() { return COLVARS_DEBUG; } /// \brief Constructor \param config_name Configuration file name /// \param restart_name (optional) Restart file name colvarmodule(colvarproxy *proxy); /// Destructor ~colvarmodule(); /// Actual function called by the destructor int reset(); /// Open a config file, load its contents, and pass it to config_string() int read_config_file(char const *config_file_name); /// \brief Parse a config string assuming it is a complete configuration /// (i.e. calling all parse functions) int read_config_string(std::string const &conf); /// \brief Parse a "clean" config string (no comments) int parse_config(std::string &conf); // Parse functions (setup internal data based on a string) /// Parse the few module's global parameters int parse_global_params(std::string const &conf); /// Parse and initialize collective variables int parse_colvars(std::string const &conf); /// Parse and initialize collective variable biases int parse_biases(std::string const &conf); /// Parse and initialize collective variable biases of a specific type template int parse_biases_type(std::string const &conf, char const *keyword, size_t &bias_count); /// Test error condition and keyword parsing /// on error, delete new bias bool check_new_bias(std::string &conf, char const *key); private: /// Useful wrapper to interrupt parsing if any error occurs int catch_input_errors(int result); public: // "Setup" functions (change internal data based on related data // from the proxy that may change during program execution) // No additional parsing is done within these functions /// (Re)initialize internal data (currently used by LAMMPS) /// Also calls setup() member functions of colvars and biases int setup(); /// (Re)initialize and (re)read the input state file calling read_restart() int setup_input(); /// (Re)initialize the output trajectory and state file (does not write it yet) int setup_output(); #ifdef NAMD_VERSION typedef ofstream_namd ofstream; #else typedef std::ofstream ofstream; #endif /// Read the input restart file std::istream & read_restart(std::istream &is); /// Write the output restart file std::ostream & write_restart(std::ostream &os); /// Open a trajectory file if requested (and leave it open) int open_traj_file(std::string const &file_name); /// Close it int close_traj_file(); /// Write in the trajectory file std::ostream & write_traj(std::ostream &os); /// Write explanatory labels in the trajectory file std::ostream & write_traj_label(std::ostream &os); /// Write all trajectory files int write_traj_files(); /// Write all restart files int write_restart_files(); /// Write all FINAL output files int write_output_files(); /// Backup a file before writing it static int backup_file(char const *filename); /// Look up a bias by name; returns NULL if not found static colvarbias * bias_by_name(std::string const &name); /// Look up a colvar by name; returns NULL if not found static colvar * colvar_by_name(std::string const &name); /// Load new configuration for the given bias - /// currently works for harmonic (force constant and/or centers) int change_configuration(std::string const &bias_name, std::string const &conf); /// Read a colvar value std::string read_colvar(std::string const &name); /// Calculate change in energy from using alt. config. for the given bias - /// currently works for harmonic (force constant and/or centers) real energy_difference(std::string const &bias_name, std::string const &conf); /// Give the total number of bins for a given bias. int bias_bin_num(std::string const &bias_name); /// Calculate the bin index for a given bias. int bias_current_bin(std::string const &bias_name); //// Give the count at a given bin index. int bias_bin_count(std::string const &bias_name, size_t bin_index); //// Share among replicas. int bias_share(std::string const &bias_name); /// Main worker function int calc(); /// Calculate collective variables int calc_colvars(); /// Calculate biases int calc_biases(); /// Integrate bias and restraint forces, send colvar forces to atoms int update_colvar_forces(); /// Perform analysis int analyze(); /// \brief Read a collective variable trajectory (post-processing /// only, not called at runtime) int read_traj(char const *traj_filename, long traj_read_begin, long traj_read_end); /// Quick conversion of an object to a string template static std::string to_str(T const &x, size_t const &width = 0, size_t const &prec = 0); /// Quick conversion of a vector of objects to a string template static std::string to_str(std::vector const &x, size_t const &width = 0, size_t const &prec = 0); /// Reduce the number of characters in a string static inline std::string wrap_string(std::string const &s, size_t const &nchars) { if (!s.size()) return std::string(nchars, ' '); else return ( (s.size() <= size_t(nchars)) ? (s+std::string(nchars-s.size(), ' ')) : (std::string(s, 0, nchars)) ); } /// Number of characters to represent a time step static size_t const it_width; /// Number of digits to represent a collective variables value(s) static size_t const cv_prec; /// Number of characters to represent a collective variables value(s) static size_t const cv_width; /// Number of digits to represent the collective variables energy static size_t const en_prec; /// Number of characters to represent the collective variables energy static size_t const en_width; /// Line separator in the log output static std::string const line_marker; // proxy functions /// \brief Value of the unit for atomic coordinates with respect to /// angstroms (used by some variables for hard-coded default values) static real unit_angstrom(); /// \brief Boltmann constant static real boltzmann(); /// \brief Temperature of the simulation (K) static real temperature(); /// \brief Time step of MD integrator (fs) static real dt(); /// Request calculation of system force from MD engine static void request_system_force(); /// Print a message to the main log static void log(std::string const &message); /// Print a message to the main log and exit with error code static void fatal_error(std::string const &message); /// Print a message to the main log and set global error code 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); // Replica exchange commands. static bool replica_enabled(); static int replica_index(); static int replica_num(); static void replica_comm_barrier(); static int replica_comm_recv(char* msg_data, int buf_len, int src_rep); static int replica_comm_send(char* msg_data, int msg_len, int dest_rep); /// \brief Get the distance between two atomic positions with pbcs handled /// correctly static rvector position_distance(atom_pos const &pos1, atom_pos const &pos2); /// \brief Get the square distance between two positions (with /// periodic boundary conditions handled transparently) /// /// Note: in the case of periodic boundary conditions, this provides /// an analytical square distance (while taking the square of /// position_distance() would produce leads to a cusp) static real position_dist2(atom_pos const &pos1, atom_pos const &pos2); /// \brief Get the closest periodic image to a reference position /// \param pos The position to look for the closest periodic image /// \param ref_pos (optional) The reference position static void select_closest_image(atom_pos &pos, atom_pos const &ref_pos); /// \brief Perform select_closest_image() on a set of atomic positions /// /// After that, distance vectors can then be calculated directly, /// without using position_distance() static void select_closest_images(std::vector &pos, atom_pos const &ref_pos); /// \brief Names of groups from a Gromacs .ndx file to be read at startup static std::list index_group_names; /// \brief Groups from a Gromacs .ndx file read at startup static std::list > index_groups; /// \brief Read a Gromacs .ndx file static int read_index_file(char const *filename); /// \brief Create atoms from a file \param filename name of the file /// (usually a PDB) \param atoms array of the atoms to be allocated /// \param pdb_field (optiona) if "filename" is a PDB file, use this /// field to determine which are the atoms to be set static int load_atoms(char const *filename, atom_group &atoms, std::string const &pdb_field, double const pdb_field_value = 0.0); /// \brief Load the coordinates for a group of atoms from a file /// (PDB or XYZ) static int load_coords(char const *filename, std::vector &pos, const std::vector &indices, std::string const &pdb_field, double const pdb_field_value = 0.0); /// \brief Load the coordinates for a group of atoms from an /// XYZ file static int load_coords_xyz(char const *filename, std::vector &pos, const std::vector &indices); /// Frequency for collective variables trajectory output static size_t cv_traj_freq; /// \brief True if only analysis is performed and not a run static bool b_analysis; /// Frequency for saving output restarts static size_t restart_out_freq; /// Output restart file name std::string restart_out_name; /// Pseudo-random number with Gaussian distribution static real rand_gaussian(void); protected: /// Configuration file std::ifstream config_s; /// Configuration file parser object colvarparse *parse; /// Name of the trajectory file std::string cv_traj_name; /// Collective variables output trajectory file colvarmodule::ofstream cv_traj_os; /// Appending to the existing trajectory file? bool cv_traj_append; /// Output restart file colvarmodule::ofstream restart_out_os; protected: /// Counter for the current depth in the object hierarchy (useg e.g. in output) static size_t depth_s; /// Thread-specific depth static std::vector depth_v; public: /// Get the current object depth in the hierarchy static size_t & depth(); /// Increase the depth (number of indentations in the output) static void increase_depth(); /// Decrease the depth (number of indentations in the output) static void decrease_depth(); static inline bool scripted_forces() { return use_scripted_forces; } /// Use scripted colvars forces? static bool use_scripted_forces; /// Wait for all biases before calculating scripted forces? static bool scripting_after_biases; /// Calculate the energy and forces of scripted biases int calc_scripted_forces(); /// \brief Pointer to the proxy object, used to retrieve atomic data /// from the hosting program; it is static in order to be accessible /// from static functions in the colvarmodule class static colvarproxy *proxy; }; /// Shorthand for the frequently used type prefix typedef colvarmodule cvm; #include "colvartypes.h" std::ostream & operator << (std::ostream &os, cvm::rvector const &v); std::istream & operator >> (std::istream &is, cvm::rvector &v); template std::string cvm::to_str(T const &x, size_t const &width, size_t const &prec) { std::ostringstream os; if (width) os.width(width); if (prec) { os.setf(std::ios::scientific, std::ios::floatfield); os.precision(prec); } os << x; return os.str(); } template std::string cvm::to_str(std::vector const &x, size_t const &width, size_t const &prec) { if (!x.size()) return std::string(""); std::ostringstream os; if (prec) { os.setf(std::ios::scientific, std::ios::floatfield); } os << "{ "; if (width) os.width(width); if (prec) os.precision(prec); os << x[0]; for (size_t i = 1; i < x.size(); i++) { os << ", "; if (width) os.width(width); if (prec) os.precision(prec); os << x[i]; } os << " }"; return os.str(); } #include "colvarproxy.h" inline cvm::real cvm::unit_angstrom() { return proxy->unit_angstrom(); } inline cvm::real cvm::boltzmann() { return proxy->boltzmann(); } inline cvm::real cvm::temperature() { return proxy->temperature(); } inline cvm::real cvm::dt() { return proxy->dt(); } // Replica exchange commands inline bool cvm::replica_enabled() { return proxy->replica_enabled(); } inline int cvm::replica_index() { return proxy->replica_index(); } inline int cvm::replica_num() { return proxy->replica_num(); } inline void cvm::replica_comm_barrier() { return proxy->replica_comm_barrier(); } inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) { return proxy->replica_comm_recv(msg_data,buf_len,src_rep); } inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) { return proxy->replica_comm_send(msg_data,msg_len,dest_rep); } inline void cvm::request_system_force() { proxy->request_system_force(true); } inline void cvm::select_closest_image(atom_pos &pos, atom_pos const &ref_pos) { proxy->select_closest_image(pos, ref_pos); } inline void cvm::select_closest_images(std::vector &pos, atom_pos const &ref_pos) { proxy->select_closest_images(pos, ref_pos); } inline cvm::rvector cvm::position_distance(atom_pos const &pos1, atom_pos const &pos2) { return proxy->position_distance(pos1, pos2); } inline cvm::real cvm::position_dist2(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) { return proxy->position_dist2(pos1, pos2); } inline cvm::real cvm::rand_gaussian(void) { return proxy->rand_gaussian(); } #endif