colvaratoms.C

00001 // -*- c++ -*-
00002 
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/Colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 #include <list>
00011 #include <vector>
00012 #include <algorithm>
00013 #include <sstream>
00014 #include <iomanip>
00015 
00016 #include "colvarmodule.h"
00017 #include "colvarproxy.h"
00018 #include "colvarparse.h"
00019 #include "colvaratoms.h"
00020 
00021 
00022 cvm::atom::atom()
00023 {
00024 index = -1;
00025 id = -1;
00026 mass = 1.0;
00027 charge = 0.0;
00028 reset_data();
00029 }
00030 
00031 
00032 cvm::atom::atom(int atom_number)
00033 {
00034 colvarproxy *p = cvm::main()->proxy;
00035 index = p->init_atom(atom_number);
00036 id = p->get_atom_id(index);
00037 update_mass();
00038 update_charge();
00039 reset_data();
00040 }
00041 
00042 
00043 cvm::atom::atom(cvm::residue_id const &residue,
00044 std::string const &atom_name,
00045 std::string const &segment_id)
00046 {
00047 colvarproxy *p = cvm::main()->proxy;
00048 index = p->init_atom(residue, atom_name, segment_id);
00049 id = p->get_atom_id(index);
00050 update_mass();
00051 update_charge();
00052 reset_data();
00053 }
00054 
00055 
00056 cvm::atom::atom(atom const &a)
00057 : index(a.index)
00058 {
00059 colvarproxy *p = cvm::main()->proxy;
00060 id = p->get_atom_id(index);
00061 p->increase_refcount(index);
00062 update_mass();
00063 update_charge();
00064 reset_data();
00065 }
00066 
00067 
00068 cvm::atom::~atom()
00069 {
00070 if (index >= 0) {
00071 (cvm::main()->proxy)->clear_atom(index);
00072 }
00073 }
00074 
00075 
00076 cvm::atom & cvm::atom::operator = (cvm::atom const &a)
00077 {
00078 index = a.index;
00079 id = (cvm::main()->proxy)->get_atom_id(index);
00080 update_mass();
00081 update_charge();
00082 reset_data();
00083 return *this;
00084 }
00085 
00086 
00087 
00088 cvm::atom_group::atom_group()
00089 {
00090 init();
00091 }
00092 
00093 
00094 cvm::atom_group::atom_group(char const *key_in)
00095 {
00096 key = key_in;
00097 init();
00098 }
00099 
00100 
00101 cvm::atom_group::atom_group(std::vector<cvm::atom> const &atoms_in)
00102 {
00103 init();
00104 atoms = atoms_in;
00105 setup();
00106 }
00107 
00108 
00109 cvm::atom_group::~atom_group()
00110 {
00111 if (is_enabled(f_ag_scalable) && !b_dummy) {
00112 (cvm::main()->proxy)->clear_atom_group(index);
00113 index = -1;
00114 }
00115 
00116 if (fitting_group) {
00117 delete fitting_group;
00118 fitting_group = NULL;
00119 }
00120 
00121 cvm::main()->unregister_named_atom_group(this);
00122 }
00123 
00124 
00125 int cvm::atom_group::add_atom(cvm::atom const &a)
00126 {
00127 if (a.id < 0) {
00128 return COLVARS_ERROR;
00129 }
00130 
00131 for (size_t i = 0; i < atoms_ids.size(); i++) {
00132 if (atoms_ids[i] == a.id) {
00133 if (cvm::debug())
00134 cvm::log("Discarding doubly counted atom with number "+
00135 cvm::to_str(a.id+1)+".\n");
00136 return COLVARS_OK;
00137 }
00138 }
00139 
00140 // for consistency with add_atom_id(), we update the list as well
00141 atoms_ids.push_back(a.id);
00142 atoms.push_back(a);
00143 total_mass += a.mass;
00144 total_charge += a.charge;
00145 
00146 return COLVARS_OK;
00147 }
00148 
00149 
00150 int cvm::atom_group::add_atom_id(int aid)
00151 {
00152 if (aid < 0) {
00153 return COLVARS_ERROR;
00154 }
00155 
00156 for (size_t i = 0; i < atoms_ids.size(); i++) {
00157 if (atoms_ids[i] == aid) {
00158 if (cvm::debug())
00159 cvm::log("Discarding doubly counted atom with number "+
00160 cvm::to_str(aid+1)+".\n");
00161 return COLVARS_OK;
00162 }
00163 }
00164 
00165 atoms_ids.push_back(aid);
00166 return COLVARS_OK;
00167 }
00168 
00169 
00170 int cvm::atom_group::remove_atom(cvm::atom_iter ai)
00171 {
00172 if (is_enabled(f_ag_scalable)) {
00173 cvm::error("Error: cannot remove atoms from a scalable group.\n", COLVARS_INPUT_ERROR);
00174 return COLVARS_ERROR;
00175 }
00176 
00177 if (!this->size()) {
00178 cvm::error("Error: trying to remove an atom from an empty group.\n", COLVARS_INPUT_ERROR);
00179 return COLVARS_ERROR;
00180 } else {
00181 total_mass -= ai->mass;
00182 total_charge -= ai->charge;
00183 atoms_ids.erase(atoms_ids.begin() + (ai - atoms.begin()));
00184 atoms.erase(ai);
00185 }
00186 
00187 return COLVARS_OK;
00188 }
00189 
00190 
00191 int cvm::atom_group::set_dummy()
00192 {
00193 if (atoms_ids.size() > 0) {
00194 return cvm::error("Error: setting group with keyword \""+key+
00195 "\" and name \""+name+"\" as dummy, but it already "
00196 "contains atoms.\n", COLVARS_INPUT_ERROR);
00197 }
00198 b_dummy = true;
00199 return COLVARS_OK;
00200 }
00201 
00202 
00203 int cvm::atom_group::set_dummy_pos(cvm::atom_pos const &pos)
00204 {
00205 if (b_dummy) {
00206 dummy_atom_pos = pos;
00207 } else {
00208 return cvm::error("Error: setting dummy position for group with keyword \""+
00209 key+"\" and name \""+name+
00210 "\", but it is not dummy.\n", COLVARS_INPUT_ERROR);
00211 }
00212 return COLVARS_OK;
00213 }
00214 
00215 
00216 int cvm::atom_group::init()
00217 {
00218 if (!key.size()) key = "unnamed";
00219 description = "atom group " + key;
00220 // These may be overwritten by parse(), if a name is provided
00221 
00222 atoms.clear();
00223 atom_group::init_dependencies();
00224 index = -1;
00225 
00226 b_dummy = false;
00227 b_user_defined_fit = false;
00228 fitting_group = NULL;
00229 
00230 noforce = false;
00231 
00232 total_mass = 0.0;
00233 total_charge = 0.0;
00234 
00235 cog.reset();
00236 com.reset();
00237 
00238 return COLVARS_OK;
00239 }
00240 
00241 
00242 int cvm::atom_group::init_dependencies() {
00243 size_t i;
00244 // Initialize static array once and for all
00245 if (features().size() == 0) {
00246 for (i = 0; i < f_ag_ntot; i++) {
00247 modify_features().push_back(new feature);
00248 }
00249 
00250 init_feature(f_ag_active, "active", f_type_dynamic);
00251 init_feature(f_ag_center, "center_to_reference", f_type_user);
00252 init_feature(f_ag_center_origin, "center_to_origin", f_type_user);
00253 init_feature(f_ag_rotate, "rotate_to_origin", f_type_user);
00254 init_feature(f_ag_fitting_group, "fitting_group", f_type_static);
00255 init_feature(f_ag_explicit_gradient, "explicit_atom_gradient", f_type_dynamic);
00256 init_feature(f_ag_fit_gradients, "fit_gradients", f_type_user);
00257 require_feature_self(f_ag_fit_gradients, f_ag_explicit_gradient);
00258 
00259 init_feature(f_ag_atom_forces, "atomic_forces", f_type_dynamic);
00260 
00261 // parallel calculation implies that we have at least a scalable center of mass,
00262 // but f_ag_scalable is kept as a separate feature to deal with future dependencies
00263 init_feature(f_ag_scalable, "scalable_group", f_type_dynamic);
00264 init_feature(f_ag_scalable_com, "scalable_group_center_of_mass", f_type_static);
00265 require_feature_self(f_ag_scalable_com, f_ag_scalable);
00266 
00267 init_feature(f_ag_collect_atom_ids, "collect_atom_ids", f_type_dynamic);
00268 exclude_feature_self(f_ag_collect_atom_ids, f_ag_scalable);
00269 
00270 // check that everything is initialized
00271 for (i = 0; i < colvardeps::f_ag_ntot; i++) {
00272 if (is_not_set(i)) {
00273 cvm::error("Uninitialized feature " + cvm::to_str(i) + " in " + description);
00274 }
00275 }
00276 }
00277 
00278 // Initialize feature_states for each instance
00279 // default as unavailable, not enabled
00280 feature_states.reserve(f_ag_ntot);
00281 for (i = 0; i < colvardeps::f_ag_ntot; i++) {
00282 feature_states.push_back(feature_state(false, false));
00283 }
00284 
00285 // Features that are implemented (or not) by all atom groups
00286 feature_states[f_ag_active].available = true;
00287 feature_states[f_ag_center].available = true;
00288 feature_states[f_ag_center_origin].available = true;
00289 feature_states[f_ag_rotate].available = true;
00290 
00291 // f_ag_scalable_com is provided by the CVC iff it is COM-based
00292 feature_states[f_ag_scalable_com].available = false;
00293 feature_states[f_ag_scalable].available = true;
00294 feature_states[f_ag_fit_gradients].available = true;
00295 feature_states[f_ag_fitting_group].available = true;
00296 feature_states[f_ag_explicit_gradient].available = true;
00297 feature_states[f_ag_collect_atom_ids].available = true;
00298 
00299 return COLVARS_OK;
00300 }
00301 
00302 
00303 int cvm::atom_group::setup()
00304 {
00305 if (atoms_ids.size() == 0) {
00306 atoms_ids.reserve(atoms.size());
00307 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00308 atoms_ids.push_back(ai->id);
00309 }
00310 }
00311 for (cvm::atom_iter ai = atoms.begin(); ai != atoms.end(); ai++) {
00312 ai->update_mass();
00313 ai->update_charge();
00314 }
00315 update_total_mass();
00316 update_total_charge();
00317 return COLVARS_OK;
00318 }
00319 
00320 
00321 void cvm::atom_group::update_total_mass()
00322 {
00323 if (b_dummy) {
00324 total_mass = 1.0;
00325 return;
00326 }
00327 
00328 if (is_enabled(f_ag_scalable)) {
00329 total_mass = (cvm::main()->proxy)->get_atom_group_mass(index);
00330 } else {
00331 total_mass = 0.0;
00332 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
00333 total_mass += ai->mass;
00334 }
00335 }
00336 if (total_mass < 1e-15) {
00337 cvm::error("ERROR: " + description + " has zero total mass.\n");
00338 }
00339 }
00340 
00341 
00342 void cvm::atom_group::update_total_charge()
00343 {
00344 if (b_dummy) {
00345 total_charge = 0.0;
00346 return;
00347 }
00348 
00349 if (is_enabled(f_ag_scalable)) {
00350 total_charge = (cvm::main()->proxy)->get_atom_group_charge(index);
00351 } else {
00352 total_charge = 0.0;
00353 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
00354 total_charge += ai->charge;
00355 }
00356 }
00357 }
00358 
00359 
00360 void cvm::atom_group::print_properties(std::string const &colvar_name,
00361 int i, int j)
00362 {
00363 if (cvm::proxy->updated_masses() && cvm::proxy->updated_charges()) {
00364 cvm::log("Re-initialized atom group for variable \""+colvar_name+"\":"+
00365 cvm::to_str(i)+"/"+
00366 cvm::to_str(j)+". "+ cvm::to_str(atoms_ids.size())+
00367 " atoms: total mass = "+cvm::to_str(total_mass)+
00368 ", total charge = "+cvm::to_str(total_charge)+".\n");
00369 }
00370 }
00371 
00372 
00373 int cvm::atom_group::parse(std::string const &group_conf)
00374 {
00375 cvm::log("Initializing atom group \""+key+"\".\n");
00376 
00377 // whether or not to include messages in the log
00378 // colvarparse::Parse_Mode mode = parse_silent;
00379 // {
00380 // bool b_verbose;
00381 // get_keyval (group_conf, "verboseOutput", b_verbose, false, parse_silent);
00382 // if (b_verbose) mode = parse_normal;
00383 // }
00384 // colvarparse::Parse_Mode mode = parse_normal;
00385 
00386 int parse_error = COLVARS_OK;
00387 
00388 // Optional group name will let other groups reuse atom definition
00389 if (get_keyval(group_conf, "name", name)) {
00390 if ((cvm::atom_group_by_name(this->name) != NULL) &&
00391 (cvm::atom_group_by_name(this->name) != this)) {
00392 cvm::error("Error: this atom group cannot have the same name, \""+this->name+
00393 "\", as another atom group.\n",
00394 COLVARS_INPUT_ERROR);
00395 return COLVARS_INPUT_ERROR;
00396 }
00397 cvm::main()->register_named_atom_group(this);
00398 description = "atom group " + name;
00399 }
00400 
00401 // We need to know about fitting to decide whether the group is scalable
00402 // and we need to know about scalability before adding atoms
00403 bool b_defined_center = get_keyval_feature(this, group_conf, "centerToOrigin", f_ag_center_origin, false);
00404 // Legacy alias
00405 b_defined_center |= get_keyval_feature(this, group_conf, "centerReference", f_ag_center, is_enabled(f_ag_center_origin), parse_deprecated);
00406 b_defined_center |= get_keyval_feature(this, group_conf, "centerToReference", f_ag_center, is_enabled(f_ag_center));
00407 
00408 if (is_enabled(f_ag_center_origin) && ! is_enabled(f_ag_center)) {
00409 return cvm::error("centerToReference may not be disabled if centerToOrigin"
00410 "is enabled.\n", COLVARS_INPUT_ERROR);
00411 }
00412 // Legacy alias
00413 bool b_defined_rotate = get_keyval_feature(this, group_conf, "rotateReference", f_ag_rotate, false, parse_deprecated);
00414 b_defined_rotate |= get_keyval_feature(this, group_conf, "rotateToReference", f_ag_rotate, is_enabled(f_ag_rotate));
00415 
00416 if (is_enabled(f_ag_rotate) || is_enabled(f_ag_center) ||
00417 is_enabled(f_ag_center_origin)) {
00418 cvm::main()->cite_feature("Moving frame of reference");
00419 }
00420 
00421 // is the user setting explicit options?
00422 b_user_defined_fit = b_defined_center || b_defined_rotate;
00423 
00424 if (is_available(f_ag_scalable_com) && !is_enabled(f_ag_rotate) && !is_enabled(f_ag_center)) {
00425 enable(f_ag_scalable_com);
00426 }
00427 
00428 {
00429 std::string atoms_of = "";
00430 if (get_keyval(group_conf, "atomsOfGroup", atoms_of)) {
00431 atom_group * ag = atom_group_by_name(atoms_of);
00432 if (ag == NULL) {
00433 cvm::error("Error: cannot find atom group with name " + atoms_of + ".\n");
00434 return COLVARS_ERROR;
00435 }
00436 parse_error |= add_atoms_of_group(ag);
00437 }
00438 }
00439 
00440 // if (get_keyval(group_conf, "copyOfGroup", source)) {
00441 // // Goal: Initialize this as a full copy
00442 // // for this we'll need an atom_group copy constructor
00443 // return COLVARS_OK;
00444 // }
00445 
00446 {
00447 std::string numbers_conf = "";
00448 size_t pos = 0;
00449 while (key_lookup(group_conf, "atomNumbers", &numbers_conf, &pos)) {
00450 parse_error |= add_atom_numbers(numbers_conf);
00451 numbers_conf = "";
00452 }
00453 }
00454 
00455 {
00456 std::string index_group_name;
00457 if (get_keyval(group_conf, "indexGroup", index_group_name)) {
00458 // use an index group from the index file read globally
00459 parse_error |= add_index_group(index_group_name);
00460 }
00461 }
00462 
00463 {
00464 std::string range_conf = "";
00465 size_t pos = 0;
00466 while (key_lookup(group_conf, "atomNumbersRange",
00467 &range_conf, &pos)) {
00468 parse_error |= add_atom_numbers_range(range_conf);
00469 range_conf = "";
00470 }
00471 }
00472 
00473 {
00474 std::vector<std::string> psf_segids;
00475 get_keyval(group_conf, "psfSegID", psf_segids, std::vector<std::string>());
00476 std::vector<std::string>::iterator psii;
00477 for (psii = psf_segids.begin(); psii < psf_segids.end(); ++psii) {
00478 if ( (psii->size() == 0) || (psii->size() > 4) ) {
00479 cvm::error("Error: invalid PSF segment identifier provided, \""+
00480 (*psii)+"\".\n", COLVARS_INPUT_ERROR);
00481 }
00482 }
00483 
00484 std::string range_conf = "";
00485 size_t pos = 0;
00486 size_t range_count = 0;
00487 psii = psf_segids.begin();
00488 while (key_lookup(group_conf, "atomNameResidueRange",
00489 &range_conf, &pos)) {
00490 range_count++;
00491 if (psf_segids.size() && (range_count > psf_segids.size())) {
00492 cvm::error("Error: more instances of \"atomNameResidueRange\" than "
00493 "values of \"psfSegID\".\n", COLVARS_INPUT_ERROR);
00494 } else {
00495 parse_error |= add_atom_name_residue_range(psf_segids.size() ?
00496 *psii : std::string(""), range_conf);
00497 if (psf_segids.size()) psii++;
00498 }
00499 range_conf = "";
00500 }
00501 }
00502 
00503 {
00504 // read the atoms from a file
00505 std::string atoms_file_name;
00506 if (get_keyval(group_conf, "atomsFile", atoms_file_name, std::string(""))) {
00507 
00508 std::string atoms_col;
00509 if (!get_keyval(group_conf, "atomsCol", atoms_col, std::string(""))) {
00510 cvm::error("Error: parameter atomsCol is required if atomsFile is set.\n",
00511 COLVARS_INPUT_ERROR);
00512 }
00513 
00514 double atoms_col_value;
00515 bool const atoms_col_value_defined = get_keyval(group_conf, "atomsColValue", atoms_col_value, 0.0);
00516 if (atoms_col_value_defined && (!atoms_col_value)) {
00517 cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", COLVARS_INPUT_ERROR);
00518 }
00519 
00520 // NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
00521 parse_error |= cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
00522 }
00523 }
00524 
00525 // Catch any errors from all the initialization steps above
00526 if (parse_error || cvm::get_error()) return (parse_error || cvm::get_error());
00527 
00528 // checks of doubly-counted atoms have been handled by add_atom() already
00529 
00530 if (get_keyval(group_conf, "dummyAtom", dummy_atom_pos, cvm::atom_pos())) {
00531 
00532 parse_error |= set_dummy();
00533 parse_error |= set_dummy_pos(dummy_atom_pos);
00534 
00535 } else {
00536 
00537 if (!(atoms_ids.size())) {
00538 parse_error |= cvm::error("Error: no atoms defined for atom group \""+
00539 key+"\".\n", COLVARS_INPUT_ERROR);
00540 }
00541 
00542 // whether these atoms will ever receive forces or not
00543 bool enable_forces = true;
00544 get_keyval(group_conf, "enableForces", enable_forces, true, colvarparse::parse_silent);
00545 noforce = !enable_forces;
00546 }
00547 
00548 // Now that atoms are defined we can parse the detailed fitting options
00549 parse_error |= parse_fitting_options(group_conf);
00550 
00551 if (is_enabled(f_ag_scalable) && !b_dummy) {
00552 cvm::log("Enabling scalable calculation for group \""+this->key+"\".\n");
00553 index = (cvm::proxy)->init_atom_group(atoms_ids);
00554 }
00555 
00556 bool b_print_atom_ids = false;
00557 get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false);
00558 
00559 // Calculate all required properties (such as total mass)
00560 setup();
00561 
00562 if (cvm::debug())
00563 cvm::log("Done initializing atom group \""+key+"\".\n");
00564 
00565 {
00566 std::string init_msg;
00567 init_msg.append("Atom group \""+key+"\" defined with "+
00568 cvm::to_str(atoms_ids.size())+" atoms requested");
00569 if ((cvm::proxy)->updated_masses()) {
00570 init_msg.append(": total mass = "+
00571 cvm::to_str(total_mass));
00572 if ((cvm::proxy)->updated_charges()) {
00573 init_msg.append(", total charge = "+
00574 cvm::to_str(total_charge));
00575 }
00576 }
00577 init_msg.append(".\n");
00578 cvm::log(init_msg);
00579 }
00580 
00581 if (b_print_atom_ids) {
00582 cvm::log("Internal definition of the atom group:\n");
00583 cvm::log(print_atom_ids());
00584 }
00585 
00586 return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
00587 }
00588 
00589 
00590 int cvm::atom_group::add_atoms_of_group(atom_group const *ag)
00591 {
00592 std::vector<int> const &source_ids = ag->atoms_ids;
00593 
00594 if (source_ids.size()) {
00595 atoms_ids.reserve(atoms_ids.size()+source_ids.size());
00596 
00597 if (is_enabled(f_ag_scalable)) {
00598 for (size_t i = 0; i < source_ids.size(); i++) {
00599 add_atom_id(source_ids[i]);
00600 }
00601 } else {
00602 atoms.reserve(atoms.size()+source_ids.size());
00603 for (size_t i = 0; i < source_ids.size(); i++) {
00604 // We could use the atom copy constructor, but only if the source
00605 // group is not scalable - whereas this works in both cases
00606 // atom constructor expects 1-based atom number
00607 add_atom(cvm::atom(source_ids[i] + 1));
00608 }
00609 }
00610 
00611 if (cvm::get_error()) return COLVARS_ERROR;
00612 } else {
00613 cvm::error("Error: source atom group contains no atoms\".\n", COLVARS_INPUT_ERROR);
00614 return COLVARS_ERROR;
00615 }
00616 
00617 return COLVARS_OK;
00618 }
00619 
00620 
00621 int cvm::atom_group::add_atom_numbers(std::string const &numbers_conf)
00622 {
00623 std::vector<int> atom_indexes;
00624 
00625 if (numbers_conf.size()) {
00626 std::istringstream is(numbers_conf);
00627 int ia;
00628 while (is >> ia) {
00629 atom_indexes.push_back(ia);
00630 }
00631 }
00632 
00633 if (atom_indexes.size()) {
00634 atoms_ids.reserve(atoms_ids.size()+atom_indexes.size());
00635 
00636 if (is_enabled(f_ag_scalable)) {
00637 for (size_t i = 0; i < atom_indexes.size(); i++) {
00638 add_atom_id((cvm::proxy)->check_atom_id(atom_indexes[i]));
00639 }
00640 } else {
00641 // if we are handling the group on rank 0, better allocate the vector in one shot
00642 atoms.reserve(atoms.size()+atom_indexes.size());
00643 for (size_t i = 0; i < atom_indexes.size(); i++) {
00644 add_atom(cvm::atom(atom_indexes[i]));
00645 }
00646 }
00647 
00648 if (cvm::get_error()) return COLVARS_ERROR;
00649 } else {
00650 cvm::error("Error: no numbers provided for \""
00651 "atomNumbers\".\n", COLVARS_INPUT_ERROR);
00652 return COLVARS_ERROR;
00653 }
00654 
00655 return COLVARS_OK;
00656 }
00657 
00658 
00659 int cvm::atom_group::add_index_group(std::string const &index_group_name)
00660 {
00661 std::vector<std::string> const &index_group_names =
00662 cvm::main()->index_group_names;
00663 std::vector<std::vector<int> *> const &index_groups =
00664 cvm::main()->index_groups;
00665 
00666 size_t i_group = 0;
00667 for ( ; i_group < index_groups.size(); i_group++) {
00668 if (index_group_names[i_group] == index_group_name)
00669 break;
00670 }
00671 
00672 if (i_group >= index_group_names.size()) {
00673 return cvm::error("Error: could not find index group "+
00674 index_group_name+" among those already provided.\n",
00675 COLVARS_INPUT_ERROR);
00676 }
00677 
00678 int error_code = COLVARS_OK;
00679 
00680 std::vector<int> const &index_group = *(index_groups[i_group]);
00681 
00682 atoms_ids.reserve(atoms_ids.size()+index_group.size());
00683 
00684 if (is_enabled(f_ag_scalable)) {
00685 for (size_t i = 0; i < index_group.size(); i++) {
00686 error_code |=
00687 add_atom_id((cvm::proxy)->check_atom_id(index_group[i]));
00688 }
00689 } else {
00690 atoms.reserve(atoms.size()+index_group.size());
00691 for (size_t i = 0; i < index_group.size(); i++) {
00692 error_code |= add_atom(cvm::atom(index_group[i]));
00693 }
00694 }
00695 
00696 return error_code;
00697 }
00698 
00699 
00700 int cvm::atom_group::add_atom_numbers_range(std::string const &range_conf)
00701 {
00702 if (range_conf.size()) {
00703 std::istringstream is(range_conf);
00704 int initial, final;
00705 char dash;
00706 if ( (is >> initial) && (initial > 0) &&
00707 (is >> dash) && (dash == '-') &&
00708 (is >> final) && (final > 0) ) {
00709 
00710 atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
00711 
00712 if (is_enabled(f_ag_scalable)) {
00713 for (int anum = initial; anum <= final; anum++) {
00714 add_atom_id((cvm::proxy)->check_atom_id(anum));
00715 }
00716 } else {
00717 atoms.reserve(atoms.size() + (final - initial + 1));
00718 for (int anum = initial; anum <= final; anum++) {
00719 add_atom(cvm::atom(anum));
00720 }
00721 }
00722 
00723 }
00724 if (cvm::get_error()) return COLVARS_ERROR;
00725 } else {
00726 cvm::error("Error: no valid definition for \"atomNumbersRange\", \""+
00727 range_conf+"\".\n", COLVARS_INPUT_ERROR);
00728 return COLVARS_ERROR;
00729 }
00730 
00731 return COLVARS_OK;
00732 }
00733 
00734 
00735 int cvm::atom_group::add_atom_name_residue_range(std::string const &psf_segid,
00736 std::string const &range_conf)
00737 {
00738 if (range_conf.size()) {
00739 std::istringstream is(range_conf);
00740 std::string atom_name;
00741 int initial, final;
00742 char dash;
00743 if ( (is >> atom_name) && (atom_name.size()) &&
00744 (is >> initial) && (initial > 0) &&
00745 (is >> dash) && (dash == '-') &&
00746 (is >> final) && (final > 0) ) {
00747 
00748 atoms_ids.reserve(atoms_ids.size() + (final - initial + 1));
00749 
00750 if (is_enabled(f_ag_scalable)) {
00751 for (int resid = initial; resid <= final; resid++) {
00752 add_atom_id((cvm::proxy)->check_atom_id(resid, atom_name, psf_segid));
00753 }
00754 } else {
00755 atoms.reserve(atoms.size() + (final - initial + 1));
00756 for (int resid = initial; resid <= final; resid++) {
00757 add_atom(cvm::atom(resid, atom_name, psf_segid));
00758 }
00759 }
00760 
00761 if (cvm::get_error()) return COLVARS_ERROR;
00762 } else {
00763 cvm::error("Error: cannot parse definition for \""
00764 "atomNameResidueRange\", \""+
00765 range_conf+"\".\n");
00766 return COLVARS_ERROR;
00767 }
00768 } else {
00769 cvm::error("Error: atomNameResidueRange with empty definition.\n");
00770 return COLVARS_ERROR;
00771 }
00772 return COLVARS_OK;
00773 }
00774 
00775 
00776 std::string const cvm::atom_group::print_atom_ids() const
00777 {
00778 size_t line_count = 0;
00779 std::ostringstream os("");
00780 for (size_t i = 0; i < atoms_ids.size(); i++) {
00781 os << " " << std::setw(9) << atoms_ids[i];
00782 if (++line_count == 7) {
00783 os << "\n";
00784 line_count = 0;
00785 }
00786 }
00787 return os.str();
00788 }
00789 
00790 
00791 int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
00792 {
00793 if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) {
00794 
00795 if (b_dummy)
00796 cvm::error("Error: centerToReference or rotateToReference "
00797 "cannot be defined for a dummy atom.\n");
00798 
00799 bool b_ref_pos_group = false;
00800 std::string fitting_group_conf;
00801 if (key_lookup(group_conf, "refPositionsGroup", &fitting_group_conf)) {
00802 b_ref_pos_group = true;
00803 cvm::log("Warning: keyword \"refPositionsGroup\" is deprecated: please use \"fittingGroup\" instead.\n");
00804 }
00805 
00806 if (b_ref_pos_group || key_lookup(group_conf, "fittingGroup", &fitting_group_conf)) {
00807 // instead of this group, define another group to compute the fit
00808 if (fitting_group) {
00809 cvm::error("Error: the atom group \""+
00810 key+"\" has already a reference group "
00811 "for the rototranslational fit, which was communicated by the "
00812 "colvar component. You should not use fittingGroup "
00813 "in this case.\n", COLVARS_INPUT_ERROR);
00814 return COLVARS_INPUT_ERROR;
00815 }
00816 cvm::log("Within atom group \""+key+"\":\n");
00817 fitting_group = new atom_group("fittingGroup");
00818 if (fitting_group->parse(fitting_group_conf) == COLVARS_OK) {
00819 fitting_group->check_keywords(fitting_group_conf, "fittingGroup");
00820 if (cvm::get_error()) {
00821 cvm::error("Error setting up atom group \"fittingGroup\".", COLVARS_INPUT_ERROR);
00822 return COLVARS_INPUT_ERROR;
00823 }
00824 }
00825 enable(f_ag_fitting_group);
00826 }
00827 
00828 atom_group *group_for_fit = fitting_group ? fitting_group : this;
00829 
00830 get_keyval(group_conf, "refPositions", ref_pos, ref_pos);
00831 
00832 std::string ref_pos_file;
00833 if (get_keyval(group_conf, "refPositionsFile", ref_pos_file, std::string(""))) {
00834 
00835 if (ref_pos.size()) {
00836 cvm::error("Error: cannot specify \"refPositionsFile\" and "
00837 "\"refPositions\" at the same time.\n");
00838 }
00839 
00840 std::string ref_pos_col;
00841 double ref_pos_col_value=0.0;
00842 
00843 if (get_keyval(group_conf, "refPositionsCol", ref_pos_col, std::string(""))) {
00844 // if provided, use PDB column to select coordinates
00845 bool found = get_keyval(group_conf, "refPositionsColValue", ref_pos_col_value, 0.0);
00846 if (found && ref_pos_col_value == 0.0) {
00847 cvm::error("Error: refPositionsColValue, "
00848 "if provided, must be non-zero.\n", COLVARS_INPUT_ERROR);
00849 return COLVARS_ERROR;
00850 }
00851 }
00852 
00853 ref_pos.resize(group_for_fit->size());
00854 cvm::load_coords(ref_pos_file.c_str(), &ref_pos, group_for_fit,
00855 ref_pos_col, ref_pos_col_value);
00856 }
00857 
00858 if (ref_pos.size()) {
00859 
00860 if (is_enabled(f_ag_rotate)) {
00861 if (ref_pos.size() != group_for_fit->size())
00862 cvm::error("Error: the number of reference positions provided("+
00863 cvm::to_str(ref_pos.size())+
00864 ") does not match the number of atoms within \""+
00865 key+
00866 "\" ("+cvm::to_str(group_for_fit->size())+
00867 "): to perform a rotational fit, "+
00868 "these numbers should be equal.\n", COLVARS_INPUT_ERROR);
00869 }
00870 
00871 // save the center of geometry of ref_pos and subtract it
00872 center_ref_pos();
00873 
00874 } else {
00875 cvm::error("Error: no reference positions provided.\n", COLVARS_INPUT_ERROR);
00876 return COLVARS_ERROR;
00877 }
00878 
00879 if (is_enabled(f_ag_rotate) && !noforce) {
00880 cvm::log("Warning: atom group \""+key+
00881 "\" will be aligned to a fixed orientation given by the reference positions provided. "
00882 "If the internal structure of the group changes too much (i.e. its RMSD is comparable "
00883 "to its radius of gyration), the optimal rotation and its gradients may become discontinuous. "
00884 "If that happens, use fittingGroup (or a different definition for it if already defined) "
00885 "to align the coordinates.\n");
00886 // initialize rot member data
00887 rot.request_group1_gradients(group_for_fit->size());
00888 }
00889 }
00890 
00891 // Enable fit gradient calculation only if necessary, and not disabled by the user
00892 // This must happen after fitting group is defined so that side-effects are performed
00893 // properly (ie. allocating fitting group gradients)
00894 {
00895 bool b_fit_gradients;
00896 get_keyval(group_conf, "enableFitGradients", b_fit_gradients, true);
00897 
00898 if (b_fit_gradients && (is_enabled(f_ag_center) || is_enabled(f_ag_rotate))) {
00899 enable(f_ag_fit_gradients);
00900 }
00901 }
00902 
00903 return COLVARS_OK;
00904 }
00905 
00906 
00907 void cvm::atom_group::do_feature_side_effects(int id)
00908 {
00909 // If enabled features are changed upstream, the features below should be refreshed
00910 switch (id) {
00911 case f_ag_fit_gradients:
00912 if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) {
00913 atom_group *group_for_fit = fitting_group ? fitting_group : this;
00914 group_for_fit->fit_gradients.assign(group_for_fit->size(), cvm::atom_pos(0.0, 0.0, 0.0));
00915 rot.request_group1_gradients(group_for_fit->size());
00916 }
00917 break;
00918 }
00919 }
00920 
00921 
00922 int cvm::atom_group::create_sorted_ids()
00923 {
00924 // Only do the work if the vector is not yet populated
00925 if (sorted_atoms_ids.size())
00926 return COLVARS_OK;
00927 
00928 // Sort the internal IDs
00929 std::list<int> sorted_atoms_ids_list;
00930 for (size_t i = 0; i < atoms_ids.size(); i++) {
00931 sorted_atoms_ids_list.push_back(atoms_ids[i]);
00932 }
00933 sorted_atoms_ids_list.sort();
00934 sorted_atoms_ids_list.unique();
00935 if (sorted_atoms_ids_list.size() != atoms_ids.size()) {
00936 return cvm::error("Error: duplicate atom IDs in atom group? (found " +
00937 cvm::to_str(sorted_atoms_ids_list.size()) +
00938 " unique atom IDs instead of " +
00939 cvm::to_str(atoms_ids.size()) + ").\n", COLVARS_BUG_ERROR);
00940 }
00941 
00942 // Compute map between sorted and unsorted elements
00943 sorted_atoms_ids.resize(atoms_ids.size());
00944 sorted_atoms_ids_map.resize(atoms_ids.size());
00945 std::list<int>::iterator lsii = sorted_atoms_ids_list.begin();
00946 size_t ii = 0;
00947 for ( ; ii < atoms_ids.size(); lsii++, ii++) {
00948 sorted_atoms_ids[ii] = *lsii;
00949 size_t const pos = std::find(atoms_ids.begin(), atoms_ids.end(), *lsii) -
00950 atoms_ids.begin();
00951 sorted_atoms_ids_map[ii] = pos;
00952 }
00953 
00954 return COLVARS_OK;
00955 }
00956 
00957 
00958 int cvm::atom_group::overlap(const atom_group &g1, const atom_group &g2){
00959 for (cvm::atom_const_iter ai1 = g1.begin(); ai1 != g1.end(); ai1++) {
00960 for (cvm::atom_const_iter ai2 = g2.begin(); ai2 != g2.end(); ai2++) {
00961 if (ai1->id == ai2->id) {
00962 return (ai1->id + 1); // 1-based index to allow boolean usage
00963 }
00964 }
00965 }
00966 return 0;
00967 }
00968 
00969 
00970 void cvm::atom_group::center_ref_pos()
00971 {
00972 ref_pos_cog = cvm::atom_pos(0.0, 0.0, 0.0);
00973 std::vector<cvm::atom_pos>::iterator pi;
00974 for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
00975 ref_pos_cog += *pi;
00976 }
00977 ref_pos_cog /= (cvm::real) ref_pos.size();
00978 for (pi = ref_pos.begin(); pi != ref_pos.end(); ++pi) {
00979 (*pi) -= ref_pos_cog;
00980 }
00981 }
00982 
00983 
00984 void cvm::atom_group::read_positions()
00985 {
00986 if (b_dummy) return;
00987 
00988 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
00989 ai->read_position();
00990 }
00991 
00992 if (fitting_group)
00993 fitting_group->read_positions();
00994 }
00995 
00996 
00997 int cvm::atom_group::calc_required_properties()
00998 {
00999 // TODO check if the com is needed?
01000 calc_center_of_mass();
01001 calc_center_of_geometry();
01002 
01003 if (!is_enabled(f_ag_scalable)) {
01004 if (is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) {
01005 if (fitting_group) {
01006 fitting_group->calc_center_of_geometry();
01007 }
01008 
01009 calc_apply_roto_translation();
01010 
01011 // update COM and COG after fitting
01012 calc_center_of_geometry();
01013 calc_center_of_mass();
01014 if (fitting_group) {
01015 fitting_group->calc_center_of_geometry();
01016 }
01017 }
01018 }
01019 
01020 // TODO calculate elements of scalable cvc's here before reduction
01021 
01022 return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
01023 }
01024 
01025 
01026 void cvm::atom_group::calc_apply_roto_translation()
01027 {
01028 // store the laborarory-frame COGs for when they are needed later
01029 cog_orig = this->center_of_geometry();
01030 if (fitting_group) {
01031 fitting_group->cog_orig = fitting_group->center_of_geometry();
01032 }
01033 
01034 if (is_enabled(f_ag_center)) {
01035 // center on the origin first
01036 cvm::atom_pos const rpg_cog = fitting_group ?
01037 fitting_group->center_of_geometry() : this->center_of_geometry();
01038 apply_translation(-1.0 * rpg_cog);
01039 if (fitting_group) {
01040 fitting_group->apply_translation(-1.0 * rpg_cog);
01041 }
01042 }
01043 
01044 if (is_enabled(f_ag_rotate)) {
01045 // rotate the group (around the center of geometry if f_ag_center is
01046 // enabled, around the origin otherwise)
01047 rot.calc_optimal_rotation(fitting_group ?
01048 fitting_group->positions() :
01049 this->positions(),
01050 ref_pos);
01051 
01052 cvm::atom_iter ai;
01053 for (ai = this->begin(); ai != this->end(); ai++) {
01054 ai->pos = rot.rotate(ai->pos);
01055 }
01056 if (fitting_group) {
01057 for (ai = fitting_group->begin(); ai != fitting_group->end(); ai++) {
01058 ai->pos = rot.rotate(ai->pos);
01059 }
01060 }
01061 }
01062 
01063 if (is_enabled(f_ag_center) && !is_enabled(f_ag_center_origin)) {
01064 // align with the center of geometry of ref_pos
01065 apply_translation(ref_pos_cog);
01066 if (fitting_group) {
01067 fitting_group->apply_translation(ref_pos_cog);
01068 }
01069 }
01070 // update of COM and COG is done from the calling routine
01071 }
01072 
01073 
01074 void cvm::atom_group::apply_translation(cvm::rvector const &t)
01075 {
01076 if (b_dummy) {
01077 cvm::error("Error: cannot translate the coordinates of a dummy atom group.\n", COLVARS_INPUT_ERROR);
01078 return;
01079 }
01080 
01081 if (is_enabled(f_ag_scalable)) {
01082 cvm::error("Error: cannot translate the coordinates of a scalable atom group.\n", COLVARS_INPUT_ERROR);
01083 return;
01084 }
01085 
01086 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01087 ai->pos += t;
01088 }
01089 }
01090 
01091 
01092 void cvm::atom_group::read_velocities()
01093 {
01094 if (b_dummy) return;
01095 
01096 if (is_enabled(f_ag_rotate)) {
01097 
01098 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01099 ai->read_velocity();
01100 ai->vel = rot.rotate(ai->vel);
01101 }
01102 
01103 } else {
01104 
01105 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01106 ai->read_velocity();
01107 }
01108 }
01109 }
01110 
01111 
01112 // TODO make this a calc function
01113 void cvm::atom_group::read_total_forces()
01114 {
01115 if (b_dummy) return;
01116 
01117 if (is_enabled(f_ag_rotate)) {
01118 
01119 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01120 ai->read_total_force();
01121 ai->total_force = rot.rotate(ai->total_force);
01122 }
01123 
01124 } else {
01125 
01126 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01127 ai->read_total_force();
01128 }
01129 }
01130 }
01131 
01132 
01133 int cvm::atom_group::calc_center_of_geometry()
01134 {
01135 if (b_dummy) {
01136 cog = dummy_atom_pos;
01137 } else {
01138 cog.reset();
01139 for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
01140 cog += ai->pos;
01141 }
01142 cog /= cvm::real(this->size());
01143 }
01144 return COLVARS_OK;
01145 }
01146 
01147 
01148 int cvm::atom_group::calc_center_of_mass()
01149 {
01150 if (b_dummy) {
01151 com = dummy_atom_pos;
01152 if (cvm::debug()) {
01153 cvm::log("Dummy atom center of mass = "+cvm::to_str(com)+"\n");
01154 }
01155 } else if (is_enabled(f_ag_scalable)) {
01156 com = (cvm::proxy)->get_atom_group_com(index);
01157 } else {
01158 com.reset();
01159 for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
01160 com += ai->mass * ai->pos;
01161 }
01162 com /= total_mass;
01163 }
01164 return COLVARS_OK;
01165 }
01166 
01167 
01168 int cvm::atom_group::calc_dipole(cvm::atom_pos const &dipole_center)
01169 {
01170 if (b_dummy) {
01171 return cvm::error("Error: trying to compute the dipole "
01172 "of a dummy group.\n", COLVARS_INPUT_ERROR);
01173 }
01174 dip.reset();
01175 for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
01176 dip += ai->charge * (ai->pos - dipole_center);
01177 }
01178 return COLVARS_OK;
01179 }
01180 
01181 
01182 void cvm::atom_group::set_weighted_gradient(cvm::rvector const &grad)
01183 {
01184 if (b_dummy) return;
01185 
01186 scalar_com_gradient = grad;
01187 
01188 if (!is_enabled(f_ag_scalable)) {
01189 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01190 ai->grad = (ai->mass/total_mass) * grad;
01191 }
01192 }
01193 }
01194 
01195 
01196 void cvm::atom_group::calc_fit_gradients()
01197 {
01198 if (b_dummy || ! is_enabled(f_ag_fit_gradients)) return;
01199 
01200 if (cvm::debug())
01201 cvm::log("Calculating fit gradients.\n");
01202 
01203 cvm::atom_group *group_for_fit = fitting_group ? fitting_group : this;
01204 
01205 if (is_enabled(f_ag_center)) {
01206 // add the center of geometry contribution to the gradients
01207 cvm::rvector atom_grad;
01208 
01209 for (size_t i = 0; i < this->size(); i++) {
01210 atom_grad += atoms[i].grad;
01211 }
01212 if (is_enabled(f_ag_rotate)) atom_grad = (rot.inverse()).rotate(atom_grad);
01213 atom_grad *= (-1.0)/(cvm::real(group_for_fit->size()));
01214 
01215 for (size_t j = 0; j < group_for_fit->size(); j++) {
01216 group_for_fit->fit_gradients[j] = atom_grad;
01217 }
01218 }
01219 
01220 if (is_enabled(f_ag_rotate)) {
01221 
01222 // add the rotation matrix contribution to the gradients
01223 cvm::rotation const rot_inv = rot.inverse();
01224 
01225 for (size_t i = 0; i < this->size(); i++) {
01226 
01227 // compute centered, unrotated position
01228 cvm::atom_pos const pos_orig =
01229 rot_inv.rotate((is_enabled(f_ag_center) ? (atoms[i].pos - ref_pos_cog) : (atoms[i].pos)));
01230 
01231 // calculate \partial(R(q) \vec{x}_i)/\partial q) \cdot \partial\xi/\partial\vec{x}_i
01232 cvm::quaternion const dxdq =
01233 rot.q.position_derivative_inner(pos_orig, atoms[i].grad);
01234 
01235 for (size_t j = 0; j < group_for_fit->size(); j++) {
01236 // multiply by {\partial q}/\partial\vec{x}_j and add it to the fit gradients
01237 for (size_t iq = 0; iq < 4; iq++) {
01238 group_for_fit->fit_gradients[j] += dxdq[iq] * rot.dQ0_1[j][iq];
01239 }
01240 }
01241 }
01242 }
01243 
01244 if (cvm::debug())
01245 cvm::log("Done calculating fit gradients.\n");
01246 }
01247 
01248 
01249 std::vector<cvm::atom_pos> cvm::atom_group::positions() const
01250 {
01251 if (b_dummy) {
01252 cvm::error("Error: positions are not available "
01253 "from a dummy atom group.\n", COLVARS_INPUT_ERROR);
01254 }
01255 
01256 if (is_enabled(f_ag_scalable)) {
01257 cvm::error("Error: atomic positions are not available "
01258 "from a scalable atom group.\n", COLVARS_INPUT_ERROR);
01259 }
01260 
01261 std::vector<cvm::atom_pos> x(this->size(), 0.0);
01262 cvm::atom_const_iter ai = this->begin();
01263 std::vector<cvm::atom_pos>::iterator xi = x.begin();
01264 for ( ; ai != this->end(); ++xi, ++ai) {
01265 *xi = ai->pos;
01266 }
01267 return x;
01268 }
01269 
01270 std::vector<cvm::atom_pos> cvm::atom_group::positions_shifted(cvm::rvector const &shift) const
01271 {
01272 if (b_dummy) {
01273 cvm::error("Error: positions are not available "
01274 "from a dummy atom group.\n", COLVARS_INPUT_ERROR);
01275 }
01276 
01277 if (is_enabled(f_ag_scalable)) {
01278 cvm::error("Error: atomic positions are not available "
01279 "from a scalable atom group.\n", COLVARS_INPUT_ERROR);
01280 }
01281 
01282 std::vector<cvm::atom_pos> x(this->size(), 0.0);
01283 cvm::atom_const_iter ai = this->begin();
01284 std::vector<cvm::atom_pos>::iterator xi = x.begin();
01285 for ( ; ai != this->end(); ++xi, ++ai) {
01286 *xi = (ai->pos + shift);
01287 }
01288 return x;
01289 }
01290 
01291 std::vector<cvm::rvector> cvm::atom_group::velocities() const
01292 {
01293 if (b_dummy) {
01294 cvm::error("Error: velocities are not available "
01295 "from a dummy atom group.\n", COLVARS_INPUT_ERROR);
01296 }
01297 
01298 if (is_enabled(f_ag_scalable)) {
01299 cvm::error("Error: atomic velocities are not available "
01300 "from a scalable atom group.\n", COLVARS_INPUT_ERROR);
01301 }
01302 
01303 std::vector<cvm::rvector> v(this->size(), 0.0);
01304 cvm::atom_const_iter ai = this->begin();
01305 std::vector<cvm::atom_pos>::iterator vi = v.begin();
01306 for ( ; ai != this->end(); vi++, ai++) {
01307 *vi = ai->vel;
01308 }
01309 return v;
01310 }
01311 
01312 std::vector<cvm::rvector> cvm::atom_group::total_forces() const
01313 {
01314 if (b_dummy) {
01315 cvm::error("Error: total forces are not available "
01316 "from a dummy atom group.\n", COLVARS_INPUT_ERROR);
01317 }
01318 
01319 if (is_enabled(f_ag_scalable)) {
01320 cvm::error("Error: atomic total forces are not available "
01321 "from a scalable atom group.\n", COLVARS_INPUT_ERROR);
01322 }
01323 
01324 std::vector<cvm::rvector> f(this->size(), 0.0);
01325 cvm::atom_const_iter ai = this->begin();
01326 std::vector<cvm::atom_pos>::iterator fi = f.begin();
01327 for ( ; ai != this->end(); ++fi, ++ai) {
01328 *fi = ai->total_force;
01329 }
01330 return f;
01331 }
01332 
01333 
01334 // TODO make this an accessor
01335 cvm::rvector cvm::atom_group::total_force() const
01336 {
01337 if (b_dummy) {
01338 cvm::error("Error: total total forces are not available "
01339 "from a dummy atom group.\n", COLVARS_INPUT_ERROR);
01340 }
01341 
01342 if (is_enabled(f_ag_scalable)) {
01343 return (cvm::proxy)->get_atom_group_total_force(index);
01344 }
01345 
01346 cvm::rvector f(0.0);
01347 for (cvm::atom_const_iter ai = this->begin(); ai != this->end(); ai++) {
01348 f += ai->total_force;
01349 }
01350 return f;
01351 }
01352 
01353 
01354 void cvm::atom_group::apply_colvar_force(cvm::real const &force)
01355 {
01356 if (cvm::debug()) {
01357 log("Communicating a colvar force from atom group to the MD engine.\n");
01358 }
01359 
01360 if (b_dummy) return;
01361 
01362 if (noforce) {
01363 cvm::error("Error: sending a force to a group that has "
01364 "\"enableForces\" set to off.\n");
01365 return;
01366 }
01367 
01368 if (is_enabled(f_ag_scalable)) {
01369 (cvm::proxy)->apply_atom_group_force(index, force * scalar_com_gradient);
01370 return;
01371 }
01372 
01373 if (is_enabled(f_ag_rotate)) {
01374 
01375 // rotate forces back to the original frame
01376 cvm::rotation const rot_inv = rot.inverse();
01377 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01378 ai->apply_force(rot_inv.rotate(force * ai->grad));
01379 }
01380 
01381 } else {
01382 
01383 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01384 ai->apply_force(force * ai->grad);
01385 }
01386 }
01387 
01388 if ((is_enabled(f_ag_center) || is_enabled(f_ag_rotate)) && is_enabled(f_ag_fit_gradients)) {
01389 
01390 atom_group *group_for_fit = fitting_group ? fitting_group : this;
01391 
01392 // Fit gradients are already calculated in "laboratory" frame
01393 for (size_t j = 0; j < group_for_fit->size(); j++) {
01394 (*group_for_fit)[j].apply_force(force * group_for_fit->fit_gradients[j]);
01395 }
01396 }
01397 }
01398 
01399 
01400 void cvm::atom_group::apply_force(cvm::rvector const &force)
01401 {
01402 if (cvm::debug()) {
01403 log("Communicating a colvar force from atom group to the MD engine.\n");
01404 }
01405 
01406 if (b_dummy) return;
01407 
01408 if (noforce) {
01409 cvm::error("Error: sending a force to a group that has "
01410 "\"enableForces\" set to off.\n");
01411 return;
01412 }
01413 
01414 if (is_enabled(f_ag_scalable)) {
01415 (cvm::proxy)->apply_atom_group_force(index, force);
01416 return;
01417 }
01418 
01419 if (is_enabled(f_ag_rotate)) {
01420 
01421 cvm::rotation const rot_inv = rot.inverse();
01422 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01423 ai->apply_force(rot_inv.rotate((ai->mass/total_mass) * force));
01424 }
01425 
01426 } else {
01427 
01428 for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
01429 ai->apply_force((ai->mass/total_mass) * force);
01430 }
01431 }
01432 }
01433 
01434 
01435 // Static members
01436 
01437 std::vector<colvardeps::feature *> cvm::atom_group::ag_features;
01438 
01439