diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index 158ff7a4f1..fa38c4d33f 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -112,22 +112,23 @@ colvar::colvar(std::string const &conf) initialize_components("distance vector", "distanceVec", distance_vec); initialize_components("Cartesian coordinates", "cartesian", cartesian); initialize_components("distance vector " - "direction", "distanceDir", distance_dir); + "direction", "distanceDir", distance_dir); initialize_components("distance projection " - "on an axis", "distanceZ", distance_z); + "on an axis", "distanceZ", distance_z); initialize_components("distance projection " - "on a plane", "distanceXY", distance_xy); + "on a plane", "distanceXY", distance_xy); initialize_components("average distance weighted by inverse power", - "distanceInv", distance_inv); + "distanceInv", distance_inv); initialize_components("N1xN2-long vector of pairwise distances", - "distancePairs", distance_pairs); + "distancePairs", distance_pairs); initialize_components("coordination " - "number", "coordNum", coordnum); + "number", "coordNum", coordnum); initialize_components("self-coordination " - "number", "selfCoordNum", selfcoordnum); + "number", "selfCoordNum", selfcoordnum); initialize_components("angle", "angle", angle); + initialize_components("dipole angle", "dipoleAngle", dipole_angle); initialize_components("dihedral", "dihedral", dihedral); initialize_components("hydrogen bond", "hBond", h_bond); @@ -136,13 +137,13 @@ colvar::colvar(std::string const &conf) initialize_components("alpha helix", "alpha", alpha_angles); initialize_components("dihedral principal " - "component", "dihedralPC", dihedPC); + "component", "dihedralPC", dihedPC); initialize_components("orientation", "orientation", orientation); initialize_components("orientation " - "angle", "orientationAngle",orientation_angle); + "angle", "orientationAngle",orientation_angle); initialize_components("orientation " - "projection", "orientationProj",orientation_proj); + "projection", "orientationProj",orientation_proj); initialize_components("tilt", "tilt", tilt); initialize_components("spin angle", "spinAngle", spin_angle); @@ -151,17 +152,17 @@ colvar::colvar(std::string const &conf) // initialize_components ("logarithm of MSD", "logmsd", logmsd); initialize_components("radius of " - "gyration", "gyration", gyration); + "gyration", "gyration", gyration); initialize_components("moment of " - "inertia", "inertia", inertia); + "inertia", "inertia", inertia); initialize_components("moment of inertia around an axis", - "inertiaZ", inertia_z); + "inertiaZ", inertia_z); initialize_components("eigenvector", "eigenvector", eigenvector); if (!cvcs.size()) { cvm::error("Error: no valid components were provided " - "for this collective variable.\n", - INPUT_ERROR); + "for this collective variable.\n", + INPUT_ERROR); return; } @@ -520,7 +521,8 @@ void colvar::build_atom_list(void) for (size_t i = 0; i < cvcs.size(); i++) { for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { - temp_id_list.push_back(cvcs[i]->atom_groups[j]->at(k).id); + cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); + temp_id_list.push_back(ag[k].id); } } } @@ -830,8 +832,9 @@ void colvar::setup() { for (size_t i = 0; i < cvcs.size(); i++) { for (size_t ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); - atoms.read_positions(); + atoms.setup(); atoms.reset_mass(name,i,ig); + atoms.read_positions(); } } } @@ -861,37 +864,40 @@ colvar::~colvar() void colvar::calc() { - size_t i, ig; if (cvm::debug()) cvm::log("Calculating colvar \""+this->name+"\".\n"); + update_cvc_flags(); + calc_cvcs(); + calc_colvar_properties(); +} + + +int colvar::calc_cvcs(int first_cvc, size_t num_cvcs) +{ + size_t i, ig; + size_t cvc_count; + size_t const cvc_max_count = num_cvcs ? num_cvcs : num_active_cvcs(); // prepare atom groups for calculation if (cvm::debug()) cvm::log("Collecting data from atom groups.\n"); - // Update the enabled/disabled status of cvcs if necessary - if (cvc_flags.size()) { - bool any = false; - for (i = 0; i < cvcs.size(); i++) { - cvcs[i]->b_enabled = cvc_flags[i]; - any = any || cvc_flags[i]; - } - if (!any) { - cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n"); - return; - } - cvc_flags.resize(0); + if (first_cvc >= cvcs.size()) { + cvm::error("Error: trying to address a component outside the " + "range defined for colvar \""+name+"\".\n", BUG_ERROR); + return BUG_ERROR; } - for (i = 0; i < cvcs.size(); i++) { + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { if (!cvcs[i]->b_enabled) continue; + cvc_count++; for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvm::atom_group &atoms = *(cvcs[i]->atom_groups[ig]); atoms.reset_atoms_data(); atoms.read_positions(); - if (atoms.b_center || atoms.b_rotate) { - atoms.calc_apply_roto_translation(); - } + atoms.calc_required_properties(); // each atom group will take care of its own ref_pos_group, if defined } } @@ -906,7 +912,11 @@ void colvar::calc() // } if (tasks[task_system_force]) { - for (i = 0; i < cvcs.size(); i++) { + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { + if (!cvcs[i]->b_enabled) continue; + cvc_count++; for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) { cvcs[i]->atom_groups[ig]->read_system_forces(); } @@ -919,9 +929,12 @@ void colvar::calc() cvm::log("Calculating colvar components.\n"); x.reset(); - // First, update component values - for (i = 0; i < cvcs.size(); i++) { + // First, calculate component values + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { if (!cvcs[i]->b_enabled) continue; + cvc_count++; cvm::increase_depth(); (cvcs[i])->calc_value(); cvm::decrease_depth(); @@ -932,22 +945,25 @@ void colvar::calc() cvm::cv_width, cvm::cv_prec)+".\n"); } - // Then combine them appropriately + // Then combine them appropriately, using either a scripted function or a polynomial if (tasks[task_scripted]) { // cvcs combined by user script int res = cvm::proxy->run_colvar_callback(scripted_function, sorted_cvc_values, x); if (res == COLVARS_NOT_IMPLEMENTED) { cvm::error("Scripted colvars are not implemented."); - return; + return COLVARS_NOT_IMPLEMENTED; } if (res != COLVARS_OK) { cvm::error("Error running scripted colvar"); - return; + return COLVARS_OK; } } else if (x.type() == colvarvalue::type_scalar) { // polynomial combination allowed - for (i = 0; i < cvcs.size(); i++) { + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { if (!cvcs[i]->b_enabled) continue; + cvc_count++; x += (cvcs[i])->sup_coeff * ( ((cvcs[i])->sup_np != 1) ? std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np) : @@ -955,8 +971,11 @@ void colvar::calc() } } else { // only linear combination allowed - for (i = 0; i < cvcs.size(); i++) { + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { if (!cvcs[i]->b_enabled) continue; + cvc_count++; x += (cvcs[i])->sup_coeff * (cvcs[i])->value(); } } @@ -970,9 +989,12 @@ void colvar::calc() if (cvm::debug()) cvm::log("Calculating gradients of colvar \""+this->name+"\".\n"); - for (i = 0; i < cvcs.size(); i++) { - // calculate the gradients of each component + // calculate the gradients of each component + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { if (!cvcs[i]->b_enabled) continue; + cvc_count++; cvm::increase_depth(); (cvcs[i])->calc_gradients(); // if requested, propagate (via chain rule) the gradients above @@ -991,40 +1013,44 @@ void colvar::calc() if (tasks[task_scripted]) { cvm::error("Collecting atomic gradients is not implemented for " - "scripted colvars."); - return; + "scripted colvars.", COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; } // Collect the atomic gradients inside colvar object for (unsigned int a = 0; a < atomic_gradients.size(); a++) { atomic_gradients[a].reset(); } - for (i = 0; i < cvcs.size(); i++) { + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { if (!cvcs[i]->b_enabled) continue; + cvc_count++; // Coefficient: d(a * x^n) = a * n * x^(n-1) * dx cvm::real coeff = (cvcs[i])->sup_coeff * cvm::real((cvcs[i])->sup_np) * std::pow((cvcs[i])->value().real_value, (cvcs[i])->sup_np-1); for (size_t j = 0; j < cvcs[i]->atom_groups.size(); j++) { + cvm::atom_group &ag = *(cvcs[i]->atom_groups[j]); + // If necessary, apply inverse rotation to get atomic // gradient in the laboratory frame if (cvcs[i]->atom_groups[j]->b_rotate) { cvm::rotation const rot_inv = cvcs[i]->atom_groups[j]->rot.inverse(); for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { - int a = std::lower_bound(atom_ids.begin(), atom_ids.end(), - cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin(); - atomic_gradients[a] += coeff * - rot_inv.rotate(cvcs[i]->atom_groups[j]->at(k).grad); + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * rot_inv.rotate(ag[k].grad); } } else { for (size_t k = 0; k < cvcs[i]->atom_groups[j]->size(); k++) { - int a = std::lower_bound(atom_ids.begin(), atom_ids.end(), - cvcs[i]->atom_groups[j]->at(k).id) - atom_ids.begin(); - atomic_gradients[a] += coeff * cvcs[i]->atom_groups[j]->at(k).grad; + size_t a = std::lower_bound(atom_ids.begin(), atom_ids.end(), + ag[k].id) - atom_ids.begin(); + atomic_gradients[a] += coeff * ag[k].grad; } } } @@ -1040,8 +1066,8 @@ void colvar::calc() // TODO see if this could reasonably be done in a generic way // from generic inverse gradients cvm::error("System force is not implemented for " - "scripted colvars."); - return; + "scripted colvars.", COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; } if (cvm::debug()) cvm::log("Calculating system force of colvar \""+this->name+"\".\n"); @@ -1054,7 +1080,11 @@ void colvar::calc() if (cvm::step_relative() > 0) { // get from the cvcs the system forces from the PREVIOUS step - for (i = 0; i < cvcs.size(); i++) { + for (i = first_cvc, cvc_count = 0; + (i < cvcs.size()) && (cvc_count < cvc_max_count); + i++) { + if (!cvcs[i]->b_enabled) continue; + cvc_count++; (cvcs[i])->calc_force_invgrads(); // linear combination is assumed cvm::increase_depth(); @@ -1072,7 +1102,12 @@ void colvar::calc() if (cvm::debug()) cvm::log("Done calculating system force of colvar \""+this->name+"\".\n"); } + return COLVARS_OK; +} + +int colvar::calc_colvar_properties() +{ if (tasks[task_fdiff_velocity]) { // calculate the velocity by finite differences if (cvm::step_relative() == 0) @@ -1109,6 +1144,8 @@ void colvar::calc() if (cvm::debug()) cvm::log("Done calculating colvar \""+this->name+"\".\n"); + + return COLVARS_OK; } @@ -1243,6 +1280,7 @@ void colvar::communicate_forces() if (tasks[task_scripted]) { std::vector > func_grads; + func_grads.reserve(cvcs.size()); for (i = 0; i < cvcs.size(); i++) { if (!cvcs[i]->b_enabled) continue; func_grads.push_back(cvm::matrix2d (x.size(), @@ -1252,7 +1290,7 @@ void colvar::communicate_forces() if (res != COLVARS_OK) { if (res == COLVARS_NOT_IMPLEMENTED) { - cvm::error("Colvar gradient scripts are not implemented."); + cvm::error("Colvar gradient scripts are not implemented.", COLVARS_NOT_IMPLEMENTED); } else { cvm::error("Error running colvar gradient script"); } @@ -1296,8 +1334,8 @@ void colvar::communicate_forces() } -int colvar::set_cvc_flags(std::vector const &flags) { - +int colvar::set_cvc_flags(std::vector const &flags) +{ if (flags.size() != cvcs.size()) { cvm::error("ERROR: Wrong number of CVC flags provided."); return COLVARS_ERROR; @@ -1309,6 +1347,36 @@ int colvar::set_cvc_flags(std::vector const &flags) { } +int colvar::update_cvc_flags() +{ + size_t i; + // Update the enabled/disabled status of cvcs if necessary + if (cvc_flags.size()) { + bool any = false; + for (i = 0; i < cvcs.size(); i++) { + cvcs[i]->b_enabled = cvc_flags[i]; + any = any || cvc_flags[i]; + } + if (!any) { + cvm::error("ERROR: All CVCs are disabled for colvar " + this->name +"\n"); + return COLVARS_ERROR; + } + cvc_flags.resize(0); + } + return COLVARS_OK; +} + + +int colvar::num_active_cvcs() const +{ + int result = 0; + for (size_t i = 0; i < cvcs.size(); i++) { + if (cvcs[i]->b_enabled) result++; + } + return result; +} + + // ******************** METRIC FUNCTIONS ******************** // Use the metrics defined by \link cvc \endlink objects diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index 8a3e40518a..7ab4467bfb 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -310,8 +310,11 @@ public: /// quantities void calc(); - /// \brief Calculate just the value, and store it in the argument - void calc_value(colvarvalue &xn); + /// \brief Calculate a given subset of colvar components (CVCs) (default: all CVCs) + int calc_cvcs(int first = 0, size_t num_cvcs = 0); + + /// \brief Calculate the quantities associated to the colvar (but not to the CVCs) + int calc_colvar_properties(); /// Set the total biasing force to zero void reset_bias_force(); @@ -332,6 +335,12 @@ public: /// \brief Enables and disables individual CVCs based on flags int set_cvc_flags(std::vector const &flags); + /// \brief Updates the flags in the CVC objects + int update_cvc_flags(); + + /// \brief Return the number of CVC objects with an active flag + int num_active_cvcs() const; + /// \brief Use the internal metrics (as from \link cvc /// \endlink objects) to calculate square distances and gradients /// @@ -486,6 +495,7 @@ public: class distance_inv; class distance_pairs; class angle; + class dipole_angle; class dihedral; class coordnum; class selfcoordnum; diff --git a/lib/colvars/colvaratoms.cpp b/lib/colvars/colvaratoms.cpp index 1f0c97516a..038506f347 100644 --- a/lib/colvars/colvaratoms.cpp +++ b/lib/colvars/colvaratoms.cpp @@ -5,56 +5,98 @@ #include "colvaratoms.h" -// member functions of the "atom" class depend tightly on the MD interface, and are -// thus defined in colvarproxy_xxx.cpp - -// in this file only atom_group functions are defined - - -// Note: "conf" is the configuration of the cvc who is using this atom group; -// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...) -cvm::atom_group::atom_group(std::string const &conf, - char const *key) - : b_center(false), b_rotate(false), b_user_defined_fit(false), - b_fit_gradients(false), - ref_pos_group(NULL), - noforce(false) +cvm::atom::atom() { - cvm::log("Defining atom group \""+ - std::string(key)+"\".\n"); - // real work is done by parse - parse(conf, key); + index = -1; + id = -1; + reset_data(); } -cvm::atom_group::atom_group(std::vector const &atoms) - : b_dummy(false), b_center(false), b_rotate(false), - b_fit_gradients(false), ref_pos_group(NULL), - noforce(false) +cvm::atom::atom(int atom_number) { - this->reserve(atoms.size()); - for (size_t i = 0; i < atoms.size(); i++) { - this->push_back(atoms[i]); + colvarproxy *p = cvm::proxy; + index = p->init_atom(atom_number); + if (cvm::debug()) { + cvm::log("The index of this atom in the colvarproxy arrays is "+ + cvm::to_str(index)+".\n"); } - total_mass = 0.0; - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - total_mass += ai->mass; + id = p->get_atom_id(index); + update_mass(); + reset_data(); +} + + +cvm::atom::atom(cvm::residue_id const &residue, + std::string const &atom_name, + std::string const &segment_id) +{ + colvarproxy *p = cvm::proxy; + index = p->init_atom(residue, atom_name, segment_id); + if (cvm::debug()) { + cvm::log("The index of this atom in the colvarproxy_namd arrays is "+ + cvm::to_str(index)+".\n"); } + id = p->get_atom_id(index); + update_mass(); + reset_data(); +} + + +cvm::atom::atom(atom const &a) + : index(a.index) +{ + id = (cvm::proxy)->get_atom_id(index); + update_mass(); + reset_data(); +} + + +cvm::atom::~atom() +{ + if (index >= 0) { + (cvm::proxy)->clear_atom(index); + } +} + + + +// TODO change this arrangement +// Note: "conf" is the configuration of the cvc who is using this atom group; +// "key" is the name of the atom group (e.g. "atoms", "group1", "group2", ...) +cvm::atom_group::atom_group(std::string const &conf, + char const *key_in) +{ + key = key_in; + cvm::log("Defining atom group \""+ + std::string(key)+"\".\n"); + init(); + // real work is done by parse + parse(conf); + setup(); +} + + +cvm::atom_group::atom_group(std::vector const &atoms_in) +{ + init(); + atoms = atoms_in; + setup(); } cvm::atom_group::atom_group() - : b_dummy(false), b_center(false), b_rotate(false), - b_user_defined_fit(false), b_fit_gradients(false), - ref_pos_group(NULL), noforce(false) { - total_mass = 0.0; + init(); } cvm::atom_group::~atom_group() { + if (index >= 0) { + (cvm::proxy)->clear_atom_group(index); + } + if (ref_pos_group) { delete ref_pos_group; ref_pos_group = NULL; @@ -62,53 +104,181 @@ cvm::atom_group::~atom_group() } -void cvm::atom_group::add_atom(cvm::atom const &a) +int cvm::atom_group::add_atom(cvm::atom const &a) +{ + if (a.id < 0) { + return COLVARS_ERROR; + } + + for (size_t i = 0; i < atoms_ids.size(); i++) { + if (atoms_ids[i] == a.id) { + if (cvm::debug()) + cvm::log("Discarding doubly counted atom with number "+ + cvm::to_str(a.id+1)+".\n"); + return COLVARS_OK; + } + } + + // for consistency with add_atom_id(), we update the list as well + atoms_ids.push_back(a.id); + atoms.push_back(a); + total_mass += a.mass; + total_charge += a.charge; + + return COLVARS_OK; +} + + +int cvm::atom_group::add_atom_id(int aid) +{ + if (aid < 0) { + return COLVARS_ERROR; + } + + for (size_t i = 0; i < atoms_ids.size(); i++) { + if (atoms_ids[i] == aid) { + if (cvm::debug()) + cvm::log("Discarding doubly counted atom with number "+ + cvm::to_str(aid+1)+".\n"); + return COLVARS_OK; + } + } + + atoms_ids.push_back(aid); + return COLVARS_OK; +} + + +int cvm::atom_group::remove_atom(cvm::atom_iter ai) +{ + if (b_scalable) { + cvm::error("Error: cannot remove atoms from a scalable group.\n", INPUT_ERROR); + return COLVARS_ERROR; + } + + if (!this->size()) { + cvm::error("Error: trying to remove an atom from an empty group.\n", INPUT_ERROR); + return COLVARS_ERROR; + } else { + total_mass -= ai->mass; + total_charge -= ai->charge; + atoms_ids.erase(atoms_ids.begin() + (ai - atoms.begin())); + atoms.erase(ai); + } + + return COLVARS_OK; +} + + +int cvm::atom_group::init() +{ + if (!key.size()) key = "atoms"; + + atoms.clear(); + + b_scalable = false; + index = -1; + + b_center = false; + b_rotate = false; + b_user_defined_fit = false; + b_fit_gradients = false; + ref_pos_group = NULL; + + noforce = false; + + total_mass = 0.0; + total_charge = 0.0; + + cog.reset(); + com.reset(); + + return COLVARS_OK; +} + + +int cvm::atom_group::setup() +{ + for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) { + ai->update_mass(); + ai->update_charge(); + } + update_total_mass(); + update_total_charge(); + return COLVARS_OK; +} + + +void cvm::atom_group::update_total_mass() { if (b_dummy) { - cvm::error("Error: cannot add atoms to a dummy group.\n", INPUT_ERROR); + total_mass = 1.0; + return; + } + + if (b_scalable) { + total_mass = (cvm::proxy)->get_atom_group_mass(index); } else { - this->push_back(a); - total_mass += a.mass; + total_mass = 0.0; + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + total_mass += ai->mass; + } } } void cvm::atom_group::reset_mass(std::string &name, int i, int j) { - total_mass = 0.0; - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - total_mass += ai->mass; - } + update_total_mass(); cvm::log("Re-initialized atom group "+name+":"+cvm::to_str(i)+"/"+ - cvm::to_str(j)+". "+ cvm::to_str(this->size())+ - " atoms: total mass = "+cvm::to_str(this->total_mass)+".\n"); + cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+ + " atoms: total mass = "+cvm::to_str(total_mass)+".\n"); } -int cvm::atom_group::parse(std::string const &conf, - char const *key) + +void cvm::atom_group::update_total_charge() +{ + if (b_dummy) { + total_charge = 0.0; + return; + } + + if (b_scalable) { + total_charge = (cvm::proxy)->get_atom_group_charge(index); + } else { + total_charge = 0.0; + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + total_charge += ai->charge; + } + } +} + + +int cvm::atom_group::parse(std::string const &conf) { std::string group_conf; + // TODO move this to the cvc class constructor/init + // save_delimiters is set to false for this call, because "conf" is // not the config string of this group, but of its parent object // (which has already taken care of the delimiters) save_delimiters = false; - key_lookup(conf, key, group_conf, dummy_pos); + key_lookup(conf, key.c_str(), group_conf, dummy_pos); // restoring the normal value, because we do want keywords checked // inside "group_conf" save_delimiters = true; if (group_conf.size() == 0) { - cvm::error("Error: atom group \""+std::string(key)+ - "\" is set, but has no definition.\n", - INPUT_ERROR); + cvm::error("Error: atom group \""+key+ + "\" is set, but has no definition.\n", + INPUT_ERROR); return COLVARS_ERROR; } cvm::increase_depth(); - cvm::log("Initializing atom group \""+std::string(key)+"\".\n"); + cvm::log("Initializing atom group \""+key+"\".\n"); // whether or not to include messages in the log // colvarparse::Parse_Mode mode = parse_silent; @@ -117,57 +287,27 @@ int cvm::atom_group::parse(std::string const &conf, // get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent); // if (b_verbose) mode = parse_normal; // } - colvarparse::Parse_Mode mode = parse_normal; + // colvarparse::Parse_Mode mode = parse_normal; + + int parse_error = COLVARS_OK; + + // if the cvc allows it, the flag has been set to true by default + get_keyval(group_conf, "scalable", b_scalable, b_scalable); { - // std::vector atom_indexes; std::string numbers_conf = ""; size_t pos = 0; - std::vector atom_indexes; while (key_lookup(group_conf, "atomNumbers", numbers_conf, pos)) { - if (numbers_conf.size()) { - std::istringstream is(numbers_conf); - int ia; - while (is >> ia) { - atom_indexes.push_back(ia); - } - } - - if (atom_indexes.size()) { - this->reserve(this->size()+atom_indexes.size()); - for (size_t i = 0; i < atom_indexes.size(); i++) { - this->push_back(cvm::atom(atom_indexes[i])); - } - if (cvm::get_error()) return COLVARS_ERROR; - } else { - cvm::error("Error: no numbers provided for \"" - "atomNumbers\".\n", INPUT_ERROR); - return COLVARS_ERROR; - } - - atom_indexes.clear(); + parse_error |= add_atom_numbers(numbers_conf); + numbers_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 - std::list::iterator names_i = cvm::index_group_names.begin(); - std::list >::iterator index_groups_i = cvm::index_groups.begin(); - for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) { - if (*names_i == index_group_name) - break; - } - if (names_i == cvm::index_group_names.end()) { - cvm::error("Error: could not find index group "+ - index_group_name+" among those provided by the index file.\n", - INPUT_ERROR); - return COLVARS_ERROR; - } - this->reserve(index_groups_i->size()); - for (size_t i = 0; i < index_groups_i->size(); i++) { - this->push_back(cvm::atom((*index_groups_i)[i])); - } - if (cvm::get_error()) return COLVARS_ERROR; + parse_error |= add_index_group(index_group_name); } } @@ -175,99 +315,55 @@ int cvm::atom_group::parse(std::string const &conf, std::string range_conf = ""; size_t pos = 0; while (key_lookup(group_conf, "atomNumbersRange", - range_conf, pos)) { - - if (range_conf.size()) { - std::istringstream is(range_conf); - int initial, final; - char dash; - if ( (is >> initial) && (initial > 0) && - (is >> dash) && (dash == '-') && - (is >> final) && (final > 0) ) { - for (int anum = initial; anum <= final; anum++) { - this->push_back(cvm::atom(anum)); - } - if (cvm::get_error()) return COLVARS_ERROR; - range_conf = ""; - continue; - } - } - - cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+ - range_conf+"\".\n", INPUT_ERROR); - } - } - - std::vector psf_segids; - get_keyval(group_conf, "psfSegID", psf_segids, std::vector (), mode); - for (std::vector::iterator psii = psf_segids.begin(); - psii < psf_segids.end(); ++psii) { - - if ( (psii->size() == 0) || (psii->size() > 4) ) { - cvm::error("Error: invalid segmend identifier provided, \""+ - (*psii)+"\".\n", INPUT_ERROR); + range_conf, pos)) { + parse_error |= add_atom_numbers_range(range_conf); + range_conf = ""; } } { + std::vector psf_segids; + get_keyval(group_conf, "psfSegID", psf_segids, std::vector()); + std::vector::iterator psii; + for (psii = psf_segids.begin(); psii < psf_segids.end(); ++psii) { + if ( (psii->size() == 0) || (psii->size() > 4) ) { + cvm::error("Error: invalid PSF segment identifier provided, \""+ + (*psii)+"\".\n", INPUT_ERROR); + } + } + std::string range_conf = ""; size_t pos = 0; size_t range_count = 0; - std::vector::iterator psii = psf_segids.begin(); + psii = psf_segids.begin(); while (key_lookup(group_conf, "atomNameResidueRange", - range_conf, pos)) { + range_conf, pos)) { range_count++; - - if (range_count > psf_segids.size()) { + if (psf_segids.size() && (range_count > psf_segids.size())) { cvm::error("Error: more instances of \"atomNameResidueRange\" than " - "values of \"psfSegID\".\n", INPUT_ERROR); - } - - std::string const &psf_segid = psf_segids.size() ? *psii : std::string(""); - - if (range_conf.size()) { - - std::istringstream is(range_conf); - std::string atom_name; - int initial, final; - char dash; - if ( (is >> atom_name) && (atom_name.size()) && - (is >> initial) && (initial > 0) && - (is >> dash) && (dash == '-') && - (is >> final) && (final > 0) ) { - for (int resid = initial; resid <= final; resid++) { - this->push_back(cvm::atom(resid, atom_name, psf_segid)); - } - if (cvm::get_error()) return COLVARS_ERROR; - range_conf = ""; - } else { - cvm::error("Error: cannot parse definition for \"" - "atomNameResidueRange\", \""+ - range_conf+"\".\n"); - } - + "values of \"psfSegID\".\n", INPUT_ERROR); } else { - cvm::error("Error: atomNameResidueRange with empty definition.\n"); + parse_error |= add_atom_name_residue_range(psf_segids.size() ? *psii : std::string(""), + range_conf); + if (psf_segids.size()) psii++; } - - if (psf_segid.size()) - ++psii; + range_conf = ""; } } { // read the atoms from a file std::string atoms_file_name; - if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""), mode)) { + if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""))) { std::string atoms_col; - if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""), mode)) { + if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""))) { cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n", - INPUT_ERROR); + INPUT_ERROR); } double atoms_col_value; - bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0, mode); + bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0); if (atoms_col_value_defined && (!atoms_col_value)) { cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", INPUT_ERROR); } @@ -277,87 +373,242 @@ int cvm::atom_group::parse(std::string const &conf, } // Catch any errors from all the initialization steps above - if (cvm::get_error()) return COLVARS_ERROR; + if (parse_error || cvm::get_error()) return COLVARS_ERROR; - for (std::vector::iterator a1 = this->begin(); - a1 != this->end(); ++a1) { - std::vector::iterator a2 = a1; - ++a2; - for ( ; a2 != this->end(); ++a2) { - if (a1->id == a2->id) { - if (cvm::debug()) - cvm::log("Discarding doubly counted atom with number "+ - cvm::to_str(a1->id+1)+".\n"); - a2 = this->erase(a2); - if (a2 == this->end()) - break; - } - } + if (b_scalable) { + index = (cvm::proxy)->init_atom_group(atoms_ids); } - if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos(), mode)) { + // checks of doubly-counted atoms have been handled by add_atom() already + + if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos())) { b_dummy = true; - this->total_mass = 1.0; - } else + // note: atoms_ids.size() is used here in lieu of atoms.size(), + // which can be empty for scalable groups + if (atoms_ids.size()) { + cvm::error("Error: cannot set up group \""+ + key+"\" as a dummy atom " + "and provide it with atom definitions.\n", INPUT_ERROR); + } + } else { b_dummy = false; - if (b_dummy && (this->size())) { - cvm::error("Error: cannot set up group \""+ - std::string(key)+"\" as a dummy atom " - "and provide it with atom definitions.\n", INPUT_ERROR); - } - -#if(! defined(COLVARS_STANDALONE)) - if ( (!b_dummy) && (!cvm::b_analysis) && (!(this->size())) ) { - cvm::error("Error: no atoms defined for atom group \""+ - std::string(key)+"\".\n"); - } -#endif - - if (!b_dummy) { - - // calculate total mass (TODO: this is the step that most needs deferred re-initialization) - this->total_mass = 0.0; - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - this->total_mass += ai->mass; + if (!(atoms_ids.size())) { + cvm::error("Error: no atoms defined for atom group \""+ + key+"\".\n", INPUT_ERROR); } // whether these atoms will ever receive forces or not bool enable_forces = true; // disableForces is deprecated - if (get_keyval(group_conf, "enableForces", enable_forces, true, mode)) { + if (get_keyval(group_conf, "enableForces", enable_forces, true)) { noforce = !enable_forces; } else { get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent); } } - // FITTING OPTIONS + parse_error |= parse_fitting_options(group_conf); - bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false, mode); - bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false, mode); + // TODO move this to colvarparse object + check_keywords(group_conf, key.c_str()); + if (cvm::get_error()) { + cvm::error("Error setting up atom group \""+key+"\"."); + return COLVARS_ERROR; + } + + // Calculate all required properties (such as total mass) + setup(); + + if (cvm::debug()) + cvm::log("Done initializing atom group \""+key+"\".\n"); + + cvm::log("Atom group \""+key+"\" defined, "+ + cvm::to_str(atoms_ids.size())+" atoms initialized: total mass = "+ + cvm::to_str(total_mass)+", total charge = "+ + cvm::to_str(total_charge)+".\n"); + + cvm::decrease_depth(); + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + + +int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf) +{ + std::vector atom_indexes; + + if (numbers_conf.size()) { + std::istringstream is(numbers_conf); + int ia; + while (is >> ia) { + atom_indexes.push_back(ia); + } + } + + if (atom_indexes.size()) { + atoms_ids.reserve(atoms_ids.size()+atom_indexes.size()); + + if (b_scalable) { + for (size_t i = 0; i < atom_indexes.size(); i++) { + add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i])); + } + } else { + // if we are handling the group on rank 0, better allocate the vector in one shot + atoms.reserve(atoms.size()+atom_indexes.size()); + for (size_t i = 0; i < atom_indexes.size(); i++) { + add_atom(cvm::atom(atom_indexes[i])); + } + } + + if (cvm::get_error()) return COLVARS_ERROR; + } else { + cvm::error("Error: no numbers provided for \"" + "atomNumbers\".\n", INPUT_ERROR); + return COLVARS_ERROR; + } + + return COLVARS_OK; +} + + +int cvm::atom_group::add_index_group(std::string const &index_group_name) +{ + std::list::iterator names_i = cvm::index_group_names.begin(); + std::list >::iterator index_groups_i = cvm::index_groups.begin(); + for ( ; names_i != cvm::index_group_names.end() ; ++names_i, ++index_groups_i) { + if (*names_i == index_group_name) + break; + } + + if (names_i == cvm::index_group_names.end()) { + cvm::error("Error: could not find index group "+ + index_group_name+" among those provided by the index file.\n", + INPUT_ERROR); + return COLVARS_ERROR; + } + + atoms_ids.reserve(atoms_ids.size()+index_groups_i->size()); + + if (b_scalable) { + for (size_t i = 0; i < index_groups_i->size(); i++) { + add_atom_id((cvm::proxy)->check_atom_id((*index_groups_i)[i])); + } + } else { + atoms.reserve(atoms.size()+index_groups_i->size()); + for (size_t i = 0; i < index_groups_i->size(); i++) { + add_atom(cvm::atom((*index_groups_i)[i])); + } + } + + if (cvm::get_error()) + return COLVARS_ERROR; + + return COLVARS_OK; +} + + +int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf) +{ + if (range_conf.size()) { + std::istringstream is(range_conf); + int initial, final; + char dash; + if ( (is >> initial) && (initial > 0) && + (is >> dash) && (dash == '-') && + (is >> final) && (final > 0) ) { + + atoms_ids.reserve(atoms_ids.size() + (final - initial + 1)); + + if (b_scalable) { + for (int anum = initial; anum <= final; anum++) { + add_atom_id((cvm::proxy)->check_atom_id(anum)); + } + } else { + atoms.reserve(atoms.size() + (final - initial + 1)); + for (int anum = initial; anum <= final; anum++) { + add_atom(cvm::atom(anum)); + } + } + + } + if (cvm::get_error()) return COLVARS_ERROR; + } else { + cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+ + range_conf+"\".\n", INPUT_ERROR); + return COLVARS_ERROR; + } + + return COLVARS_OK; +} + + +int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid, + std::string const &range_conf) +{ + if (range_conf.size()) { + std::istringstream is(range_conf); + std::string atom_name; + int initial, final; + char dash; + if ( (is >> atom_name) && (atom_name.size()) && + (is >> initial) && (initial > 0) && + (is >> dash) && (dash == '-') && + (is >> final) && (final > 0) ) { + + atoms_ids.reserve(atoms_ids.size() + (final - initial + 1)); + + if (b_scalable) { + for (int resid = initial; resid <= final; resid++) { + add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid)); + } + } else { + atoms.reserve(atoms.size() + (final - initial + 1)); + for (int resid = initial; resid <= final; resid++) { + add_atom(cvm::atom(resid, atom_name, psf_segid)); + } + } + + if (cvm::get_error()) return COLVARS_ERROR; + } else { + cvm::error("Error: cannot parse definition for \"" + "atomNameResidueRange\", \""+ + range_conf+"\".\n"); + return COLVARS_ERROR; + } + } else { + cvm::error("Error: atomNameResidueRange with empty definition.\n"); + return COLVARS_ERROR; + } + return COLVARS_OK; +} + + +int cvm::atom_group::parse_fitting_options(std::string const &group_conf) +{ + bool b_defined_center = get_keyval(group_conf, "centerReference", b_center, false); + bool b_defined_rotate = get_keyval(group_conf, "rotateReference", b_rotate, false); // is the user setting explicit options? b_user_defined_fit = b_defined_center || b_defined_rotate; - get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true, mode); + get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true); if (b_center || b_rotate) { if (b_dummy) cvm::error("Error: centerReference or rotateReference " - "cannot be defined for a dummy atom.\n"); + "cannot be defined for a dummy atom.\n"); if (key_lookup(group_conf, "refPositionsGroup")) { // instead of this group, define another group to compute the fit if (ref_pos_group) { cvm::error("Error: the atom group \""+ - std::string(key)+"\" has already a reference group " - "for the rototranslational fit, which was communicated by the " - "colvar component. You should not use refPositionsGroup " - "in this case.\n"); + key+"\" has already a reference group " + "for the rototranslational fit, which was communicated by the " + "colvar component. You should not use refPositionsGroup " + "in this case.\n"); } - cvm::log("Within atom group \""+std::string(key)+"\":\n"); + cvm::log("Within atom group \""+key+"\":\n"); ref_pos_group = new atom_group(group_conf, "refPositionsGroup"); // regardless of the configuration, fit gradients must be calculated by refPositionsGroup @@ -367,25 +618,25 @@ int cvm::atom_group::parse(std::string const &conf, atom_group *group_for_fit = ref_pos_group ? ref_pos_group : this; - get_keyval(group_conf, "refPositions", ref_pos, ref_pos, mode); + get_keyval(group_conf, "refPositions", ref_pos, ref_pos); std::string ref_pos_file; - if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""), mode)) { + if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""))) { if (ref_pos.size()) { cvm::error("Error: cannot specify \"refPositionsFile\" and " - "\"refPositions\" at the same time.\n"); + "\"refPositions\" at the same time.\n"); } std::string ref_pos_col; double ref_pos_col_value=0.0; - if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""), mode)) { + if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""))) { // if provided, use PDB column to select coordinates - bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0, mode); + bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0); if (found && ref_pos_col_value == 0.0) cvm::error("Error: refPositionsColValue, " - "if provided, must be non-zero.\n"); + "if provided, must be non-zero.\n", INPUT_ERROR); } else { // if not, rely on existing atom indices for the group group_for_fit->create_sorted_ids(); @@ -393,7 +644,7 @@ int cvm::atom_group::parse(std::string const &conf, } cvm::load_coords(ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids, - ref_pos_col, ref_pos_col_value); + ref_pos_col, ref_pos_col_value); } if (ref_pos.size()) { @@ -401,21 +652,20 @@ int cvm::atom_group::parse(std::string const &conf, if (b_rotate) { if (ref_pos.size() != group_for_fit->size()) cvm::error("Error: the number of reference positions provided("+ - cvm::to_str(ref_pos.size())+ - ") does not match the number of atoms within \""+ - std::string(key)+ - "\" ("+cvm::to_str(group_for_fit->size())+ - "): to perform a rotational fit, "+ - "these numbers should be equal.\n", INPUT_ERROR); + cvm::to_str(ref_pos.size())+ + ") does not match the number of atoms within \""+ + key+ + "\" ("+cvm::to_str(group_for_fit->size())+ + "): to perform a rotational fit, "+ + "these numbers should be equal.\n", INPUT_ERROR); } // save the center of geometry of ref_pos and subtract it center_ref_pos(); } else { -#if(! defined(COLVARS_STANDALONE)) - cvm::error("Error: no reference positions provided.\n"); -#endif + cvm::error("Error: no reference positions provided.\n", INPUT_ERROR); + return COLVARS_ERROR; } if (b_fit_gradients) { @@ -424,35 +674,18 @@ int cvm::atom_group::parse(std::string const &conf, } if (b_rotate && !noforce) { - cvm::log("Warning: atom group \""+std::string(key)+ - "\" will be aligned to a fixed orientation given by the reference positions provided. " - "If the internal structure of the group changes too much (i.e. its RMSD is comparable " - "to its radius of gyration), the optimal rotation and its gradients may become discontinuous. " - "If that happens, use refPositionsGroup (or a different definition for it if already defined) " - "to align the coordinates.\n"); + cvm::log("Warning: atom group \""+key+ + "\" will be aligned to a fixed orientation given by the reference positions provided. " + "If the internal structure of the group changes too much (i.e. its RMSD is comparable " + "to its radius of gyration), the optimal rotation and its gradients may become discontinuous. " + "If that happens, use refPositionsGroup (or a different definition for it if already defined) " + "to align the coordinates.\n"); // initialize rot member data rot.request_group1_gradients(this->size()); } - } - if (cvm::debug()) - cvm::log("Done initializing atom group with name \""+ - std::string(key)+"\".\n"); - - this->check_keywords(group_conf, key); - if (cvm::get_error()) { - cvm::error("Error setting up atom group \""+std::string(key)+"\"."); - return COLVARS_ERROR; - } - - cvm::log("Atom group \""+std::string(key)+"\" defined, "+ - cvm::to_str(this->size())+" atoms initialized: total mass = "+ - cvm::to_str(this->total_mass)+".\n"); - - cvm::decrease_depth(); - - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return COLVARS_OK; } @@ -463,17 +696,16 @@ int cvm::atom_group::create_sorted_ids(void) return COLVARS_OK; std::list temp_id_list; - cvm::atom_iter ai; - for (ai = this->begin(); ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { temp_id_list.push_back(ai->id); } temp_id_list.sort(); temp_id_list.unique(); if (temp_id_list.size() != this->size()) { cvm::error("Error: duplicate atom IDs in atom group? (found " + - cvm::to_str(temp_id_list.size()) + - " unique atom IDs instead of" + - cvm::to_str(this->size()) + ").\n"); + cvm::to_str(temp_id_list.size()) + + " unique atom IDs instead of" + + cvm::to_str(this->size()) + ").\n"); return COLVARS_ERROR; } sorted_ids = std::vector (temp_id_list.size()); @@ -486,13 +718,9 @@ int cvm::atom_group::create_sorted_ids(void) return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); } + void cvm::atom_group::center_ref_pos() { - // save the center of geometry of ref_pos and then subtract it from - // them; in this way it will be possible to use ref_pos also for - // the rotational fit - // This is called either by atom_group::parse or by CVCs that set - // reference positions (eg. RMSD, eigenvector) ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0); std::vector::iterator pi; for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) { @@ -504,12 +732,12 @@ void cvm::atom_group::center_ref_pos() } } + void cvm::atom_group::read_positions() { if (b_dummy) return; - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_position(); } @@ -517,6 +745,26 @@ void cvm::atom_group::read_positions() ref_pos_group->read_positions(); } + +int cvm::atom_group::calc_required_properties() +{ + if (b_dummy) return COLVARS_OK; + + // TODO check if the com is needed? + calc_center_of_mass(); + if (!b_scalable) { + // TODO check if calc_center_of_geometry() is needed without a fit? + calc_center_of_geometry(); + if (b_center || b_rotate) { + calc_apply_roto_translation(); + } + } + + // TODO calculate elements of scalable cvc's here before reduction + + return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); +} + void cvm::atom_group::calc_apply_roto_translation() { atom_group *fit_group = ref_pos_group ? ref_pos_group : this; @@ -524,10 +772,7 @@ void cvm::atom_group::calc_apply_roto_translation() if (b_center) { // center on the origin first cvm::atom_pos const cog = fit_group->center_of_geometry(); - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - ai->pos -= cog; - } + apply_translation(-1.0 * cog); } if (b_rotate) { @@ -535,38 +780,53 @@ void cvm::atom_group::calc_apply_roto_translation() // true, around the origin otherwise) rot.calc_optimal_rotation(fit_group->positions(), ref_pos); - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos = rot.rotate(ai->pos); } } if (b_center) { // align with the center of geometry of ref_pos - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - ai->pos += ref_pos_cog; - } + apply_translation(ref_pos_cog); + // update the center of geometry for external use + cog = ref_pos_cog; } + + // recalculate the center of mass + calc_center_of_mass(); } void cvm::atom_group::apply_translation(cvm::rvector const &t) { - if (b_dummy) return; + if (b_dummy) { + cvm::error("Error: cannot translate the coordinates of a dummy atom group.\n", INPUT_ERROR); + return; + } - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + if (b_scalable) { + cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", INPUT_ERROR); + return; + } + + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->pos += t; } } void cvm::atom_group::apply_rotation(cvm::rotation const &rot) { - if (b_dummy) return; + if (b_dummy) { + cvm::error("Error: cannot rotate the coordinates of a dummy atom group.\n", INPUT_ERROR); + return; + } - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - ai->pos = rot.rotate(ai->pos); + if (b_scalable) { + cvm::error("Error: cannot rotate the coordinates of a scalable atom group.\n", INPUT_ERROR); + return; + } + + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + ai->pos = rot.rotate(ai->pos - center_of_geometry()) + center_of_geometry(); } } @@ -577,16 +837,14 @@ void cvm::atom_group::read_velocities() if (b_rotate) { - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_velocity(); ai->vel = rot.rotate(ai->vel); } } else { - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_velocity(); } } @@ -598,47 +856,63 @@ void cvm::atom_group::read_system_forces() if (b_rotate) { - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_system_force(); ai->system_force = rot.rotate(ai->system_force); } } else { - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->read_system_force(); } } } -cvm::atom_pos cvm::atom_group::center_of_geometry() const -{ - if (b_dummy) - return dummy_atom_pos; - cvm::atom_pos cog(0.0, 0.0, 0.0); - for (cvm::atom_const_iter ai = this->begin(); - ai != this->end(); ai++) { - cog += ai->pos; +int cvm::atom_group::calc_center_of_geometry() +{ + if (b_dummy) { + cog = dummy_atom_pos; + } else { + cog.reset(); + for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { + cog += ai->pos; + } + cog /= this->size(); } - cog /= this->size(); - return cog; + return COLVARS_OK; } -cvm::atom_pos cvm::atom_group::center_of_mass() const -{ - if (b_dummy) - return dummy_atom_pos; - cvm::atom_pos com(0.0, 0.0, 0.0); - for (cvm::atom_const_iter ai = this->begin(); - ai != this->end(); ai++) { - com += ai->mass * ai->pos; +int cvm::atom_group::calc_center_of_mass() +{ + if (b_dummy) { + com = dummy_atom_pos; + } else if (b_scalable) { + com = (cvm::proxy)->get_atom_group_com(index); + } else { + com.reset(); + for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { + com += ai->mass * ai->pos; + } + com /= total_mass; } - com /= this->total_mass; - return com; + return COLVARS_OK; +} + + +int cvm::atom_group::calc_dipole(cvm::atom_pos const &com) +{ + if (b_dummy) { + cvm::error("Error: trying to compute the dipole of an empty group.\n", INPUT_ERROR); + return COLVARS_ERROR; + } + dip.reset(); + for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) { + dip += ai->charge * (ai->pos - com); + } + return COLVARS_OK; } @@ -646,9 +920,13 @@ void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad) { if (b_dummy) return; - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - ai->grad = (ai->mass/this->total_mass) * grad; + if (b_scalable) { + scalar_com_gradient = grad; + return; + } + + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + ai->grad = (ai->mass/total_mass) * grad; } } @@ -670,8 +948,8 @@ void cvm::atom_group::calc_fit_gradients() for (size_t i = 0; i < this->size(); i++) { // need to bring the gradients in original frame first cvm::rvector const atom_grad = b_rotate ? - (rot.inverse()).rotate((*this)[i].grad) : - (*this)[i].grad; + (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())) * @@ -684,16 +962,16 @@ 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 = this->center_of_geometry(); + 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 ? ((*this)[i].pos - cog) : ((*this)[i].pos))); + cvm::atom_pos const pos_orig = rot_inv.rotate((b_center ? (atoms[i].pos - cog) : (atoms[i].pos))); for (size_t j = 0; j < group_for_fit->size(); j++) { // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i cvm::quaternion const dxdq = - rot.q.position_derivative_inner(pos_orig, (*this)[i].grad); + 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]; @@ -710,7 +988,12 @@ std::vector cvm::atom_group::positions() const { if (b_dummy) { cvm::error("Error: positions are not available " - "from a dummy atom group.\n"); + "from a dummy atom group.\n", INPUT_ERROR); + } + + if (b_scalable) { + cvm::error("Error: atomic positions are not available " + "from a scalable atom group.\n", INPUT_ERROR); } std::vector x(this->size(), 0.0); @@ -726,7 +1009,12 @@ std::vector cvm::atom_group::positions_shifted(cvm::rvector const { if (b_dummy) { cvm::error("Error: positions are not available " - "from a dummy atom group.\n"); + "from a dummy atom group.\n", INPUT_ERROR); + } + + if (b_scalable) { + cvm::error("Error: atomic positions are not available " + "from a scalable atom group.\n", INPUT_ERROR); } std::vector x(this->size(), 0.0); @@ -742,7 +1030,12 @@ std::vector cvm::atom_group::velocities() const { if (b_dummy) { cvm::error("Error: velocities are not available " - "from a dummy atom group.\n"); + "from a dummy atom group.\n", INPUT_ERROR); + } + + if (b_scalable) { + cvm::error("Error: atomic velocities are not available " + "from a scalable atom group.\n", INPUT_ERROR); } std::vector v(this->size(), 0.0); @@ -758,7 +1051,12 @@ std::vector cvm::atom_group::system_forces() const { if (b_dummy) { cvm::error("Error: system forces are not available " - "from a dummy atom group.\n"); + "from a dummy atom group.\n", INPUT_ERROR); + } + + if (b_scalable) { + cvm::error("Error: atomic system forces are not available " + "from a scalable atom group.\n", INPUT_ERROR); } std::vector f(this->size(), 0.0); @@ -773,8 +1071,12 @@ std::vector cvm::atom_group::system_forces() const cvm::rvector cvm::atom_group::system_force() const { if (b_dummy) { - cvm::error("Error: system forces are not available " - "from a dummy atom group.\n"); + cvm::error("Error: total system forces are not available " + "from a dummy atom group.\n", INPUT_ERROR); + } + + if (b_scalable) { + return (cvm::proxy)->get_atom_group_system_force(index); } cvm::rvector f(0.0); @@ -787,12 +1089,21 @@ cvm::rvector cvm::atom_group::system_force() const void cvm::atom_group::apply_colvar_force(cvm::real const &force) { + if (cvm::debug()) { + log("Communicating a colvar force from atom group to the MD engine.\n"); + } + if (b_dummy) return; if (noforce) { cvm::error("Error: sending a force to a group that has " - "\"enableForces\" set to off.\n"); + "\"enableForces\" set to off.\n"); + return; + } + + if (b_scalable) { + (cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient); return; } @@ -800,15 +1111,13 @@ void cvm::atom_group::apply_colvar_force(cvm::real const &force) // rotate forces back to the original frame cvm::rotation const rot_inv = rot.inverse(); - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force(rot_inv.rotate(force * ai->grad)); } } else { - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { ai->apply_force(force * ai->grad); } } @@ -841,23 +1150,25 @@ void cvm::atom_group::apply_force(cvm::rvector const &force) if (noforce) { cvm::error("Error: sending a force to a group that has " - "\"disableForces\" defined.\n"); + "\"disableForces\" defined.\n"); return; } + if (b_scalable) { + (cvm::proxy)->apply_atom_group_force(index, force); + } + if (b_rotate) { cvm::rotation const rot_inv = rot.inverse(); - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - ai->apply_force(rot_inv.rotate((ai->mass/this->total_mass) * force)); + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + ai->apply_force(rot_inv.rotate((ai->mass/total_mass) * force)); } } else { - for (cvm::atom_iter ai = this->begin(); - ai != this->end(); ai++) { - ai->apply_force((ai->mass/this->total_mass) * force); + for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) { + ai->apply_force((ai->mass/total_mass) * force); } } } @@ -870,11 +1181,11 @@ void cvm::atom_group::apply_forces(std::vector const &forces) if (noforce) cvm::error("Error: sending a force to a group that has " - "\"disableForces\" defined.\n"); + "\"disableForces\" defined.\n"); if (forces.size() != this->size()) { cvm::error("Error: trying to apply an array of forces to an atom " - "group which does not have the same length.\n"); + "group which does not have the same length.\n"); } if (b_rotate) { diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index ef3bbcb39a..746d4c2918 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -6,43 +6,44 @@ #include "colvarmodule.h" #include "colvarparse.h" + /// \brief Stores numeric id, mass and all mutable data for an atom, /// mostly used by a \link cvc \endlink /// -/// This class may be used (although not necessarily) to keep atomic -/// data (id, mass, position and collective variable derivatives) -/// altogether. There may be multiple instances with identical +/// This class may be used to keep atomic data such as id, mass, +/// position and collective variable derivatives) altogether. +/// There may be multiple instances with identical /// numeric id, all acting independently: forces communicated through /// these instances will be summed together. -/// -/// Read/write operations depend on the underlying code: hence, some -/// member functions are defined in colvarproxy_xxx.h. + class colvarmodule::atom { protected: - /// \brief Index in the list of atoms involved by the colvars (\b - /// NOT in the global topology!) - int index; + /// Index in the colvarproxy arrays (\b NOT in the global topology!) + int index; public: - /// Internal identifier (zero-based) - int id; + /// Identifier for the MD program (0-based) + int id; /// Mass - cvm::real mass; + cvm::real mass; + + /// Charge + cvm::real charge; /// \brief Current position (copied from the program, can be - /// manipulated) + /// modified if necessary) cvm::atom_pos pos; /// \brief Current velocity (copied from the program, can be - /// manipulated) + /// modified if necessary) cvm::rvector vel; /// \brief System force at the previous step (copied from the - /// program, can be manipulated) + /// program, can be modified if necessary) cvm::rvector system_force; /// \brief Gradient of a scalar collective variable with respect @@ -57,13 +58,13 @@ public: /// implementation cvm::rvector grad; - /// \brief Default constructor, setting index and id to invalid numbers - atom() : index(-1), id(-1) { reset_data(); } + /// \brief Default constructor (sets index and id both to -1) + atom(); /// \brief Initialize an atom for collective variable calculation /// and get its internal identifier \param atom_number Atom index in /// the system topology (starting from 1) - atom(int const &atom_number); + atom(int atom_number); /// \brief Initialize an atom for collective variable calculation /// and get its internal identifier \param residue Residue number @@ -71,8 +72,8 @@ public: /// segment_id For PSF topologies, the segment identifier; for other /// type of topologies, may not be required atom(cvm::residue_id const &residue, - std::string const &atom_name, - std::string const &segment_id = std::string("")); + std::string const &atom_name, + std::string const &segment_id); /// Copy constructor atom(atom const &a); @@ -80,57 +81,176 @@ public: /// Destructor ~atom(); - /// Set non-constant data (everything except id and mass) to zero - inline void reset_data() { - pos = atom_pos(0.0); - vel = grad = system_force = rvector(0.0); + /// Set mutable data (everything except id and mass) to zero; update mass + inline void reset_data() + { + pos = cvm::atom_pos(0.0); + vel = grad = system_force = cvm::rvector(0.0); + } + + /// Get the latest value of the mass + inline void update_mass() + { + mass = (cvm::proxy)->get_atom_mass(index); + } + + /// Get the latest value of the charge + inline void update_charge() + { + charge = (cvm::proxy)->get_atom_charge(index); } /// Get the current position - void read_position(); + inline void read_position() + { + pos = (cvm::proxy)->get_atom_position(index); + } /// Get the current velocity - void read_velocity(); + inline void read_velocity() + { + vel = (cvm::proxy)->get_atom_velocity(index); + } /// Get the system force - void read_system_force(); + inline void read_system_force() + { + system_force = (cvm::proxy)->get_atom_system_force(index); + } /// \brief Apply a force to the atom /// - /// The force will be used later by the MD integrator, the - /// collective variables module does not integrate equations of - /// motion. Multiple calls to this function by either the same + /// Note: the force is not applied instantly, but will be used later + /// by the MD integrator (the colvars module does not integrate + /// equations of motion. + /// + /// Multiple calls to this function by either the same /// \link atom \endlink object or different objects with identical - /// \link id \endlink, will all add to the existing MD force. - void apply_force(cvm::rvector const &new_force); + /// \link id \endlink will all be added together. + inline void apply_force(cvm::rvector const &new_force) const + { + (cvm::proxy)->apply_atom_force(index, new_force); + } }; - /// \brief Group of \link atom \endlink objects, mostly used by a -/// \link cvc \endlink -/// -/// This class inherits from \link colvarparse \endlink and from -/// std::vector, and hence all functions and -/// operators (including the bracket operator, group[i]) can be used -/// on an \link atom_group \endlink object. It can be initialized as -/// a vector, or by parsing a keyword in the configuration. +/// \link cvc \endlink object to gather all atomic data class colvarmodule::atom_group - : public std::vector, - public colvarparse + : public colvarparse { public: - // Note: all members here are kept public, to allow any - // object accessing and manipulating them + /// \brief Initialize the group by looking up its configuration + /// string in conf and parsing it; this is actually done by parse(), + /// which is a member function so that a group can be initialized + /// also after construction + atom_group(std::string const &conf, + char const *key); + + /// \brief Keyword used to define the group + // TODO Make this field part of the data structures that link a group to a CVC + std::string key; + + /// \brief Set default values for common flags + int init(); + + /// \brief Update data required to calculate cvc's + int setup(); + + /// \brief Initialize the group by looking up its configuration + /// string in conf and parsing it + int parse(std::string const &conf); + + int add_atom_numbers(std::string const &numbers_conf); + int add_index_group(std::string const &index_group_name); + int add_atom_numbers_range(std::string const &range_conf); + int add_atom_name_residue_range(std::string const &psf_segid, + std::string const &range_conf); + int parse_fitting_options(std::string const &group_conf); + + /// \brief Initialize the group after a (temporary) vector of atoms + atom_group(std::vector const &atoms_in); + + /// \brief Add an atom object to this group + int add_atom(cvm::atom const &a); + + /// \brief Add an atom ID to this group (the actual atomicdata will be not be handled by the group) + int add_atom_id(int aid); + + /// \brief Remove an atom object from this group + int remove_atom(cvm::atom_iter ai); + + /// \brief Re-initialize the total mass of a group. + /// This is needed in case the hosting MD code has an option to + /// change atom masses after their initialization. + void reset_mass(std::string &name, int i, int j); + + /// \brief Whether or not the properties of this group will be computed in parallel + bool b_scalable; + + /// \brief Default constructor + atom_group(); + + /// \brief Destructor + ~atom_group(); + +protected: + + /// \brief Array of atom objects + std::vector atoms; + + /// \brief Array of atom identifiers for the MD program (0-based) + std::vector atoms_ids; + + /// \brief Dummy atom position + cvm::atom_pos dummy_atom_pos; + + /// \brief Index in the colvarproxy arrays (if the group is scalable) + int index; + +public: + + inline cvm::atom & operator [] (size_t const i) + { + return atoms[i]; + } + + inline cvm::atom const & operator [] (size_t const i) const + { + return atoms[i]; + } + + inline cvm::atom_iter begin() + { + return atoms.begin(); + } + + inline cvm::atom_const_iter begin() const + { + return atoms.begin(); + } + + inline cvm::atom_iter end() + { + return atoms.end(); + } + + inline cvm::atom_const_iter end() const + { + return atoms.end(); + } + + inline size_t size() const + { + return atoms.size(); + } /// \brief If this option is on, this group merely acts as a wrapper /// for a fixed position; any calls to atoms within or to /// functions that return disaggregated data will fail bool b_dummy; - /// \brief dummy atom position - cvm::atom_pos dummy_atom_pos; /// Sorted list of zero-based (internal) atom ids /// (populated on-demand by create_sorted_ids) @@ -177,56 +297,35 @@ public: /// Total mass of the atom group cvm::real total_mass; + void update_total_mass(); + + /// Total charge of the atom group + cvm::real total_charge; + void update_total_charge(); /// \brief Don't apply any force on this group (use its coordinates /// only to calculate a colvar) bool noforce; - - /// \brief Initialize the group by looking up its configuration - /// string in conf and parsing it; this is actually done by parse(), - /// which is a member function so that a group can be initialized - /// also after construction - atom_group(std::string const &conf, - char const *key); - - /// \brief Initialize the group by looking up its configuration - /// string in conf and parsing it - int parse(std::string const &conf, - char const *key); - - /// \brief Initialize the group after a temporary vector of atoms - atom_group(std::vector const &atoms); - - /// \brief Add an atom to this group - void add_atom(cvm::atom const &a); - - /// \brief Re-initialize the total mass of a group. - /// This is needed in case the hosting MD code has an option to - /// change atom masses after their initialization. - void reset_mass(std::string &name, int i, int j); - - /// \brief Default constructor - atom_group(); - - /// \brief Destructor - ~atom_group(); - - /// \brief Get the current positions; if b_center or b_rotate are - /// true, calc_apply_roto_translation() will be called too + /// \brief Get the current positions void read_positions(); /// \brief (Re)calculate the optimal roto-translation void calc_apply_roto_translation(); - /// \brief Save center of geometry fo ref positions, - /// then subtract it + /// \brief Save aside the center of geometry of the reference positions, + /// then subtract it from them + /// + /// In this way it will be possible to use ref_pos also for the + /// rotational fit. + /// This is called either by atom_group::parse or by CVCs that assign + /// reference positions (eg. RMSD, eigenvector). void center_ref_pos(); /// \brief Move all positions void apply_translation(cvm::rvector const &t); - /// \brief Rotate all positions + /// \brief Rotate all positions around the center of geometry void apply_rotation(cvm::rotation const &q); @@ -243,36 +342,68 @@ public: /// Call reset_data() for each atom inline void reset_atoms_data() { - for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) + for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) ai->reset_data(); if (ref_pos_group) ref_pos_group->reset_atoms_data(); } + /// \brief Recompute all mutable quantities that are required to compute CVCs + int calc_required_properties(); + /// \brief Return a copy of the current atom positions std::vector positions() const; + /// \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; +public: + /// \brief Return the center of geometry of the atomic positions + inline cvm::atom_pos center_of_geometry() const + { + return cog; + } + + /// \brief Calculate the center of mass of the atomic positions, assuming that + /// they are already pbc-wrapped + int calc_center_of_mass(); +private: + /// \brief Center of mass + cvm::atom_pos com; + /// \brief The derivative of a scalar variable with respect to the COM + // TODO for scalable calculations of more complex variables (e.g. rotation), + // use a colvarvalue of vectors to hold the entire derivative + cvm::rvector scalar_com_gradient; +public: + /// \brief Return the center of mass of the atomic positions + inline cvm::atom_pos center_of_mass() const + { + return com; + } + /// \brief Return a copy of the current atom positions, shifted by a constant vector std::vector positions_shifted(cvm::rvector const &shift) const; - /// \brief Return the center of geometry of the positions, assuming - /// that coordinates are already pbc-wrapped - cvm::atom_pos center_of_geometry() const; - - /// \brief Return the center of mass of the positions, assuming that - /// coordinates are already pbc-wrapped - cvm::atom_pos center_of_mass() const; - - /// \brief Atom positions at the previous step - std::vector old_pos; - - /// \brief Return a copy of the current atom velocities std::vector velocities() const; + ///\brief Calculate the dipole of the atom group around the specified center + int calc_dipole(cvm::atom_pos const &com); +private: + cvm::rvector dip; +public: + ///\brief Return the (previously calculated) dipole of the atom group + inline cvm::rvector dipole() const + { + return dip; + } /// \brief Return a copy of the system forces std::vector system_forces() const; + /// \brief Return a copy of the aggregated total force on the group cvm::rvector system_force() const; diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index c05ab3f98b..88e44ddd86 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -46,12 +46,14 @@ colvar::cvc::cvc(std::string const &conf) void colvar::cvc::parse_group(std::string const &conf, - char const *group_key, - cvm::atom_group &group, - bool optional) + char const *group_key, + cvm::atom_group &group, + bool optional) { if (key_lookup(conf, group_key)) { - if (group.parse(conf, group_key) != COLVARS_OK) { + // TODO turn on scalable flag for group objects in cvc init function + group.key = group_key; + if (group.parse(conf) != COLVARS_OK) { cvm::error("Error parsing definition for atom group \""+ std::string(group_key)+"\".\n"); return; @@ -130,10 +132,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) group.read_positions(); // change one coordinate group[ia].pos[id] += cvm::debug_gradients_step_size; - // (re)do the fit (if defined) - if (group.b_center || group.b_rotate) { - group.calc_apply_roto_translation(); - } + 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]; @@ -167,7 +166,7 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) // ref.read_positions(); // // change one coordinate // ref[ia].pos[id] += cvm::debug_gradients_step_size; -// group.calc_apply_roto_translation(); +// 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"); diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 0a66253b1d..2f0e4b43b7 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -670,6 +670,49 @@ public: colvarvalue const &x2) const; }; +/// \brief Colvar component: angle between the dipole of a molecule and an axis +/// formed by two groups of atoms(colvarvalue::type_scalar type, range [0:PI]) +class colvar::dipole_angle + : public colvar::cvc +{ +protected: + + /// Dipole atom group + cvm::atom_group group1; + /// Atom group + cvm::atom_group group2; + /// Atom group + cvm::atom_group group3; + + /// Inter site vectors + cvm::rvector r21, r23; + /// Inter site vector norms + cvm::real r21l, r23l; + /// Derivatives wrt group centers of mass + cvm::rvector dxdr1, dxdr3; + + /// Compute system force on first site only to avoid unwanted + /// coupling to other colvars (see e.g. Ciccotti et al., 2005) + /// (or to allow dummy atoms) + bool b_1site_force; +public: + + /// Initialize by parsing the configuration + dipole_angle (std::string const &conf); + /// \brief Initialize the three groups after three atoms + dipole_angle (cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3); + dipole_angle(); + virtual inline ~dipole_angle() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force (colvarvalue const &force); + virtual cvm::real dist2 (colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_lgrad (colvarvalue const &x1, + colvarvalue const &x2) const; + virtual colvarvalue dist2_rgrad (colvarvalue const &x1, + colvarvalue const &x2) const; +}; /// \brief Colvar component: dihedral between the centers of mass of /// four groups (colvarvalue::type_scalar type, range [-PI:PI]) @@ -1190,6 +1233,7 @@ simple_scalar_dist_functions(distance) simple_scalar_dist_functions(distance_xy) simple_scalar_dist_functions(distance_inv) simple_scalar_dist_functions(angle) +simple_scalar_dist_functions(dipole_angle) simple_scalar_dist_functions(coordnum) simple_scalar_dist_functions(selfcoordnum) simple_scalar_dist_functions(h_bond) diff --git a/lib/colvars/colvarcomp_angles.cpp b/lib/colvars/colvarcomp_angles.cpp index 726b6befcc..35b98c0e46 100644 --- a/lib/colvars/colvarcomp_angles.cpp +++ b/lib/colvars/colvarcomp_angles.cpp @@ -29,9 +29,9 @@ colvar::angle::angle(std::string const &conf) colvar::angle::angle(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3) - : group1(std::vector (1, a1)), - group2(std::vector (1, a2)), - group3(std::vector (1, a3)) + : group1(std::vector(1, a1)), + group2(std::vector(1, a2)), + group3(std::vector(1, a3)) { function_type = "angle"; b_inverse_gradients = true; @@ -95,6 +95,12 @@ void colvar::angle::calc_gradients() group3[i].grad = (group3[i].mass/group3.total_mass) * (dxdr3); } + + if (b_debug_gradients) { + debug_gradients(group1); + debug_gradients(group2); + debug_gradients(group3); + } } void colvar::angle::calc_force_invgrads() @@ -143,6 +149,116 @@ void colvar::angle::apply_force(colvarvalue const &force) +colvar::dipole_angle::dipole_angle(std::string const &conf) + : cvc(conf) +{ + function_type = "dipole_angle"; + parse_group(conf, "group1", group1); + parse_group(conf, "group2", group2); + parse_group(conf, "group3", group3); + + atom_groups.push_back(&group1); + atom_groups.push_back(&group2); + atom_groups.push_back(&group3); + if (get_keyval(conf, "oneSiteSystemForce", b_1site_force, false)) { + cvm::log("Computing system force on group 1 only"); + } + x.type(colvarvalue::type_scalar); +} + + +colvar::dipole_angle::dipole_angle(cvm::atom const &a1, + cvm::atom const &a2, + cvm::atom const &a3) + : group1(std::vector(1, a1)), + group2(std::vector(1, a2)), + group3(std::vector(1, a3)) +{ + function_type = "dipole_angle"; + b_1site_force = false; + atom_groups.push_back(&group1); + atom_groups.push_back(&group2); + atom_groups.push_back(&group3); + + x.type(colvarvalue::type_scalar); +} + + +colvar::dipole_angle::dipole_angle() +{ + function_type = "dipole_angle"; + x.type(colvarvalue::type_scalar); +} + + +void colvar::dipole_angle::calc_value() +{ + cvm::atom_pos const g1_pos = group1.center_of_mass(); + cvm::atom_pos const g2_pos = group2.center_of_mass(); + cvm::atom_pos const g3_pos = group3.center_of_mass(); + + group1.calc_dipole(g1_pos); + + r21 = group1.dipole(); + r21l = r21.norm(); + r23 = cvm::position_distance(g2_pos, g3_pos); + r23l = r23.norm(); + + cvm::real const cos_theta = (r21*r23)/(r21l*r23l); + + x.real_value = (180.0/PI) * std::acos(cos_theta); +} + +//to be implemented +//void colvar::dipole_angle::calc_force_invgrads(){} +//void colvar::dipole_angle::calc_Jacobian_derivative(){} + +void colvar::dipole_angle::calc_gradients() +{ + cvm::real const cos_theta = (r21*r23)/(r21l*r23l); + cvm::real const dxdcos = -1.0 / std::sqrt(1.0 - cos_theta*cos_theta); + + dxdr1 = (180.0/PI) * dxdcos * + (1.0/r21l)* (r23/r23l + (-1.0) * cos_theta * r21/r21l ); + + dxdr3 = (180.0/PI) * dxdcos * + (1.0/r23l) * ( r21/r21l + (-1.0) * cos_theta * r23/r23l ); + + //this auxiliar variables are to avoid numerical errors inside "for" + double aux1 = group1.total_charge/group1.total_mass; + // double aux2 = group2.total_charge/group2.total_mass; + // double aux3 = group3.total_charge/group3.total_mass; + + size_t i; + for (i = 0; i < group1.size(); i++) { + group1[i].grad =(group1[i].charge + (-1)* group1[i].mass * aux1) * (dxdr1); + } + + for (i = 0; i < group2.size(); i++) { + group2[i].grad = (group2[i].mass/group2.total_mass)* dxdr3 * (-1.0); + } + + for (i = 0; i < group3.size(); i++) { + group3[i].grad =(group3[i].mass/group3.total_mass) * (dxdr3); + } +} + + +void colvar::dipole_angle::apply_force(colvarvalue const &force) +{ + if (!group1.noforce) + group1.apply_colvar_force(force.real_value); + + if (!group2.noforce) + group2.apply_colvar_force(force.real_value); + + if (!group3.noforce) + group3.apply_colvar_force(force.real_value); +} + + + + colvar::dihedral::dihedral(std::string const &conf) : cvc(conf) { @@ -171,10 +287,10 @@ colvar::dihedral::dihedral(cvm::atom const &a1, cvm::atom const &a2, cvm::atom const &a3, cvm::atom const &a4) - : group1(std::vector (1, a1)), - group2(std::vector (1, a2)), - group3(std::vector (1, a3)), - group4(std::vector (1, a4)) + : group1(std::vector(1, a1)), + group2(std::vector(1, a2)), + group3(std::vector(1, a3)), + group4(std::vector(1, a4)) { if (cvm::debug()) cvm::log("Initializing dihedral object from atom groups.\n"); diff --git a/lib/colvars/colvarcomp_rotations.cpp b/lib/colvars/colvarcomp_rotations.cpp index 4bc0dcff0e..7e0cf75879 100644 --- a/lib/colvars/colvarcomp_rotations.cpp +++ b/lib/colvars/colvarcomp_rotations.cpp @@ -88,9 +88,6 @@ colvar::orientation::orientation() void colvar::orientation::calc_value() { - // atoms.reset_atoms_data(); - // atoms.read_positions(); - atoms_cog = atoms.center_of_geometry(); rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); @@ -145,9 +142,6 @@ colvar::orientation_angle::orientation_angle() void colvar::orientation_angle::calc_value() { - // atoms.reset_atoms_data(); - // atoms.read_positions(); - atoms_cog = atoms.center_of_geometry(); rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); @@ -261,9 +255,6 @@ colvar::tilt::tilt() void colvar::tilt::calc_value() { - // atoms.reset_atoms_data(); - // atoms.read_positions(); - atoms_cog = atoms.center_of_geometry(); rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); @@ -331,9 +322,6 @@ colvar::spin_angle::spin_angle() void colvar::spin_angle::calc_value() { - // atoms.reset_atoms_data(); - // atoms.read_positions(); - atoms_cog = atoms.center_of_geometry(); rot.calc_optimal_rotation(ref_pos, atoms.positions_shifted(-1.0 * atoms_cog)); diff --git a/lib/colvars/colvargrid.cpp b/lib/colvars/colvargrid.cpp index e2e660b984..522d866409 100644 --- a/lib/colvars/colvargrid.cpp +++ b/lib/colvars/colvargrid.cpp @@ -127,7 +127,6 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os) integral = 0.0; int_vals.push_back( 0.0 ); - bin = 0.0; min = 0.0; // correction for periodic colvars, so that the PMF is periodic @@ -138,7 +137,7 @@ void colvar_grid_gradient::write_1D_integral(std::ostream &os) corr = 0.0; } - for (std::vector ix = new_index(); index_ok(ix); incr(ix), bin += 1.0 ) { + for (std::vector ix = new_index(); index_ok(ix); incr(ix)) { if (samples) { size_t const samples_here = samples->value(ix); diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index 2d62cace56..f7c33c5e93 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -100,38 +100,24 @@ int colvarmodule::read_config_string(std::string const &config_str) int colvarmodule::parse_config(std::string &conf) { - int error_code = 0; - // parse global options - error_code |= parse_global_params(conf); - - if (error_code != COLVARS_OK || cvm::get_error()) { - set_error_bits(INPUT_ERROR); - return COLVARS_ERROR; + if (catch_input_errors(parse_global_params(conf))) { + return get_error(); } // parse the options for collective variables - error_code |= parse_colvars(conf); - - if (error_code != COLVARS_OK || cvm::get_error()) { - set_error_bits(INPUT_ERROR); - return COLVARS_ERROR; + if (catch_input_errors(parse_colvars(conf))) { + return get_error(); } // parse the options for biases - error_code |= parse_biases(conf); - - if (error_code != COLVARS_OK || cvm::get_error()) { - set_error_bits(INPUT_ERROR); - return COLVARS_ERROR; + if (catch_input_errors(parse_biases(conf))) { + return get_error(); } // done parsing known keywords, check that all keywords found were valid ones - error_code |= parse->check_keywords(conf, "colvarmodule"); - - if (error_code != COLVARS_OK || cvm::get_error()) { - set_error_bits(INPUT_ERROR); - return COLVARS_ERROR; + if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) { + return get_error(); } cvm::log(cvm::line_marker); @@ -143,7 +129,7 @@ int colvarmodule::parse_config(std::string &conf) write_traj_label(cv_traj_os); } - return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); + return get_error(); } @@ -181,7 +167,7 @@ int colvarmodule::parse_global_params(std::string const &conf) colvarparse::parse_silent); if (use_scripted_forces && !proxy->force_script_defined) { - cvm::fatal_error("User script for scripted colvar forces not found."); + cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); @@ -315,6 +301,17 @@ 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); + 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(); @@ -652,7 +649,6 @@ int colvarmodule::setup() (*cvi)->setup(); } return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK); - } colvarmodule::~colvarmodule() @@ -664,6 +660,8 @@ colvarmodule::~colvarmodule() 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 cvi = colvars.rbegin(); @@ -1072,6 +1070,8 @@ void cvm::error(std::string const &message, int code) 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); proxy->fatal_error(message); } @@ -1143,7 +1143,7 @@ int cvm::read_index_file(char const *filename) } int cvm::load_atoms(char const *file_name, - std::vector &atoms, + cvm::atom_group &atoms, std::string const &pdb_field, double const pdb_field_value) { diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index fbb2525983..ba3be84532 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -4,7 +4,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2015-10-22" +#define COLVARS_VERSION "2016-02-28" #endif #ifndef COLVARS_DEBUG @@ -31,6 +31,7 @@ #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 #include @@ -109,6 +110,7 @@ public: static inline void set_error_bits(int code) { errorCode |= code; + errorCode |= COLVARS_ERROR; } static inline int get_error() { @@ -227,6 +229,12 @@ public: /// 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 @@ -420,17 +428,17 @@ public: /// \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, - std::vector &atoms, - std::string const &pdb_field, - double const pdb_field_value = 0.0); + 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); + 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 diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index a97d95e36e..42d7b6e01d 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -110,8 +110,8 @@ size_t colvarparse::dummy_pos = 0; data_count++; \ } \ if (data_count == 0) \ - cvm::fatal_error("Error: in parsing \""+ \ - std::string(key)+"\".\n"); \ + cvm::error("Error: in parsing \""+ \ + std::string(key)+"\".\n", INPUT_ERROR); \ if (data_count > 1) { \ cvm::error("Error: multiple values " \ "are not allowed for keyword \""+ \ @@ -283,8 +283,8 @@ bool colvarparse::get_keyval(std::string const &conf, (data == std::string("false")) ) { value = false; } else - cvm::fatal_error("Error: boolean values only are allowed " - "for \""+std::string(key)+"\".\n"); + cvm::error("Error: boolean values only are allowed " + "for \""+std::string(key)+"\".\n", INPUT_ERROR); if (parse_mode != parse_silent) { cvm::log("# "+std::string(key)+" = "+ (value ? "on" : "off")+"\n"); @@ -398,10 +398,9 @@ int colvarparse::check_keywords(std::string &conf, char const *key) } } if (!found_keyword) { - cvm::log("Error: keyword \""+uk+"\" is not supported, " - "or not recognized in this context.\n"); - cvm::set_error_bits(INPUT_ERROR); - return COLVARS_ERROR; + cvm::error("Error: keyword \""+uk+"\" is not supported, " + "or not recognized in this context.\n", INPUT_ERROR); + return INPUT_ERROR; } } @@ -527,10 +526,10 @@ bool colvarparse::key_lookup(std::string const &conf, while (brace == std::string::npos) { // add a new line if (line_end >= conf.size()) { - cvm::fatal_error("Parse error: reached the end while " - "looking for closing brace; until now " - "the following was parsed: \"\n"+ - line+"\".\n"); + cvm::error("Parse error: reached the end while " + "looking for closing brace; until now " + "the following was parsed: \"\n"+ + line+"\".\n", INPUT_ERROR); return false; } size_t const old_end = line.size(); diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index 9e3811b660..b635711f9c 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -147,6 +147,7 @@ public: /// \brief Use this after parsing a config string (note that check_keywords() calls it already) void clear_keyword_registry(); + inline void reset() { clear_keyword_registry(); } /// \brief Return a lowercased copy of the string static inline std::string to_lower_cppstr(std::string const &in) diff --git a/lib/colvars/colvarproxy.h b/lib/colvars/colvarproxy.h index 47806990b7..d18ec3715a 100644 --- a/lib/colvars/colvarproxy.h +++ b/lib/colvars/colvarproxy.h @@ -9,16 +9,13 @@ #include "colvarmodule.h" #include "colvarvalue.h" -// return values for the frame() routine -#define COLVARS_NO_SUCH_FRAME -1 - // forward declarations class colvarscript; /// \brief Interface between the collective variables module and -/// the simulation or analysis program. -/// This is the base class: each program is supported by a derived class. -/// Only pure virtual functions ("= 0") must be reimplemented in a new interface. +/// the simulation or analysis program (NAMD, VMD, LAMMPS...). +/// This is the base class: each interfaced program is supported by a derived class. +/// Only pure virtual functions ("= 0") must be reimplemented to ensure baseline functionality. class colvarproxy { @@ -31,13 +28,29 @@ public: inline colvarproxy() : script(NULL) {} /// Default destructor - virtual inline ~colvarproxy() {} + virtual ~colvarproxy() {} - /// (Re)initialize member data after construction - virtual void setup() {} + /// (Re)initialize required member data after construction + virtual int setup() + { + return COLVARS_OK; + } + /// \brief Update data required by the colvars module (e.g. cache atom positions) + /// + /// TODO Break up colvarproxy_namd and colvarproxy_lammps function into these + virtual int update_input() + { + return COLVARS_OK; + } - // **************** SYSTEM-WIDE PHYSICAL QUANTITIES **************** + /// \brief Update data based from the results of a module update (e.g. send forces) + virtual int update_output() + { + return COLVARS_OK; + } + + // **************** SIMULATION PARAMETERS **************** /// \brief Value of the unit for atomic coordinates with respect to /// angstroms (used by some variables for hard-coded default values) @@ -62,6 +75,45 @@ public: // return 0 on success, -1 on failure virtual int frame(int) { return COLVARS_NOT_IMPLEMENTED; } + /// \brief Prefix to be used for input files (restarts, not + /// configuration) + std::string input_prefix_str, output_prefix_str, restart_output_prefix_str; + + inline std::string input_prefix() + { + return input_prefix_str; + } + + /// \brief Prefix to be used for output restart files + inline std::string restart_output_prefix() + { + return restart_output_prefix_str; + } + + /// \brief Prefix to be used for output files (final system + /// configuration) + inline std::string output_prefix() + { + return output_prefix_str; + } + + /// \brief Restarts will be written each time this number of steps has passed + virtual size_t restart_frequency() + { + return 0; + } + +protected: + + /// \brief Currently opened output files: by default, these are ofstream objects. + /// Allows redefinition to implement different output mechanisms + std::list output_files; + /// \brief Identifiers for output_stream objects: by default, these are the names of the files + std::list output_stream_names; + +public: + + // **************** MULTIPLE REPLICAS COMMUNICATION **************** // Replica exchange commands: @@ -87,68 +139,6 @@ public: return COLVARS_NOT_IMPLEMENTED; } - // **************** SIMULATION PARAMETERS **************** - - /// \brief Prefix to be used for input files (restarts, not - /// configuration) - std::string input_prefix_str, output_prefix_str, restart_output_prefix_str; - - inline std::string input_prefix() - { - return input_prefix_str; - } - - /// \brief Prefix to be used for output restart files - inline std::string restart_output_prefix() - { - return restart_output_prefix_str; - } - - /// \brief Prefix to be used for output files (final system - /// configuration) - inline std::string output_prefix() - { - return output_prefix_str; - } - - /// \brief Restarts will be written each time this number of steps has passed - virtual size_t restart_frequency() = 0; - - - - // **************** ACCESS ATOMIC DATA **************** - - /// Pass restraint energy value for current timestep to MD engine - virtual void add_energy(cvm::real energy) = 0; - - /// Tell the proxy whether system forces are needed (may not always be available) - virtual void request_system_force(bool yesno) = 0; - - - - // **************** PERIODIC BOUNDARY CONDITIONS **************** - - /// \brief Get the PBC-aware distance vector between two positions - virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) = 0; - - /// \brief Get the PBC-aware square distance between two positions; - /// may be implemented independently from position_distance() for optimization purposes - virtual cvm::real position_dist2(cvm::atom_pos const &pos1, - cvm::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 The reference position - virtual void select_closest_image(cvm::atom_pos &pos, - cvm::atom_pos const &ref_pos) = 0; - - /// \brief Perform select_closest_image() on a set of atomic positions - /// - /// After that, distance vectors can then be calculated directly, - /// without using position_distance() - void select_closest_images(std::vector &pos, - cvm::atom_pos const &ref_pos); // **************** SCRIPTING INTERFACE **************** @@ -175,6 +165,7 @@ public: std::vector > &gradient) { return COLVARS_NOT_IMPLEMENTED; } + // **************** INPUT/OUTPUT **************** /// Print a message to the main log @@ -187,36 +178,11 @@ public: virtual void fatal_error(std::string const &message) = 0; /// Print a message to the main log and exit normally - virtual void exit(std::string const &message) = 0; - - /// \brief Read atom identifiers from a file \param filename name of - /// the file (usually a PDB) \param atoms array to which atoms read - /// from "filename" will be appended \param pdb_field (optiona) if - /// "filename" is a PDB file, use this field to determine which are - /// the atoms to be set - virtual int load_atoms(char const *filename, - std::vector &atoms, - std::string const &pdb_field, - double const pdb_field_value = 0.0) = 0; - - /// \brief Load the coordinates for a group of atoms from a file - /// (usually a PDB); if "pos" is already allocated, the number of its - /// elements must match the number of atoms in "filename" - virtual 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) = 0; - -protected: - - /// \brief Open output files: by default, these are ofstream objects. - /// Allows redefinition to implement different output mechanisms - std::list output_files; - /// \brief Identifiers for output_stream objects: by default, these are the names of the files - std::list output_stream_names; - -public: + virtual void exit(std::string const &message) + { + cvm::error("Error: exiting without error is not implemented, returning error code.\n", + COLVARS_NOT_IMPLEMENTED); + } // TODO the following definitions may be moved to a .cpp file @@ -254,6 +220,8 @@ public: return COLVARS_OK; } } + cvm::error("Error: trying to close an output file or stream that wasn't open.\n", + BUG_ERROR); return COLVARS_ERROR; } @@ -262,22 +230,331 @@ public: { return COLVARS_NOT_IMPLEMENTED; } + + + + // **************** ACCESS SYSTEM DATA **************** + + /// Pass restraint energy value for current timestep to MD engine + virtual void add_energy(cvm::real energy) = 0; + + /// Tell the proxy whether system forces are needed (may not always be available) + virtual void request_system_force(bool yesno) + { + if (yesno == true) + cvm::error("Error: system forces are currently not implemented.\n", + COLVARS_NOT_IMPLEMENTED); + } + + /// \brief Get the PBC-aware distance vector between two positions + virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) = 0; + + /// \brief Get the PBC-aware square distance between two positions; + /// may need to be reimplemented independently from position_distance() for optimization purposes + virtual cvm::real position_dist2(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) + { + return (position_distance(pos1, pos2)).norm2(); + } + + /// \brief Get the closest periodic image to a reference position + /// \param pos The position to look for the closest periodic image + /// \param ref_pos The reference position + virtual void select_closest_image(cvm::atom_pos &pos, + cvm::atom_pos const &ref_pos) + { + pos = position_distance(ref_pos, pos) + 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() + void select_closest_images(std::vector &pos, + cvm::atom_pos const &ref_pos) + { + for (std::vector::iterator pi = pos.begin(); + pi != pos.end(); ++pi) { + select_closest_image(*pi, ref_pos); + } + } + + + // **************** ACCESS ATOMIC DATA **************** +protected: + + /// \brief Array of 0-based integers used to uniquely associate atoms + /// within the host program + std::vector atoms_ids; + /// \brief Keep track of how many times each atom is used by a separate colvar object + std::vector atoms_ncopies; + /// \brief Masses of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) + std::vector atoms_masses; + /// \brief Charges of the atoms (allow redefinition during a run, as done e.g. in LAMMPS) + std::vector atoms_charges; + /// \brief Current three-dimensional positions of the atoms + std::vector atoms_positions; + /// \brief Most recent total forces on each atom + std::vector atoms_total_forces; + /// \brief Most recent forces applied by external potentials onto each atom + std::vector atoms_applied_forces; + /// \brief Forces applied from colvars, to be communicated to the MD integrator + std::vector atoms_new_colvar_forces; + + /// Used by all init_atom() functions: create a slot for an atom not requested yet + inline int add_atom_slot(int atom_id) + { + atoms_ids.push_back(atom_id); + atoms_ncopies.push_back(1); + atoms_masses.push_back(1.0); + atoms_charges.push_back(0.0); + atoms_positions.push_back(cvm::rvector(0.0, 0.0, 0.0)); + atoms_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); + atoms_applied_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); + atoms_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); + return (atoms_ids.size() - 1); + } + +public: + + /// Prepare this atom for collective variables calculation, selecting it by numeric index (1-based) + virtual int init_atom(int atom_number) = 0; + + /// Check that this atom number is valid, but do not initialize the corresponding atom yet + virtual int check_atom_id(int atom_number) = 0; + + /// Select this atom for collective variables calculation, using name and residue number. + /// Not all programs support this: leave this function as is in those cases. + virtual int init_atom(cvm::residue_id const &residue, + std::string const &atom_name, + std::string const &segment_id) + { + cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n", + COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; + } + + /// Check that this atom is valid, but do not initialize it yet + virtual int check_atom_id(cvm::residue_id const &residue, + std::string const &atom_name, + std::string const &segment_id) + { + cvm::error("Error: initializing an atom by name and residue number is currently not supported.\n", + COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; + } + + /// \brief Used by the atom class destructor: rather than deleting the array slot + /// (costly) set the corresponding atoms_ncopies to zero + virtual void clear_atom(int index) + { + if (((size_t) index) >= atoms_ids.size()) { + cvm::error("Error: trying to disable an atom that was not previously requested.\n", + INPUT_ERROR); + } + if (atoms_ncopies[index] > 0) { + atoms_ncopies[index] -= 1; + } + } + + /// Get the numeric ID of the given atom (for the program) + inline int get_atom_id(int index) const + { + return atoms_ids[index]; + } + + /// Get the mass of the given atom + inline cvm::real get_atom_mass(int index) const + { + return atoms_masses[index]; + } + + /// Get the charge of the given atom + inline cvm::real get_atom_charge(int index) const + { + return atoms_charges[index]; + } + + /// Read the current position of the given atom + inline cvm::rvector get_atom_position(int index) const + { + return atoms_positions[index]; + } + + /// Read the current total system force of the given atom + inline cvm::rvector get_atom_system_force(int index) const + { + return atoms_total_forces[index] - atoms_applied_forces[index]; + } + + /// Request that this force is applied to the given atom + inline void apply_atom_force(int index, cvm::rvector const &new_force) + { + atoms_new_colvar_forces[index] += new_force; + } + + /// Read the current velocity of the given atom + virtual cvm::rvector get_atom_velocity(int index) + { + cvm::error("Error: reading the current velocity of an atom is not yet implemented.\n", + COLVARS_NOT_IMPLEMENTED); + return cvm::rvector(0.0); + } + + // useful functions for data management outside this class + inline std::vector *modify_atom_ids() { return &atoms_ids; } + inline std::vector *modify_atom_masses() { return &atoms_masses; } + inline std::vector *modify_atom_charges() { return &atoms_charges; } + inline std::vector *modify_atom_positions() { return &atoms_positions; } + inline std::vector *modify_atom_total_forces() { return &atoms_total_forces; } + inline std::vector *modify_atom_applied_forces() { return &atoms_applied_forces; } + inline std::vector *modify_atom_new_colvar_forces() { return &atoms_new_colvar_forces; } + + /// \brief Read atom identifiers from a file \param filename name of + /// the file (usually a PDB) \param atoms array to which atoms read + /// from "filename" will be appended \param pdb_field (optiona) if + /// "filename" is a PDB file, use this field to determine which are + /// the atoms to be set + virtual int load_atoms(char const *filename, + cvm::atom_group &atoms, + std::string const &pdb_field, + double const pdb_field_value = 0.0) + { + cvm::error("Error: loading atom identifiers from a file is currently not implemented.\n", + COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; + } + + /// \brief Load the coordinates for a group of atoms from a file + /// (usually a PDB); if "pos" is already allocated, the number of its + /// elements must match the number of atoms in "filename" + virtual 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) + { + cvm::error("Error: loading atomic coordinates from a file is currently not implemented.\n"); + return COLVARS_NOT_IMPLEMENTED; + } + + // **************** ACCESS GROUP DATA **************** + +protected: + + /// \brief Array of 0-based integers used to uniquely associate atom groups + /// within the host program + std::vector atom_groups_ids; + /// \brief Keep track of how many times each group is used by a separate cvc + std::vector atom_groups_ncopies; + /// \brief Total masses of the atom groups + std::vector atom_groups_masses; + /// \brief Total charges of the atom groups (allow redefinition during a run, as done e.g. in LAMMPS) + std::vector atom_groups_charges; + /// \brief Current centers of mass of the atom groups + std::vector atom_groups_coms; + /// \brief Most recently updated total forces on the com of each group + std::vector atom_groups_total_forces; + /// \brief Most recent forces applied by external potentials onto each group + std::vector atom_groups_applied_forces; + /// \brief Forces applied from colvars, to be communicated to the MD integrator + std::vector atom_groups_new_colvar_forces; + + /// TODO Add here containers of handles to cvc objects that are computed in parallel + +public: + + /// \brief Whether this proxy implementation has capability for scalable groups + virtual bool has_scalable_groups() const + { + return false; + } + + /// Used by all init_atom_group() functions: create a slot for an atom group not requested yet + // TODO Add a handle to cvc objects + inline int add_atom_group_slot(int atom_group_id) + { + atom_groups_ids.push_back(atom_group_id); + atom_groups_masses.push_back(1.0); + atom_groups_charges.push_back(0.0); + atom_groups_coms.push_back(cvm::rvector(0.0, 0.0, 0.0)); + atom_groups_total_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); + atom_groups_applied_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); + atom_groups_new_colvar_forces.push_back(cvm::rvector(0.0, 0.0, 0.0)); + return (atom_groups_ids.size() - 1); + } + + /// Prepare this group for collective variables calculation, selecting atoms by internal ids (0-based) + virtual int init_atom_group(std::vector const &atoms_ids) // TODO Add a handle to cvc objects + { + cvm::error("Error: initializing a group outside of the colvars module is currently not supported.\n", + COLVARS_NOT_IMPLEMENTED); + return COLVARS_NOT_IMPLEMENTED; + } + + /// \brief Used by the atom_group class destructor + virtual void clear_atom_group(int index) + { + if (((size_t) index) >= atom_groups_ids.size()) { + cvm::error("Error: trying to disable an atom group that was not previously requested.\n", + INPUT_ERROR); + } + + atom_groups_ids.erase(atom_groups_ids.begin()+index); + atom_groups_masses.erase(atom_groups_masses.begin()+index); + atom_groups_charges.erase(atom_groups_charges.begin()+index); + atom_groups_coms.erase(atom_groups_coms.begin()+index); + atom_groups_total_forces.erase(atom_groups_total_forces.begin()+index); + atom_groups_applied_forces.erase(atom_groups_applied_forces.begin()+index); + atom_groups_new_colvar_forces.erase(atom_groups_new_colvar_forces.begin()+index); + } + + /// Get the numeric ID of the given atom group (for the MD program) + inline cvm::real get_atom_group_id(int index) const + { + return atom_groups_ids[index]; + } + + /// Get the mass of the given atom group + inline cvm::real get_atom_group_mass(int index) const + { + return atom_groups_masses[index]; + } + + /// Get the charge of the given atom group + inline cvm::real get_atom_group_charge(int index) const + { + return atom_groups_charges[index]; + } + + /// Read the current position of the center of mass given atom group + inline cvm::rvector get_atom_group_com(int index) const + { + return atom_groups_coms[index]; + } + + /// Read the current total system force of the given atom group + inline cvm::rvector get_atom_group_system_force(int index) const + { + return atom_groups_total_forces[index] - atom_groups_applied_forces[index]; + } + + /// Request that this force is applied to the given atom group + inline void apply_atom_group_force(int index, cvm::rvector const &new_force) + { + atom_groups_new_colvar_forces[index] += new_force; + } + + /// Read the current velocity of the given atom group + virtual cvm::rvector get_atom_group_velocity(int index) + { + cvm::error("Error: reading the current velocity of an atom group is not yet implemented.\n", + COLVARS_NOT_IMPLEMENTED); + return cvm::rvector(0.0); + } + }; -inline void colvarproxy::select_closest_images(std::vector &pos, - cvm::atom_pos const &ref_pos) -{ - for (std::vector::iterator pi = pos.begin(); - pi != pos.end(); ++pi) { - select_closest_image(*pi, ref_pos); - } -} - -inline cvm::real colvarproxy::position_dist2(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) -{ - return (position_distance(pos1, pos2)).norm2(); -} - #endif diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index 4c999eee84..efc6731ef7 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -26,11 +26,13 @@ int colvarscript::run(int argc, char const *argv[]) { if (argc < 2) { result = help_string(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } std::string cmd = argv[1]; + int error_code = COLVARS_OK; + if (cmd == "colvar") { return proc_colvar(argc-1, &(argv[1])); } @@ -41,13 +43,13 @@ int colvarscript::run(int argc, char const *argv[]) { if (cmd == "version") { result = COLVARS_VERSION; - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (cmd == "reset") { /// Delete every child object colvars->reset(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (cmd == "delete") { @@ -55,12 +57,17 @@ int colvarscript::run(int argc, char const *argv[]) { // Note: the delete bit may be ignored by some backends // it is mostly useful in VMD colvars->set_error_bits(DELETE_COLVARS); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (cmd == "update") { - colvars->calc(); - return COLVARSCRIPT_OK; + error_code |= proxy->update_input(); + error_code |= colvars->calc(); + error_code |= proxy->update_output(); + if (error_code) { + result += "Error updating the colvars module.\n"; + } + return error_code; } if (cmd == "list") { @@ -70,14 +77,14 @@ int colvarscript::run(int argc, char const *argv[]) { ++cvi) { result += (cvi == colvars->colvars.begin() ? "" : " ") + (*cvi)->name; } - return COLVARSCRIPT_OK; + return COLVARS_OK; } else if (argc == 3 && !strcmp(argv[2], "biases")) { for (std::vector::iterator bi = colvars->biases.begin(); bi != colvars->biases.end(); ++bi) { result += (bi == colvars->biases.begin() ? "" : " ") + (*bi)->name; } - return COLVARSCRIPT_OK; + return COLVARS_OK; } else { result = "Wrong arguments to command \"list\"\n" + help_string(); return COLVARSCRIPT_ERROR; @@ -91,7 +98,7 @@ int colvarscript::run(int argc, char const *argv[]) { return COLVARSCRIPT_ERROR; } if (colvars->read_config_file(argv[2]) == COLVARS_OK) { - return COLVARSCRIPT_OK; + return COLVARS_OK; } else { result = "Error parsing configuration file"; return COLVARSCRIPT_ERROR; @@ -106,7 +113,7 @@ int colvarscript::run(int argc, char const *argv[]) { } std::string conf = argv[2]; if (colvars->read_config_string(conf) == COLVARS_OK) { - return COLVARSCRIPT_OK; + return COLVARS_OK; } else { result = "Error parsing configuration string"; return COLVARSCRIPT_ERROR; @@ -121,7 +128,7 @@ int colvarscript::run(int argc, char const *argv[]) { } proxy->input_prefix_str = argv[2]; if (colvars->setup_input() == COLVARS_OK) { - return COLVARSCRIPT_OK; + return COLVARS_OK; } else { result = "Error loading state file"; return COLVARSCRIPT_ERROR; @@ -138,7 +145,7 @@ int colvarscript::run(int argc, char const *argv[]) { int error = 0; error |= colvars->setup_output(); error |= colvars->write_output_files(); - return error ? COLVARSCRIPT_ERROR : COLVARSCRIPT_OK; + return error ? COLVARSCRIPT_ERROR : COLVARS_OK; } /// Print the values that would go on colvars.traj @@ -146,13 +153,13 @@ int colvarscript::run(int argc, char const *argv[]) { std::ostringstream os; colvars->write_traj_label(os); result = os.str(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (cmd == "printframe") { std::ostringstream os; colvars->write_traj(os); result = os.str(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (cmd == "frame") { @@ -160,7 +167,7 @@ int colvarscript::run(int argc, char const *argv[]) { int f = proxy->frame(); if (f >= 0) { result = cvm::to_str(f); - return COLVARSCRIPT_OK; + return COLVARS_OK; } else { result = "Frame number is not available"; return COLVARSCRIPT_ERROR; @@ -171,7 +178,7 @@ int colvarscript::run(int argc, char const *argv[]) { long int f = proxy->frame(strtol(argv[2], NULL, 10)); colvars->it = proxy->frame(); result = cvm::to_str(f); - return COLVARSCRIPT_OK; + return COLVARS_OK; } else { result = "Wrong arguments to command \"frame\"\n" + help_string(); return COLVARSCRIPT_ERROR; @@ -199,24 +206,24 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { if (subcmd == "value") { result = (cv->value()).to_simple_string(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "width") { result = cvm::to_str(cv->width, 0, cvm::cv_prec); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "type") { result = cv->value().type_desc(cv->value().value_type); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "update") { cv->calc(); cv->update(); result = (cv->value()).to_simple_string(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "delete") { @@ -228,12 +235,12 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { delete cv; // TODO this could be done by the destructors colvars->write_traj_label(colvars->cv_traj_os); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "getconfig") { result = cv->get_config(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "addforce") { @@ -253,7 +260,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { } cv->add_bias_force(force); result = force.to_simple_string(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "cvcflags") { @@ -276,7 +283,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { return COLVARSCRIPT_ERROR; } result = "0"; - return COLVARSCRIPT_OK; + return COLVARS_OK; } result = "Syntax error\n" + help_string(); @@ -301,25 +308,25 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { if (subcmd == "energy") { result = cvm::to_str(b->get_energy()); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "update") { b->update(); result = cvm::to_str(b->get_energy()); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "getconfig") { result = b->get_config(); - return COLVARSCRIPT_OK; + return COLVARS_OK; } // Subcommands for MW ABF if (subcmd == "bin") { int r = b->current_bin(); result = cvm::to_str(r); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "binnum") { @@ -329,7 +336,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { return COLVARSCRIPT_ERROR; } result = cvm::to_str(r); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (subcmd == "share") { @@ -339,7 +346,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { return COLVARSCRIPT_ERROR; } result = cvm::to_str(r); - return COLVARSCRIPT_OK; + return COLVARS_OK; } // End commands for MW ABF @@ -348,7 +355,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { delete b; // TODO this could be done by the destructors colvars->write_traj_label(colvars->cv_traj_os); - return COLVARSCRIPT_OK; + return COLVARS_OK; } if (argc >= 4) { @@ -360,7 +367,7 @@ int colvarscript::proc_bias(int argc, char const *argv[]) { return COLVARSCRIPT_ERROR; } result = cvm::to_str(b->bin_count(index)); - return COLVARSCRIPT_OK; + return COLVARS_OK; } result = "Syntax error\n" + help_string(); diff --git a/lib/colvars/colvarscript.h b/lib/colvars/colvarscript.h index 745914c76b..25138387ca 100644 --- a/lib/colvars/colvarscript.h +++ b/lib/colvars/colvarscript.h @@ -9,6 +9,7 @@ #include "colvarbias.h" #include "colvarproxy.h" +// TODO merge these into colvarmodule.h #define COLVARSCRIPT_ERROR -1 #define COLVARSCRIPT_OK 0