Main Page   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

colvaratoms.C

Go to the documentation of this file.
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 

Generated on Thu Mar 28 02:42:44 2024 for VMD (current) by doxygen1.2.14 written by Dimitri van Heesch, © 1997-2002